C C ________________________________________________________ C | | C | COMPUTE EIGENVECTOR CORRESPONDING TO GIVEN REAL | C | EIGENVALUE FOR A REAL HESSENBERG MATRIX | C | | C | INPUT: | C | | C | E --ESTIMATE OF EIGENVALUE | C | | C | A --HESSENBERG MATRIX PACKED AT START OF | C | ARRAY WITH LENGTH AT LEAST N(N+7)/2 - 2| C | (ALGORITHM DESTROYS ELEMENTS) | C | | C | N --DIMENSION OF MATRIX STORED IN A | C | | C | OUTPUT: | C | | C | E --IMPROVED ESTIMATE FOR EIGENVALUE | C | | C | X --ARRAY CONTAINING EIGENVECTOR | C | | C | BUILTIN FUNCTIONS: ABS | C |________________________________________________________| C SUBROUTINE EVECT(E,X,A,N) REAL A(1),X(1),E,R,S,T,U,V INTEGER F,G,H,I,J,K,L,M,N,O,P,Q IF ( N .GT. 1 ) GOTO 10 E = A(1) X(1) = 1. RETURN 10 G = (N*(N+1))/2 F = G H = G + N - 2 M = 1 L = 2 J = 0 K = 1 C ------------------------------------------------------ C |*** MAKE ROOM FOR PIVOTS AND SUBTRACT EIGENVALUE ***| C ------------------------------------------------------ 20 IF ( K .GT. N ) GOTO 40 P = L - 1 A(P) = A(P) - E IF ( K .EQ. N ) L = P DO 30 I = M,L 30 A(I-J) = A(I) J = K K = K + 1 A(H+K) = A(L) M = L + 1 L = M + K GOTO 20 40 I = G + 1 J = H + 2 K = H + N 50 IF ( I .GT. K ) GOTO 60 A(I) = A(J) I = I + 2 J = J + 1 GOTO 50 60 K = N M = G L = G - N + 1 G = H + N - 1 R = ABS(A(M)) + ABS(A(G)) 70 IF ( K .EQ. 1 ) GOTO 140 J = K K = K - 1 S = A(M) T = A(G) IF ( ABS(S) .GE. ABS(T) ) GOTO 100 C -------------------------------------------- C |*** PIVOT AND PERFORM ELIMINATION STEP ***| C -------------------------------------------- A(M) = T IF ( R .LT. ABS(T) ) GOTO 80 R = ABS(T) P = J O = M 80 T = S/T A(G) = T A(G+1) = 1. G = G - 2 L = L - K M = M - J DO 90 I = L,M Q = I + K S = A(I) A(I) = A(Q) - T*S 90 A(Q) = S GOTO 70 C -------------------------- C |*** ELIMINATION STEP ***| C -------------------------- 100 IF ( R .LT. ABS(T) ) GOTO 110 R = ABS(T) P = J O = M 110 IF ( S .EQ. 0. ) GOTO 130 T = T/S L = L - K M = M - J DO 120 I = L,M 120 A(I) = A(I) - T*A(I+K) 130 A(G) = T A(G+1) = 0. G = G - 2 GOTO 70 140 T = ABS(A(1)) IF ( T .GT. R ) GOTO 150 R = T P = K C ------------------------------------------------------ C |*** COMPUTE INITIAL APPROXIMATION TO EIGENVECTOR ***| C ------------------------------------------------------ 150 DO 160 I = 1,N 160 X(I) = 0. T = 1. L = F J = O - P K = P GOTO 180 170 T = X(K)/A(J+K) 180 X(K) = T IF ( K .EQ. 1 ) GOTO 200 K = K - 1 DO 190 I = 1,K 190 X(I) = X(I) - T*A(I+J) J = J - K GOTO 170 200 IF ( K .EQ. N ) GOTO 210 J = K K = K + 1 L = L + 2 X(K) = X(K) - A(L-1)*X(J) IF ( A(L) .EQ. 0. ) GOTO 200 T = X(K) X(K) = X(J) X(J) = T GOTO 200 210 IF ( R .EQ. 0. ) GOTO 310 C ----------------------------------------------------- C |*** APPLY ONE ITERATION OF INVERSE POWER METHOD ***| C ----------------------------------------------------- S = 0. DO 220 I = 1,N T = ABS(X(I)) 220 IF ( T .GT. S ) S = T P = H + N R = 0. S = 1./S DO 230 I = 1,N T = S*X(I) R = R + T*T X(I) = T 230 A(I+P) = T V = 0. L = F J = L - N K = N 240 T = X(K)/A(J+K) X(K) = T IF ( K .EQ. 1 ) GOTO 260 K = K - 1 DO 250 I = 1,K 250 X(I) = X(I) - T*A(I+J) J = J - K GOTO 240 260 U = ABS(X(1)) 270 IF ( K .EQ. N ) GOTO 290 J = K K = K + 1 L = L + 2 T = X(J) S = X(K) - A(L-1)*T X(K) = S IF ( ABS(S) .GT. U ) U = ABS(S) IF ( A(L) .EQ. 1. ) GOTO 280 V = V + A(P+J)*T GOTO 270 280 X(J) = S X(K) = T V = V + A(P+J)*X(J) GOTO 270 290 V = V + A(P+N)*X(N) IF ( V .NE. 0. ) V = R/V S = 0. T = 1./U DO 300 I = 1,N S = S + (V*X(I))**2 300 X(I) = T*X(I) C -------------------------------------------------------------- C |*** USE RAYLEIGH QUOTIENT TO IMPROVE EIGENVALUE ESTIMATE ***| C -------------------------------------------------------------- IF ( R+R .GE. S ) E = E + V RETURN 310 U = 0. DO 320 I = 1,N T = ABS(X(I)) 320 IF ( T .GT. U ) U = T T = 1./U DO 330 I = 1,N 330 X(I) = T*X(I) RETURN END