! Last change: DOS 27 Jul 2000 8:50 pm ! *** copyright 2000 *** ! *** filename csvd0.f95 *** John F. Monahan ** ! ********************** PROGRAM PCSVD0 ! TEST PROGRAM FOR CSVD0 -- SINGULAR VALUE DECOMPOSITION ! COMPUTES SINGULAR VALUES AND RIGHT AND LEFT SINGULAR VECTORS ! IMPLICIT NONE COMPLEX, DIMENSION(10,10) :: A COMPLEX, DIMENSION(10) :: D,F REAL, DIMENSION(10) :: PIU,PIV INTEGER I,J,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 CSVD0(D,F,N,MIT) IMPLICIT NONE INTEGER, INTENT(IN) :: N,MIT COMPLEX, DIMENSION(:), INTENT(IN OUT) :: D,F END SUBROUTINE CSVD0 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 ') ! THREE EXAMPLES IN 'csvdex.dat' OPEN( UNIT=5, FILE='csvdex.dat') OPEN( UNIT=6, FILE='csvd0.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) 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) ! NOW FIND SINGULAR VALUES CALL CSVD0(D,F,N,20) WRITE(6,27) WRITE(6,22) (D(I),I=1,N) WRITE(6,22) (F(I),I=1,N) ! END DO ! LOOP ON ICASE STOP END PROGRAM PCSVD0 SUBROUTINE CSVD0(D,F,N,MIT) ! COMPUTES SINGULAR VALUES * ONLY * 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 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) :: N,MIT COMPLEX, DIMENSION(:), INTENT(IN OUT) :: D,F COMPLEX A,B,X,Y,GAMA,SIGMA,CNU REAL C,DD,SHIFT REAL, PARAMETER :: UNITM = 1.E-6 ! MACHINE UNIT ! INTEGER 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 D(NCUR) = CABS( D(NCUR) ) ! 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 D(NCUR) = CABS( D(NCUR) ) 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 D(L) = CABS( D(L) ) ! 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 D(L) = CABS( D(L) ) 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 ! 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 ! 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) ! 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 CSVD0 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 ! *** end of filename csvd0.f95 *********************