! Last change: DOS 27 Jul 2000 9:05 pm ! *** copyright 2000 *** ! *** filename rsvd1.f95 *** John F. Monahan ** ! ********************** PROGRAM PRSVD1 ! TEST PROGRAM FOR RSVD1 -- SINGULAR VALUE DECOMPOSITION ! COMPUTES SINGULAR VALUES AND RIGHT AND LEFT SINGULAR VECTORS ! IMPLICIT NONE REAL, DIMENSION(10,10) :: A,ACOPY,U,V REAL, DIMENSION(10) :: D,F,PIU,PIV REAL S INTEGER I,J,K,L,M,N,ICASE ! INTERFACE BLOCK INTERFACE SUBROUTINE BIDIAG(A,M,N,D,F,PIU,PIV) REAL, DIMENSION(:,:), INTENT(IN OUT) :: A REAL, DIMENSION(:), INTENT(OUT) :: D,F,PIU REAL, DIMENSION(:), INTENT(OUT) :: PIV INTEGER, INTENT(IN) :: M,N END SUBROUTINE BIDIAG SUBROUTINE EXPNDB(A,M,N,PIU,PIV,U,V) REAL, DIMENSION(:,:), INTENT(IN) :: A REAL, DIMENSION(:), INTENT(IN) :: PIU REAL, DIMENSION(:), INTENT(IN) :: PIV REAL, DIMENSION(:,:), INTENT(OUT) :: U REAL, DIMENSION(:,:), INTENT(OUT) :: V INTEGER, INTENT(IN) :: M,N END SUBROUTINE EXPNDB SUBROUTINE RSVD1(D,F,M,N,MIT,U,V) INTEGER, INTENT(IN) :: M,N,MIT REAL, DIMENSION(:), INTENT(IN OUT) :: D,F REAL, DIMENSION(:,:) :: U REAL, DIMENSION(:,:) :: V END SUBROUTINE RSVD1 END INTERFACE ! 20 FORMAT(2I4) 21 FORMAT(5F6.2) 22 FORMAT(1X,6F12.6) 25 FORMAT(//' ORIGINAL MATRIX, M=',I2,' N=',I2) 26 FORMAT(' BIDIAGONALIZED, DIAGONAL, THEN SUPERDIAGONAL') 27 FORMAT(' SINGULAR VALUES, THEN OFF-DIAGONAL ') 28 FORMAT(' IS V''A U DIAGONAL ?') 29 FORMAT(' IS V ORTHOGONAL, V''V = I ?') 30 FORMAT(' IS U ORTHOGONAL, U''U = I ?') 31 FORMAT(' ORTHOGONAL MATRIX V') 32 FORMAT(' ORTHOGONAL MATRIX U') ! THREE EXAMPLES IN 'svdex.dat' OPEN( UNIT=5, FILE='svdex.dat') OPEN( UNIT=6, FILE='rsvd1.out') ! DO ICASE = 1,3 ! READ IN MATRIX AND MAKE A COPY OF IT TO CHECK READ(5,20) M,N WRITE(6,25) M,N DO I = 1,M READ(5,21) (A(I,J),J=1,N) WRITE(6,21) (A(I,J),J=1,N) DO J = 1,N ACOPY(I,J) = A(I,J) END DO ! LOOP ON J END DO ! LOOP ON I ! FIRST BIDIAGONALIZE MATRIX CALL BIDIAG(A,M,N,D,F,PIU,PIV) ! WRITE(6,26) WRITE(6,22) (D(I),I=1,N) WRITE(6,22) (F(I),I=1,N) ! EXPAND U AND V STORED IN COMPRESSED FORM CALL EXPNDB(A,M,N,PIU,PIV,U,V) ! NOW HAVE U,V SO V'A U IS BIDIAGONAL ! NOW FIND SINGULAR VALUES AND VECTORS CALL RSVD1(D,F,M,N,20,U,V) WRITE(6,27) WRITE(6,22) (D(I),I=1,N) WRITE(6,22) (F(I),I=1,N) ! ! TEST BY COMPUTING V'A U WRITE(6,28) DO I = 1,M DO J = 1,N S = 0. DO K = 1,M DO L = 1,N S = S + V(K,I)*ACOPY(K,L)*U(L,J) END DO ! LOOP ON L END DO ! LOOP ON K A(I,J) = S END DO ! LOOP ON J WRITE(6,22) (A(I,J),J=1,N) END DO ! LOOP ON I ! NOW CHECK TO SEE IF U AND V ARE ORTHOGONAL WRITE(6,29) DO I = 1,M DO J = 1,M S = 0. DO K = 1,M S = S + V(I,K)*V(J,K) END DO ! LOOP ON K D(J) = S END DO ! LOOP ON J WRITE(6,22) (D(J),J=1,M) END DO ! LOOP ON I ! IS U ORTHOGONAL? WRITE(6,30) DO I = 1,N DO J = 1,N S = 0. DO K = 1,N S = S + U(I,K)*U(J,K) END DO ! LOOP ON K D(J) = S END DO ! LOOP ON J WRITE(6,22) (D(J),J=1,N) END DO ! LOOP ON I ! WRITE OUT THE TWO ORTHOGONAL MATRICES WRITE(6,31) DO I = 1,M WRITE(6,22) (V(I,J),J=1,M) END DO ! LOOP ON I WRITE(6,32) DO I = 1,N WRITE(6,22) (U(I,J),J=1,N) END DO ! LOOP ON I ! END DO ! LOOP ON ICASE STOP END PROGRAM PRSVD1 SUBROUTINE RSVD1(D,F,M,N,MIT,U,V) ! COMPUTES SINGULAR VALUES AND VECTORS OF BIDIAGONAL MATRIX OF SIZE N ! D(1) ... D(N) HOLDS DIAGONAL ELEMENTS ! F(1) ... F(N-1) HOLDS SUPERDIAGONAL ELEMENTS ! MIT IS MAXIMUM NUMBER OF ITERATIONS ! ! *** ON INPUT *** MATRICES U AND V MUST HOLD EITHER OF THESE ! 1) IDENTITY MATRIX -- OR ! 2) ORTHOGONAL MATRIX THAT FORMED BIDIAGONAL FROM BIDIAG AND EXPNDB ! ! *** ON OUTPUT *** MATRICES U AND V WILL HOLD EITHER ! 1) RIGHT AND LEFT SINGULAR VECTORS OF BIDIAGONAL MATRIX ! 2) RIGHT AND LEFT SINGULAR VECTORS OF ORIGINAL MATRIX (BEFORE BIDIAG) ! -- STORED AS COLUMNS ! ! ON OUTPUT, F(N) HOLDS THE NUMBER OF SINGULAR VALUES FOUND BEFORE MIT ! SINGULAR VALUES FOUND IN D(1) ... D(N) ON SUCCESSFUL EXIT ! IF UNSUCCESSFUL EXIT (F(N).NE.FLOAT(N)) THEN SINGULAR VALUES FOUND ! WILL BE AT END OF D D(N-FOUND+1) ... D(N) ! ! METHOD IS SIMPLE IMPLEMENTATION OF GOLUB-REINSCH SVD ALGORITHM ! ! J F MONAHAN (JULY 1986) ! REVISED NOVEMBER 1992, JUNE 1993 ! REVISED OCTOBER 1999 FOR FORTRAN 95, corrected March 2000 ! IMPLICIT NONE INTEGER, INTENT(IN) :: M,N,MIT REAL, DIMENSION(:), INTENT(IN OUT) :: D,F REAL, DIMENSION(:,:) :: U REAL, DIMENSION(:,:) :: V REAL A,B,C,DD,X,Y,GAMA,SIGMA,RNU,UI1,UI2,SHIFT REAL, PARAMETER :: UNITM = 1.E-6 ! MACHINE UNIT ! INTEGER I,KIT,NCUR,NCURM1,KT,KTP1,KTP2,L,LP1 ! ! COUNTS NUMBER OF ITERATIONS KIT = 0 ! COUNTS NUMBER OF SINGULAR VALUES FOUND F(N) = 0. ! COUNTS CURRENT SIZE OF PROBLEM LEFT TO BE DONE NCUR = N NCURM1 = NCUR - 1 ! COUNT BEGINING POINT - MAY CHANGE AS TOP ONES ! CONVERGE SLOWLY L = 1 LP1 = L + 1 ! DO WHILE( KIT .LT. MIT ) ! TEST FIRST BEFORE DOING ANY WORK IF(ABS(F(NCURM1)) .LE. UNITM*(ABS(D(NCUR))+ABS(D(NCURM1)))) THEN ! ! DEFLATE PROBLEM BY SHORTENING BY ONE F(N) = F(N) + 1. NCUR = NCURM1 IF(NCUR.EQ.1) THEN ! DEFLATED TO DONE F(N) = FLOAT(N) RETURN END IF NCURM1 = NCUR - 1 ELSE ! ! START THE ITERATION KIT = KIT + 1 ! FIRST COMPUTE THE SHIFT ! CREATE LAST 2 BY 2 SUBMATRIX OF A'A C = D(NCUR)*D(NCUR) + F(NCURM1)*F(NCURM1) B = D(NCURM1)*F(NCURM1) A = D(NCURM1)*D(NCURM1) IF( NCUR .GT. 2 ) A = A + F(NCUR-2)*F(NCUR-2) DD = (A - C) / 2. SHIFT = C + DD - SIGN( SQRT( DD*DD + B*B ), DD ) ! ! TEST BEGINNING BEFORE SHIFTING DO WHILE( ABS(F(L)).LE.UNITM*(ABS(D(L))+ABS(D(LP1))) ) ! INCREMENT L AS UPPER OFF DIAGONAL HAS CONVERGED L = LP1 LP1 = L + 1 F(N) = F(N) + 1 IF( L .GE. NCUR ) THEN ! DEFLATED TO DONE F(N) = FLOAT(N) RETURN END IF END DO ! WHILE LOOP ! ! HAVE SHIFT NOW DO THE FIRST TRANSFORMATION ! A = D(L)*D(L) - SHIFT B = D(L)*F(L) ! COMPUTE GIVENS ROTATION TO PRODUCE ZERO IN A'A CALL ROT734(A,B,GAMA,SIGMA,RNU) ! MULTIPLY ROTATION ON BIDIAGONAL MATRIX A = D(L)*GAMA + F(L)*SIGMA B = F(L)*GAMA - D(L)*SIGMA D(L) = A F(L) = B ! NEW ELEMENT INTRODUCED IS X IN (2,1) OR (L,LP1) POSITION X = D(LP1)*SIGMA D(LP1) = D(LP1)*GAMA ! APPLY ROTATION ON RIGHT VECTOR MATRIX U DO I = 1,N UI1 = U(I,L) UI2 = U(I,LP1) U(I,L) = UI1*GAMA + UI2*SIGMA U(I,LP1) = UI2*GAMA - UI1*SIGMA END DO ! LOOP ON I ! NOW DO SEQUENCE OF TRANSFORMATIONS TO CHASE NONZERO OUT DO KT = L,NCURM1 ! TRANSFORM ON LEFT TO ZERO OUT X IN (KT+1,KT) ELEMENT KTP1 = KT + 1 A = D(KT) CALL ROT734(A,X,GAMA,SIGMA,RNU) D(KT) = RNU A = GAMA*F(KT) + SIGMA*D(KTP1) B = GAMA*D(KTP1) - SIGMA*F(KT) F(KT) = A D(KTP1) = B ! APPLY ROTATION TO LEFT VECTOR V IN TRANSPOSED FORM DO I = 1,M UI1 = V(I,KT) UI2 = V(I,KTP1) V(I,KT) = UI1*GAMA + UI2*SIGMA V(I,KTP1) = UI2*GAMA - UI1*SIGMA END DO ! LOOP ON I ! IS THE CHASING OVER? IF( KT .LT. NCURM1 ) THEN Y = SIGMA*F(KTP1) F(KTP1) = GAMA*F(KTP1) ! TRANSFORM ON RIGHT TO ZERO OUT Y IN (KT,KT+2) ELEMENT A = F(KT) CALL ROT734(A,Y,GAMA,SIGMA,RNU) F(KT) = RNU A = GAMA*D(KTP1) + SIGMA*F(KTP1) B = GAMA*F(KTP1) - SIGMA*D(KTP1) D(KTP1) = A F(KTP1) = B KTP2 = KT + 2 X = SIGMA*D(KTP2) D(KTP2) = GAMA*D(KTP2) ! APPLY ROTATION TO RIGHT VECTORS IN U DO I = 1,N UI1 = U(I,KTP1) UI2 = U(I,KTP2) U(I,KTP1) = UI1*GAMA + UI2*SIGMA U(I,KTP2) = UI2*GAMA - UI1*SIGMA END DO ! LOOP ON I ! END IF ! ( KT .EQ. NCURM1 ) END DO ! LOOP ON KT ! END OF ITERATION ENDIF ! END DO ! WHILE ( KIT .LT. MIT ) ! FINISH HERE WHEN TOO MANY ITERATIONS RETURN END SUBROUTINE RSVD1 SUBROUTINE BIDIAG(A,M,N,D,F,PIU,PIV) ! BIDIAGONALIZES MATRIX A (M BY N) USING HOUSEHOLDER TRANSFORMATIONS ! U (N BY N) AND V (M BY M) ORTHOGONAL MATRICES SO THAT ! T ! V A U = B *** A IS M BY N WITH M .GE. N *** ! ! WHERE B IS BIDIAGONAL WITH DIAGONAL ELEMENTS STORED IN D AND ! SUPERDIAGONAL ELEMENTS STORED IN F ! ! MATRICES U AND V STORED IN COMPACT FORM OVERWRITING A ! PI VALUES OF HOUSEHOLDER TRANSFORMATIONS STORED IN PIU, PIV ! CALL EXPNDB TO EXPAND FROM COMPACT FORM TO FULL MATRICES U AND V ! ! ADAPTED FROM GOLUB - REINSCH ALGORITHM ! J F MONAHAN (JULY,1986) DEPT OF STATISTICS, N C S U, RALEIGH NC USA ! CORRECTED JULY, 1987 ! RECODED SEPTEMBER 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN OUT) :: A REAL, DIMENSION(:), INTENT(OUT) :: D,F,PIU REAL, DIMENSION(:), INTENT(OUT) :: PIV INTEGER, INTENT(IN) :: M,N REAL ETA,S INTEGER I,J,K,KP1 ! ! DO N STEPS ON LEFT, N-2 STEPS ON RIGHT ! *** BECAUSE M .GT. N *** DO K = 1,N KP1 = K + 1 ! DO V TRANSFORMATION, PUT ZEROS IN COLUMN KC ETA = 0 DO I = K,M ETA = MAX(ETA, ABS(A(I,K)) ) END DO ! LOOP ON I IF( (ETA .EQ. 0.) .OR. (K .EQ. M) ) THEN ! IF ALL ZERO THEN DO NO TRANSFORMATION PIV(K) = 0. D(K) = A(K,K) ELSE S = 0. DO I = K,M A(I,K) = A(I,K)/ETA S = S + A(I,K)*A(I,K) END DO ! LOOP ON I S = SIGN( SQRT(S), A(K,K) ) A(K,K) = A(K,K) + S D(K) = -S*ETA PIV(K) = S*A(K,K) ! PI VALUE FOR TFM STORED, VECTOR IN COLUMN K (K,M) ! IF( K .LT. N ) THEN ! APPLY TRANSFORM TO REMAINDER OF MATRIX DO J = KP1,N S = 0. DO I = K,M S = S + A(I,K)*A(I,J) END DO ! LOOP ON I S = S / PIV(K) DO I = K,M A(I,J) = A(I,J) - S*A(I,K) END DO ! LOOP ON I END DO ! LOOP ON J ! ! DONE WITH THIS TRANSFORMATION ON COLUMN K ! END IF ! ( KC .LE. N ) END IF ! ( (ETA .EQ. 0.) ! OR (K .EQ. M) ) IF( K .LE. N-2 ) THEN ! NOW TRANSFORM ON RIGHT AS U -- ZEROS IN ROWS ETA = 0. DO I = KP1,N ETA = MAX( ETA, ABS(A(K,I)) ) END DO ! LOOP ON I IF( ETA .EQ. 0. ) THEN ! IF RELEVANT PART OF ROW IS ALL ZERO, THEN NO TFM PIU(K) = 0. F(K) = 0. ELSE S = 0. DO I = KP1,N A(K,I) = A(K,I) / ETA S = S +A(K,I)*A(K,I) END DO ! LOOP ON I S = SIGN( SQRT(S), A(K,KP1) ) A(K,KP1) = A(K,KP1) + S PIU(K) = S*A(K,KP1) F(K) = - S * ETA ! NOW APPLY TRANSFORMATION TO REST OF MATRIX DO J = KP1,M S = 0. DO I = KP1,N S = S + A(K,I)*A(J,I) END DO ! LOOP ON I S = S / PIU(K) DO I = KP1,N A(J,I) = A(J,I) - S * A(K,I) END DO ! LOOP ON I END DO ! LOOP ON J ! END IF ! ( ETA .EQ. 0. ) END IF ! ( K .LE. N-2 ) END DO ! LOOP ON K F(N-1) = A(N-1,N) RETURN END SUBROUTINE BIDIAG SUBROUTINE EXPNDB(A,M,N,PIU,PIV,U,V) ! TAKES MATRIX A HOLDING HOUSEHOLDER TRANSFORMATIONS FROM BIDIAG ! IN COMPACT FORM, WITH CONSTANT PI IN PIU AND PIV ! AND CREATES MATRICES U AND V ! IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN) :: A REAL, DIMENSION(:), INTENT(IN) :: PIU REAL, DIMENSION(:), INTENT(IN) :: PIV REAL, DIMENSION(:,:), INTENT(OUT) :: U REAL, DIMENSION(:,:), INTENT(OUT) :: V INTEGER, INTENT(IN) :: M,N REAL S INTEGER I,J,K,KP1,KVMAX,KUMAX ! CREATE IDENTITY IN V DO I = 1,M DO J = 1,M V(I,J) = 0. END DO ! LOOP ON J V(I,I) = 1. END DO ! LOOP ON I ! CREATE IDENTITY IN U DO I = 1,N DO J = 1,N U(I,J) = 0. END DO ! LOOP ON J U(I,I) = 1. END DO ! LOOP ON I ! HOW MANY TRANSFORMATIONS ? KVMAX = MIN(M-1,N) KUMAX = MIN(M,N-2) ! DO K = 1,KVMAX ! IF NO TRANSFORMATION SKIP IT ALL IF( PIV(K) .NE. 0. ) THEN DO I = 1,M S = 0. DO J = K,M S = S + V(I,J)*A(J,K) END DO ! LOOP ON J S = S / PIV(K) DO J = K,M V(I,J) = V(I,J) - S * A(J,K) END DO ! LOOP ON J END DO ! LOOP ON I END IF ! ( PIV(K) .NE. 0. ) END DO ! LOOP ON K -- TRANSFORMATIONS ! NOW DO THE U TRANSFORMATIONS IF( KUMAX .LE. 0 ) RETURN DO K = 1,KUMAX IF(PIU(K) .NE. 0. ) THEN KP1 = K + 1 DO I = 1,N S = 0. DO J = KP1,N S = S + U(I,J)*A(K,J) END DO ! LOOP ON J S = S / PIU(K) DO J = KP1,N U(I,J) = U(I,J) - S * A(K,J) END DO ! LOOP ON J END DO ! LOOP ON I END IF ! (PIU(K) .NE. 0.) END DO ! LOOP ON K -- TRANSFORMATIONS RETURN END SUBROUTINE EXPNDB 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 rsvd1.f95 *********************