! Last change: DOS 27 Jul 2000 8:59 pm ! *** copyright 2000 *** ! *** filename qreig1.f95 *** John F. Monahan ** ! ********************** PROGRAM PQREIG1 ! DEMONSTRATION PROGRAM QREIG1 -- 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,XLI REAL, DIMENSION(6,6) :: X,ACOPY,AX REAL S INTEGER L,I,J,K,N ! INTERFACE BLOCK INTERFACE SUBROUTINE EXPNDH(A,N,Q,ID) REAL, DIMENSION(:), INTENT(IN) :: A INTEGER, INTENT(IN) :: N REAL, DIMENSION(:,:), INTENT(IN OUT) :: Q INTEGER, INTENT(IN) :: ID END SUBROUTINE EXPNDH SUBROUTINE QREIG1(A,B,N,MIT,Q) REAL, DIMENSION(:), INTENT(IN OUT) :: A,B REAL, DIMENSION(:,:), INTENT(IN OUT) :: Q INTEGER, INTENT(IN) :: N,MIT END SUBROUTINE QREIG1 END INTERFACE ! 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='qreig1.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 ! MAKE A COPY OF THE MATRIX A DO J = 1,N L = (I*(I-1))/2 + J ! LOC(I,J) ACOPY(I,J) = A(L) ACOPY(J,I) = A(L) END DO ! LOOP ON J 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) ! EXPAND TO FULL FORM FOR NEXT STEP CALL EXPNDH(A,N,X,1) ! CALL EIGENVALUE AND VECTOR ROUTINE CALL QREIG1(D,G,N,15,X) WRITE(6,27) WRITE(6,21) (D(J),J=1,N) WRITE(6,21) (G(J),J=1,N) ! MULTIPLY A * X AND SEE IF IT MATCHES X * LAMDA ! WRITE(6,28) DO I = 1,N DO J = 1,N S = 0. DO K = 1,N S = S + ACOPY(I,K) * X(K,J) END DO ! LOOP ON K AX(I,J) = S XLI(J) = X(I,J) * D(J) END DO ! LOOP ON J WRITE(6,22) I WRITE(6,21) (AX(I,J),J=1,N) WRITE(6,21) (XLI(J),J=1,N) END DO ! LOOP ON I ! LASTLY, JUST PRINT OUT EIGENVECTORS WRITE(6,29) DO I = 1,N WRITE(6,21) (X(I,J),J=1,N) END DO ! LOOP ON I ! ! FINISH BY CHECKING X FOR ORTHOGONALITY WRITE(6,30) DO I = 1,N DO J = 1,N S = 0. DO K = 1,N S = S + X(K,I) * X(K,J) END DO ! LOOP ON K XLI(J) = S END DO ! LOOP ON J WRITE(6,21) (XLI(J),J=1,N) END DO ! LOOP ON I ! END OF TEST FOR THIS SIZE MATRIX END DO ! LOOP ON N STOP END PROGRAM PQREIG1 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, DIMENSION(:), INTENT(IN OUT) :: A,B REAL, DIMENSION(:,:), INTENT(IN OUT) :: Q INTEGER, INTENT(IN) :: N,MIT ! 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,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 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 EXPNDH(A,N,Q,ID) ! SUBROUTINE TO EXPAND A PRODUCT OF HOUSEHOLDER TFMS STORED ! IN COMPACT FORM IN A, MULTIPLYING IN BACK OF MATRIX Q ! ! A IN MATRIX WITH HOUSEHOLDER TFMS CREATED BY TRIDIG ! N IN DIMENSION OF MATRIX ! Q OUT ORTHOGONAL MATRIX THAT TRIDIAGONALIZED ORIGINAL A ! ID IN FLAG FOR WHETHER TO USE Q OR CREATE IDENTITY IN IT ! IF ID = 0 THEN COMPUTE Q U U U U ... U ! IF ID NE 0 THEN JUST PRODUCT OF U'S IN ASCENDING ORDER (Q=I) ! ! J F MONAHAN (JUNE, 1986) ! RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE REAL, DIMENSION(:), INTENT(IN) :: A INTEGER, INTENT(IN) :: N REAL, DIMENSION(:,:), INTENT(IN OUT) :: Q INTEGER, INTENT(IN) :: ID REAL R,PIK INTEGER NM2,I,J,K,KP1,L ! ! A IS SYMMETRIC, SO TO ACCESS A(I,J) USE LOCATION GIVEN BY ! LOC(I,J) = ( I * (I-1) ) /2 + J ! ! IF IDENT IS NOT ZERO, CREATE IDENTITY IF( ID .NE. 0 ) THEN ! DO J = 1,N DO I = 1,N Q(I,J) = 0. END DO ! LOOP ON I Q(J,J) = 1. END DO ! LOOP ON J END IF ! ( ID .NE. 0 ) ! ! NOW DO THE WORK (IF THERE IS ANY) IF( N .LE. 2 ) RETURN ! NM2 = N - 2 DO K = 1,NM2 L = (K*(K+1))/2 ! LOC(K,K) PIK = A(L) IF( PIK .NE. 0. ) THEN ! PIK .EQ. 0. MEANS NO TRANSFORMATION KP1 = K + 1 ! GO THROUGH THE ROWS OF Q (BY I) DO I = 1,N R = 0. DO J = KP1,N L = (J*(J-1))/2 + K ! LOC(J,K) R = R + Q(I,J) * A(L) END DO ! LOOP ON J R = R/PIK DO J = KP1,N L = (J*(J-1))/2 + K ! LOC(J,K) Q(I,J) = Q(I,J) - R * A(L) END DO ! LOOP ON J END DO ! LOOP ON I END IF ! ( PIK .NE. 0. ) ! END OF LOOP ON ROWS OF Q END DO ! LOOP ON K RETURN END SUBROUTINE EXPNDH 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 qreig1.f95 *********************