! Last change: DOS 27 Jul 2000 9:01 pm ! *** copyright 2000 *** ! *** filename rsvd0.f95 *** John F. Monahan ** ! ********************** PROGRAM PRSVD0 ! TEST PROGRAM FOR RSVD0 -- SINGULAR VALUE DECOMPOSITION ! COMPUTES SINGULAR VALUES * ONLY * ! IMPLICIT NONE REAL, DIMENSION(10,10) :: A REAL, DIMENSION(10) :: D,F,PIU,PIV INTEGER I,J,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 RSVD0(D,F,N,MIT) INTEGER, INTENT(IN) :: N,MIT REAL, DIMENSION(:), INTENT(IN OUT) :: D,F END SUBROUTINE RSVD0 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 ') ! THREE EXAMPLES IN 'svdex.dat' OPEN( UNIT=5, FILE='svdex.dat') OPEN( UNIT=6, FILE='rsvd0.out') ! DO ICASE = 1,3 ! READ IN MATRIX 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 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) ! NOW FIND SINGULAR VALUES CALL RSVD0(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 PRSVD0 SUBROUTINE RSVD0(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 ! ! J F MONAHAN (JULY 1986) ! REVISED NOVEMBER 1992, JUNE 1993 ! REVISED OCTOBER 1999 FOR FORTRAN 95, corrected March 2000 ! IMPLICIT NONE INTEGER, INTENT(IN) :: N,MIT REAL, DIMENSION(:), INTENT(IN OUT) :: D,F REAL A,B,C,DD,X,Y,GAMA,SIGMA,RNU,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. ! 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 ! 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 ! 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) ! 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 RSVD0 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 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 rsvd0.f95 *********************