C C ________________________________________________________ C | | C | COMPUTE EIGENVECTOR CORRESPONDING TO GIVEN REAL | C | EIGENVALUE FOR A REAL SYMMETRIC BAND MATRIX | C | | C | INPUT: | C | | C | E --EIGENVALUE | C | | C | A --ARRAY CONTAINING BANDS ON DIAGONAL AND | C | BELOW DIAGONAL | C | | C | LA --LEADING (ROW) DIMENSION OF ARRAY A | C | | C | N --MATRIX DIMENSION | C | | C | H --HALF BANDWIDTH | C | | C | W --WORK ARRAY (LENGTH AT LEAST (3H+2)N) | C | | C | OUTPUT: | C | | C | E --IMPROVED ESTIMATE FOR EIGENVALUE | C | | C | X --EIGENVECTOR | C | | C | BUILTIN FUNCTIONS: ABS,MAX0,MIN0,SQRT | C |________________________________________________________| C SUBROUTINE HVECT(E,X,A,LA,N,H,W) INTEGER B,C,D,F,G,H,I,J,K,L,LA,M,N,O,P,Q,Y,Z REAL A(LA,1),W(1),X(1),E,R,S,T,U,V IF ( H .GT. 0 ) GOTO 20 S = ABS(A(1,1)-E) DO 10 J = 1,N X(J) = 0. T = ABS(A(1,J)-E) IF ( T .GT. S ) GOTO 10 S = T K = J 10 CONTINUE X(K) = 1. E = A(1,K) RETURN C --------------------------------------------------------- C |*** STORE MATRIX IN W ARRAY AND SUBTRACT EIGENVALUE ***| C --------------------------------------------------------- 20 C = 1 + 3*H D = H + H P = H + 1 M = 1 - H DO 70 J = 1,N X(J) = 0. L = M + D M = M + P DO 30 I = M,L 30 W(I) = 0. M = L + P W(M) = A(1,J) - E F = M L = N - J IF ( L .GE. H ) GOTO 40 IF ( L .EQ. 0 ) GOTO 70 L = L + M GOTO 50 40 L = M + H 50 M = M + 1 K = 2 - M DO 60 I = M,L T = A(K+I,J) W(I) = T F = F + C 60 W(F) = T 70 CONTINUE R = 0. Q = P + P L = C + 1 DO 80 I = Q,L T = ABS(W(I)) IF ( T .GT. R ) R = T 80 CONTINUE I = -H Z = P + H C --------------------------- C |*** FACTOR THE MATRIX ***| C --------------------------- K = 0 90 K = K + 1 I = I + L IF ( K .EQ. N ) GOTO 160 M = I + 1 Q = I O = MIN0(H,N-K) P = I + O C --------------------------------------- C |*** FIND PIVOT AND START ROW SWAP ***| C --------------------------------------- S = ABS(W(I)) DO 100 J = M,P T = ABS(W(J)) IF ( S .GT. T ) GOTO 100 Q = J S = T 100 CONTINUE IF ( R .LT. S ) GOTO 110 R = S Y = K IF ( R .EQ. 0 ) GOTO 170 110 J = I - Z B = Q - I W(J) = K + B T = W(Q) W(Q) = W(I) W(I) = T C ----------------------------- C |*** COMPUTE MULTIPLIERS ***| C ----------------------------- DO 120 J = M,P 120 W(J) = W(J)/T F = I + C*MIN0(D,N-K) G = C - O 130 M = P + G P = M + B T = W(P) W(P) = W(M) W(M) = T P = M + O IF ( T .EQ. 0. ) GOTO 150 Q = I - M M = M + 1 C ------------------------------ C |*** ELIMINATE BY COLUMNS ***| C ------------------------------ DO 140 J = M,P 140 W(J) = W(J) - T*W(J+Q) 150 IF ( P .LT. F ) GOTO 130 GOTO 90 160 IF ( ABS(W(I)) .GE. R ) GOTO 170 R = W(I) Y = N C --------------------------------------------------- C |*** COMPUTE INITIAL EIGENVECTOR APPROXIMATION ***| C --------------------------------------------------- 170 K = Y J = C*K - H T = 1. GOTO 190 180 T = X(K)/W(J+K) 190 X(K) = T IF ( K .EQ. 1 ) GOTO 210 Q = MAX0(1,K-D) K = K - 1 DO 200 I = Q,K 200 X(I) = X(I) - T*W(I+J) J = J - C GOTO 180 210 IF ( R .EQ. 0. ) GOTO 360 S = 0. DO 220 I = 1,N T = ABS(X(I)) 220 IF ( T .GT. S ) S = T R = 0. S = 1./S DO 230 I = 1,N T = S*X(I) R = R + T*T 230 X(I) = T J = -H Q = 1 P = H + 1 DO 240 I = 1,P 240 W(I+1) = X(I) C ----------------------------------------------------- C |*** APPLY ONE ITERATION OF INVERSE POWER METHOD ***| C ----------------------------------------------------- 250 J = J + C K = Q IF ( K .EQ. N ) GOTO 270 Q = K + 1 I = W(J+K-Z) T = X(I) X(I) = X(K) X(K) = T O = MIN0(K+H,N) DO 260 I = Q,O 260 X(I) = X(I) - T*W(I+J) IF ( O .EQ. N ) GOTO 250 W(Q+J) = X(O+1) GOTO 250 C --------------------------- C |*** BACK SUBSTITUTION ***| C --------------------------- 270 U = 0. 280 T = X(K)/W(J+K) IF ( ABS(T) .GT. U ) U = ABS(T) X(K) = T IF ( K .EQ. 1 ) GOTO 300 Q = MAX0(1,K-D) K = K - 1 DO 290 I = Q,K 290 X(I) = X(I) - T*W(I+J) J = J - C GOTO 280 300 V = 0. DO 310 I = 1,P 310 V = V + X(I)*W(I+1) IF ( P .EQ. N ) GOTO 330 J = Z + 1 K = P + 1 DO 320 I = K,N V = V + X(I)*W(J) 320 J = J + L 330 IF ( V .NE. 0. ) V = R/V S = 0. T = 1./U U = 0. DO 340 I = 1,N S = S + (X(I)*V)**2 340 U = U + (T*X(I))**2 T = T/SQRT(U) DO 350 I = 1,N 350 X(I) = T*X(I) C -------------------------------------------------------------- C |*** USE RAYLEIGH QUOTIENT TO IMPROVE EIGENVALUE ESTIMATE ***| C -------------------------------------------------------------- IF ( R+R .GE. S ) E = E + V RETURN 360 U = 0. DO 370 I = 1,N T = ABS(X(I)) 370 IF ( T .GT. U ) U = T T = 1./U U = 0. DO 380 I = 1,N 380 U = U + (T*X(I))**2 T = T/SQRT(U) DO 390 I = 1,N 390 X(I) = T*X(I) RETURN END