SUBROUTINE QR2NOZ(NM, N, LOW, IGH, H, WR, WI, IERR) QR2 10 C PURPOSE C THE FORTRAN SUBROUTINE QR2NOZ COMPUTES THE EIGENVALUES OF A REAL C UPPER HESSENBERG MATRIX USING THE QR METHOD AND REDUCES THE C MATRIX TO A STANDARDIZED QUASI-TRIANGULAR FORM. COMPUTATIONS C ARE DONE IN REAL ARITHMETIC. C THE SUBROUTINE STATEMENT IS C SUBROUTINE QR2NOZ(NM,N,LOW,IGH,H,WR,WI,IERR). C NM IS AN INTEGER INPUT VARIABLE SET EQUAL TO THE ROW DIMENSION C OF THE ARRAY H AS SPECIFIED IN THE CALLING PROGRAM. C N IS AN INTEGER INPUT VARIABLE SET EQUAL TO THE ORDER OF THE C MATRIX H. N .LE. NM C LOW,IGH ARE INTEGER INPUT VARIABLES INDICATING THE BOUNDARY INDICES C FOR THE BALANCED MATRIX. IF THE MATRIX IS NOT BALANCED SET C LOW TO 1 AND IGH TO N. C H IS A REAL TWO-DIMENSIONAL ARRAY WITH ROW DIMENSION NM AND C COLUMN DIMENSION AT LEAST N. ON INPUT IT CONTAINS THE C UPPER HESSENBERG MATRIX OF ORDER N. ON OUTPUT IT CONTAINS C THE STANDARDIZED QUASI-TRIANGULAR MATRIX. C WR,WI ARE REAL OUTPUT ONE-DIMENSIONAL VARIABLES OF DIMENSION AT C LEAST N CONTAINING THE REAL AND IMAGINARY PARTS, C RESPECTIVELY, OF THE EIGENVALUES OF THE HESSENBERG MATRIX. C THE EIGENVALUES ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE C PAIRS OF EIGENVALUES APPEAR CONSECUTIVELY WITH THE C EIGENVALUE HAVING THE POSITIVE IMAGINARY PART FIRST. C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE. IF MORE THAN 30 ITERATIONS ARE REQUIRED C TO DETERMINE AN EIGENVALUE, THIS SUBROUTINE TERMINATES WITH C IERR SET EQUAL TO THE INDEX OF THE EIGENVALUE FOR WHICH C FAILURE OCCURS. THE EIGENVALUES IN THE WR AND WI ARRAYS C SHOULD BE CORRECT FOR INDICES IERR+1,IERR+2,...,N. IF ALL C THE EIGENVALUES ARE DETERMINED WITHIN 30 ITERATIONS, IERR C IS SET TO ZERO. DIMENSION H(NM,N), WR(N), WI(N) REAL MACHEP INTEGER EN, ENM2 LOGICAL NOTLAS C MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING THE RELATIVE C PRECISION OF FLOATING POINT ARITHMETIC. THE VALUE BELOW IS C EQUAL TO 2.**(-46), WHICH IS APPROPRIATE FOR CDC 6000-SERIES C MACHINES. DATA MACHEP /16424000000000000000B/ IERR = 0 C STORE ROOTS ISOLATED BY BALANC DO 10 I=1,N IF (I.GE.LOW .AND. I.LE.IGH) GO TO 10 WR(I) = H(I,I) WI(I) = 0.0 10 CONTINUE EN = IGH T = 0.0 C SEARCH FOR NEXT EIGENVALUES 20 IF (EN.LT.LOW) RETURN ITS = 0 NA = EN - 1 ENM2 = NA - 1 C LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT C FOR L=EN STEP -1 UNTIL LOW DO 30 IF (EN.EQ.LOW) GO TO 50 DO 40 LL=LOW,NA L = EN + LOW - LL IF (ABS(H(L,L-1)).LE.MACHEP*(ABS(H(L-1,L-1))+ABS(H(L,L)))) * GO TO 60 40 CONTINUE 50 L = LOW C FORM SHIFT 60 X = H(EN,EN) IF (L.EQ.EN) GO TO 200 Y = H(NA,NA) W = H(EN,NA)*H(NA,EN) IF (L.EQ.NA) GO TO 210 IF (ITS.EQ.30) GO TO 270 IF (ITS.NE.10 .AND. ITS.NE.20) GO TO 80 C FORM EXCEPTIONAL SHIFT T = T + X DO 70 I=LOW,EN H(I,I) = H(I,I) - X 70 CONTINUE S = ABS(H(EN,NA)) + ABS(H(NA,ENM2)) X = 0.75*S Y = X W = -0.4275*S*S 80 ITS = ITS + 1 C LOOK FOR TWO CONSECUTIVE SMALL SUB-DIAGONAL C ELEMENTS. FOR M=EN-2 STEP -1 UNTIL L DO DO 90 MM=L,ENM2 M = ENM2 + L - MM ZZ = H(M,M) R = X - ZZ S = Y - ZZ P = (R*S-W)/H(M+1,M) + H(M,M+1) Q = H(M+1,M+1) - ZZ - R - S R = H(M+2,M+1) S = ABS(P) + ABS(Q) + ABS(R) P = P/S Q = Q/S R = R/S IF (M.EQ.L) GO TO 100 IF (ABS(H(M,M-1))*(ABS(Q)+ABS(R)).LE.MACHEP*ABS(P)* * (ABS(H(M-1,M-1))+ABS(ZZ)+ABS(H(M+1,M+1)))) GO TO 100 90 CONTINUE 100 MP2 = M + 2 DO 110 I=MP2,EN H(I,I-2) = 0.0 IF (I.EQ.MP2) GO TO 110 H(I,I-3) = 0.0 110 CONTINUE C DOUBLE QR STEP INVOLVING ROWS L TO EN C AND COLUMNS M TO EN. DO 190 K=M,NA NOTLAS = K.NE.NA IF (K.EQ.M) GO TO 120 P = H(K,K-1) Q = H(K+1,K-1) R = 0.0 IF (NOTLAS) R = H(K+2,K-1) X = ABS(P) + ABS(Q) + ABS(R) IF (X.EQ.0.0) GO TO 190 P = P/X Q = Q/X R = R/X 120 S = SIGN(SQRT(P*P+Q*Q+R*R),P) IF (K.EQ.M) GO TO 130 H(K,K-1) = -S*X GO TO 140 130 IF (L.NE.M) H(K,K-1) = -H(K,K-1) 140 P = P + S X = P/S Y = Q/S ZZ = R/S Q = Q/P R = R/P C ROW MODIFICATION DO 160 J=K,N P = H(K,J) + Q*H(K+1,J) IF (.NOT.NOTLAS) GO TO 150 P = P + R*H(K+2,J) H(K+2,J) = H(K+2,J) - P*ZZ 150 H(K+1,J) = H(K+1,J) - P*Y H(K,J) = H(K,J) - P*X 160 CONTINUE J = MIN0(EN,K+3) C COLUMN MODIFICATION DO 180 I=1,J P = X*H(I,K) + Y*H(I,K+1) IF (.NOT.NOTLAS) GO TO 170 P = P + ZZ*H(I,K+2) H(I,K+2) = H(I,K+2) - P*R 170 H(I,K+1) = H(I,K+1) - P*Q H(I,K) = H(I,K) - P 180 CONTINUE 190 CONTINUE GO TO 30 C ONE ROOT FOUND 200 H(EN,EN) = X + T WR(EN) = H(EN,EN) WI(EN) = 0.0 EN = NA GO TO 20 C TWO ROOTS FOUND 210 P = (Y-X)/2.0 Q = P*P + W ZZ = SQRT(ABS(Q)) H(EN,EN) = X + T X = H(EN,EN) H(NA,NA) = Y + T IF (Q.LT.0.0) GO TO 220 ZZ = P + SIGN(ZZ,P) C REAL PAIR WR(NA) = X + ZZ WR(EN) = WR(NA) IF (ZZ.NE.0.0) WR(EN) = X - W/ZZ WI(NA) = 0.0 WI(EN) = 0.0 X = H(EN,NA) R = SQRT(X*X+ZZ*ZZ) P = X/R Q = ZZ/R GO TO 230 C COMPLEX PAIR 220 WR(NA) = X + P WR(EN) = X + P WI(NA) = ZZ WI(EN) = -ZZ C MAKE DIAGONAL ELEMENTS EQUAL IF (P.EQ.0.0) GO TO 260 BPC = H(EN,NA) + H(NA,EN) TX = SQRT(BPC*BPC+4.0*P*P) Q = SQRT(.5*(1.0+ABS(BPC)/TX)) P = SIGN(P/(Q*TX),-BPC*P) C ROW MODIFICATION 230 DO 240 J=NA,N ZZ = H(NA,J) H(NA,J) = Q*ZZ + P*H(EN,J) H(EN,J) = Q*H(EN,J) - P*ZZ 240 CONTINUE C COLUMN MODIFICATION DO 250 I=1,N ZZ = H(I,NA) H(I,NA) = Q*ZZ + P*H(I,EN) H(I,EN) = Q*H(I,EN) - P*ZZ 250 CONTINUE 260 EN = ENM2 GO TO 20 270 IERR = EN RETURN END SUBROUTINE CONDIT(NM, N, A, V1, V2, WI, COND) CON 10 C CONDIT COMPUTES THE CONDITION NUMBERS OF THE EIGENVALUES OF A C STANDARDIZED QUASI-TRIANGULAR MATRIX. C THE SUBROUTINE STATEMENT IS C SUBROUTINE CONDIT(NM,N,A,V1,V2,WI,COND). C ON INPUT C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO DIMENSIONAL C ARRAY AS DECLARED IN THE CALLING PROGRAM. C N IS THE ORDER OF THE MATRIX. N.LE.NM C A CONTAINS THE STANDARDIZED QUASI-TRIANGULAR MATRIX PRODUCED C BY QR2NOZ. C WI CONTAINS THE IMAGINARY PARTS OF THE EIGENVALUES. THE C EIGENVALUES ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE C PAIRS APPEAR CONSECUTIVELY. C V1,V2 ARE FOR TEMPORARY STORAGE. C ON OUTPUT C A IS UNALTERED. C COND CONTAINS THE CONDITION NUMBERS CORRESPONDING TO THE C EIGENVALUES IN (V2,WI). COND = 1./TOL IF THE USUAL C FORMULA WOULD CAUSE OVERFLOW OR YIELD A VALUE EXCEEDING C 1/TOL. TOL NEED NOT DEPEND ON THE COMPUTER. C V2 CONTAINS THE REAL PARTS OF THE EIGENVALUES. C TYPICAL USAGE C DIMENSION A(50,50),WR(50),WI(50),COND(50),ORT(50) C *******************ENTER MATRIX A AND DIMENSIONS N,NM***************** C LOW=1 C IGH=N C CALL ORTHES(NM,N,LOW,IGH,A,ORT) C CALL QR2NOZ(NM,N,LOW,IGH,A,WR,WI,IERR) C CALL CONDIT(NM,N,A,ORT,WR,WI,COND) C ********************************************************************** C NOTE THE USE OF ORT AND WR IN CONDIT DIMENSION A(NM,NM), V1(NM), V2(NM), WI(NM), COND(NM) DIMENSION R1(2), R2(2) DATA TOL /1.E-30/ C ---------------------------------------------------------------------- I = 1 10 IF (I.GT.N) GO TO 190 VALR = A(I,I) VALI = WI(I) VALI2 = VALI*VALI C NJ GIVES EIGENVECTOR TYPE, 0 FOR COLUMN, 1 FOR ROW NJ = 0 C INITIALIZE NONZERO ELEMENTS OF EIGENVECTOR (V1,V2) V1(I) = 1.0 V2(I) = 0.0 20 J = I - 1 + 2*NJ IF (VALI.EQ.0.0) GO TO 30 V2(I+1) = VALI/A(I,I+1) V1(I+1) = 0.0 IF (NJ.EQ.1) V2(I+1) = 1.0/V2(I+1) J = I - 1 + 3*NJ C FIND THE INDICES OF ELEMENTS COMPUTED SO FAR 30 KS = J + 1 + NJ*(I-J-1) KF = I + 1 + NJ*(J-I-2) IF (VALI.EQ.0.0 .AND. NJ.EQ.0) KF = KF - 1 C TEST FOR COMPLETION OF EIGENVECTOR IF ((J+NJ.LT.1) .OR. (J+NJ.GT.N+1)) GO TO 130 C ********************************************************************** C *SOLVE -D*V + V*E = R FOR V = (V1,V2). D IS A DIAGONAL BLOCK IN ROWS C * J1,J2, AND E IS THE REAL CANONICAL FORM OF THE ITH EIGENVALUE. C * EITHER D OR E OR BOTH CAN BE 1 BY 1 C ********************************************************************** C FIND J1 AND J2 (J1.LE.J2) FOR ALL CASES JJ = J IF (WI(J).NE.0.0) JJ = J - 1 + 2*NJ J0 = NJ*(J-JJ) J1 = JJ + J0 J2 = J - J0 D1 = VALR - A(J,J) C CALCULATE RIGHT HAND SIDE R DO 70 L=J1,J2 LJ = L - J1 + 1 R1(LJ) = 0.0 R2(LJ) = 0.0 IF (VALI.NE.0.0) GO TO 50 DO 40 K=KS,KF LK = NJ*(K-L) AA = A(L+LK,K-LK) R1(LJ) = R1(LJ) + AA*V1(K) 40 CONTINUE GO TO 70 50 DO 60 K=KS,KF LK = NJ*(K-L) AA = A(L+LK,K-LK) R1(LJ) = R1(LJ) + AA*V1(K) R2(LJ) = R2(LJ) + AA*V2(K) 60 CONTINUE 70 CONTINUE IF (JJ.NE.J) GO TO 100 C *****************************D IS 1 BY 1 ***************************** IF (VALI.NE.0.0) GO TO 80 C E IS 1 BY 1 (D IS 1 BY 1) IF (ABS(D1).LT.TOL*ABS(R1(1))) GO TO 180 V1(J) = 0.0 V2(J) = 0.0 IF (D1.NE.0.0) V1(J) = R1(1)/D1 GO TO 90 C E IS 2 BY 2 ( D IS 1 BY 1 ) 80 DEN = D1*D1 + VALI2 VAL = VALI*(-1.0)**NJ V1(J) = R1(1)*D1 + R2(1)*VAL V2(J) = R2(1)*D1 - R1(1)*VAL VMAX = AMAX1(ABS(V1(J)),ABS(V2(J))) IF (DEN.LT.TOL*VMAX) GO TO 180 V1(J) = V1(J)/DEN V2(J) = V2(J)/DEN C NEXT J 90 J = J - 1 + 2*NJ GO TO 30 C ********************************D IS 2 BY 2*************************** 100 IF (VALI.NE.0.0) GO TO 110 C E IS 1 BY 1 (D IS 2 BY 2) DEN = D1*D1 + WI(J)**2 V2(J1) = 0.0 V2(J2) = 0.0 V1(J1) = R1(1)*D1 + R1(2)*A(JJ,J) V1(J2) = R1(1)*A(J,JJ) + R1(2)*D1 VMAX = AMAX1(ABS(V1(J1)),ABS(V1(J2))) IF (DEN.LT.TOL*VMAX) GO TO 180 V1(J1) = V1(J1)/DEN V1(J2) = V1(J2)/DEN GO TO 120 C E IS 2 BY 2 (D IS 2 BY 2). CLOSED FORM SOLUTION 110 B = A(JJ,J) C = A(J,JJ) VAL = VALI*(-1.0)**NJ BXC = B*C H = D1*D1 + VALI2 - BXC E = D1*H F = VAL*(H+2.0*BXC) G = H - 2.0*VALI2 H = 2.0*D1*VAL V1(J1) = R1(1)*E + R2(1)*F + R1(2)*B*G + R2(2)*B*H V2(J1) = -R1(1)*F + R2(1)*E - R1(2)*B*H + R2(2)*B*G V1(J2) = R1(1)*C*G + R2(1)*C*H + R1(2)*E + R2(2)*F V2(J2) = -R1(1)*C*H + R2(1)*C*G - R1(2)*F + R2(2)*E VMAX = AMAX1(ABS(V1(J1)),ABS(V2(J1)),ABS(V1(J2)),ABS(V2(J2))) DEN = G*G + H*H IF (DEN.LT.TOL*VMAX) GO TO 180 IF (DEN.EQ.0.0) GO TO 120 V1(J1) = V1(J1)/DEN V2(J1) = V2(J1)/DEN V1(J2) = V1(J2)/DEN V2(J2) = V2(J2)/DEN C NEXT J 120 J = J - 2 + 4*NJ GO TO 30 C COMPUTE EIGENVECTOR NORM 130 VMAX = 0.0 DO 140 K=KS,KF VMAX = VMAX + V1(K)**2 + V2(K)**2 140 CONTINUE IF (NJ.EQ.1) GO TO 150 C PREPARE TO COMPUTE ROW EIGENVECTOR NJ = 1 CNORM2 = VMAX GO TO 20 C COMPUTE CONDITION NUMBER 150 COND(I) = SQRT(CNORM2*VMAX) 160 IF (VALI.EQ.0.0) GO TO 170 COND(I) = COND(I)/2.0 COND(I+1) = COND(I) I = I + 1 C NEXT I 170 I = I + 1 GO TO 10 C DEFECTIVE CASE 180 COND(I) = 1.0/TOL GO TO 160 C PLACE REAL PART OF EIGENVALUE IN V2 190 DO 200 I=1,N V2(I) = A(I,I) 200 CONTINUE RETURN END