C C ________________________________________________________ C | | C | COMPUTE ALL EIGENVALUES OF A COMPLEX HESSENBERG MATRIX | C | | C | INPUT: | C | | C | A --COEFFICIENTS OF HESSENBERG MATRIX | C | PACKED AT START OF COMPLEX ARRAY | C | | C | N --DIMENSION OF MATRIX STORED IN A | C | | C | W --COMPLEX WORK ARRAY WITH AT LEAST | C | N ELEMENTS | C | | C | D,P --REAL WORK ARRAYS WITH AT LEAST N | C | ELEMENTS IN EACH ARRAY | C | | C | OUTPUT: | C | | C | E --COMPLEX ARRAY OF EIGENVALUES | C | | C | BUILTIN FUNCTIONS: AMAX1,CONJG,SQRT | C | PACKAGE FUNCTIONS: MAG,SQR | C | PACKAGE SUBROUTINES: CEIG | C |________________________________________________________| C SUBROUTINE VLS(E,A,N,D,P,W) COMPLEX A(1),E(1),W(1),X,Y,Z,Z1,Z2,Z3,Z4 REAL D(1),P(1),B,Q,R,S,T,T0,T1,T2,U,MAG,SQR INTEGER F,F0,G,G0,H,I,J,K,KL,KM,KS,K0,L,LL,L0,L1,M,N,NS,O DO 10 I = 1,N 10 D(I) = 1. IF ( N .EQ. 1 ) GOTO 330 B = 65536.**(-3) C --------------------------------- C |*** COMPUTE MACHINE EPSILON ***| C --------------------------------- T = 1. 20 T = .5*T S = 1. + T IF ( S .GT. 1. ) GOTO 20 T0 = T + T T2 = T0*T0 NS = 50*N LL = 0 KL = (N*(N+3))/2 - 1 M = N KM = KL - M 30 I = KM J = M C -------------------------------------------- C |*** CHECK FOR ZERO SUBDIAGONAL ELEMENT ***| C -------------------------------------------- 40 IF ( MAG(A(I)) .EQ. 0. ) GOTO 50 I = I - J J = J - 1 IF ( J .GT. 1 ) GOTO 40 50 K0 = J L = J + 1 G = I + L F = G - 1 F0 = F G0 = G L0 = L IF ( J .EQ. M ) GOTO 320 DO 60 I = J,M 60 P(I) = SQRT(D(I)) S = 0. K = J C -------------------------------- C |*** COMPUTE NORM OF MATRIX ***| C -------------------------------- 70 H = F - J T = P(K) DO 80 I = F,G 80 S = S + T*MAG(A(I))*P(I-H) K = L L = L + 1 F = F + K G = G + L IF ( K .LT. M ) GOTO 70 IF ( K .GT. M ) GOTO 90 G = G - 1 GOTO 70 90 IF ( S .EQ. 0. ) S = 1. T1 = 1./(S*T0) L1 = 0 GOTO 300 100 LL = LL + 1 IF ( LL .GT. NS ) GOTO 340 T = 1. DO 110 I = K0,M 110 IF ( D(I) .LT. T ) T = D(I) IF ( T .GT. B ) GOTO 150 C ---------------------------- C |*** RESCALE THE MATRIX ***| C ---------------------------- DO 120 I = K0,M P(I) = SQRT(D(I)) 120 D(I) = 1. F = F0 G = G0 L = L0 K = K0 130 H = F - K0 T = P(K) DO 140 I = F,G 140 A(I) = A(I)*T*P(I-H) K = L L = L + 1 F = F + K G = G + L IF ( K .LT. M ) GOTO 130 IF ( K .GT. M ) GOTO 150 G = G - 1 GOTO 130 C ----------------------- C |*** COMPUTE SHIFT ***| C ----------------------- 150 S = D(M-1) T = D(M) Z1 = A(KM-1)*S Z2 = A(KM)*S Z3 = A(KL-1)*T Z4 = A(KL)*T CALL CEIG(Z,Z,Z1,Z2,Z3,Z4) 160 K = K0 L = L0 F = F0 G = G0 KS = K Z1 = D(K)*A(F) - Z Z2 = D(K)*A(G) C --------------------------------- C |*** COMPUTE GIVENS ROTATION ***| C --------------------------------- 170 T = MAG(Z1) + MAG(Z2) IF ( T .EQ. 0. ) GOTO 190 T = 1./T Q = SQR(Z1,T) R = SQR(Z2,T) Q = D(K)*Q R = D(L)*R IF ( Q .GT. R ) GOTO 180 Z4 = Z1/Z2 Z3 = (D(K)/D(L))*CONJG(Z4) E(K) = Z3 W(K) = Z4 P(K) = 1. S = R/(Q+R) R = D(K) D(K) = D(L)*S Q = D(L) D(L) = R*S GOTO 200 180 Z4 = Z2/Z1 Z3 = (D(L)/D(K))*CONJG(Z4) E(K) = Z3 W(K) = Z4 P(K) = 0. S = Q/(Q+R) D(K) = D(K)*S Q = D(L) D(L) = D(L)*S GOTO 200 190 P(K) = 0. Z3 = (0.,0.) Z4 = (0.,0.) E(K) = Z3 W(K) = Z4 Q = D(L) 200 IF ( K .GT. KS ) GOTO 220 C ------------------------------------------ C |*** PROCESS INITIAL COLUMN FROM LEFT ***| C ------------------------------------------ IF ( P(K) .EQ. 0. ) GOTO 210 Y = A(G) + Z3*A(F) A(G) = Z4*A(G) - A(F) A(F) = Y GOTO 240 210 Y = A(F) + Z3*A(G) A(G) = A(G) - Z4*A(F) A(F) = Y GOTO 240 C ----------------------------------- C |*** PROCESS COLUMNS FROM LEFT ***| C ----------------------------------- 220 I = G - 1 IF ( P(K) .EQ. 0. ) GOTO 230 X = A(H)*Z3 + X A(H) = X Y = A(G) + Z3*A(I) A(G) = Z4*A(G) - A(I) A(I) = Y IF ( MAG(X) .NE. 0. ) GOTO 240 K0 = K L0 = L F0 = I G0 = G GOTO 240 230 X = X*Z3 + A(H) A(H) = X Y = A(I) + Z3*A(G) A(G) = A(G) - Z4*A(I) A(I) = Y IF ( MAG(X) .NE. 0. ) GOTO 240 K0 = K L0 = L F0 = I G0 = G 240 J = K K = L L = L + 1 DO 260 I = KS,J O = G + I H = O + 1 IF ( P(I) .EQ. 0. ) GOTO 250 Y = A(H) + E(I)*A(O) A(H) = W(I)*A(H) - A(O) A(O) = Y GOTO 260 250 Y = A(O) + E(I)*A(H) A(H) = A(H) - W(I)*A(O) A(O) = Y 260 CONTINUE C ------------------------------------ C |*** PROCESS COLUMNS FROM RIGHT ***| C ------------------------------------ Z3 = CONJG(Z3) Z4 = CONJG(Z4) IF ( P(J) .EQ. 0. ) GOTO 280 Z1 = Q*A(G+K) - Z*W(J) DO 270 I = F,G H = I + K Y = A(H) + Z3*A(I) A(H) = Z4*A(H) - A(I) 270 A(I) = Y IF ( L .GT. M ) GOTO 300 F = F + K H = G G = G + L Z2 = Q*A(G) X = A(G) A(G) = Z4*A(G) GOTO 170 280 Z1 = Q*A(G+K) - Z DO 290 I = F,G H = I + K Y = A(I) + Z3*A(H) A(H) = A(H) - Z4*A(I) 290 A(I) = Y IF ( L .GT. M ) GOTO 300 F = F + K H = G G = G + L Z2 = Q*A(G) X = Z3*A(G) GOTO 170 300 T = D(M) S = D(M-1) Q = MAG(A(KM)) C ------------------------------ C |*** TEST FOR CONVERGENCE ***| C ------------------------------ IF ( ((T*Q)*T1)*((S*Q)*T1) .GT. 1. ) GOTO 100 L1 = L1 + 1 IF ( L1 .GT. 30 ) GOTO 310 R = MAG(A(KL)) U = AMAX1(Q,R) IF ( U .EQ. 0. ) GOTO 310 IF ( (S*Q)*(Q/U) .GT. T2*(T*R)*(R/U) ) GOTO 100 310 L1 = 0 E(M) = A(KL)*D(M) KL = KM - 1 KM = KM - M M = M - 1 IF ( M .GT. K0 ) GOTO 300 320 E(M) = A(KL)*D(M) KL = KM - 1 KM = KM - M M = M - 1 IF ( M .GT. 1 ) GOTO 30 330 E(1) = A(1)*D(1) P(1) = N RETURN 340 M = N - M + 1 WRITE(6,*) 'SINCE THE STOPPING CRITERION NOT SATISFIED' WRITE(6,*) 'AFTER',NS,'ITERATIONS, WE STOP WHILE COMPUTING' WRITE(6,*) 'EIGENVALUE NUMBER',M P(1) = M RETURN END