C C ________________________________________________________ C | | C | COMPUTE SINGULAR VALUE DECOMPOSITION OF A BIDIAGONAL | C | MATRIX: A = Q TIMES DIAGONAL MATRIX TIMES P TRANSPOSE | C | | C | INPUT: | C | | C | D --DIAGONAL | C | | C | N --DIMENSION OF BIDIAGONAL MATRIX | C | | C | U --OFF DIAGONAL | C | | C | IU --= 0 IF U IS SUPERDIAGONAL AND = 1 IF | C | U IS SUBDIAGONAL | C | | C | Q --EITHER A SEGMENT OF AN IDENTITY MATRIX | C | OR A SEGMENT OF THE MATRIX USED TO | C | BIDIAGONALIZE THE COEFFICIENT MATRIX | C | | C | LQ --LEADING (ROW) DIMENSION OF ARRAY Q | C | | C | MQ --NUMBER OF MATRIX ELEMENTS CONTAINED IN | C | EACH COLUMN OF Q | C | | C | IQ --AN INTEGER WHICH INDICATES WHETHER | C | COLUMNS OF Q TO BE PROCESSED (= 0 MEANS| C | NO AND = 1 MEANS YES) | C | | C | LP --LEADING (ROW) DIMENSION OF ARRAY P | C | | C | MP --NUMBER OF MATRIX ELEMENTS CONTAINED IN | C | EACH COLUMN OF P | C | | C | IP --AN INTEGER DEFINED LIKE IQ | C | | C | E --WORK ARRAY WITH AT LEAST N ELEMENTS | C | | C | F --WORK ARRAY WITH AT LEAST N ELEMENTS | C | | C | OUTPUT: | C | | C | Q --Q FACTOR IN THE SINGULAR VALUE DECOMP. | C | | C | D --SINGULAR VALUES IN DECREASING ORDER | C | | C | P --P FACTOR IN THE SINGULAR VALUE DECOMP. | C | | C | BUILTIN FUNCTIONS: ABS,AMAX1,SIGN,SQRT | C | PACKAGE SUBROUTINES: EIG3,FGV,HSR1-HSR5,SCL,SFT, | C | SNG0,SNG1,SORT2 | C |________________________________________________________| C SUBROUTINE SINGB(D,N,U,IU,Q,LQ,MQ,IQ,P,LP,MP,IP,E,F) INTEGER G,H,I,ID,IP,IQ,IU,J,JL,JR,J0,K,K2,L,LL,LP,LQ,L1 INTEGER M,MP,MQ,N,NS REAL D(1),E(1),F(1),Q(LQ,1),P(LP,1),U(1) REAL B,C,R,S,T,T0,T1,T2,T3,V,W,X,Y,Z IF ( N .GT. 1 ) GOTO 10 IF ( IQ .EQ. 1 ) Q(1,1) = 1. IF ( IP .EQ. 1 ) P(1,1) = 1. IF ( D(1) .GE. 0. ) RETURN D(1) = -D(1) IF ( IQ .EQ. 1 ) Q(1,1) = -1. RETURN 10 JL = IQ IF ( JL .EQ. 0 ) JL = 3 IF ( JL .NE. 3 ) JL = 1 JR = IP IF ( JR .EQ. 0 ) JR = 3 IF ( JR .NE. 3 ) JR = 2 IF ( IU .EQ. 0 ) GOTO 20 I = JL JL = JR JR = I 20 J = 0 L = N - 1 K2 = N - 2 DO 30 I = 1,L E(I) = 1. F(I) = 1. IF ( U(I) .EQ. 0. ) J = I 30 IF ( D(I) .EQ. 0. ) J = I E(N) = 1. F(N) = 1. B = 65536.**(-3) T = 1. 40 T = .5*T S = 1. + T IF ( S .GT. 1. ) GOTO 40 T0 = 1./(T+T) T2 = (T+T)**2 NS = 50*N L1 = 0 K = N LL = 0 GOTO 70 50 J = 0 DO 60 I = 1,L IF ( U(I) .EQ. 0. ) J = I 60 IF ( D(I) .EQ. 0. ) J = I 70 IF ( J .EQ. 0 ) GOTO 140 IF ( U(J) .EQ. 0. ) GOTO 140 C ------------------------------- C |*** ZERO DIAGONAL ELEMENT ***| C ------------------------------- I = J V = U(J) U(J) = 0. S = -V*ABS(E(J)) C --------------------------- C |*** PROCESS LEFT SIDE ***| C --------------------------- 80 H = I I = I + 1 R = ABS(E(I))*D(I) CALL FGV(X,Y,T,R,S,E(J),E(I)) IF ( T .EQ. 1. ) GOTO 90 D(I) = D(I) - Y*V CALL SNG0(Q,LQ,MQ,P,LP,MP,JL,J,I,X,Y) IF ( I .EQ. K ) GOTO 100 V = X*U(I) S = -V*ABS(E(J)) GOTO 80 90 D(I) = D(I)*Y - V CALL SNG1(Q,LQ,MQ,P,LP,MP,JL,J,I,X,Y) IF ( I .EQ. K ) GOTO 100 V = U(I) U(I) = V*Y S = -V*ABS(E(J)) GOTO 80 100 IF ( J .EQ. 1 ) GOTO 130 I = J - 1 S = U(I) U(I) = 0. C ---------------------------- C |*** PROCESS RIGHT SIDE ***| C ---------------------------- 110 H = I I = I - 1 R = D(H) CALL FGV(X,Y,T,R,S,F(H),F(J)) IF ( T .EQ. 1. ) GOTO 120 D(H) = R + X*S CALL SNG0(Q,LQ,MQ,P,LP,MP,JR,H,J,X,Y) IF ( H .EQ. 1 ) GOTO 130 S = -Y*U(I) GOTO 110 120 D(H) = X*R + S CALL SNG1(Q,LQ,MQ,P,LP,MP,JR,H,J,X,Y) IF ( H .EQ. 1 ) GOTO 130 S = -U(I) U(I) = X*U(I) GOTO 110 130 CALL SCL(D,U,N,Q,LQ,MQ,P,LP,MP,E,F,B,1,K,JL,JR) 140 J = J + 1 IF ( J .EQ. K ) GOTO 320 S = 0. T = 0. C ----------------------------- C |*** SET ERROR TOLERANCE ***| C ----------------------------- DO 150 I = J,L X = ABS((D(I)*E(I))*(D(I)*F(I))) Y = ABS((U(I)*E(I))*(U(I)*F(I+1))) S = AMAX1(S,X,Y) 150 T = T + X + Y X = ABS((E(K)*D(K))*(F(K)*D(K))) S = AMAX1(S,X) T3 = S*T2 T = T + X IF ( T .EQ. 0. ) T = 1. T1 = T0/T GOTO 280 160 LL = LL + 1 IF ( LL .GT. NS ) GOTO 530 IF ( L .GT. J ) GOTO 170 S = 0. T = 0. GOTO 180 170 S = U(K2) T = E(K2) C ----------------------- C |*** COMPUTE SHIFT ***| C ----------------------- 180 CALL SFT(C,D(K),D(L),U(L),S,E(K),E(L),T,F(K),F(L)) V = ABS(E(J)) W = ABS(F(J)) T = D(J)*W R = T*(D(J)*V) - C S = T*(U(J)*V) ID = 1 J0 = J I = J H = J - 1 190 G = H H = I I = I + 1 Z = ABS(F(I)) C ---------------------------- C |*** PROCESS RIGHT SIDE ***| C ---------------------------- CALL FGV(X,Y,T,R,S,F(H),F(I)) V = D(H) IF ( T .EQ. 1 ) GOTO 210 IF ( H .EQ. J ) GOTO 200 T = U(G) + X*S U(G) = T IF ( ABS((T*E(G))*(T*F(H))) .GT. T3 ) GOTO 200 J = H ID = 1 200 R = V + X*U(H) S = X*D(I) U(H) = U(H) - Y*V CALL SNG0(Q,LQ,MQ,P,LP,MP,JR,H,I,X,Y) GOTO 230 210 IF ( H .EQ. J ) GOTO 220 T = X*U(G) + S U(G) = T IF ( ABS((T*E(G))*(T*F(H))) .GT. T3 ) GOTO 220 J = H ID = 1 220 R = X*V + U(H) S = D(I) U(H) = Y*U(H) - V D(I) = Y*S CALL SNG1(Q,LQ,MQ,P,LP,MP,JR,H,I,X,Y) C --------------------------- C |*** PROCESS LEFT SIDE ***| C --------------------------- 230 CALL FGV(X,Y,T,R,S,E(H),E(I)) IF ( T .EQ. 1. ) GOTO 250 T = R + X*S D(H) = T IF ( ABS((T*E(H))*(T*F(H))) .GT. T3 ) GOTO 240 ID = 0 J = H 240 R = U(H) + X*D(I) D(I) = D(I) - Y*U(H) U(H) = R CALL SNG0(Q,LQ,MQ,P,LP,MP,JL,H,I,X,Y) IF ( I .EQ. K ) GOTO 270 S = X*U(I) GOTO 190 250 T = S + X*R D(H) = T IF ( ABS((T*E(H))*(T*F(H))) .GT. T3 ) GOTO 260 ID = 0 J = H 260 R = D(I) + X*U(H) D(I) = Y*D(I) - U(H) U(H) = R CALL SNG1(Q,LQ,MQ,P,LP,MP,JL,H,I,X,Y) IF ( I .EQ. K ) GOTO 270 S = U(I) U(I) = S*Y GOTO 190 270 CALL SCL(D,U,N,Q,LQ,MQ,P,LP,MP,E,F,B,J0,K,JL,JR) IF ( ID .EQ. 0 ) GOTO 70 280 W = E(L) X = E(K) Y = F(L) Z = F(K) R = ABS((X*D(K))*(Z*D(K))) S = ABS((W*D(L))*(Y*D(L))) T = ABS((X*U(L))*(Y*U(L))) C ------------------------------ C |*** TEST FOR CONVERGENCE ***| C ------------------------------ IF ( (S*T1)*(T*T1) .GT. 1. ) GOTO 160 L1 = L1 + 1 IF ( L1 .GT. 40 ) GOTO 290 IF ( T .EQ. 0. ) GOTO 290 R = R + T IF ( (S/R)*(T/R) .GT. T2 ) GOTO 160 290 L1 = 0 IF ( S .GT. R ) GOTO 310 R = -D(K)*ABS(X) S = U(L)*ABS(W) CALL FGV(W,Y,T,R,S,E(L),E(K)) X = E(K) IF ( T .EQ. 1. ) GOTO 300 D(K) = D(K) - Y*U(L) CALL SNG0(Q,LQ,MQ,P,LP,MP,JL,L,K,W,Y) GOTO 310 300 D(L) = D(L)*W D(K) = Y*D(K) - U(L) CALL SNG1(Q,LQ,MQ,P,LP,MP,JL,L,K,W,Y) 310 T = SIGN(SQRT(ABS(X)),X)*D(K)*SIGN(SQRT(ABS(Z)),Z) IF ( T .LT. 0. ) E(K) = -E(K) D(K) = ABS(T) K = L L = K2 K2 = K2 - 1 IF ( K .GT. J ) GOTO 280 320 X = E(K) Z = F(K) T = SIGN(SQRT(ABS(X)),X)*D(K)*SIGN(SQRT(ABS(Z)),Z) IF ( T .LT. 0. ) E(K) = -E(K) D(K) = ABS(T) K = L L = K2 K2 = K2 - 1 IF ( K .GT. 1 ) GOTO 50 IF ( K .EQ. 0 ) GOTO 330 GOTO 320 C ------------------------- C |*** RESCALE P AND Q ***| C ------------------------- 330 GOTO (340,360,380),JL 340 DO 350 J = 1,N T = E(J) T = SIGN(SQRT(ABS(T)),T) DO 350 I = 1,MQ 350 Q(I,J) = Q(I,J)*T GOTO 380 360 DO 370 J = 1,N T = E(J) T = SIGN(SQRT(ABS(T)),T) DO 370 I = 1,MP 370 P(I,J) = P(I,J)*T 380 GOTO (390,410,430),JR 390 DO 400 J = 1,N T = F(J) T = SIGN(SQRT(ABS(T)),T) DO 400 I = 1,MQ 400 Q(I,J) = Q(I,J)*T GOTO 430 410 DO 420 J = 1,N T = F(J) T = SIGN(SQRT(ABS(T)),T) DO 420 I = 1,MP 420 P(I,J) = P(I,J)*T C -------------------------------- C |*** REORDER THE EIGENPAIRS ***| C -------------------------------- 430 CALL SORT2(D,E,F,N) DO 440 I = 1,N J = E(I) 440 F(J) = I M = N + 1 DO 520 J = 1,N L = M - J K = E(L) I = F(J) E(I) = K F(K) = I T = D(J) D(J) = D(K) D(K) = T IF ( JL .EQ. 1 ) GOTO 450 IF ( JR .NE. 1 ) GOTO 480 450 S = 0. DO 460 I = 1,MQ T = Q(I,K) Q(I,K) = Q(I,J) Q(I,J) = T 460 S = S + T*T S = 1./SQRT(S) DO 470 I = 1,MQ 470 Q(I,J) = S*Q(I,J) 480 IF ( JR .EQ. 2 ) GOTO 490 IF ( JL .NE. 2 ) GOTO 520 490 S = 0. DO 500 I = 1,MP T = P(I,K) P(I,K) = P(I,J) P(I,J) = T 500 S = S + T*T S = 1./SQRT(S) DO 510 I = 1,MP 510 P(I,J) = S*P(I,J) 520 CONTINUE E(1) = N RETURN 530 K = N - K + 1 WRITE(6,*) 'SINCE THE STOPPING CRITERION NOT SATISFIED' WRITE(6,*) 'AFTER',NS,'ITERATIONS, WE STOP WHILE COMPUTING' WRITE(6,*) 'EIGENVALUE NUMBER',K E(1) = K RETURN END C% SUBROUTINE FGV(X,Y,S,P,Q,A,B) REAL A,B,C,P,Q,R,S,T,X,Y IF ( ABS(P) .GT. ABS(Q) ) GOTO 10 IF ( Q .EQ. 0. ) GOTO 110 R = A/B S = P/Q T = ABS(R)*S*S IF ( T .LT. 1 ) GOTO 70 T = T/(1.+T) R = B/A S = Q/P GOTO 20 10 R = B/A S = Q/P T = ABS(R)*S*S IF ( T .GT. 1. ) GOTO 60 T = 1./(1.+T) 20 A = SIGN(A*T,P) IF ( SIGN(R,P) .EQ. R ) GOTO 40 B = -ABS(B*T) GOTO 50 40 B = ABS(B*T) 50 Y = S X = ABS(R)*S S = 0. RETURN 60 T = T/(1.+T) R = A/B S = P/Q GOTO 80 70 T = 1./(1.+T) 80 C = A A = SIGN(B*T,Q) IF ( SIGN(R,Q) .EQ. R ) GOTO 90 B = -ABS(C*T) GOTO 100 90 B = ABS(C*T) 100 Y = S X = ABS(R)*S S = 1. RETURN 110 X = 0. Y = 0. S = 0. RETURN END C% SUBROUTINE SNG0(Q,LQ,M,P,LP,N,L,J,K,X,Y) INTEGER I,J,K,L,LP,LQ,M,N REAL Q(LQ,1),P(LP,1),S,T,X,Y GOTO (10,30,50),L 10 DO 20 I = 1,M T = Q(I,J) S = Q(I,K) Q(I,J) = T + X*S 20 Q(I,K) = S - Y*T RETURN 30 DO 40 I = 1,N T = P(I,J) S = P(I,K) P(I,J) = T + X*S 40 P(I,K) = S - Y*T 50 RETURN END C% SUBROUTINE SNG1(Q,LQ,M,P,LP,N,L,J,K,X,Y) INTEGER I,J,K,L,LP,LQ,M,N REAL Q(LQ,1),P(LP,1),S,T,X,Y GOTO (10,30,50),L 10 DO 20 I = 1,M T = Q(I,J) S = Q(I,K) Q(I,J) = X*T + S 20 Q(I,K) = Y*S - T RETURN 30 DO 40 I = 1,N T = P(I,J) S = P(I,K) P(I,J) = X*T + S 40 P(I,K) = Y*S - T 50 RETURN END C% SUBROUTINE SFT(S,A,B,C,D,E2,E1,E0,F2,F1) REAL A,B,C,D,E0,E1,E2,F1,F2,G0,G1,G2,H1,H2,S,W,X,Y,Z G0 = ABS(E0) G1 = ABS(E1) G2 = ABS(E2) H1 = ABS(F1) H2 = ABS(F2) W = (A*G2)*(A*H2) + (C*G1)*(C*H2) X = (B*G1)*(B*H1) + (D*G0)*(D*H1) Y = (B*G1)*(C*H1) Z = (B*G1)*(C*H2) CALL EIG3(S,S,X,W,Y,Z) RETURN END C% SUBROUTINE SCL(D,U,N,Q,LQ,MQ,P,LP,MP,E,F,B,J,K,JL,JR) INTEGER I,J,JL,JR,K,L,LP,LQ,M,MP,MQ,N REAL D(1),E(1),F(1),P(LP,1),Q(LQ,1),U(1),B,R,S,T,V INTEGER H T = 1. DO 10 I = J,K IF ( ABS(E(I)) .LT. T ) T = ABS(E(I)) 10 IF ( ABS(F(I)) .LT. T ) T = ABS(F(I)) IF ( T .GT. B ) RETURN C ---------------------------- C |*** RESCALE THE MATRIX ***| C ---------------------------- R = E(J) R = SIGN(SQRT(ABS(R)),R) E(J) = 1. S = F(J) S = SIGN(SQRT(ABS(S)),S) F(J) = 1. D(J) = R*D(J)*S T = R IF ( JL .EQ. 1 ) GOTO 20 IF ( JR .NE. 1 ) GOTO 40 T = S 20 DO 30 I = 1,MQ 30 Q(I,J) = Q(I,J)*T 40 T = S IF ( JR .EQ. 2 ) GOTO 50 IF ( JL .NE. 2 ) GOTO 70 T = R 50 DO 60 I = 1,MP 60 P(I,J) = P(I,J)*T 70 L = J + 1 H = J DO 130 M = L,K T = E(M) T = SIGN(SQRT(ABS(T)),T) E(M) = 1. S = F(M) S = SIGN(SQRT(ABS(S)),S) F(M) = 1. D(M) = S*D(M)*T U(H) = R*U(H)*S H = M R = T V = T IF ( JL .EQ. 1 ) GOTO 80 IF ( JR .NE. 1 ) GOTO 100 V = S 80 DO 90 I = 1,MQ 90 Q(I,M) = Q(I,M)*V 100 V = S IF ( JR .EQ. 2 ) GOTO 110 IF ( JL .NE. 2 ) GOTO 130 V = T 110 DO 120 I = 1,MP 120 P(I,M) = P(I,M)*V 130 CONTINUE RETURN END