! Last change: DOS 31 Jul 2000 7:47 pm ! *** copyright 2000 *** ! *** filename qgaustb.f95 *** John F. Monahan ** ! ********************** program qgaustb ! ! compute gauss quadrature tables of abscissas and weights ! for gauss-hermite ! gauss-hermite (modified by rescale) ! gauss-legendre (modified by shift) ! ! implicit none real(KIND=8), parameter :: sqrpi = 1.772453850905516d0 real(KIND=8), dimension(20) :: d,g real(KIND=8), dimension(20,20) :: q integer k,j,n,ktype ! INTERFACE BLOCK INTERFACE SUBROUTINE QREIG1(A,B,N,MIT,Q) REAL(KIND=8), DIMENSION(:), INTENT(IN OUT) :: A,B REAL(KIND=8), DIMENSION(:,:), INTENT(IN OUT) :: Q INTEGER, INTENT(IN) :: N,MIT END SUBROUTINE QREIG1 END INTERFACE ! 20 format(//' size of problem',i4,4x,'type',i2,4x,'found',f4.0) 21 format(5f16.8) 22 format(i4,2f16.8) 29 format('eigenvectors') ! output unit open( unit=6, file='qgaustb.out' ) ! make some tables open( unit=8, file='qgaushm.tab' ) open( unit=9, file='qgauslg.tab' ) ! do a few values of n do n = 3,20 ! do three types ! ktype = 1 is the regular gauss-hermite ! ktype = 2 is the rescaled gauss-hermite ! ktype = 3 is the shifted gauss-legendre do ktype = 1,3 ! do k = 1,n ! diagonal and off diagonal d(k) = 0.d0 if( ktype .eq. 3 ) d(k) = .5d0 g(k) = sqrt( real(k,8) /2.d0 ) !KIND if( ktype .eq. 2 ) g(k) = sqrt( real(k,8) ) !KIND if( ktype .eq. 3 ) & & g(k) = (real(k,8)/2.d0)/sqrt(real(4*k*k-1,8)) !KIND ! no need to tridiagonalize, but get identity do j = 1,n q(k,j) = 0.d0 end do ! loop on j q(k,k) = 1.d0 end do ! loop on k ! call eigenvalue and eigenvector routine ! already tridiagonalized, q has identity matrix ! max 100 iterations call qreig1(d,g,n,100,q) write(6,20) n,ktype,g(n) ! now get the weights from the eigenvectors do k = 1,n g(k) = q(1,k)*q(1,k) ! renormalize only for usual gauss-hermite: sqrt(pi) if( ktype .eq. 1 ) g(k) = g(k)*sqrpi write(6,22) k,d(k),g(k) ! write again in table files if( ktype .eq. 2 ) write(8,22) k,d(k),g(k) if( ktype .eq. 3 ) write(9,22) k,d(k),g(k) end do ! loop on k ! end do ! loop on ktype ! end do ! loop on n stop end program qgaustb SUBROUTINE QREIG1(A,B,N,MIT,Q) ! COMPUTES EIGENVALUES AND VECTORS OF TRIDIAGONAL MATRIX OF ORDER N ! A(1) ... A(N) HOLDS THE DIAGONAL ELEMENTS ! B(1) ... B(N-1) HOLDS THE OFF DIAGONAL ELEMENTS ! MIT IS THE MAXIMUM NUMBER OF ITERATIONS (2N TO 3N IS GOOD) ! ! *** ON INPUT *** MATRIX Q MUST HOLD EITHER OF THE FOLLOWING ! 1) IDENTITY MATRIX OR ! 2) ORTHOGONAL MATRIX THAT FORMED TRIDIAGONAL FROM TRIDIG AND EXPNDH ! *** ON OUTPUT *** MATRIX Q WILL HOLD EITHER ! 1) EIGENVECTORS OF TRIDIAGONAL MATRIX (STORED AS COLUMNS) ! 2) EIGENVECTORS OF ORIGINAL MATRIX (BEFORE TRIDIAGONALIZATION) ! ! ! ON OUTPUT, B(N) COUNTS NUMBER OF EIGENVALUES FOUND BEFORE MIT ! EIGENVALUES FOUND IN A(1) ... A(N) ON SUCCESSFUL EXIT ! IF UNSUCCESSFUL EXIT ( B(N) NOT FLOAT(N) ) THEN EIGENVALUES FOUND ! WILL BE AT END OF A A(N-FOUND+1) ... A(N) ! ! METHOD IS IMPLICIT SYMMETRIC QR ALGORITHM WITH WILKINSON SHIFT ! SEE STEWART (7.4.3) OR GOLUB AND VAN LOAN (8.2-2) ! ! J F MONAHAN (JUNE 1986) Edited July, 1996 ! RECODED OCTOBER, 1999 FOR FORTRAN 95, correction March 2000 ! IMPLICIT NONE REAL(KIND=8), DIMENSION(:), INTENT(IN OUT) :: A,B REAL(KIND=8), DIMENSION(:,:), INTENT(IN OUT) :: Q INTEGER, INTENT(IN) :: N,MIT ! UNITM IS MACHINE UNIT FOR THIS COMPUTER REAL(KIND=8), PARAMETER :: UNITM = 1.d-15 REAL(KIND=8) D,SHIFT,PI,RHO,TAU,OMEGA,RNU,SIGMA,GAMA,SIG2, & & GAM2,GAMSIG,QJK,QJKP1 INTEGER NCUR,NCURM1,J,K,KP1,KIT ! ! COUNT THE NUMBER OF ITERATIONS KIT = 0 ! COUNT THE NUMBER OF EIGENVALUES FOUND B(N) = 0. ! NCUR = N NCURM1 = NCUR - 1 ! ! RESTART HERE DO WHILE( KIT .LE. MIT ) ! TEST FIRST BEFORE DOING ANY WORK IF( ABS(B(NCURM1)) .LE. & & UNITM*( ABS(A(NCURM1)) + ABS(A(NCUR)) ) ) THEN ! ! DEFLATE PROBLEM BY SHORTENING BY ONE B(N) = B(N) + 1. NCUR = NCURM1 IF( NCUR .EQ. 1 ) THEN ! SUCCESSFUL ** AND ALL DONE B(N) = N RETURN END IF NCURM1 = NCUR - 1 ELSE ! ! START THE ITERATION KIT = KIT + 1 ! HAVE WE DONE ALL WE SHOULD DO IF( KIT .EQ. MIT ) RETURN ! ! COMPUTE THE SHIFT D = ( A(NCURM1) - A(NCUR) ) / 2. ! USE WILKINSON SHIFT SHIFT = A(NCUR) + D & & - SIGN( SQRT( D*D + B(NCURM1)*B(NCURM1) ), D ) ! INITIAL VALUES DEFINE FIRST TRANSFORMATION PI = A(1) - SHIFT RHO = B(1) TAU = A(1) OMEGA = B(1) ! DO GIVENS ROTATIONS TO GET QR FACTORIZATION DO K = 1,NCURM1 KP1 = K + 1 ! COMPUTE THE ROTATION CALL ROT734(PI,RHO,GAMA,SIGMA,RNU) ! NEW OFF DIAGONAL ELEMENT IF( K .NE. 1 ) B(K-1) = RNU ! GET SOME COMMON FACTORS FIRST GAM2 = GAMA * GAMA SIG2 = SIGMA * SIGMA GAMSIG = GAMA * SIGMA ! NEW DIAGONAL ELEMENT A(K) = GAM2*TAU + 2.*GAMSIG*OMEGA + SIG2*A(KP1) ! NEW VALUES FOR NEXT TRANSFORMATION PI = (GAM2-SIG2)*OMEGA + GAMSIG*( A(KP1) - TAU ) TAU = SIG2*TAU - 2.*GAMSIG*OMEGA + GAM2*A(KP1) RHO = SIGMA*B(KP1) OMEGA = GAMA*B(KP1) ! HAVE DONE FRONT AND BACK FOR K AND KP1 ! NOW DO TFM TO Q DO J = 1,N ! SAVE THESE TWO FIRST QJK = Q(J,K) QJKP1 = Q(J,KP1) ! REPLACE WITH ROTATION OF THESE Q(J,K) = GAMA * QJK + SIGMA * QJKP1 Q(J,KP1) = GAMA * QJKP1 - SIGMA * QJK END DO ! LOOP ON J ! END DO ! LOOP ON K ! LAST PART OF ITERATION B(NCURM1) = PI A(NCUR) = TAU END IF ! END OF ITERATION END DO ! WHILE( KIT .LE. MIT) ! FINISH HERE WHEN TOO MANY ITERATIONS RETURN END SUBROUTINE QREIG1 SUBROUTINE ROT734(ALPHA,BETA,GAMA,SIGMA,RNU) ! GIVEN ALPHA AND BETA, PRODUCES GIVENS TRANSFORMATION VALUES ! ! ( GAMA SIGMA ) ( ALPHA ) ( RNU ) ! ( ) ( ) = ( ) ! ( -SIGMA GAMA ) ( BETA ) ( ZERO ) ! ! ALGORITHM 7.3.4 FROM STEWART ! REAL(KIND=8), INTENT(IN) :: ALPHA,BETA REAL(KIND=8), INTENT(OUT) :: GAMA,SIGMA,RNU REAL(KIND=8) ETA,DELTA ! ETA = MAX( ABS(ALPHA), ABS(BETA) ) ! IF BOTH ALPHA AND BETA ARE ZERO THEN ! USE IDENTITY TRANSFORMATION IF( ETA .EQ. 0. ) then ! SET UP IDENTITY TRANSFORMATION IF NO WORK DONE RNU = 0. GAMA = 1. SIGMA = 0. RETURN else DELTA = (ALPHA/ETA)*(ALPHA/ETA) + (BETA/ETA)*(BETA/ETA) DELTA = SQRT(DELTA) RNU = ETA * DELTA GAMA = ALPHA / RNU SIGMA = BETA / RNU END if RETURN END subroutine rot734 ! *** end of filename qgaustb.f95 *********************