! Last change: DOS 27 Jul 2000 9:08 pm ! *** copyright 2000 *** ! *** filename tridig.f95 *** John F. Monahan ** ! ********************** PROGRAM PTRIDIG ! TEST PROGRAM FOR TRIDIG WHICH TRIDIAGONALIZES A SYMMETRIC MATRIX ! AND EXPNDH WHICH EXPANDS HOUSEHOLDER TFMS TO FULL ORTHOGONAL 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 ! A STORED IN SYMMETRIC FORMAT REAL, DIMENSION(21) :: A ! BIGGEST IS 6, 6*7/2 = 21 REAL, DIMENSION(6) :: D,G REAL, DIMENSION(6,6) :: ACOPY,Q,QAQ INTEGER :: ID = 1 ! JUST DO ID NE 0 INTEGER N,I,J,K,L ! 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 END INTERFACE ! 20 FORMAT(6F5.1) 21 FORMAT(6F12.6) 25 FORMAT(///' THE ORIGINAL MATRIX, N =',I2) 26 FORMAT(/' MATRIX AFTER WITH PACKED TRANSFORMATIONS ') 27 FORMAT(/' DIAGONAL, THEN SUPER/SUB-DIAGONAL ') 28 FORMAT(/' ORTHOGONAL MATRIX THAT DOES TRIDIAGONALIZATION') 29 FORMAT(/' DOES IT TRIDIAGONALIZE? Q''A Q IS') 30 FORMAT(/' IS IT ORTHOGONAL? Q''Q IS') ! JUST DO IDENT NE 0 ! 4 EXAMPLES IN THIS FILE OPEN( UNIT=5, FILE='eigen.dat') OPEN( UNIT=6, FILE='tridig.out') ! FOUR MATRICES OF DIMENSION 3, 4, 5, 6 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) N K = 0 DO I = 1,N WRITE(6,21) (A((I*(I-1))/2 + J),J=1,I) ! MAKE COPY DO J = 1,I K = K + 1 ACOPY(I,J) = A(K) ACOPY(J,I) = A(K) END DO ! LOOP ON J END DO ! LOOP ON I ! PUT INTO TRIDIAGONAL FORM CALL TRIDIG(A,N,D,G) ! WRITE OUT MATRIX WITH TFMS STORED WRITE(6,26) DO I = 1,N WRITE(6,21) (A((I*(I-1))/2 + J),J=1,I) END DO ! LOOP ON I ! WRITE DIAGONAL AND SUPER/SUB WRITE(6,27) WRITE(6,21) (D(J),J=1,N) WRITE(6,21) (G(J),J=1,N) ! EXPAND INTO ORTHOGONAL MATRIX CALL EXPNDH(A,N,Q,ID) ! WRITE IT OUT WRITE(6,28) DO I = 1,N WRITE(6,21) (Q(I,J),J=1,N) END DO ! LOOP ON I ! CHECK IT OUT -- FIRST DOES IT TRIDIAGONALIZE ? WRITE(6,29) DO I = 1,N DO J = 1,N QAQ(I,J) = 0. DO K = 1,N DO L = 1,N QAQ(I,J) = QAQ(I,J) + Q(K,I)*ACOPY(K,L)*Q(L,J) END DO ! LOOP ON L END DO ! LOOP ON K END DO ! LOOP ON J WRITE(6,21) (QAQ(I,J),J=1,N) END DO ! LOOP ON I ! IS IT ORTHOGONAL ? WRITE(6,30) DO I = 1,N DO J = 1,N QAQ(I,J) = 0. DO K = 1,N QAQ(I,J) = QAQ(I,J) + Q(I,K)*Q(J,K) END DO ! LOOP ON K END DO ! LOOP ON J WRITE(6,21) (QAQ(I,J),J=1,N) END DO ! LOOP ON I ! DONE WITH THIS MATRIX END DO ! LOOP ON N STOP END PROGRAM PTRIDIG 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 ! *** end of filename tridig.f95 *********************