SUBROUTINE LZHES(N, A, NA, B, NB, X, NX, WANTX) 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*16 Y, W, Z, A(NA,N), B(NB,N), X(NX,N) REAL*8 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, * 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*16 A(NA,N), B(NB,N), EIGA(N), EIGB(N) COMPLEX*16 S, W, Y, Z COMPLEX*16 X(NX,N) INTEGER ITER(N) COMPLEX*16 ANNM1, ALFM, BETM, D, SL, DEN, NUM, ANM1M1 REAL*8 EPSA, EPSB, SS, R, ANORM, BNORM, ANI, BNI, C REAL*8 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 = SQRT(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