*DECK HQR SUBROUTINE HQR (NM, N, LOW, IGH, H, WR, WI, IERR) C***BEGIN PROLOGUE HQR C***PURPOSE Compute the eigenvalues of a real upper Hessenberg matrix C using the QR method. C***LIBRARY SLATEC (EISPACK) C***CATEGORY D4C2B C***TYPE SINGLE PRECISION (HQR-S, COMQR-C) C***KEYWORDS EIGENVALUES, EIGENVECTORS, EISPACK C***AUTHOR Smith, B. T., et al. C***DESCRIPTION C C This subroutine is a translation of the ALGOL procedure HQR, C NUM. MATH. 14, 219-231(1970) by Martin, Peters, and Wilkinson. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971). C C This subroutine finds the eigenvalues of a REAL C UPPER Hessenberg matrix by the QR method. C C On INPUT C C NM must be set to the row dimension of the two-dimensional C array parameter, H, as declared in the calling program C dimension statement. NM is an INTEGER variable. C C N is the order of the matrix H. N is an INTEGER variable. C N must be less than or equal to NM. C C LOW and IGH are two INTEGER variables determined by the C balancing subroutine BALANC. If BALANC has not been C used, set LOW=1 and IGH equal to the order of the matrix, N. C C H contains the upper Hessenberg matrix. Information about C the transformations used in the reduction to Hessenberg C form by ELMHES or ORTHES, if performed, is stored C in the remaining triangle under the Hessenberg matrix. C H is a two-dimensional REAL array, dimensioned H(NM,N). C C On OUTPUT C C H has been destroyed. Therefore, it must be saved before C calling HQR if subsequent calculation and back C transformation of eigenvectors is to be performed. C C WR and WI contain the real and imaginary parts, respectively, C of the eigenvalues. The eigenvalues are unordered except C that complex conjugate pairs of values appear consecutively C with the eigenvalue having the positive imaginary part first. C If an error exit is made, the eigenvalues should be correct C for indices IERR+1, IERR+2, ..., N. WR and WI are one- C dimensional REAL arrays, dimensioned WR(N) and WI(N). C C IERR is an INTEGER flag set to C Zero for normal return, C J if the J-th eigenvalue has not been C determined after a total of 30*N iterations. C The eigenvalues should be correct for indices C IERR+1, IERR+2, ..., N. C C Questions and comments should be directed to B. S. Garbow, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C ------------------------------------------------------------------ C C***REFERENCES B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, C Y. Ikebe, V. C. Klema and C. B. Moler, Matrix Eigen- C system Routines - EISPACK Guide, Springer-Verlag, C 1976. C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 760101 DATE WRITTEN C 890531 Changed all specific intrinsics to generic. (WRB) C 890831 Modified array declarations. (WRB) C 890831 REVISION DATE from Version 3.2 C 891214 Prologue converted to Version 4.0 format. (BAB) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE HQR C INTEGER I,J,K,L,M,N,EN,LL,MM,NA,NM,IGH,ITN,ITS,LOW,MP2,ENM2,IERR REAL H(NM,*),WR(*),WI(*) REAL P,Q,R,S,T,W,X,Y,ZZ,NORM,S1,S2 LOGICAL NOTLAS C C***FIRST EXECUTABLE STATEMENT HQR IERR = 0 NORM = 0.0E0 K = 1 C .......... STORE ROOTS ISOLATED BY BALANC C AND COMPUTE MATRIX NORM .......... DO 50 I = 1, N C DO 40 J = K, N 40 NORM = NORM + ABS(H(I,J)) C K = I IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50 WR(I) = H(I,I) WI(I) = 0.0E0 50 CONTINUE C EN = IGH T = 0.0E0 ITN = 30*N C .......... SEARCH FOR NEXT EIGENVALUES .......... 60 IF (EN .LT. LOW) GO TO 1001 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 -- .......... 70 DO 80 LL = LOW, EN L = EN + LOW - LL IF (L .EQ. LOW) GO TO 100 S = ABS(H(L-1,L-1)) + ABS(H(L,L)) IF (S .EQ. 0.0E0) S = NORM S2 = S + ABS(H(L,L-1)) IF (S2 .EQ. S) GO TO 100 80 CONTINUE C .......... FORM SHIFT .......... 100 X = H(EN,EN) IF (L .EQ. EN) GO TO 270 Y = H(NA,NA) W = H(EN,NA) * H(NA,EN) IF (L .EQ. NA) GO TO 280 IF (ITN .EQ. 0) GO TO 1000 IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130 C .......... FORM EXCEPTIONAL SHIFT .......... T = T + X C DO 120 I = LOW, EN 120 H(I,I) = H(I,I) - X C S = ABS(H(EN,NA)) + ABS(H(NA,ENM2)) X = 0.75E0 * S Y = X W = -0.4375E0 * S * S 130 ITS = ITS + 1 ITN = ITN - 1 C .......... LOOK FOR TWO CONSECUTIVE SMALL C SUB-DIAGONAL ELEMENTS. C FOR M=EN-2 STEP -1 UNTIL L DO -- .......... DO 140 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 150 S1 = ABS(P) * (ABS(H(M-1,M-1)) + ABS(ZZ) + ABS(H(M+1,M+1))) S2 = S1 + ABS(H(M,M-1)) * (ABS(Q) + ABS(R)) IF (S2 .EQ. S1) GO TO 150 140 CONTINUE C 150 MP2 = M + 2 C DO 160 I = MP2, EN H(I,I-2) = 0.0E0 IF (I .EQ. MP2) GO TO 160 H(I,I-3) = 0.0E0 160 CONTINUE C .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND C COLUMNS M TO EN .......... DO 260 K = M, NA NOTLAS = K .NE. NA IF (K .EQ. M) GO TO 170 P = H(K,K-1) Q = H(K+1,K-1) R = 0.0E0 IF (NOTLAS) R = H(K+2,K-1) X = ABS(P) + ABS(Q) + ABS(R) IF (X .EQ. 0.0E0) GO TO 260 P = P / X Q = Q / X R = R / X 170 S = SIGN(SQRT(P*P+Q*Q+R*R),P) IF (K .EQ. M) GO TO 180 H(K,K-1) = -S * X GO TO 190 180 IF (L .NE. M) H(K,K-1) = -H(K,K-1) 190 P = P + S X = P / S Y = Q / S ZZ = R / S Q = Q / P R = R / P C .......... ROW MODIFICATION .......... DO 210 J = K, EN P = H(K,J) + Q * H(K+1,J) IF (.NOT. NOTLAS) GO TO 200 P = P + R * H(K+2,J) H(K+2,J) = H(K+2,J) - P * ZZ 200 H(K+1,J) = H(K+1,J) - P * Y H(K,J) = H(K,J) - P * X 210 CONTINUE C J = MIN(EN,K+3) C .......... COLUMN MODIFICATION .......... DO 230 I = L, J P = X * H(I,K) + Y * H(I,K+1) IF (.NOT. NOTLAS) GO TO 220 P = P + ZZ * H(I,K+2) H(I,K+2) = H(I,K+2) - P * R 220 H(I,K+1) = H(I,K+1) - P * Q H(I,K) = H(I,K) - P 230 CONTINUE C 260 CONTINUE C GO TO 70 C .......... ONE ROOT FOUND .......... 270 WR(EN) = X + T WI(EN) = 0.0E0 EN = NA GO TO 60 C .......... TWO ROOTS FOUND .......... 280 P = (Y - X) / 2.0E0 Q = P * P + W ZZ = SQRT(ABS(Q)) X = X + T IF (Q .LT. 0.0E0) GO TO 320 C .......... REAL PAIR .......... ZZ = P + SIGN(ZZ,P) WR(NA) = X + ZZ WR(EN) = WR(NA) IF (ZZ .NE. 0.0E0) WR(EN) = X - W / ZZ WI(NA) = 0.0E0 WI(EN) = 0.0E0 GO TO 330 C .......... COMPLEX PAIR .......... 320 WR(NA) = X + P WR(EN) = X + P WI(NA) = ZZ WI(EN) = -ZZ 330 EN = ENM2 GO TO 60 C .......... SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30*N ITERATIONS .......... 1000 IERR = EN 1001 RETURN END