SUBROUTINE LZHES(N, A, NA, B, NB, X, NX, WANTX) LZH 10 C THIS SUBROUTINE REDUCES THE COMPLEX MATRIX A TO UPPER C HESSENBERG FORM AND REDUCES THE COMPLEX MATRIX B TO C TRIANGULAR FORM C INPUT PARAMETERS.. C N THE ORDER OF THE A AND B MATRICES C A A COMPLEX MATRIX C NA THE ROW DIMENSION OF THE A MATRIX C B A COMPLEX MATRIX C NB THE ROW DIMENSION OF THE B MATRIX C NX THE ROW DIMENSION OF THE X MATRIX C WANTX A LOGICAL VARIABLE WHICH IS SET TO .TRUE. IF C THE EIGENVECTORS ARE WANTED. OTHERWISE IT SHOULD C BE SET TO .FALSE. C OUTPUT PARAMETERS.. C A ON OUTPUT A IS AN UPPER HESSENBERG MATRIX, THE C ORIGINAL MATRIX HAS BEEN DESTROYED C B AN UPPER TRIANGULAR MATRIX, THE ORIGINAL MATRIX C HAS BEEN DESTROYED C X CONTAINS THE TRANSFORMATIONS NEEDED TO COMPUTE C THE EIGENVECTORS OF THE ORIGINAL SYSTEM COMPLEX Y, W, Z, A(NA,N), B(NB,N), X(NX,N) REAL C, D LOGICAL WANTX NM1 = N - 1 C REDUCE B TO TRIANGULAR FORM USING ELEMENTARY C TRANSFORMATIONS DO 80 I=1,NM1 D = 0.00 IP1 = I + 1 DO 10 K=IP1,N Y = B(K,I) C = ABS(REAL(Y)) + ABS(AIMAG(Y)) IF (C.LE.D) GO TO 10 D = C II = K 10 CONTINUE IF (D.EQ.0.0) GO TO 80 Y = B(I,I) IF (D.LE.ABS(REAL(Y))+ABS(AIMAG(Y))) GO TO 40 C MUST INTERCHANGE DO 20 J=1,N Y = A(I,J) A(I,J) = A(II,J) A(II,J) = Y 20 CONTINUE DO 30 J=I,N Y = B(I,J) B(I,J) = B(II,J) B(II,J) = Y 30 CONTINUE 40 DO 70 J=IP1,N Y = B(J,I)/B(I,I) IF (REAL(Y).EQ.0.0 .AND. AIMAG(Y).EQ.0.0) GO TO 70 DO 50 K=1,N A(J,K) = A(J,K) - Y*A(I,K) 50 CONTINUE DO 60 K=IP1,N B(J,K) = B(J,K) - Y*B(I,K) 60 CONTINUE 70 CONTINUE B(IP1,I) = (0.0,0.0) 80 CONTINUE C INITIALIZE X IF (.NOT.WANTX) GO TO 110 DO 100 I=1,N DO 90 J=1,N X(I,J) = (0.0,0.0) 90 CONTINUE X(I,I) = (1.0,0.00) 100 CONTINUE C REDUCE A TO UPPER HESSENBERG FORM 110 NM2 = N - 2 IF (NM2.LT.1) GO TO 270 DO 260 J=1,NM2 JM2 = NM1 - J JP1 = J + 1 DO 250 II=1,JM2 I = N + 1 - II IM1 = I - 1 IMJ = I - J W = A(I,J) Z = A(IM1,J) IF (ABS(REAL(W))+ABS(AIMAG(W)).LE.ABS(REAL(Z)) * +ABS(AIMAG(Z))) GO TO 140 C MUST INTERCHANGE ROWS DO 120 K=J,N Y = A(I,K) A(I,K) = A(IM1,K) A(IM1,K) = Y 120 CONTINUE DO 130 K=IM1,N Y = B(I,K) B(I,K) = B(IM1,K) B(IM1,K) = Y 130 CONTINUE 140 Z = A(I,J) IF (REAL(Z).EQ.0.0 .AND. AIMAG(Z).EQ.0.0) GO TO 170 Y = Z/A(IM1,J) DO 150 K=JP1,N A(I,K) = A(I,K) - Y*A(IM1,K) 150 CONTINUE DO 160 K=IM1,N B(I,K) = B(I,K) - Y*B(IM1,K) 160 CONTINUE C TRANSFORMATION FROM THE RIGHT 170 W = B(I,IM1) Z = B(I,I) IF (ABS(REAL(W))+ABS(AIMAG(W)).LE.ABS(REAL(Z)) * +ABS(AIMAG(Z))) GO TO 210 C MUST INTERCHANGE COLUMNS DO 180 K=1,I Y = B(K,I) B(K,I) = B(K,IM1) B(K,IM1) = Y 180 CONTINUE DO 190 K=1,N Y = A(K,I) A(K,I) = A(K,IM1) A(K,IM1) = Y 190 CONTINUE IF (.NOT.WANTX) GO TO 210 DO 200 K=IMJ,N Y = X(K,I) X(K,I) = X(K,IM1) X(K,IM1) = Y 200 CONTINUE 210 Z = B(I,IM1) IF (REAL(Z).EQ.0.0 .AND. AIMAG(Z).EQ.0.0) GO TO 250 Y = Z/B(I,I) DO 220 K=1,IM1 B(K,IM1) = B(K,IM1) - Y*B(K,I) 220 CONTINUE B(I,IM1) = (0.0,0.0) DO 230 K=1,N A(K,IM1) = A(K,IM1) - Y*A(K,I) 230 CONTINUE IF (.NOT.WANTX) GO TO 250 DO 240 K=IMJ,N X(K,IM1) = X(K,IM1) - Y*X(K,I) 240 CONTINUE 250 CONTINUE A(JP1+1,J) = (0.0,0.0) 260 CONTINUE 270 RETURN END SUBROUTINE LZIT(N, A, NA, B, NB, X, NX, WANTX, ITER, EIGA, LZI 10 * EIGB) C THIS SUBROUTINE SOLVES THE GENERALIZED EIGENVALUE PROBLEM C A X = LAMBDA B X C WHERE A IS A COMPLEX UPPER HESSENBERG MATRIX OF C ORDER N AND B IS A COMPLEX UPPER TRIANGULAR MATRIX OF ORDER N C INPUT PARAMETERS C N ORDER OF A AND B C A AN N X N UPPER HESSENBERG COMPLEX MATRIX C NA THE ROW DIMENSION OF THE A MATRIX C B AN N X N UPPER TRIANGULAR COMPLEX MATRIX C NB THE ROW DIMENSION OF THE B MATRIX C X CONTAINS TRANSFORMATIONS TO OBTAIN EIGENVECTORS OF C ORIGINAL SYSTEM. IF EIGENVECTORS ARE REQUESTED AND QZHES C IS NOT CALLED, X SHOULD BE SET TO THE IDENTITY MATRIX C NX THE ROW DIMENSION OF THE X MATRIX C WANTX LOGICAL VARIABLE WHICH SHOULD BE SET TO .TRUE. C IF EIGENVECTORS ARE WANTED. OTHERWISE IT C SHOULD BE SET TO .FALSE. C OUTPUT PARAMETERS C X THE ITH COLUMN CONTAINS THE ITH EIGENVECTOR C IF EIGENVECTORS ARE REQUESTED C ITER AN INTEGER ARRAY OF LENGTH N WHOSE ITH ENTRY C CONTAINS THE NUMBER OF ITERATIONS NEEDED TO FIND C THE ITH EIGENVALUE. FOR ANY I IF ITER(I) =-1 THEN C AFTER 30 ITERATIONS THERE HAS NOT BEEN A SUFFICIENT C DECREASE IN THE LAST SUBDIAGONAL ELEMENT OF A C TO CONTINUE ITERATING. C EIGA A COMPLEX ARRAY OF LENGTH N CONTAINING THE DIAGONAL OF A C EIGB A COMPLEX ARRAY OF LENGTH N CONTAINING THE DIAGONAL OF B C THE ITH EIGENVALUE CAN BE FOUND BY DIVIDING EIGA(I) BY C EIGB(I). WATCH OUT FOR EIGB(I) BEING ZERO COMPLEX A(NA,N), B(NB,N), EIGA(N), EIGB(N) COMPLEX S, W, Y, Z, CSQRT COMPLEX X(NX,N) INTEGER ITER(N) COMPLEX ANNM1, ALFM, BETM, D, SL, DEN, NUM, ANM1M1 REAL EPSA, EPSB, SS, R, ANORM, BNORM, ANI, BNI, C REAL D0, D1, D2, E0, E1 LOGICAL WANTX NN = N C COMPUTE THE MACHINE PRECISION TIMES THE NORM OF A AND B ANORM = 0. BNORM = 0. DO 30 I=1,N ANI = 0. IF (I.EQ.1) GO TO 10 Y = A(I,I-1) ANI = ANI + ABS(REAL(Y)) + ABS(AIMAG(Y)) 10 BNI = 0. DO 20 J=I,N ANI = ANI + ABS(REAL(A(I,J))) + ABS(AIMAG(A(I,J))) BNI = BNI + ABS(REAL(B(I,J))) + ABS(AIMAG(B(I,J))) 20 CONTINUE IF (ANI.GT.ANORM) ANORM = ANI IF (BNI.GT.BNORM) BNORM = BNI 30 CONTINUE IF (ANORM.EQ.0.) ANORM = 1.0 IF (BNORM.EQ.0.) BNORM = 1.0 EPSB = BNORM EPSA = ANORM 40 EPSA = EPSA/2.0 EPSB = EPSB/2.0 C = ANORM + EPSA IF (C.GT.ANORM) GO TO 40 IF (N.LE.1) GO TO 320 50 ITS = 0 NM1 = NN - 1 C CHECK FOR NEGLIGIBLE SUBDIAGONAL ELEMENTS 60 D2 = ABS(REAL(A(NN,NN))) + ABS(AIMAG(A(NN,NN))) DO 70 LB=2,NN L = NN + 2 - LB SS = D2 Y = A(L-1,L-1) D2 = ABS(REAL(Y)) + ABS(AIMAG(Y)) SS = SS + D2 Y = A(L,L-1) R = SS + ABS(REAL(Y)) + ABS(AIMAG(Y)) IF (R.EQ.SS) GO TO 80 70 CONTINUE L = 1 80 IF (L.EQ.NN) GO TO 320 IF (ITS.LT.30) GO TO 90 ITER(NN) = -1 IF (ABS(REAL(A(NN,NM1)))+ABS(AIMAG(A(NN,NM1))).GT.0.8* * ABS(REAL(ANNM1))+ABS(AIMAG(ANNM1))) RETURN 90 IF (ITS.EQ.10 .OR. ITS.EQ.20) GO TO 110 C COMPUTE SHIFT AS EIGENVALUE OF LOWER 2 BY 2 ANNM1 = A(NN,NM1) ANM1M1 = A(NM1,NM1) S = A(NN,NN)*B(NM1,NM1) - ANNM1*B(NM1,NN) W = ANNM1*B(NN,NN)*(A(NM1,NN)*B(NM1,NM1)-B(NM1,NN)*ANM1M1) Y = (ANM1M1*B(NN,NN)-S)/2. Z = CSQRT(Y*Y+W) IF (REAL(Z).EQ.0.0 .AND. AIMAG(Z).EQ.0.0) GO TO 100 D0 = REAL(Y/Z) IF (D0.LT.0.0) Z = -Z 100 DEN = (Y+Z)*B(NM1,NM1)*B(NN,NN) IF (REAL(DEN).EQ.0.0 .AND. AIMAG(DEN).EQ.0.0) DEN = * CMPLX(EPSA,0.0) NUM = (Y+Z)*S - W GO TO 120 C AD-HOC SHIFT 110 Y = A(NM1,NN-2) NUM = CMPLX(ABS(REAL(ANNM1))+ABS(AIMAG(ANNM1)),ABS(REAL(Y)) * +ABS(AIMAG(Y))) DEN = (1.0,0.0) C CHECK FOR 2 CONSECUTIVE SMALL SUBDIAGONAL ELEMENTS 120 IF (NN.EQ.L+1) GO TO 140 D2 = ABS(REAL(A(NM1,NM1))) + ABS(AIMAG(A(NM1,NM1))) E1 = ABS(REAL(ANNM1)) + ABS(AIMAG(ANNM1)) D1 = ABS(REAL(A(NN,NN))) + ABS(AIMAG(A(NN,NN))) NL = NN - (L+1) DO 130 MB=1,NL M = NN - MB E0 = E1 Y = A(M,M-1) E1 = ABS(REAL(Y)) + ABS(AIMAG(Y)) D0 = D1 D1 = D2 Y = A(M-1,M-1) D2 = ABS(REAL(Y)) + ABS(AIMAG(Y)) Y = A(M,M)*DEN - B(M,M)*NUM D0 = (D0+D1+D2)*(ABS(REAL(Y))+ABS(AIMAG(Y))) E0 = E0*E1*(ABS(REAL(DEN))+ABS(AIMAG(DEN))) + D0 IF (E0.EQ.D0) GO TO 150 130 CONTINUE 140 M = L 150 CONTINUE ITS = ITS + 1 W = A(M,M)*DEN - B(M,M)*NUM Z = A(M+1,M)*DEN D1 = ABS(REAL(Z)) + ABS(AIMAG(Z)) D2 = ABS(REAL(W)) + ABS(AIMAG(W)) C FIND L AND M AND SET A=LAM AND B=LBM C NP1 = N + 1 LOR1 = L NNORN = NN IF (.NOT.WANTX) GO TO 160 LOR1 = 1 NNORN = N 160 DO 310 I=M,NM1 J = I + 1 C FIND ROW TRANSFORMATIONS TO RESTORE A TO C UPPER HESSENBERG FORM. APPLY TRANSFORMATIONS C TO A AND B IF (I.EQ.M) GO TO 170 W = A(I,I-1) Z = A(J,I-1) D1 = ABS(REAL(Z)) + ABS(AIMAG(Z)) D2 = ABS(REAL(W)) + ABS(AIMAG(W)) IF (D1.EQ.0.0) GO TO 60 170 IF (D2.GT.D1) GO TO 190 C MUST INTERCHANGE ROWS DO 180 K=I,NNORN Y = A(I,K) A(I,K) = A(J,K) A(J,K) = Y Y = B(I,K) B(I,K) = B(J,K) B(J,K) = Y 180 CONTINUE IF (I.GT.M) A(I,I-1) = A(J,I-1) IF (D2.EQ.0.0) GO TO 220 C THE SCALING OF W AND Z IS DESIGNED TO AVOID A DIVISION BY ZERO C WHEN THE DENOMINATOR IS SMALL Y = CMPLX(REAL(W)/D1,AIMAG(W)/D1)/CMPLX(REAL(Z)/D1,AIMAG(Z)/ * D1) GO TO 200 190 Y = CMPLX(REAL(Z)/D2,AIMAG(Z)/D2)/CMPLX(REAL(W)/D2,AIMAG(W)/ * D2) 200 DO 210 K=I,NNORN A(J,K) = A(J,K) - Y*A(I,K) B(J,K) = B(J,K) - Y*B(I,K) 210 CONTINUE 220 IF (I.GT.M) A(J,I-1) = (0.0,0.0) C PERFORM TRANSFORMATIONS FROM RIGHT TO RESTORE B TO C TRIANGLULAR FORM C APPLY TRANSFORMATIONS TO A Z = B(J,I) W = B(J,J) D2 = ABS(REAL(W)) + ABS(AIMAG(W)) D1 = ABS(REAL(Z)) + ABS(AIMAG(Z)) IF (D1.EQ.0.0) GO TO 60 IF (D2.GT.D1) GO TO 270 C MUST INTERCHANGE COLUMNS DO 230 K=LOR1,J Y = A(K,J) A(K,J) = A(K,I) A(K,I) = Y Y = B(K,J) B(K,J) = B(K,I) B(K,I) = Y 230 CONTINUE IF (I.EQ.NM1) GO TO 240 Y = A(J+1,J) A(J+1,J) = A(J+1,I) A(J+1,I) = Y 240 IF (.NOT.WANTX) GO TO 260 DO 250 K=1,N Y = X(K,J) X(K,J) = X(K,I) X(K,I) = Y 250 CONTINUE 260 IF (D2.EQ.0.0) GO TO 310 Z = CMPLX(REAL(W)/D1,AIMAG(W)/D1)/CMPLX(REAL(Z)/D1,AIMAG(Z)/ * D1) GO TO 280 270 Z = CMPLX(REAL(Z)/D2,AIMAG(Z)/D2)/CMPLX(REAL(W)/D2,AIMAG(W)/ * D2) 280 DO 290 K=LOR1,J A(K,I) = A(K,I) - Z*A(K,J) B(K,I) = B(K,I) - Z*B(K,J) 290 CONTINUE B(J,I) = (0.0,0.0) IF (I.LT.NM1) A(I+2,I) = A(I+2,I) - Z*A(I+2,J) IF (.NOT.WANTX) GO TO 310 DO 300 K=1,N X(K,I) = X(K,I) - Z*X(K,J) 300 CONTINUE 310 CONTINUE GO TO 60 320 CONTINUE EIGA(NN) = A(NN,NN) EIGB(NN) = B(NN,NN) IF (NN.EQ.1) GO TO 330 ITER(NN) = ITS NN = NM1 IF (NN.GT.1) GO TO 50 ITER(1) = 0 GO TO 320 C FIND EIGENVECTORS USING B FOR INTERMEDIATE STORAGE 330 IF (.NOT.WANTX) RETURN M = N 340 CONTINUE ALFM = A(M,M) BETM = B(M,M) B(M,M) = (1.0,0.0) L = M - 1 IF (L.EQ.0) GO TO 370 350 CONTINUE L1 = L + 1 SL = (0.0,0.0) DO 360 J=L1,M SL = SL + (BETM*A(L,J)-ALFM*B(L,J))*B(J,M) 360 CONTINUE Y = BETM*A(L,L) - ALFM*B(L,L) IF (REAL(Y).EQ.0.0 .AND. AIMAG(Y).EQ.0.0) Y = * CMPLX((EPSA+EPSB)/2.0,0.0) B(L,M) = -SL/Y L = L - 1 370 IF (L.GT.0) GO TO 350 M = M - 1 IF (M.GT.0) GO TO 340 C TRANSFORM TO ORIGINAL COORDINATE SYSTEM M = N 380 CONTINUE DO 400 I=1,N S = (0.0,0.0) DO 390 J=1,M S = S + X(I,J)*B(J,M) 390 CONTINUE X(I,M) = S 400 CONTINUE M = M - 1 IF (M.GT.0) GO TO 380 C NORMALIZE SO THAT LARGEST COMPONENT = 1. M = N 410 CONTINUE SS = 0. DO 420 I=1,N R = ABS(REAL(X(I,M))) + ABS(AIMAG(X(I,M))) IF (R.LT.SS) GO TO 420 SS = R D = X(I,M) 420 CONTINUE IF (SS.EQ.0.0) GO TO 440 DO 430 I=1,N X(I,M) = X(I,M)/D 430 CONTINUE 440 M = M - 1 IF (M.GT.0) GO TO 410 RETURN END