SUBROUTINE HQR3(A, V, N, NLOW, NUP, EPS, ER, EI, TYPE, NA, NV) HQR 10 INTEGER N, NA, NLOW, NUP, NV, TYPE(N) REAL A(NA,N), EI(N), ER(N), EPS, V(NV,N) C HQR3 REDUCES THE UPPER HESSENBERG MATRIX A TO QUASI- C TRIANGULAR FORM BY UNITARY SIMILARITY TRANSFORMATIONS. C THE EIGENVALUES OF A, WHICH ARE CONTAINED IN THE 1X1 C AND 2X2 DIAGONAL BLOCKS OF THE REDUCED MATRIX, ARE C ORDERED IN DESCENDING ORDER OF MAGNITUDE ALONG THE C DIAGONAL. THE TRANSFORMATIONS ARE ACCUMULATED IN THE C ARRAY V. HQR3 REQUIRES THE SUBROUTINES EXCHNG, C QRSTEP, AND SPLIT. THE PARAMETERS IN THE CALLING C SEQUENCE ARE (STARRED PARAMETERS ARE ALTERED BY THE C SUBROUTINE) C *A AN ARRAY THAT INITIALLY CONTAINS THE N X N C UPPER HESSENBERG MATRIX TO BE REDUCED. ON C RETURN A CONTAINS THE REDUCED, QUASI- C TRIANGULAR MATRIX. C *V AN ARRAY THAT CONTAINS A MATRIX INTO WHICH C THE REDUCING TRANSFORMATIONS ARE TO BE C MULTIPLIED. C N THE ORDER OF THE MATRICES A AND V. C NLOW A(NLOW,NLOW-1) AND A(NUP,+1,NUP) ARE C NUP ASSUMED TO BE ZERO, AND ONLY ROWS NLOW C THROUGH NUP AND COLUMNS NLOW THROUGH C NUP ARE TRANSFORMED, RESULTING IN THE C CALCULATION OF EIGENVALUES NLOW C THROUGH NUP. C EPS A CONVERGENCE CRITERION. C *ER AN ARRAY THAT ON RETURN CONTAINS THE REAL C PARTS OF THE EIGENVALUES. C *EI AN ARRAY THAT ON RETURN CONTAINS THE C IMAGINARY PARTS OF THE EIGENVALUES. C *TYPE AN INTEGER ARRAY WHOSE I-TH ENTRY IS C 0 IF THE I-TH EIGENVALUE IS REAL, C 1 IF THE I-TH EIGENVALUE IS COMPLEX C WITH POSITIVE IMAGINARY PART. C 2 IF THE I-TH EIGENVALUE IS COMPLEX C WITH NEGATIVE IMAGINARY PART, C -1 IF THE I-TH EIGENVALUE WAS NOT C CALCULATED SUCCESSFULLY. C NA THE FIRST DIMENSION OF THE ARRAY A. C NV THE FIRST DIMENSION OF THE ARRAY V. C THE CONVERGENCE CRITERION EPS IS USED TO DETERMINE C WHEN A SUBDIAGONAL ELEMENT OF A IS NEGLIGIBLE. C SPECIFICALLY A(I+1,I) IS REGARDED AS NEGLIGIBLE C IF C ABS(A(I+1),I)) .LE. EPS*(ABS(A(I,I))+ABS(A(I+1,I+1))). C THIS MEANS THAT THE FINAL MATRIX RETURNED BY THE C PROGRAM WILL BE EXACTLY SIMILAR TO A + E WHERE E IS C OF ORDER EPS*NORM(A), FOR ANY REASONABLY BALANCED NORM C SUCH AS THE ROW-SUM NORM. C INTERNAL VARIABLES INTEGER I, IT, L, MU, NL, NU REAL E1, E2, P, Q, R, S, T, W, X, Y, Z LOGICAL FAIL C INITIALIZE. DO 10 I=NLOW,NUP TYPE(I) = -1 10 CONTINUE T = 0. C MAIN LOOP. FIND AND ORDER EIGENVALUES. NU = NUP 20 IF (NU.LT.NLOW) GO TO 240 IT = 0 C QR LOOP. FIND NEGLIGIBLE ELEMENTS AND PERFORM C QR STEPS. 30 CONTINUE C SEARCH BACK FOR NEGLIGIBLE ELEMENTS. L = NU 40 CONTINUE IF (L.EQ.NLOW) GO TO 50 IF (ABS(A(L,L-1)).LE.EPS*(ABS(A(L-1,L-1))+ABS(A(L,L)))) GO TO * 50 L = L - 1 GO TO 40 50 CONTINUE C TEST TO SEE IF AN EIGENVALUE OR A 2X2 BLOCK C HAS BEEN FOUND. X = A(NU,NU) IF (L.EQ.NU) GO TO 160 Y = A(NU-1,NU-1) W = A(NU,NU-1)*A(NU-1,NU) IF (L.EQ.NU-1) GO TO 100 C TEST ITERATION COUNT. IF IT IS 30 QUIT. IF C IT IS 10 OR 20 SET UP AN AD-HOC SHIFT. IF (IT.EQ.30) GO TO 240 IF (IT.NE.10 .AND. IT.NE.20) GO TO 70 C AD-HOC SHIFT. T = T + X DO 60 I=NLOW,NU A(I,I) = A(I,I) - X 60 CONTINUE S = ABS(A(NU,NU-1)) + ABS(A(NU-1,NU-2)) X = 0.75*S Y = X W = -0.4375*S**2 70 CONTINUE IT = IT + 1 C LOOK FOR TWO CONSECUTIVE SMALL SUB-DIAGONAL C ELEMENTS. NL = NU - 2 80 CONTINUE Z = A(NL,NL) R = X - Z S = Y - Z P = (R*S-W)/A(NL+1,NL) + A(NL,NL+1) Q = A(NL+1,NL+1) - Z - R - S R = A(NL+2,NL+1) S = ABS(P) + ABS(Q) + ABS(R) P = P/S Q = Q/S R = R/S IF (NL.EQ.L) GO TO 90 IF (ABS(A(NL,NL-1))*(ABS(Q)+ABS(R)).LE.EPS*ABS(P)*(ABS(A(NL-1, * NL-1))+ABS(Z)+ABS(A(NL+1,NL+1)))) GO TO 90 NL = NL - 1 GO TO 80 90 CONTINUE C PERFORM A QR STEP BETWEEN NL AND NU. CALL QRSTEP(A, V, P, Q, R, NL, NU, N, NA, NV) GO TO 30 C 2X2 BLOCK FOUND. 100 IF (NU.NE.NLOW+1) A(NU-1,NU-2) = 0. A(NU,NU) = A(NU,NU) + T A(NU-1,NU-1) = A(NU-1,NU-1) + T TYPE(NU) = 0 TYPE(NU-1) = 0 MU = NU C LOOP TO POSITION 2X2 BLOCK. 110 CONTINUE NL = MU - 1 C ATTEMPT TO SPLIT THE BLOCK INTO TWO REAL C EIGENVALUES. CALL SPLIT(A, V, N, NL, E1, E2, NA, NV) C IF THE SPLIT WAS SUCCESSFUL, GO AND ORDER THE C REAL EIGENVALUES. IF (A(MU,MU-1).EQ.0.) GO TO 170 C TEST TO SEE IF THE BLOCK IS PROPERLY POSITIONED, C AND IF NOT EXCHANGE IT IF (MU.EQ.NUP) GO TO 230 IF (MU.EQ.NUP-1) GO TO 130 IF (A(MU+2,MU+1).EQ.0.) GO TO 130 C THE NEXT BLOCK IS 2X2. IF (A(MU-1,MU-1)*A(MU,MU)-A(MU-1,MU)*A(MU,MU-1).GE.A(MU+1, * MU+1)*A(MU+2,MU+2)-A(MU+1,MU+2)*A(MU+2,MU+1)) GO TO 230 CALL EXCHNG(A, V, N, NL, 2, 2, EPS, FAIL, NA, NV) IF (.NOT.FAIL) GO TO 120 TYPE(NL) = -1 TYPE(NL+1) = -1 TYPE(NL+2) = -1 TYPE(NL+3) = -1 GO TO 240 120 CONTINUE MU = MU + 2 GO TO 150 130 CONTINUE C THE NEXT BLOCK IS 1X1. IF (A(MU-1,MU-1)*A(MU,MU)-A(MU-1,MU)*A(MU,MU-1).GE.A(MU+1, * MU+1)**2) GO TO 230 CALL EXCHNG(A, V, N, NL, 2, 1, EPS, FAIL, NA, NV) IF (.NOT.FAIL) GO TO 140 TYPE(NL) = -1 TYPE(NL+1) = -1 TYPE(NL+2) = -1 GO TO 240 140 CONTINUE MU = MU + 1 150 CONTINUE GO TO 110 C SINGLE EIGENVALUE FOUND. 160 NL = 0 A(NU,NU) = A(NU,NU) + T IF (NU.NE.NLOW) A(NU,NU-1) = 0. TYPE(NU) = 0 MU = NU C LOOP TO POSITION ONE OR TWO REAL EIGENVALUES. 170 CONTINUE C POSITION THE EIGENVALUE LOCATED AT A(NL,NL). 180 CONTINUE IF (MU.EQ.NUP) GO TO 220 IF (MU.EQ.NUP-1) GO TO 200 IF (A(MU+2,MU+1).EQ.0.) GO TO 200 C THE NEXT BLOCK IS 2X2. IF (A(MU,MU)**2.GE.A(MU+1,MU+1)*A(MU+2,MU+2)-A(MU+1,MU+2)* * A(MU+2,MU+1)) GO TO 220 CALL EXCHNG(A, V, N, MU, 1, 2, EPS, FAIL, NA, NV) IF (.NOT.FAIL) GO TO 190 TYPE(MU) = -1 TYPE(MU+1) = -1 TYPE(MU+2) = -1 GO TO 240 190 CONTINUE MU = MU + 2 GO TO 210 200 CONTINUE C THE NEXT BLOCK IS 1X1. IF (ABS(A(MU,MU)).GE.ABS(A(MU+1,MU+1))) GO TO 220 CALL EXCHNG(A, V, N, MU, 1, 1, EPS, FAIL, NA, NV) MU = MU + 1 210 CONTINUE GO TO 180 220 CONTINUE MU = NL NL = 0 IF (MU.NE.0) GO TO 170 C GO BACK AND GET THE NEXT EIGENVALUE. 230 CONTINUE NU = L - 1 GO TO 20 C ALL THE EIGENVALUES HAVE BEEN FOUND AND ORDERED. C COMPUTE THEIR VALUES AND TYPE. 240 IF (NU.LT.NLOW) GO TO 260 DO 250 I=NLOW,NU A(I,I) = A(I,I) + T 250 CONTINUE 260 CONTINUE NU = NUP 270 CONTINUE IF (TYPE(NU).NE.-1) GO TO 280 NU = NU - 1 GO TO 310 280 CONTINUE IF (NU.EQ.NLOW) GO TO 290 IF (A(NU,NU-1).EQ.0.) GO TO 290 C 2X2 BLOCK. CALL SPLIT(A, V, N, NU-1, E1, E2, NA, NV) IF (A(NU,NU-1).EQ.0.) GO TO 290 ER(NU) = E1 EI(NU-1) = E2 ER(NU-1) = ER(NU) EI(NU) = -EI(NU-1) TYPE(NU-1) = 1 TYPE(NU) = 2 NU = NU - 2 GO TO 300 290 CONTINUE C SINGLE ROOT. ER(NU) = A(NU,NU) EI(NU) = 0. NU = NU - 1 300 CONTINUE 310 CONTINUE IF (NU.GE.NLOW) GO TO 270 RETURN END SUBROUTINE SPLIT(A, V, N, L, E1, E2, NA, NV) SPL 10 INTEGER L, N, NA, NV REAL A(NA,N), V(NV,N) C GIVEN THE UPPER HESSENBERG MATRIX A WITH A 2X2 BLOCK C STARTING AT A(L,L), SPLIT DETERMINES IF THE C CORRESPONDING EIGENVALUES ARE REAL OR COMPLEX. IF THEY C ARE REAL, A ROTATION IS DETERMINED THAT REDUCES THE C BLOCK TO UPPER TRIANGULAR FORM WITH THE EIGENVALUE C OF LARGEST ABSOLUTE VALUE APPEARING FIRST. THE C ROTATION IS ACCUMULATED IN V. THE EIGENVALUES (REAL C OR COMPLEX) ARE RETURNED IN E1 AND E2. THE PARAMETERS C IN THE CALLING SEQUENCE ARE (STARRED PARAMETERS ARE C ALTERED BY THE SUBROUTINE) C *A THE UPPER HESSENVERG MATRIX WHOSE 2X2 C BLOCK IS TO BE SPLIT. C *V THE ARRAY IN WHICH THE SPLITTING TRANS- C FORMATION IS TO BE ACCUMULATED. C N THE ORDER OF THE MATRIX A. C L THE POSITION OF THE 2X2 BLOCK. C *E1 ON RETURN IF THE EIGENVALUES ARE COMPLEX C *E2 E1 CONTAINS THEIR COMMON REAL PART AND C E2 CONTAINS THE POSITIVE IMAGINARY PART. C IF THE EIGENVALUES ARE REAL, E1 CONTAINS C THE ONE LARGEST IN ABSOLUTE VALUE AND E2 C CONTAINS THE OTHER ONE. C NA THE FIRST DIMENSION OF THE ARRAY A. C NV THE FIRST DIMENSION OF THE ARRAY V. C INTERNAL VARIABLES INTEGER I, J, L1 REAL P, Q, R, T, U, W, X, Y, Z X = A(L+1,L+1) Y = A(L,L) W = A(L,L+1)*A(L+1,L) P = (Y-X)/2. Q = P**2 + W IF (Q.GE.0.) GO TO 10 C COMPLEX EIGENVALUE. E1 = P + X E2 = SQRT(-Q) RETURN 10 CONTINUE C TWO REAL EIGENVALUES. SET UP TRANSFORMATION. Z = SQRT(Q) IF (P.LT.0.) GO TO 20 Z = P + Z GO TO 30 20 CONTINUE Z = P - Z 30 CONTINUE IF (Z.EQ.0.) GO TO 40 R = -W/Z GO TO 50 40 CONTINUE R = 0. 50 CONTINUE IF (ABS(X+Z).GE.ABS(X+R)) Z = R Y = Y - X - Z X = -Z T = A(L,L+1) U = A(L+1,L) IF (ABS(Y)+ABS(U).LE.ABS(T)+ABS(X)) GO TO 60 Q = U P = Y GO TO 70 60 CONTINUE Q = X P = T 70 CONTINUE R = SQRT(P**2+Q**2) IF (R.GT.0.) GO TO 80 E1 = A(L,L) E2 = A(L+1,L+1) A(L+1,L) = 0. RETURN 80 CONTINUE P = P/R Q = Q/R C PREMULTIPLY. DO 90 J=L,N Z = A(L,J) A(L,J) = P*Z + Q*A(L+1,J) A(L+1,J) = P*A(L+1,J) - Q*Z 90 CONTINUE C POSTMULTIPLY. L1 = L + 1 DO 100 I=1,L1 Z = A(I,L) A(I,L) = P*Z + Q*A(I,L+1) A(I,L+1) = P*A(I,L+1) - Q*Z 100 CONTINUE C ACCUMULATE THE TRANSFORMATION IN V. DO 110 I=1,N Z = V(I,L) V(I,L) = P*Z + Q*V(I,L+1) V(I,L+1) = P*V(I,L+1) - Q*Z 110 CONTINUE A(L+1,L) = 0. E1 = A(L,L) E2 = A(L+1,L+1) RETURN END SUBROUTINE EXCHNG(A, V, N, L, B1, B2, EPS, FAIL, NA, NV) EXC 10 INTEGER B1, B2, L, NA, NV REAL A(NA,N), EPS, V(NV,N) LOGICAL FAIL C GIVEN THE UPPER HESSENBERG MATRIX A WITH CONSECUTIVE C B1XB1 AND B2XB2 DIAGONAL BLOCKS (B1,B2 .LE. 2) C STARTING AT A(L,L), EXCHNG PRODUCES A UNITARY C SIMILARITY TRANSFORMATION THAT EXCHANGES THE BLOCKS C ALONG WITH THEIR EIGENVALUES. THE TRANSFORMATION C IS ACCUMULATED IN V. EXCHNG REQUIRES THE SUBROUTINE C QRSTEP. THE PARAMETERS IN THE CALLING SEQUENCE ARE C %STARRED PARAMETERS ARE ALTERED BY THE SUBROUTINE) C *A THE MATRIX WHOSE BLOCKS ARE TO BE C INTERCHANGED. C *V THE ARRAY INTO WHICH THE TRANSFORMATIONS C ARE TO BE ACCUMULATED. C N THE ORDER OF THE MATRIX A. C L THE POSITION OF THE BLOCKS. C B1 AN INTEGER CONTAINING THE SIZE OF THE C FIRST BLOCK. C B2 AN INTEGER CONTAINING THE SIZE OF THE C SECOND BLOCK. C EPS A CONVERGENCE CRITERION (CF. HQR3). C *FAIL A LOGICAL VARIABLE WHICH IS FALSE ON A C NORMAL RETURN. IF THIRTY ITERATIONS WERE C PERFORMED WITHOUT CONVERGENCE, FAIL IS SET C TO TRUE AND THE ELEMENT C A(L+B2,L+B2-1) CANNOT BE ASSUMED ZERO. C NA THE FIRST DIMENSION OF THE ARRAY A. C NV THE FIRST DIMENSION OF THE ARRAY V. C INTERNAL VARIABLES. INTEGER I, IT, J, L1, M REAL P, Q, R, S, W, X, Y, Z FAIL = .FALSE. IF (B1.EQ.2) GO TO 70 IF (B2.EQ.2) GO TO 40 C INTERCHANGE 1X1 AND 1X1 BLOCKS. L1 = L + 1 Q = A(L+1,L+1) - A(L,L) P = A(L,L+1) R = AMAX1(ABS(P),ABS(Q)) IF (R.EQ.0.) RETURN P = P/R Q = Q/R R = SQRT(P**2+Q**2) P = P/R Q = Q/R DO 10 J=L,N S = P*A(L,J) + Q*A(L+1,J) A(L+1,J) = P*A(L+1,J) - Q*A(L,J) A(L,J) = S 10 CONTINUE DO 20 I=1,L1 S = P*A(I,L) + Q*A(I,L+1) A(I,L+1) = P*A(I,L+1) - Q*A(I,L) A(I,L) = S 20 CONTINUE DO 30 I=1,N S = P*V(I,L) + Q*V(I,L+1) V(I,L+1) = P*V(I,L+1) - Q*V(I,L) V(I,L) = S 30 CONTINUE A(L+1,L) = 0. RETURN 40 CONTINUE C INTERCHANGE 1X1 AND 2X2 BLOCKS. X = A(L,L) P = 1. Q = 1. R = 1. CALL QRSTEP(A, V, P, Q, R, L, L+2, N, NA, NV) IT = 0 50 IT = IT + 1 IF (IT.LE.30) GO TO 60 FAIL = .TRUE. RETURN 60 CONTINUE P = A(L,L) - X Q = A(L+1,L) R = 0. CALL QRSTEP(A, V, P, Q, R, L, L+2, N, NA, NV) IF (ABS(A(L+2,L+1)).GT.EPS*(ABS(A(L+1,L+1))+ABS(A(L+2,L+2)))) * GO TO 50 A(L+2,L+1) = 0. RETURN 70 CONTINUE C INTERCHANGE 2X2 AND B2XB2 BLOCKS. M = L + 2 IF (B2.EQ.2) M = M + 1 X = A(L+1,L+1) Y = A(L,L) W = A(L+1,L)*A(L,L+1) P = 1. Q = 1. R = 1. CALL QRSTEP(A, V, P, Q, R, L, M, N, NA, NV) IT = 0 80 IT = IT + 1 IF (IT.LE.30) GO TO 90 FAIL = .TRUE. RETURN 90 CONTINUE Z = A(L,L) R = X - Z S = Y - Z P = (R*S-W)/A(L+1,L) + A(L,L+1) Q = A(L+1,L+1) - Z - R - S R = A(L+2,L+1) S = ABS(P) + ABS(Q) + ABS(R) P = P/S Q = Q/S R = R/S CALL QRSTEP(A, V, P, Q, R, L, M, N, NA, NV) IF (ABS(A(M-1,M-2)).GT.EPS*(ABS(A(M-1,M-1))+ABS(A(M-2,M-2)))) * GO TO 80 A(M-1,M-2) = 0. RETURN END SUBROUTINE QRSTEP(A, V, P, Q, R, NL, NU, N, NA, NV) QRS 10 INTEGER N, NA, NL, NU, NV REAL A(NA,N), P, Q, R, V(NV,N) C QRSTEP PERFORMS ONE IMPLICIT QR STEP ON THE C UPPER HESSENBERG MATRIX A. THE SHIFT IS DETERMINED C BY THE NUMBERS P,Q, AND R, AND THE STEP IS APPLIED TO C ROWS AND COLUMNS NL THROUGH NU. THE TRANSFORMATIONS C ARE ACCUMULATED IN V. THE PARAMETERS IN THE CALLING C SEQUENCE ARE (STARRED APRAMETERS ARE ALTERED BY THE C SUBROUTINE) C *A THE UPPER HESSENBERG MATRIX ON WHICH THE C QR STEP IS TO BE PERFORMED. C *V THE ARRAY IN WHICH THE TRANSFORMATIONS C ARE TO BE ACCUMULATED C *P PARAMETERS THAT DETERMINE THE SHIFT. C *Q C *R C NL THE LOWER LIMIT OF THE STEP. C NU THE UPPER LIMIT OF THE STEP. C N THE ORDER OF THE MATRIX A. C NA THE FIRST DIMENSION OF THE ARRAY A. C NV THE FIRST DIMENSION OF THE ARRAY V. C INTERNAL VARIABLES. INTEGER I, J, K, NL2, NL3, NUM1 REAL S, X, Y, Z LOGICAL LAST NL2 = NL + 2 DO 10 I=NL2,NU A(I,I-2) = 0. 10 CONTINUE IF (NL2.EQ.NU) GO TO 30 NL3 = NL + 3 DO 20 I=NL3,NU A(I,I-3) = 0. 20 CONTINUE 30 CONTINUE NUM1 = NU - 1 DO 130 K=NL,NUM1 C DETERMINE THE TRANSFORMATION. LAST = K.EQ.NUM1 IF (K.EQ.NL) GO TO 40 P = A(K,K-1) Q = A(K+1,K-1) R = 0. IF (.NOT.LAST) R = A(K+2,K-1) X = ABS(P) + ABS(Q) + ABS(R) IF (X.EQ.0.) GO TO 130 P = P/X Q = Q/X R = R/X 40 CONTINUE S = SQRT(P**2+Q**2+R**2) IF (P.LT.0.) S = -S IF (K.EQ.NL) GO TO 50 A(K,K-1) = -S*X GO TO 60 50 CONTINUE IF (NL.NE.1) A(K,K-1) = -A(K,K-1) 60 CONTINUE P = P + S X = P/S Y = Q/S Z = R/S Q = Q/P R = R/P C PREMULTIPLY. DO 80 J=K,N P = A(K,J) + Q*A(K+1,J) IF (LAST) GO TO 70 P = P + R*A(K+2,J) A(K+2,J) = A(K+2,J) - P*Z 70 CONTINUE A(K+1,J) = A(K+1,J) - P*Y A(K,J) = A(K,J) - P*X 80 CONTINUE C POSTMULTIPLY. J = MIN0(K+3,NU) DO 100 I=1,J P = X*A(I,K) + Y*A(I,K+1) IF (LAST) GO TO 90 P = P + Z*A(I,K+2) A(I,K+2) = A(I,K+2) - P*R 90 CONTINUE A(I,K+1) = A(I,K+1) - P*Q A(I,K) = A(I,K) - P 100 CONTINUE C ACCUMULATE THE TRANSFORMATION IN V. DO 120 I=1,N P = X*V(I,K) + Y*V(I,K+1) IF (LAST) GO TO 110 P = P + Z*V(I,K+2) V(I,K+2) = V(I,K+2) - P*R 110 CONTINUE V(I,K+1) = V(I,K+1) - P*Q V(I,K) = V(I,K) - P 120 CONTINUE 130 CONTINUE RETURN END