C C ________________________________________________________ C | | C | COMPUTE DIAGONALIZATION OF A COMPLEX HESSENBERG MATRIX | C | | C | INPUT: | C | | C | V --EITHER THE IDENTITY MATRIX OR THE SIM- | C | ILARITY TRANSFORMATION USED TO OBTAIN | C | HESSENBERG FORM | C | | C | LV --LEADING (ROW) DIMENSION OF COMPLEX | C | ARRAY V | 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 | V --COMPLEX ARRAY OF EIGENVECTORS | C | | C | BUILTIN FUNCTIONS: AMAX1,CONJG,SQRT | C | PACKAGE FUNCTIONS: MAG,SQR | C | PACKAGE SUBROUTINES: CEIG | C |________________________________________________________| C SUBROUTINE DAG(E,V,LV,A,N,D,P,W) INTEGER F,F0,G,G0,H,I,IG,II,J,K,KL,KM,KS,K0,K1 INTEGER L,LL,LV,L0,L1,M,N,NS,O COMPLEX A(1),E(1),V(LV,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 DO 10 I = 1,N 10 D(I) = 1. IF ( N .EQ. 1 ) GOTO 440 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 L1 = 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 430 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) GOTO 410 100 LL = LL + 1 IF ( LL .GT. NS ) GOTO 550 T = 1. DO 110 I = K0,M 110 IF ( D(I) .LT. T ) T = D(I) IF ( T .GT. B ) GOTO 210 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 T = P(K) IF ( K0 .EQ. 1 ) GOTO 150 II = F - K0 + 1 H = F - 1 DO 140 I = II,H 140 A(I) = A(I)*T 150 H = F - K0 DO 160 I = F,G 160 A(I) = A(I)*T*P(I-H) DO 170 I = 1,N 170 V(I,K) = V(I,K)*T K = L L = L + 1 F = F + K G = G + L IF ( K .LT. M ) GOTO 130 IF ( K .GT. M ) GOTO 180 G = G - 1 GOTO 130 180 G = G - 1 190 IF ( K .GT. N ) GOTO 210 H = F - K0 DO 200 I = F,G 200 A(I) = A(I)*P(I-H) K = K + 1 F = F + K G = G + K GOTO 190 C ----------------------- C |*** COMPUTE SHIFT ***| C ----------------------- 210 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) K = K0 L = L0 F = F0 G = G0 KS = K K1 = K - 1 Z1 = D(K)*A(F) - Z Z2 = D(K)*A(G) C --------------------------------- C |*** COMPUTE GIVENS ROTATION ***| C --------------------------------- 220 T = MAG(Z1) + MAG(Z2) IF ( T .EQ. 0. ) GOTO 240 T = 1./T Q = SQR(Z1,T) R = SQR(Z2,T) Q = D(K)*Q R = D(L)*R IF ( Q .GT. R ) GOTO 230 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 250 230 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 250 240 P(K) = 0. Z3 = (0.,0.) Z4 = (0.,0.) E(K) = Z3 W(K) = Z4 Q = D(L) 250 IF ( K .GT. KS ) GOTO 270 C ------------------------------------------ C |*** PROCESS INITIAL COLUMN FROM LEFT ***| C ------------------------------------------ IF ( P(K) .EQ. 0. ) GOTO 260 Y = A(G) + Z3*A(F) A(G) = Z4*A(G) - A(F) A(F) = Y GOTO 290 260 Y = A(F) + Z3*A(G) A(G) = A(G) - Z4*A(F) A(F) = Y GOTO 290 C ----------------------------------- C |*** PROCESS COLUMNS FROM LEFT ***| C ----------------------------------- 270 I = G - 1 IF ( P(K) .EQ. 0. ) GOTO 280 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 290 K0 = K L0 = L F0 = I G0 = G GOTO 290 280 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 290 K0 = K L0 = L F0 = I G0 = G 290 J = K K = L L = L + 1 DO 310 I = KS,J O = G + I H = O + 1 IF ( P(I) .EQ. 0. ) GOTO 300 Y = A(H) + E(I)*A(O) A(H) = W(I)*A(H) - A(O) A(O) = Y GOTO 310 300 Y = A(O) + E(I)*A(H) A(H) = A(H) - W(I)*A(O) A(O) = Y 310 CONTINUE C ------------------------------------ C |*** PROCESS COLUMNS FROM RIGHT ***| C ------------------------------------ Z3 = CONJG(Z3) Z4 = CONJG(Z4) O = F - K1 IF ( P(J) .EQ. 0. ) GOTO 340 Z1 = Q*A(G+K) - Z*W(J) DO 320 I = O,G H = I + K Y = A(H) + Z3*A(I) A(H) = Z4*A(H) - A(I) 320 A(I) = Y DO 330 I = 1,N Y = V(I,K) + Z3*V(I,J) V(I,K) = Z4*V(I,K) - V(I,J) 330 V(I,J) = Y F = F + K H = G G = G + L IF ( L .GT. M ) GOTO 370 Z2 = Q*A(G) X = A(G) A(G) = Z4*A(G) GOTO 220 340 Z1 = Q*A(G+K) - Z DO 350 I = O,G H = I + K Y = A(I) + Z3*A(H) A(H) = A(H) - Z4*A(I) 350 A(I) = Y DO 360 I = 1,N Y = V(I,J) + Z3*V(I,K) V(I,K) = V(I,K) - Z4*V(I,J) 360 V(I,J) = Y F = F + K H = G G = G + L IF ( L .GT. M ) GOTO 370 Z2 = Q*A(G) X = Z3*A(G) GOTO 220 C --------------------------------------- C |*** PROCESS MATRIX TAIL FROM LEFT ***| C --------------------------------------- 370 IG = G II = L 380 IF ( II .GT. N ) GOTO 410 II = II + 1 DO 400 I = KS,J O = IG + I H = O + 1 IF ( P(I) .EQ. 0. ) GOTO 390 Y = A(H) + E(I)*A(O) A(H) = W(I)*A(H) - A(O) A(O) = Y GOTO 400 390 Y = A(O) + E(I)*A(H) A(H) = A(H) - W(I)*A(O) A(O) = Y 400 CONTINUE IG = IG + II GOTO 380 410 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 420 R = MAG(A(KL)) U = AMAX1(Q,R) IF ( U .EQ. 0. ) GOTO 420 IF ( (S*Q)*(Q/U) .GT. T2*(T*R)*(R/U) ) GOTO 100 420 L1 = 0 E(M) = A(KL)*D(M) KL = KM - 1 KM = KM - M M = M - 1 IF ( M .GT. K0 ) GOTO 410 430 E(M) = A(KL)*D(M) KL = KM - 1 KM = KM - M M = M - 1 IF ( M .GT. 1 ) GOTO 30 440 E(1) = A(1)*D(1) P(1) = N C ------------------------------ C |*** COMPUTE EIGENVECTORS ***| C ------------------------------ KM = (N*(N+1))/2 - 1 K = N 450 X = E(K) W(K) = (1.,0.) L = K K = K - 1 IF ( K .EQ. 0 ) GOTO 520 DO 460 I = 1,K 460 W(I) = -A(I+KM) KM = KM - L F = KM J = K 470 Y = A(F+J) - X/D(J) Z = (0.,0.) IF ( MAG(Y) .NE. 0. ) Z = W(J)/Y W(J) = Z IF ( J .EQ. 1 ) GOTO 490 G = F F = F - J J = J - 1 DO 480 I = 1,J 480 W(I) = W(I) - Z*A(I+G) GOTO 470 490 X = W(L) DO 500 I = 1,N 500 V(I,L) = X*V(I,L) DO 510 J = 1,K X = W(J) DO 510 I = 1,N 510 V(I,L) = V(I,L) + X*V(I,J) C -------------------------------- C |*** NORMALIZE EIGENVECTORS ***| C -------------------------------- 520 T = 0. DO 530 I = 1,N S = MAG(V(I,L)) 530 IF ( T .LT. S ) T = S IF ( T .NE. 0. ) T = 1./T DO 540 I = 1,N 540 V(I,L) = V(I,L)*T IF ( K .GT. 0 ) GOTO 450 RETURN 550 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