C C ________________________________________________________ C | | C | COMPUTE THE DIAGONALIZATION OF A | C | SYMMETRIC TRIDIAGONAL MATRIX | C | | C | INPUT: | C | | C | V --EITHER THE IDENTITY MATRIX OR THE SIM- | C | ILARITY TRANSFORMATION USED TO TRIDIAG-| C | ONALIZE THE COEFFICIENT MATRIX | C | | C | LV --LEADING (ROW) DIMENSION OF ARRAY V | C | | C | D --DIAGONAL | C | | C | U --SUPERDIAGONAL | C | | C | N --MATRIX DIMENSION | C | | C | OUTPUT: | C | | C | E --ARRAY OF EIGENVALUES | C | | C | V --ARRAY OF EIGENVECTORS | C | | C | NOTE: THE ALGORITHM DESTROYS THE STARTING D AND U | C | | C | BUILTIN FUNCTIONS: ABS,AMAX1,SQRT | C | PACKAGE SUBROUTINES: EIG3,SORT2 | C |________________________________________________________| C SUBROUTINE TDG(E,V,LV,D,U,N) INTEGER G,H,I,J,K,L,LL,LV,L1,M,N,NS REAL D(1),E(1),U(1),V(LV,1),A,B,F,P,Q,R,S,T,T0,T1,T2,W,X,Y,Z E(1) = D(1) IF ( N .LE. 1 ) RETURN J = 1 H = 1 D(1) = 1. T = 0. DO 10 I = 2,N E(I) = D(I) D(I) = 1. T = T + ABS(E(H)) + ABS(U(H)) IF ( U(H) .EQ. 0. ) J = I 10 H = I B = 65536.**(-3) R = 1. 20 R = .5*R S = 1. + R IF ( S .GT. 1. ) GOTO 20 T0 = R + R T2 = T0*T0 NS = 50*N K = N LL = 0 L = N - 1 L1 = 0 IF ( K .EQ. J ) GOTO 250 GOTO 60 30 K = L L = L - 1 J = 0 DO 40 I = 1,L 40 IF ( U(I) .EQ. 0. ) J = I J = J + 1 IF ( K .EQ. J ) GOTO 250 T = 0. R = SQRT(D(J)) C ----------------------------- C |*** SET ERROR TOLERANCE ***| C ----------------------------- DO 50 I = J,L S = SQRT(D(I+1)) T = T + D(I)*ABS(E(I)) + R*ABS(U(I))*S 50 R = S 60 T = T + ABS(E(K)) IF ( T .EQ. 0. ) T = 1. T1 = 1./(T*T0) GOTO 230 70 LL = LL + 1 IF ( LL .GT. NS ) GOTO 330 T = 1. DO 80 I = J,K 80 IF ( D(I) .LT. T ) T = D(I) IF ( T .GT. B ) GOTO 120 C ---------------------------- C |*** RESCALE THE MATRIX ***| C ---------------------------- E(J) = E(J)*D(J) S = SQRT(D(J)) D(J) = 1. DO 90 I = 1,N 90 V(I,J) = V(I,J)*S IF ( J .GT. 1 ) U(J-1) = U(J-1)*S H = J G = J + 1 DO 110 M = G,K E(M) = E(M)*D(M) T = SQRT(D(M)) D(M) = 1. U(H) = S*U(H)*T S = T DO 100 I = 1,N 100 V(I,M) = V(I,M)*T 110 H = M 120 T = D(K) S = D(L) C ----------------------- C |*** COMPUTE SHIFT ***| C ----------------------- CALL EIG3(F,F,S*E(L),T*E(K),S*U(L),T*U(L)) H = J - 1 I = J P = D(J)*E(J) - F Q = D(J)*U(J) 130 G = H H = I I = I + 1 Y = D(H) Z = D(I) C ------------------------- C |*** GIVENS ROTATION ***| C ------------------------- IF ( ABS(P) .GT. ABS(Q) ) GOTO 140 R = Y/Z S = P/Q T = R*S*S IF ( T .LT. 1. ) GOTO 190 T = T/(1.+T) R = Z/Y S = Q/P GOTO 150 140 R = Z/Y S = Q/P T = R*S*S IF ( T .GT. 1. ) GOTO 180 T = 1./(1.+T) 150 D(H) = Y*T D(I) = Z*T Y = S X = R*S C ------------------------------ C |*** PROCESS EIGENVECTORS ***| C ------------------------------ DO 160 M = 1,N T = V(M,H) S = V(M,I) V(M,H) = T + X*S 160 V(M,I) = S - Y*T IF ( H .EQ. J ) GOTO 170 T = U(G) + X*A U(G) = T IF ( T .EQ. 0. ) J = H 170 T = E(H) P = U(H) Q = T + X*P R = P - Y*T W = E(I) - Y*P S = P + X*E(I) U(H) = S - Y*Q E(H) = Q + X*S E(I) = W - R*Y IF ( I .EQ. K ) GOTO 230 Q = U(I) A = X*Q P = Z*W - F Q = Q*Z GOTO 130 180 T = T/(1.+T) R = Y/Z S = P/Q GOTO 200 190 T = 1./(1.+T) 200 D(H) = Z*T D(I) = Y*T Y = S X = R*S C ------------------------------ C |*** PROCESS EIGENVECTORS ***| C ------------------------------ DO 210 M = 1,N T = V(M,H) S = V(M,I) V(M,H) = X*T + S 210 V(M,I) = Y*S - T IF ( H .EQ. J ) GOTO 220 T = X*U(G) + A U(G) = T IF ( T .EQ. 0. ) J = H 220 T = E(H) P = U(H) Q = X*T + P R = Y*P - T W = Y*E(I) - P S = X*P + E(I) U(H) = Y*S - Q E(H) = X*Q + S E(I) = Y*W - R IF ( I .EQ. K ) GOTO 230 A = U(I) Q = A U(I) = Y*A P = Z*W - Y*F Q = Q*Z GOTO 130 230 T = D(K) S = D(L) Q = E(K) P = U(L) C ------------------------------ C |*** TEST FOR CONVERGENCE ***| C ------------------------------ IF ( ((T*P)*T1)*((S*P)*T1) .GT. 1. ) GOTO 70 L1 = L1 + 1 IF ( L1 .GT. 30 ) GOTO 240 R = AMAX1(ABS(P),ABS(Q)) IF ( R .EQ. 0. ) GOTO 240 IF ( ABS((S*P)*(P/R)) .GT. T2*ABS((T*Q)*(Q/R)) ) GOTO 70 240 L1 = 0 E(K) = Q*T K = L L = L - 1 IF ( K .GT. J ) GOTO 230 250 E(J) = E(J)*D(J) IF ( J .GT. 2 ) GOTO 30 IF ( J .EQ. 1 ) GOTO 270 E(1) = E(1)*D(1) C ---------------------------------- C |*** RESCALE THE EIGENVECTORS ***| C ---------------------------------- 270 DO 280 J = 1,N T = SQRT(D(J)) DO 280 I = 1,N 280 V(I,J) = V(I,J)*T C -------------------------------- C |*** REORDER THE EIGENPAIRS ***| C -------------------------------- CALL SORT2(E,D,U,N) DO 290 I = 1,N J = D(I) 290 U(J) = I DO 320 J = 1,N K = D(J) I = U(J) D(I) = K U(K) = I T = E(J) E(J) = E(K) E(K) = T S = 0. DO 300 I = 1,N T = V(I,K) V(I,K) = V(I,J) V(I,J) = T 300 S = S + T*T S = 1./SQRT(S) DO 310 I = 1,N 310 V(I,J) = S*V(I,J) 320 CONTINUE D(1) = N RETURN 330 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 D(1) = K RETURN END