! Last change: DOS 27 Jul 2000 8:56 pm ! *** copyright 2000 *** ! *** filename qreig0.f95 *** John F. Monahan ** ! ********************** PROGRAM PQREIG0 ! DEMONSTRATION PROGRAM QREIG0 -- COMPUTES EIGENVALUES AND VECTORS ! FOR A SYMMETRIC MATRIX ! ! N = 3 EXAMPLE 8.2-1 FROM GOLUB AND VAN LOAN ! N = 4 EXAMPLE FROM JENNINGS, TABLE 8.2, P 260 ! N = 5,6 EXAMPLES FROM MARTIN, REINCH, AND WILKINSON, HACLA, P 223 ! IMPLICIT NONE REAL, DIMENSION(21) :: A ! BIGGEST IS 6, 21 = 6*7/2 REAL, DIMENSION(6) :: D,G INTEGER I,J,N ! 20 FORMAT(6F5.1) 21 FORMAT(6F12.6) 22 FORMAT(I3) 25 FORMAT(/' ORIGINAL MATRIX') 26 FORMAT(' TRIDIAGONALIZED, FIRST DIAGONAL, THEN SUPER/SUB') 27 FORMAT(' DONE, FIRST EIGENVALUES, THEN OFF-DIAGONAL') 28 FORMAT(' DOES A*X = X*S ? ROWS COME IN PAIRS, L,R ') 29 FORMAT(' ORTHOGONAL MATRIX X IS NEXT') 30 FORMAT(' CHECKING ORTHOGONALITY -- IS X''X = I ?') ! ! SYMMETRIC MATRIX WITH (I,J) ELEMENT STORED IN LOCATION ! LOC(I,J) = ( I * (I-1) ) /2 + J ! ! FOUR EXAMPLES IN 'eigen.dat' OPEN( UNIT=5, FILE='eigen.dat') OPEN( UNIT=6, FILE='qreig0.out') ! DO N = 3,6 ! READ IN THE MATRIX DO I = 1,N READ(5,20) (A((I*(I-1))/2 + J),J=1,I) END DO ! LOOP ON I ! WRITE IT BACK OUT WRITE(6,25) DO I = 1,N WRITE(6,21) (A((I*(I-1))/2 + J),J=1,I) END DO ! LOOP ON I ! TRIDIAGONALIZE CALL TRIDIG(A,N,D,G) WRITE(6,26) WRITE(6,21) (D(J),J=1,N) WRITE(6,21) (G(J),J=1,N) ! CALL EIGENVALUE ROUTINE CALL QREIG0(D,G,N,15) WRITE(6,27) WRITE(6,21) (D(J),J=1,N) WRITE(6,21) (G(J),J=1,N) ! END OF TEST FOR THIS SIZE MATRIX END DO ! LOOP ON N STOP END PROGRAM PQREIG0 SUBROUTINE QREIG0(A,B,N,MIT) ! COMPUTES EIGENVALUES * ONLY * 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 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, corrections March, April 2000 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N,MIT REAL, DIMENSION( N ), INTENT(IN OUT) :: A,B ! UNITM IS MACHINE UNIT FOR THIS COMPUTER REAL, PARAMETER :: UNITM = 1.E-6 REAL D,SHIFT,PI,RHO,TAU,OMEGA,RNU,SIGMA,GAMA,SIG2, & & GAM2,GAMSIG INTEGER NCUR,NCURM1,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 ! 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 QREIG0 SUBROUTINE TRIDIG(A,N,D,G) ! TRANSFORMS SYMMETRIC MATRIX A (N BY N) TO A TRIDIAGONAL MATRIX ! BY ORTHOGONAL SIMILARITY TRANSFORMATIONS ! ! ALGORITHM 1.2 FROM CHAPTER 7 OF G. W. STEWART, INTO TO MX COMP ! ! A in Matrix in symmetric storage ! out Householder transformations ! N in Dimension of matrix ! D out Diagonal part of tridiagonal matrix ! G out Super/subdiagonal part of tridiagonal matrix ! ! J F MONAHAN (JUNE, 1986) ! recoded October 1999, April 2000 for Fortran 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION( (N*(N+1))/2 ), INTENT(IN OUT) :: A REAL, DIMENSION( N ), INTENT(OUT) :: D,G REAL ETA,R,S INTEGER I,IP1,J,K,KP1,L,LIJ,LIK,LJK,NM2 ! ! A IS SYMMETRIC, SO TO ACCESS A(I,J) USE LOCATION GIVEN BY ! LOC(I,J) = ( I * (I-1) ) /2 + J ! IF( N .GT. 2 ) THEN ! IF N = 1 OR 2, ALREADY TRIDIAGONAL ! NM2 = N-2 DO K = 1,NM2 L = (K*(K+1))/2 ! LOC(K,K) ! COPY DIAGONAL ELEMENT D(K) = A(L) ! DETERMINE TRANSFORMATION KP1 = K + 1 ! FIND LARGEST ELEMENT IN COLUMN - AVOID OVERFLOW ETA = 0. DO I = KP1,N L = (I*(I-1))/2 + K ! LOC(I,K) ETA = MAX( ETA, ABS(A(L)) ) END DO ! LOOP ON I IF( ETA .EQ. 0. ) THEN ! IF ALL ZERO NO TRANSFORMATION NEEDED L = (K*(K+1))/2 ! LOC(K,K) A(L) = 0. G(K) = 0. ELSE ! CONTINUE WITH DETERMINING TRANSFORMATION S = 0. DO I = KP1,N L = (I*(I-1))/2 + K ! LOC(I,K) R = A(L) / ETA A(L) = R S = S + R*R END DO ! LOOP ON I ! L = (KP1*(KP1-1))/2 + K ! LOC(KP1,K) S = SIGN( SQRT( S ), A(L) ) A(L) = A(L) + S R = A(L) L = (K*(K+1))/2 ! LOC(K,K) A(L) = S * R G(K) = - S * ETA ! HAVE PI(K) STORED IN A(K,K) ! U STORED IN COLUMN K, DIAG ELEMENT STORED ! NOW APPLY TFM TO REST OF MATRIX S = 0. DO I = KP1,N R = 0. DO J = KP1,I LIJ = (I*(I-1))/2 + J ! LOC(I,J) LJK = (J*(J-1))/2 + K ! LOC(J,K) R = R + A(LIJ) * A(LJK) END DO ! LOOP ON J IF( I .NE. N ) THEN IP1 = I + 1 DO J = IP1,N L = (J*(J-1))/2 + I ! LOC(J,I) LJK = (J*(J-1))/2 + K ! LOC(J,K) R = R + A(L) * A(LJK) END DO ! LOOP ON J END IF ! (I .NE. N ) ! FINISH WITH R(K) G(I) = R L = (I*(I-1))/2 + K ! LOC(I,K) S = S + R * A(L) END DO ! LOOP ON I ! DONE WITH R'S NOW MAKE T'S DO I = KP1,N L = (K*(K+1))/2 ! LOC(K,K) R = A(L) ! R IS PI, S IS CORRECTED L = (I*(I-1))/2 + K ! LOC(I,K) G(I) = ( G(I) - S*A(L) / (2. * R ) ) / R END DO ! LOOP ON I ! FINISH WITH UPDATE OF REST OF A DO I = KP1,N LIK = (I*(I-1))/2 + K ! LOC(I,K) DO J = KP1,I LIJ = (I*(I-1))/2 + J ! LOC(I,J) LJK = (J*(J-1))/2 + K ! LOC(J,K) A(LIJ) = A(LIJ) - A(LIK)*G(J) - A(LJK)*G(I) END DO ! LOOP ON J END DO ! LOOP ON I ! END IF ! ( ETA .EQ. 0. ) END DO ! LOOP ON K ! END OF LOOP ON K END IF ! ( N .GT. 2 ) ! FINISH WITH LAST CASES (NO TFM) L = ((N-1)*N)/2 ! LOC(N-1,N-1) D(N-1) = A(L) L = (N*(N+1))/2 ! LOC(N,N) D(N) = A(L) L = L - 1 ! LOC(N,N-1) G(N-1) = A(L) ! RETURN END SUBROUTINE TRIDIG 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 ! IMPLICIT NONE REAL, INTENT(IN) :: ALPHA,BETA REAL, INTENT(OUT) :: GAMA,SIGMA,RNU REAL 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 qreig0.f95 *********************