C C ________________________________________________________ C | | C |COMPUTE SINGULAR VALUE DECOMPOSITION OF A GENERAL MATRIX| C | A = Q TIMES DIAGONAL MATRIX TIMES P TRANSPOSE | C | | C | INPUT: | C | | C | S --ARRAY WITH AT LEAST N ELEMENTS | C | | C | LQ --LEADING (ROW) DIMENSION OF ARRAY Q | C | | C | IQ --AN INTEGER WHICH INDICATES WHICH COL- | C | UMNS OF Q TO COMPUTE (= 0 MEANS NONE, | C | = 1 MEANS FIRST L, = 2 MEANS LAST M-L, | C | = 3 MEANS ALL M WHERE L = MIN(M,N)) | C | | C | LP --LEADING (ROW) DIMENSION OF ARRAY P | C | | C | IP --AN INTEGER (LIKE IQ) WHICH INDICATES | C | WHICH COLUMNS OF P TO COMPUTE | C | | C | A --ARRAY CONTAINING COEFFICIENT MATRIX | C | | C | LA --LEADING (ROW) DIMENSION OF ARRAY A | C | | C | M --ROW DIMENSION OF MATRIX STORED IN A | C | | C | N --COLUMN DIMENSION OF MATRIX STORED IN A | C | | C | W --WORK ARRAY(LENGTH AT LEAST MAX(M,3L-1))| C | | C | OUTPUT: | C | | C | Q --Q FACTOR IN THE SINGULAR VALUE DECOMP. | C | | C | S --SINGULAR VALUES IN DECREASING ORDER | C | | C | P --P FACTOR IN THE SINGULAR VALUE DECOMP. | C | | C | A --THE HOUSEHOLDER VECTORS USED IN THE | C | REDUCTION PROCESS | C | | C | NOTE: | C | | C | EITHER P OR Q CAN BE IDENTIFIED WITH A BUT NOT | C | BOTH. WHEN EITHER P OR Q ARE IDENTIFIED WITH A,| C | THE HOUSEHOLDER VECTORS IN A ARE DESTROYED. IF | C | IQ = 2, Q MUST HAVE M COLUMNS EVEN THOUGH THE | C | OUTPUT ARRAY HAS JUST M-L COLUMNS. SIMILARLY IF| C | IP = 2, P MUST HAVE N COLUMNS EVEN THOUGH THE | C | OUTPUT ARRAY HAS JUST N-L COLUMNS. | C | | C | BUILTIN FUNCTIONS: MIN0 | C | PACKAGE SUBROUTINES: BIDAG2,EIG3,FGV,HSR1-HSR5,SCL, | C | SFT,SINGB,SNG0,SNG1,SORT2 | C |________________________________________________________| C SUBROUTINE SING(Q,LQ,IQ,S,P,LP,IP,A,LA,M,N,W) INTEGER I,IP,IQ,IU,JL,JR,L,LA,LP,LQ,M,N REAL A(LA,1),S(1),Q(LQ,1),P(LP,1),W(1) IF ( IQ .GE. 0 ) GOTO 20 10 WRITE(6,*) 'ERROR: INPUT PARAMETER IQ FOR SUBROUTINE SING' WRITE(6,*) 'EITHER LESS THAN 0 OR GREATER THAN 3' STOP 20 IF ( IQ .GT. 3 ) GOTO 10 JL = 0 IF ( IQ .EQ. 0 ) GOTO 30 IF ( IQ .EQ. 2 ) GOTO 30 JL = 1 30 IF ( IP .GE. 0 ) GOTO 50 40 WRITE(6,*) 'ERROR: INPUT PARAMETER IP FOR SUBROUTINE SING' WRITE(6,*) 'EITHER LESS THAN 0 OR GREATER THAN 3' STOP 50 IF ( IP .GT. 3 ) GOTO 40 JR = 0 IF ( IP .EQ. 0 ) GOTO 60 IF ( IP .EQ. 2 ) GOTO 60 JR = 1 60 CALL BIDAG2(S,W,Q,LQ,IQ,P,LP,IP,A,LA,M,N) L = MIN0(M,N) IF ( L .GT. 1 ) GOTO 80 IF ( S(1) .GE. 0. ) RETURN S(1) = -S(1) IF ( JL .EQ. 0. ) RETURN DO 70 I = 1,M 70 Q(I,1) = -Q(I,1) RETURN 80 IU = 0 IF ( M .GE. N ) GOTO 90 IU = 1 90 CALL SINGB(S,L,W,IU,Q,LQ,M,JL,P,LP,N,JR,W(L),W(L+L)) RETURN END C C ________________________________________________________ C | | C | REDUCE A GENERAL MATRIX A TO BIDIAGONAL FORM | C | A = Q TIMES BIDIAGONAL MATRIX TIMES P TRANSPOSE | C | | C | INPUT: | C | | C | D --ARRAY WITH AT LEAST N ELEMENTS | C | | C | B --ARRAY WITH AT LEAST M ELEMENTS | C | | C | LQ --LEADING (ROW) DIMENSION OF ARRAY Q | C | | C | IQ --AN INTEGER WHICH INDICATES WHICH COL- | C | UMNS OF Q TO COMPUTE (= 0 MEANS NONE, | C | = 1 MEANS FIRST L, = 2 MEANS LAST M-L, | C | = 3 MEANS ALL M WHERE L = MIN(M,N)) | C | | C | LP --LEADING (ROW) DIMENSION OF ARRAY P | C | | C | IP --AN INTEGER (LIKE IQ) WHICH INDICATES | C | WHICH COLUMNS OF P TO COMPUTE | C | | C | A --ARRAY CONTAINING COEFFICIENT MATRIX | C | | C | LA --LEADING (ROW) DIMENSION OF ARRAY A | C | | C | M --ROW DIMENSION OF MATRIX STORED IN A | C | | C | N --COLUMN DIMENSION OF MATRIX STORED IN A | C | | C | OUTPUT: | C | | C | D --DIAGONAL OF BIDIAGONAL FORM | C | | C | B --SUPERDIAGONAL (IF M .GE. N) OR SUBDIAG-| C | ONAL (IF M .LT. N) OF BIDIAGONAL FORM | C | | C | A --THE HOUSEHOLDER VECTORS USED IN THE | C | REDUCTION PROCESS | C | | C | Q --THE Q FACTOR IN THE BIDIAGONALIZATION | C | | C | P --THE P FACTOR IN THE BIDIAGONALIZATION | C | | C | NOTE: | C | | C | EITHER P OR Q CAN BE IDENTIFIED WITH A BUT NOT | C | BOTH. WHEN EITHER P OR Q ARE IDENTIFIED WITH A,| C | THEN THE HOUSEHOLDER VECTORS IN A ARE DESTROYED| C | | C | BUILTIN FUNCTIONS: ABS,MIN0,SQRT | C | PACKAGE SUBROUTINES: HSR1-HSR5 | C |________________________________________________________| C SUBROUTINE BIDAG2(D,B,Q,LQ,IQ,P,LP,IP,A,LA,M,N) INTEGER H,I,IP,IQ,J,JP,JQ,K,L,LA,LP,LQ,M,N REAL A(LA,1),B(1),D(1),Q(LQ,1),P(LP,1),R,S,T,U L = MIN0(M,N) IF ( IQ .GE. 0 ) GOTO 20 10 WRITE(6,*) 'ERROR: INPUT PARAMETER IQ FOR SUBROUTINE BIDAG2' WRITE(6,*) 'EITHER LESS THAN 0 OR GREATER THAN 3' STOP 20 IF ( IQ .GT. 3 ) GOTO 10 JQ = IQ IF ( IQ .LE. 1 ) GOTO 30 IF ( IQ .EQ. 3 ) GOTO 30 IF ( M .EQ. L ) JQ = 0 30 IF ( IP .GE. 0 ) GOTO 50 40 WRITE(6,*) 'ERROR: INPUT PARAMETER IP FOR SUBROUTINE BIDAG2' WRITE(6,*) 'EITHER LESS THAN 0 OR GREATER THAN 3' STOP 50 IF ( IP .GT. 3 ) GOTO 40 JP = IP IF ( IP .LE. 1 ) GOTO 60 IF ( IP .EQ. 3 ) GOTO 60 IF ( N .EQ. L ) JP = 0 60 K = 1 H = 2 IF ( M .LT. N ) GOTO 330 IF ( M .GT. 1 ) GOTO 70 D(1) = A(1,1) IF ( IQ .GT. 0 ) Q(1,1) = 1. IF ( IP .GT. 0 ) P(1,1) = 1. RETURN 70 J = K K = H DO 80 I = K,M 80 IF ( A(I,J) .NE. 0. ) GOTO 110 D(J) = A(J,J) A(J,J) = 0. DO 90 I = K,N 90 D(I) = A(J,I) IF ( JQ .EQ. 0 ) GOTO 200 DO 100 I = J,M 100 Q(I,J) = 0. GOTO 200 110 T = ABS(A(J,J)) IF ( T .NE. 0. ) U = 1./T R = 1. DO 130 L = I,M S = ABS(A(L,J)) IF ( S .LE. T ) GOTO 120 U = 1./S R = 1. + R*(T*U)**2 T = S GOTO 130 120 R = R + (S*U)**2 130 CONTINUE S = T*SQRT(R) R = A(J,J) U = 1./SQRT(S*(S+ABS(R))) IF ( R .LT. 0. ) S = -S D(J) = -S A(J,J) = U*(R+S) DO 140 I = K,M 140 A(I,J) = A(I,J)*U IF ( JQ .EQ. 0 ) GOTO 160 DO 150 I = J,M 150 Q(I,J) = A(I,J) 160 IF ( K .GT. N ) GOTO 620 DO 190 L = K,N T = 0. DO 170 I = J,M 170 T = T + A(I,J)*A(I,L) A(J,L) = A(J,L) - T*A(J,J) D(L) = A(J,L) DO 180 I = K,M 180 A(I,L) = A(I,L) - T*A(I,J) 190 CONTINUE 200 H = K + 1 IF ( K .LT. N ) GOTO 210 IF ( K .GT. N ) GOTO 620 IF ( M .EQ. N ) GOTO 610 B(J) = A(J,N) GOTO 70 210 DO 220 I = H,N 220 IF ( D(I) .NE. 0. ) GOTO 240 B(J) = D(K) A(J,K) = 0. IF ( IP .EQ. 0 ) GOTO 70 DO 230 I = K,N 230 P(I,J) = 0. GOTO 70 240 T = ABS(D(K)) IF ( T .NE. 0. ) U = 1./T R = 1. DO 260 L = I,N S = ABS(D(L)) IF ( S .LE. T ) GOTO 250 U = 1./S R = 1. + R*(T*U)**2 T = S GOTO 260 250 R = R + (S*U)**2 260 CONTINUE S = T*SQRT(R) R = D(K) U = 1./SQRT(S*(S+ABS(R))) IF ( R .LT. 0. ) S = -S D(K) = U*(R+S) DO 270 I = H,N 270 D(I) = D(I)*U IF ( IP .EQ. 0 ) GOTO 290 DO 280 I = K,N 280 P(I,J) = D(I) 290 B(J) = -S DO 300 I = K,M 300 B(I) = 0. DO 310 L = K,N T = D(L) A(J,L) = T DO 310 I = K,M 310 B(I) = B(I) + T*A(I,L) DO 320 L = K,N T = D(L) DO 320 I = K,M 320 A(I,L) = A(I,L) - T*B(I) GOTO 70 330 DO 340 I = K,N 340 D(I) = A(K,I) 350 J = K K = H DO 360 I = K,N 360 IF ( D(I) .NE. 0. ) GOTO 370 U = D(J) D(J) = 0. GOTO 440 370 T = ABS(D(J)) IF ( T .NE. 0. ) U = 1./T R = 1. DO 390 L = I,N S = ABS(D(L)) IF ( S .LE. T ) GOTO 380 U = 1./S R = 1. + R*(T*U)**2 T = S GOTO 390 380 R = R + (S*U)**2 390 CONTINUE S = T*SQRT(R) R = D(J) U = 1./SQRT(S*(S+ABS(R))) IF ( R .LT. 0. ) S = -S D(J) = U*(R+S) DO 400 I = K,N 400 D(I) = D(I)*U U = -S IF ( K .GT. M ) GOTO 470 DO 410 I = K,M 410 B(I) = 0. DO 420 L = J,N T = D(L) A(J,L) = T DO 420 I = K,M 420 B(I) = B(I) + T*A(I,L) DO 430 L = J,N T = D(L) DO 430 I = K,M 430 A(I,L) = A(I,L) - T*B(I) 440 H = K + 1 IF ( IP .EQ. 0 ) GOTO 460 DO 450 I = J,N 450 P(I,J) = D(I) 460 D(J) = U IF ( K .LT. M ) GOTO 490 IF ( K .GT. M ) GOTO 620 B(J) = A(M,J) GOTO 330 470 DO 480 I = J,N 480 A(J,I) = D(I) GOTO 440 490 DO 500 I = H,M 500 IF ( A(I,J) .NE. 0. ) GOTO 520 B(J) = A(K,J) A(K,J) = 0. IF ( IQ .EQ. 0 ) GOTO 330 DO 510 I = K,M 510 Q(I,J) = 0. GOTO 330 520 T = ABS(A(K,J)) IF ( T .NE. 0. ) U = 1./T R = 1. DO 540 L = I,M S = ABS(A(L,J)) IF ( S .LE. T ) GOTO 530 U = 1./S R = 1. + R*(T*U)**2 T = S GOTO 540 530 R = R + (S*U)**2 540 CONTINUE S = T*SQRT(R) R = A(K,J) U = 1./SQRT(S*(S+ABS(R))) IF ( R .LT. 0. ) S = -S B(J) = -S A(K,J) = U*(R+S) DO 550 I = H,M 550 A(I,J) = A(I,J)*U IF ( IQ .EQ. 0 ) GOTO 570 DO 560 I = K,M 560 Q(I,J) = A(I,J) 570 DO 600 L = K,N T = 0. DO 580 I = K,M 580 T = T + A(I,J)*A(I,L) A(K,L) = A(K,L) - T*A(K,J) D(L) = A(K,L) DO 590 I = H,M 590 A(I,L) = A(I,L) - T*A(I,J) 600 CONTINUE GOTO 350 610 D(N) = A(N,N) B(N-1) = A(N-1,N) 620 IF ( JQ .EQ. 0 ) GOTO 650 IF ( N .GT. M ) GOTO 640 IF ( N .EQ. M ) GOTO 630 IF ( JQ .EQ. 1 ) CALL HSR3(Q,LQ,M,N) IF ( JQ .EQ. 2 ) CALL HSR4(Q,LQ,M,N) IF ( JQ .EQ. 3 ) CALL HSR5(Q,LQ,M,N) GOTO 650 630 CALL HSR2(Q,LQ,M) GOTO 650 640 CALL HSR1(Q,LQ,M) 650 IF ( JP .EQ. 0 ) RETURN IF ( N .LE. M ) GOTO 660 IF ( JP .EQ. 1 ) CALL HSR3(P,LP,N,M) IF ( JP .EQ. 2 ) CALL HSR4(P,LP,N,M) IF ( JP .EQ. 3 ) CALL HSR5(P,LP,N,M) RETURN 660 CALL HSR1(P,LP,N) RETURN END C SUBROUTINE HSR1(A,LA,N) INTEGER I,J,K,L,LA,M,N REAL A(LA,1),S A(1,1) = 1. IF ( N .EQ. 1 ) RETURN A(1,2) = 0. IF ( N .GT. 2 ) GOTO 10 A(2,1) = 0. A(2,2) = 1. RETURN 10 L = N - 2 M = N - 1 K = N S = A(K,L) A(N,N) = 1. - S*A(N,L) A(M,N) = -S*A(M,L) 20 J = M M = L L = L - 1 IF ( L .EQ. 0 ) GOTO 50 S = 0. DO 30 I = J,N 30 S = S + A(I,L)*A(I,K) A(M,K) = -S*A(M,L) DO 40 I = J,N 40 A(I,K) = A(I,K) - S*A(I,L) GOTO 20 50 A(1,K) = 0. K = K - 1 M = K L = K - 1 S = -A(M,L) DO 60 I = K,N 60 A(I,K) = S*A(I,L) A(K,K) = 1. + A(K,K) IF ( L .GT. 1 ) GOTO 20 DO 70 I = 2,N 70 A(I,1) = 0. RETURN END SUBROUTINE HSR2(A,LA,N) INTEGER I,J,K,L,LA,M,N REAL A(LA,1),S IF ( N .GT. 1 ) GOTO 10 A(1,1) = 1. RETURN 10 L = N - 2 M = N - 1 K = N S = A(N,M) A(N,N) = 1. - S*A(N,M) A(M,N) = -S*A(M,M) 20 J = M M = M - 1 IF ( M .EQ. 0 ) GOTO 50 S = 0. DO 30 I = J,N 30 S = S + A(I,M)*A(I,K) A(M,K) = -S*A(M,M) DO 40 I = J,N 40 A(I,K) = A(I,K) - S*A(I,M) GOTO 20 50 K = K - 1 M = K S = -A(K,K) DO 60 I = K,N 60 A(I,K) = S*A(I,K) A(K,K) = 1. + A(K,K) IF ( K .GT. 1 ) GOTO 20 RETURN END SUBROUTINE HSR3(A,LA,M,N) INTEGER I,J,K,L,LA,M,N REAL A(LA,1),S K = N IF ( M .GE. N ) GOTO 10 WRITE(6,*) 'ERROR: ARGUMENT M MUST BE .GE. N IN SUBROUTINE HSR3' STOP 10 S = -A(K,K) DO 20 I = K,M 20 A(I,K) = S*A(I,K) A(K,K) = 1. + A(K,K) IF ( K .EQ. 1 ) RETURN L = K 30 J = L L = L - 1 IF ( L .EQ. 0 ) GOTO 60 S = 0. DO 40 I = J,M 40 S = S + A(I,L)*A(I,K) A(L,K) = -S*A(L,L) DO 50 I = J,M 50 A(I,K) = A(I,K) - S*A(I,L) GOTO 30 60 K = K - 1 GOTO 10 END SUBROUTINE HSR4(A,LA,M,N) INTEGER I,J,K,L,LA,M,N REAL A(LA,1),S K = M IF ( M .GT. N ) GOTO 10 WRITE(6,*) 'ERROR: ARGUMENT M MUST BE .GT. N IN SUBROUTINE HSR4' STOP 10 S = -A(K,N) DO 20 I = N,M 20 A(I,K) = S*A(I,N) A(K,K) = 1. + A(K,K) L = N 30 J = L L = L - 1 IF ( L .EQ. 0 ) GOTO 60 S = 0. DO 40 I = J,M 40 S = S + A(I,L)*A(I,K) A(L,K) = -S*A(L,L) DO 50 I = J,M 50 A(I,K) = A(I,K) - S*A(I,L) GOTO 30 60 K = K - 1 IF ( K .GT. N ) GOTO 10 K = M - N DO 70 J = 1,K DO 70 I = 1,M 70 A(I,J) = A(I,J+N) RETURN END C ________________________________________________________ C | | C | PACKAGE SUBROUTINES: HSR3 | C |________________________________________________________| C SUBROUTINE HSR5(A,LA,M,N) INTEGER I,J,K,L,LA,M,N REAL A(LA,1),S IF ( M .GE. N ) GOTO 10 WRITE(6,*) 'ERROR: ARGUMENT M MUST BE .GE. N IN SUBROUTINE HSR5' STOP 10 IF ( M .EQ. N ) GOTO 90 K = M 20 S = -A(K,N) DO 30 I = N,M 30 A(I,K) = S*A(I,N) A(K,K) = 1. + A(K,K) L = N 40 J = L L = L - 1 IF ( L .EQ. 0 ) GOTO 70 S = 0. DO 50 I = J,M 50 S = S + A(I,L)*A(I,K) A(L,K) = -S*A(L,L) DO 60 I = J,M 60 A(I,K) = A(I,K) - S*A(I,L) GOTO 40 70 K = K - 1 IF ( K .GT. N ) GOTO 20 90 CALL HSR3(A,LA,M,N) RETURN END