C C ________________________________________________________ C | | C | COMPUTE EIGENVECTOR CORRESPONDING TO GIVEN COMPLEX | C | EIGENVALUE FOR A COMPLEX HESSENBERG MATRIX | C | | C | INPUT: | C | | C | E --EIGENVALUE ESTIMATE (COMPLEX VARIABLE) | C | | C | A --HESSENBERG MATRIX PACKED AT START OF | C | COMPLEX ARRAY WITH AT LEAST N(N+7)/2 -2| C | ELEMENTS (ALGORITHM DESTROYS ELEMENTS) | C | | C | OUTPUT: | C | | C | E --IMPROVED ESTIMATE FOR EIGENVALUE | C | | C | X --ARRAY CONTAINING EIGENVECTOR | C | | C | BUILTIN FUNCTIONS: CABS,CONJG,REAL | C | PACKAGE FUNCTIONS: MAG,SQ | C |________________________________________________________| C SUBROUTINE CEVECT(E,X,A,N) COMPLEX A(1),X(1),E,Y,Z REAL R,S,T,MAG,SQ 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.,0.) 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 = CABS(A(M)) + CABS(A(G)) 70 IF ( K .EQ. 1 ) GOTO 140 J = K K = K - 1 Y = A(M) Z = A(G) S = CABS(Y) T = CABS(Z) IF ( S .GE. T ) GOTO 100 C -------------------------------------------- C |*** PIVOT AND PERFORM ELIMINATION STEP ***| C -------------------------------------------- A(M) = Z IF ( R .LT. T ) GOTO 80 R = T P = J O = M 80 Z = Y/Z A(G) = Z A(G+1) = (1.,0.) G = G - 2 L = L - K M = M - J DO 90 I = L,M Q = I + K Y = A(I) A(I) = A(Q) - Y*Z 90 A(Q) = Y GOTO 70 C -------------------------- C |*** ELIMINATION STEP ***| C -------------------------- 100 IF ( R .LT. T ) GOTO 110 R = T P = J O = M 110 IF ( S .EQ. 0. ) GOTO 130 Z = Z/Y L = L - K M = M - J DO 120 I = L,M 120 A(I) = A(I) - Z*A(I+K) 130 A(G) = Z A(G+1) = (0.,0.) G = G - 2 GOTO 70 140 T = CABS(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.,0.) Z = (1.,0.) L = F J = O - P K = P GOTO 180 170 Z = X(K)/A(J+K) 180 X(K) = Z IF ( K .EQ. 1 ) GOTO 200 K = K - 1 DO 190 I = 1,K 190 X(I) = X(I) - Z*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 ( REAL(A(L)) .EQ. 0. ) GOTO 200 Z = X(K) X(K) = X(J) X(J) = Z 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 = MAG(X(I)) 220 IF ( T .GT. S ) S = T P = H + N R = 0. S = 1./S DO 230 I = 1,N Z = S*X(I) R = R + SQ(Z) X(I) = Z 230 A(I+P) = Z Y = (0.,0.) L = F J = L - N K = N 240 Z = X(K)/A(J+K) X(K) = Z IF ( K .EQ. 1 ) GOTO 260 K = K - 1 DO 250 I = 1,K 250 X(I) = X(I) - Z*A(I+J) J = J - K GOTO 240 260 T = MAG(X(1)) 270 IF ( K .EQ. N ) GOTO 290 J = K K = K + 1 L = L + 2 Z = X(J) X(K) = X(K) - A(L-1)*Z S = MAG(X(K)) IF ( S .GT. T ) T = S IF ( REAL(A(L)) .EQ. 1. ) GOTO 280 Y = Y + CONJG(A(P+J))*Z GOTO 270 280 X(J) = X(K) X(K) = Z Y = Y + CONJG(A(P+J))*X(J) GOTO 270 290 Y = Y + CONJG(A(P+N))*X(N) IF ( MAG(Y) .NE. 0. ) Y = R/Y S = 0. T = 1./T DO 300 I = 1,N S = S + SQ(Y*X(I)) 300 X(I) = T*X(I) C -------------------------------------------------------------- C |*** USE RAYLEIGH QUOTIENT TO IMPROVE EIGENVALUE ESTIMATE ***| C -------------------------------------------------------------- IF ( R+R .GE. S ) E = E + Y RETURN 310 T = 0. DO 320 I = 1,N S = MAG(X(I)) 320 IF ( S .GT. T ) T = S T = 1./T DO 330 I = 1,N 330 X(I) = T*X(I) RETURN END FUNCTION SQ(A) COMPLEX A,B REAL C(2) EQUIVALENCE (B,C) B = A SQ = C(1)**2 + C(2)**2 RETURN END