C C ________________________________________________________ C | | C | COMPUTE EIGENVECTOR CORRESPONDING TO GIVEN REAL | C | EIGENVALUE FOR A REAL TRIDIAGONAL MATRIX | C | | C | INPUT: | C | | C | E --EIGENVALUE | C | | C | L --SUBDIAGONAL | C | | C | D --DIAGONAL | C | | C | U --SUPERDIAGONAL | C | | C | N --MATRIX DIMENSION | C | | C | W --WORK ARRAY (LENGTH AT LEAST 4N) | C | | C | OUTPUT: | C | | C | E --IMPROVED ESTIMATE FOR EIGENVALUE | C | | C | X --EIGENVECTOR | C | | C | BUILTIN FUNCTIONS: ABS,SQRT | C |________________________________________________________| C SUBROUTINE TVECT(E,X,L,D,U,N,W) REAL D(1),L(1),U(1),W(1),X(1),E,O,P,Q,R,S,T,V,Y,Z INTEGER F,G,H,I,J,K,M,N IF ( N .GT. 1 ) GOTO 10 E = D(1) X(1) = 1. RETURN 10 M = N - 1 J = 2 C --------------------------------------------------------- C |*** STORE MATRIX IN W ARRAY AND SUBTRACT EIGENVALUE ***| C --------------------------------------------------------- DO 20 I = 1,M X(I) = 0. W(J) = U(I) W(J+1) = D(I) - E W(J+2) = L(I) 20 J = J + 4 W(J) = 0. W(J+1) = D(N) - E O = 65536.**(-4) T = .5*O S = T 30 T = .5*T P = S S = S + T IF ( S .GE. O ) GOTO 40 IF ( S+T .GT. S ) GOTO 30 40 R = W(3) V = ABS(R) + ABS(W(4)) F = 4 M = 4*N - 3 G = -2 C --------------------------- C |*** FACTOR THE MATRIX ***| C --------------------------- 50 G = G + 4 IF ( G .GT. M ) GOTO 110 H = G - 1 I = G + 2 J = G + 5 C -------------------- C |*** FIND PIVOT ***| C -------------------- Q = W(I) Y = ABS(Q) Z = ABS(R) IF ( Z .GE. Y ) GOTO 80 C ------------------- C |*** SWAP ROWS ***| C ------------------- IF ( V .GE. Y ) GOTO 60 V = Y F = I 60 T = W(G) W(G) = W(J) W(J) = T T = R/Q K = G + 1 W(K) = Q K = J - 1 S = W(K) IF ( S .EQ. 0. ) GOTO 70 IF ( S .EQ. O ) S = P W(K) = -S*T W(H) = S GOTO 100 70 W(K) = S W(H) = O GOTO 100 80 W(H) = 0. IF ( V .GE. Z ) GOTO 90 V = Z F = I 90 IF ( R .EQ. 0. ) GOTO 120 T = Q/R C ------------------- C |*** ELIMINATE ***| C ------------------- 100 R = W(J) - T*W(G) W(J) = R W(I) = T GOTO 50 110 IF ( ABS(R) .GE. V ) GOTO 120 V = R F = J + 1 C --------------------------------------------------- C |*** COMPUTE INITIAL EIGENVECTOR APPROXIMATION ***| C --------------------------------------------------- 120 J = F/4 X(J) = 1. IF ( J .EQ. 1 ) GOTO 140 K = F - 5 J = J - 1 X(J) = (X(J)-W(K-1)*X(J+1))/W(K) 130 IF ( J .EQ. 1 ) GOTO 140 J = J - 1 K = K - 4 T = W(K-2) IF ( T .EQ. O ) T = 0. X(J) = (X(J)-W(K-1)*X(J+1)-T*X(J+2))/W(K) GOTO 130 140 IF ( V .EQ. 0. ) GOTO 230 S = 0. DO 150 I = 1,N T = ABS(X(I)) 150 IF ( T .GT. S ) S = T R = 0. S = 1./S DO 160 I = 1,N T = S*X(I) R = R + T*T 160 X(I) = T K = 0 J = 1 Y = X(1) C ----------------------------------------------------- C |*** APPLY ONE ITERATION OF INVERSE POWER METHOD ***| C ----------------------------------------------------- 170 K = K + 4 I = J J = J + 1 S = W(K) W(K) = Y Y = X(J) IF ( W(K-3) .EQ. 0. ) GOTO 180 T = X(J) X(J) = X(I) X(I) = T 180 X(J) = X(J) - S*X(I) IF ( J .LT. N ) GOTO 170 C --------------------------- C |*** BACK SUBSTITUTION ***| C --------------------------- S = X(J)/W(K+3) X(J) = S T = ABS(S) V = S*Y J = J - 1 K = K - 1 S = (X(J)-W(K-1)*S)/W(K) X(J) = S IF ( ABS(S) .GT. T ) T = ABS(S) V = V + S*W(K+1) 190 IF ( J .EQ. 1 ) GOTO 200 J = J - 1 K = K - 4 Z = W(K-2) IF ( Z .EQ. O ) Z = 0. S = (X(J)-W(K-1)*S-Z*X(J+2))/W(K) X(J) = S IF ( ABS(S) .GT. T ) T = ABS(S) V = V + S*W(K+1) GOTO 190 200 IF ( V .NE. 0. ) V = R/V S = 0. T = 1./T Z = 0. DO 210 I = 1,N S = S + (X(I)*V)**2 210 Z = Z + (T*X(I))**2 T = T/SQRT(Z) DO 220 I = 1,N 220 X(I) = T*X(I) C -------------------------------------------------------------- C |*** USE RAYLEIGH QUOTIENT TO IMPROVE EIGENVALUE ESTIMATE ***| C -------------------------------------------------------------- IF ( R+R .GE. S ) E = E + V RETURN 230 T = 0. DO 240 I = 1,N S = ABS(X(I)) 240 IF ( S .GT. T ) T = S T = 1./T Z = 0. DO 250 I = 1,N 250 Z = Z + (T*X(I))**2 T = T/SQRT(Z) DO 260 I = 1,N 260 X(I) = T*X(I) RETURN END