! Last change: DOS 27 Jul 2000 8:52 pm ! *** copyright 2000 *** ! *** filename csvd1.f95 *** John F. Monahan ** ! ********************** PROGRAM PCSVD1 ! TEST PROGRAM FOR CSVD1 -- SINGULAR VALUE DECOMPOSITION ! COMPUTES SINGULAR VALUES AND RIGHT AND LEFT SINGULAR VECTORS ! IMPLICIT NONE COMPLEX, DIMENSION(10,10) :: A,ACOPY,U,V COMPLEX, DIMENSION(10) :: D,F REAL, DIMENSION(10) :: PIU,PIV COMPLEX S INTEGER I,J,K,L,M,N,ICASE ! INTERFACE BLOCK INTERFACE SUBROUTINE BIDIAC(A,M,N,D,F,PIU,PIV) complex, DIMENSION(:,:), INTENT(IN OUT) :: A complex, DIMENSION(:), INTENT(OUT) :: D,F REAL, DIMENSION(:), INTENT(OUT) :: PIU REAL, DIMENSION(:), INTENT(OUT) :: PIV INTEGER, INTENT(IN) :: M,N END SUBROUTINE BIDIAC SUBROUTINE EXPNDC(A,M,N,PIU,PIV,U,V) complex, DIMENSION(:,:), INTENT(IN) :: A REAL, DIMENSION(:), INTENT(IN) :: PIU REAL, DIMENSION(:), INTENT(IN) :: PIV complex, DIMENSION(:,:), INTENT(OUT) :: U complex, DIMENSION(:,:), INTENT(OUT) :: V INTEGER, INTENT(IN) :: M,N END SUBROUTINE EXPNDC SUBROUTINE CSVD1(D,F,M,N,MIT,U,V) IMPLICIT NONE INTEGER, INTENT(IN) :: M,N,MIT COMPLEX, DIMENSION(:), INTENT(IN OUT) :: D,F COMPLEX, DIMENSION(:,:), INTENT(IN OUT) :: U COMPLEX, DIMENSION(:,:), INTENT(IN OUT) :: V END SUBROUTINE CSVD1 END INTERFACE ! 20 FORMAT(2I4) 21 FORMAT(10F6.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 VH A U DIAGONAL ?') 29 FORMAT(' IS V UNITARY, VH V = I ?') 30 FORMAT(' IS U UNITARY, UH U = I ?') 31 FORMAT(' UNITARY MATRIX V') 32 FORMAT(' UNITARY MATRIX U') ! THREE EXAMPLES IN 'csvdex.dat' OPEN( UNIT=5, FILE='csvdex.dat') OPEN( UNIT=6, FILE='csvd1.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 BIDIAC(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 EXPNDC(A,M,N,PIU,PIV,U,V) ! NOW HAVE U,V SO V'A U IS BIDIAGONAL ! NOW FIND SINGULAR VALUES AND VECTORS CALL CSVD1(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., 0. ) DO K = 1,M DO L = 1,N S = S + conjg( 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., 0. ) DO K = 1,M S = S + conjg( 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., 0. ) DO K = 1,N S = S + conjg( 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 PCSVD1 SUBROUTINE CSVD1(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) UNITARY 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 ! MODIFIED FOR COMPLEX MATRICES ! ! J F MONAHAN (JULY 1986 real version) ! REVISED NOVEMBER 1992, JUNE 1993 ! REVISED MARCH 2000 FOR FORTRAN 95 ! IMPLICIT NONE INTEGER, INTENT(IN) :: M,N,MIT COMPLEX, DIMENSION(:), INTENT(IN OUT) :: D,F COMPLEX, DIMENSION(:,:), INTENT(IN OUT) :: U COMPLEX, DIMENSION(:,:), INTENT(IN OUT) :: V COMPLEX A,B,X,Y,GAMA,SIGMA,CNU,UI1,UI2 REAL C,DD,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., 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(CABS(F(NCURM1)).LE.UNITM*(CABS(D(NCUR))+CABS(D(NCURM1)))) THEN ! ! ZERO OFF-DIAGONAL, MAKE DIAGONAL REAL DD = CABS( D(NCUR) ) B = conjg( D(NCUR) )/ DD D(NCUR) = DD DO I = 1,N U(I,NCUR) = U(I,NCUR)*B ! FIX U ACCORDINGLY END DO ! LOOP ON I ! DEFLATE PROBLEM BY SHORTENING BY ONE F(N) = F(N) + ( 1., 0. ) NCUR = NCURM1 IF(NCUR.EQ.1) THEN ! DEFLATED TO DONE ! ZERO OFF-DIAGONAL, MAKE DIAGONAL REAL DD = CABS( D(NCUR) ) B = conjg( D(NCUR) )/ DD D(NCUR) = DD DO I = 1,N U(I,NCUR) = U(I,NCUR)*B ! FIX U ACCORDINGLY END DO ! LOOP ON I 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)*conjg(D(NCUR)) + F(NCURM1)*conjg(F(NCURM1)) B = D(NCURM1)*conjg( F(NCURM1) ) A = D(NCURM1)*conjg( D(NCURM1) ) IF( NCUR .GT. 2 ) A = A + F(NCUR-2)*conjg( F(NCUR-2) ) DD = real( (A - C) / 2. ) SHIFT = C + DD - SIGN( SQRT( real( DD*DD + B*conjg(B) ) ), DD ) ! ! TEST BEGINNING BEFORE SHIFTING DO WHILE( CABS(F(L)).LE.UNITM*(CABS(D(L))+CABS(D(LP1))) ) ! ZERO OFF DIAGONAL, MAKE DIAGONAL REAL DD = CABS( D(L) ) B = conjg( D(L) )/ DD D(L) = DD DO I = 1,N U(I,L) = U(I,L)*B ! FIX U ACCORDINGLY END DO ! LOOP ON I ! INCREMENT L AS UPPER OFF DIAGONAL HAS CONVERGED L = LP1 LP1 = L + 1 F(N) = F(N) + ( 1., 0. ) IF( L .GE. NCUR ) THEN ! DEFLATED TO DONE ! ZERO OFF-DIAGONAL, MAKE DIAGONAL REAL DD = CABS( D(L) ) B = conjg( D(L) )/ DD D(L) = DD DO I = 1,N U(I,L) = U(I,L)*B ! FIX U ACCORDINGLY END DO ! LOOP ON I F(N) = FLOAT(N) RETURN END IF END DO ! WHILE LOOP ! ! HAVE SHIFT NOW DO THE FIRST TRANSFORMATION ! A = conjg( D(L) )*D(L) - SHIFT B = D(L)*conjg( F(L) ) ! COMPUTE GIVENS ROTATION TO PRODUCE ZERO IN A'A CALL RCT734(A,B,GAMA,SIGMA,CNU) ! MULTIPLY ROTATION ON BIDIAGONAL MATRIX A = D(L)*conjg( GAMA ) + F(L)*conjg( 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)*conjg( 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*conjg( GAMA ) + UI2*conjg( 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 RCT734(A,X,GAMA,SIGMA,CNU) D(KT) = CNU A = conjg( GAMA )*F(KT) + conjg( 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*conjg( GAMA ) - UI1*conjg( SIGMA ) END DO ! LOOP ON I ! IS THE CHASING OVER? IF( KT .LT. NCURM1 ) THEN Y = conjg( SIGMA )*F(KTP1) F(KTP1) = GAMA*F(KTP1) ! TRANSFORM ON RIGHT TO ZERO OUT Y IN (KT,KT+2) ELEMENT A = F(KT) CALL RCT734(A,Y,GAMA,SIGMA,CNU) F(KT) = CNU A = conjg( GAMA )*D(KTP1) + conjg( SIGMA )*F(KTP1) B = GAMA*F(KTP1) - SIGMA*D(KTP1) D(KTP1) = A F(KTP1) = B KTP2 = KT + 2 X = conjg( 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*conjg( GAMA ) + UI2*conjg( 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 CSVD1 SUBROUTINE RCT734(ALPHA,BETA,GAMA,SIGMA,CNU) ! GIVEN ALPHA AND BETA, PRODUCES GIVENS TRANSFORMATION VALUES ! ! ( conjg(GAMA) conjg(SIGMA) ) ( ALPHA ) ( CNU ) ! ( ) ( ) = ( ) ! ( -SIGMA GAMA ) ( BETA ) ( ZERO ) ! ! COMPLEX VERSION OF ALGORITHM 7.3.4 FROM STEWART ! IMPLICIT NONE COMPLEX, INTENT(IN) :: ALPHA,BETA COMPLEX, INTENT(OUT) :: GAMA,SIGMA,CNU REAL RNU,ETA,DELTA ! ETA = MAX( ABS(REAL(ALPHA)), ABS(AIMAG(ALPHA)), & & ABS(REAL(BETA)), ABS(AIMAG(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 CNU = ( 0., 0. ) GAMA = ( 1., 0. ) SIGMA = ( 0., 0. ) RETURN else DELTA = real( (ALPHA/ETA)*conjg(ALPHA/ETA) & & + (BETA/ETA)*conjg(BETA/ETA) ) DELTA = SQRT(DELTA) RNU = ETA * DELTA GAMA = ALPHA / RNU SIGMA = BETA / RNU CNU = RNU END if RETURN END subroutine rct734 SUBROUTINE BIDIAC(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) unitary MATRICES SO THAT ! H ! 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 EXPNDC 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 MARCH 2000 FOR FORTRAN 95 ! complex form of BIDAIG ! IMPLICIT NONE complex, DIMENSION(:,:), INTENT(IN OUT) :: A complex, DIMENSION(:), INTENT(OUT) :: D,F REAL, DIMENSION(:), INTENT(OUT) :: PIU REAL, DIMENSION(:), INTENT(OUT) :: PIV INTEGER, INTENT(IN) :: M,N REAL ETA,AS,C complex 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( real(A(I,K))), abs(aimag(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 AS = ( 0., 0. ) DO I = K,M A(I,K) = A(I,K)/ETA AS = AS + A(I,K)*conjg( A(I,K) ) END DO ! LOOP ON I AS = SQRT(AS) c = Cabs( A(K,K) ) s = AS if( c .GT. 0. ) s = A(K,K)*AS/C A(K,K) = A(K,K) + S D(K) = -S*ETA PIV(K) = AS * (AS + C ) ! 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 + conjg( 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(real(A(K,I))), abs(aimag(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 AS = ( 0., 0. ) DO I = KP1,N A(K,I) = A(K,I) / ETA AS = AS +A(K,I)*conjg( A(K,I) ) END DO ! LOOP ON I AS = SQRT(AS) C = cabs( A(K,KP1) ) S = AS if( c .gt. 0. ) S = A(K,KP1)*AS/C A(K,KP1) = A(K,KP1) + S PIU(K) = AS * (AS + C) F(K) = - S * ETA ! NOW APPLY TRANSFORMATION TO REST OF MATRIX DO J = KP1,M S = ( 0., 0. ) DO I = KP1,N S = S + conjg( 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 BIDIAC SUBROUTINE EXPNDC(A,M,N,PIU,PIV,U,V) ! TAKES MATRIX A HOLDING HOUSEHOLDER TRANSFORMATIONS FROM BIDIAC ! IN COMPACT FORM, WITH CONSTANT PI IN PIU AND PIV ! AND CREATES MATRICES U AND V ! IMPLICIT NONE complex, DIMENSION(:,:), INTENT(IN) :: A REAL, DIMENSION(:), INTENT(IN) :: PIU REAL, DIMENSION(:), INTENT(IN) :: PIV complex, DIMENSION(:,:), INTENT(OUT) :: U complex, DIMENSION(:,:), INTENT(OUT) :: V INTEGER, INTENT(IN) :: M,N complex S INTEGER I,J,K,KP1,KVMAX,KUMAX ! CREATE IDENTITY IN V DO I = 1,M DO J = 1,M V(I,J) = ( 0., 0. ) END DO ! LOOP ON J V(I,I) = ( 1., 0. ) END DO ! LOOP ON I ! CREATE IDENTITY IN U DO I = 1,N DO J = 1,N U(I,J) = ( 0., 0. ) END DO ! LOOP ON J U(I,I) = ( 1., 0. ) END DO ! LOOP ON I ! HOW MANY TRANSFORMATIONS ? KVMAX = MIN0(M-1,N) KUMAX = MIN0(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., 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 * conjg( 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., 0.) DO J = KP1,N S = S + U(I,J)*conjg( 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 EXPNDC ! *** end of filename csvd1.f95 *********************