* * $Id: hqr2.F,v 1.1.1.1 1996/04/01 15:02:35 mclareni Exp $ * * $Log: hqr2.F,v $ * Revision 1.1.1.1 1996/04/01 15:02:35 mclareni * Mathlib gen * * #include "gen/pilot.h" SUBROUTINE HQR2(NM,N,LOW,IGH,H,WR,WI,Z,IERR) INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NA,NM,NN, X IGH,ITS,LOW,MP2,ENM2,IERR REAL H(NM,N),WR(N),WI(N),Z(NM,N) REAL P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,MACHEP LOGICAL NOTLAS COMPLEX Z3 REAL T3(2) EQUIVALENCE (Z3,T3(1)) #if defined(CERNLIB_CDC) MACHEP=2.**(-47) #endif #if !defined(CERNLIB_CDC) MACHEP=2.**(-23) #endif IERR = 0 DO 50 I = 1, N IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50 WR(I) = H(I,I) WI(I) = 0.0 50 CONTINUE EN = IGH T = 0.0 60 IF (EN .LT. LOW) GO TO 340 ITS = 0 NA = EN - 1 ENM2 = NA - 1 70 DO 80 LL = LOW, EN L = EN + LOW - LL IF (L .EQ. LOW) GO TO 100 IF (ABS(H(L,L-1)) .LE. MACHEP * (ABS(H(L-1,L-1)) X + ABS(H(L,L)))) GO TO 100 80 CONTINUE 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 (ITS .EQ. 30) GO TO 1000 IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130 T = T + X DO 120 I = LOW, EN 120 H(I,I) = H(I,I) - X S = ABS(H(EN,NA)) + ABS(H(NA,ENM2)) X = 0.75 * S Y = X W = -0.4375 * S * S 130 ITS = ITS + 1 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 IF (ABS(H(M,M-1)) * (ABS(Q) + ABS(R)) .LE. MACHEP * ABS(P) X * (ABS(H(M-1,M-1)) + ABS(ZZ) + ABS(H(M+1,M+1)))) GO TO 150 140 CONTINUE 150 MP2 = M + 2 DO 160 I = MP2, EN H(I,I-2) = 0.0 IF (I .EQ. MP2) GO TO 160 H(I,I-3) = 0.0 160 CONTINUE 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.0 IF (NOTLAS) R = H(K+2,K-1) X = ABS(P) + ABS(Q) + ABS(R) IF (X .EQ. 0.0) 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 DO 210 J = K, N 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 J = MIN(EN,K+3) DO 230 I = 1, 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 DO 250 I = LOW, IGH P = X * Z(I,K) + Y * Z(I,K+1) IF (.NOT. NOTLAS) GO TO 240 P = P + ZZ * Z(I,K+2) Z(I,K+2) = Z(I,K+2) - P * R 240 Z(I,K+1) = Z(I,K+1) - P * Q Z(I,K) = Z(I,K) - P 250 CONTINUE 260 CONTINUE GO TO 70 270 H(EN,EN) = X + T WR(EN) = H(EN,EN) WI(EN) = 0.0 EN = NA GO TO 60 280 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 320 ZZ = P + SIGN(ZZ,P) 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 DO 290 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 290 CONTINUE DO 300 I = 1, EN ZZ = H(I,NA) H(I,NA) = Q * ZZ + P * H(I,EN) H(I,EN) = Q * H(I,EN) - P * ZZ 300 CONTINUE DO 310 I = LOW, IGH ZZ = Z(I,NA) Z(I,NA) = Q * ZZ + P * Z(I,EN) Z(I,EN) = Q * Z(I,EN) - P * ZZ 310 CONTINUE GO TO 330 320 WR(NA) = X + P WR(EN) = X + P WI(NA) = ZZ WI(EN) = -ZZ 330 EN = ENM2 GO TO 60 340 NORM = 0.0 K = 1 DO 360 I = 1, N DO 350 J = K, N 350 NORM = NORM + ABS(H(I,J)) K = I 360 CONTINUE IF (NORM .EQ. 0.0) GO TO 1001 DO 800 NN = 1, N EN = N + 1 - NN P = WR(EN) Q = WI(EN) NA = EN - 1 IF (Q) 710, 600, 800 600 M = EN H(EN,EN) = 1.0 IF (NA .EQ. 0) GO TO 800 DO 700 II = 1, NA I = EN - II W = H(I,I) - P R = H(I,EN) IF (M .GT. NA) GO TO 620 DO 610 J = M, NA 610 R = R + H(I,J) * H(J,EN) 620 IF (WI(I) .GE. 0.0) GO TO 630 ZZ = W S = R GO TO 700 630 M = I IF (WI(I) .NE. 0.0) GO TO 640 T = W IF (W .EQ. 0.0) T = MACHEP * NORM H(I,EN) = -R / T GO TO 700 640 X = H(I,I+1) Y = H(I+1,I) Q = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) T = (X * S - ZZ * R) / Q H(I,EN) = T IF (ABS(X) .LE. ABS(ZZ)) GO TO 650 H(I+1,EN) = (-R - W * T) / X GO TO 700 650 H(I+1,EN) = (-S - Y * T) / ZZ 700 CONTINUE GO TO 800 710 M = NA IF (ABS(H(EN,NA)) .LE. ABS(H(NA,EN))) GO TO 720 H(NA,NA) = Q / H(EN,NA) H(NA,EN) = -(H(EN,EN) - P) / H(EN,NA) GO TO 730 720 Z3 = CMPLX(0.0,-H(NA,EN)) / CMPLX(H(NA,NA)-P,Q) H(NA,NA) = T3(1) H(NA,EN) = T3(2) 730 H(EN,NA) = 0.0 H(EN,EN) = 1.0 ENM2 = NA - 1 IF (ENM2 .EQ. 0) GO TO 800 DO 790 II = 1, ENM2 I = NA - II W = H(I,I) - P RA = 0.0 SA = H(I,EN) DO 760 J = M, NA RA = RA + H(I,J) * H(J,NA) SA = SA + H(I,J) * H(J,EN) 760 CONTINUE IF (WI(I) .GE. 0.0) GO TO 770 ZZ = W R = RA S = SA GO TO 790 770 M = I IF (WI(I) .NE. 0.0) GO TO 780 Z3 = CMPLX(-RA,-SA) / CMPLX(W,Q) H(I,NA) = T3(1) H(I,EN) = T3(2) GO TO 790 780 X = H(I,I+1) Y = H(I+1,I) VR = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) - Q * Q VI = (WR(I) - P) * 2.0 * Q IF (VR .EQ. 0.0 .AND. VI .EQ. 0.0) VR = MACHEP * NORM X * (ABS(W) + ABS(Q) + ABS(X) + ABS(Y) + ABS(ZZ)) Z3 = CMPLX(X*R-ZZ*RA+Q*SA,X*S-ZZ*SA-Q*RA) / CMPLX(VR,VI) H(I,NA) = T3(1) H(I,EN) = T3(2) IF (ABS(X) .LE. ABS(ZZ) + ABS(Q)) GO TO 785 H(I+1,NA) = (-RA - W * H(I,NA) + Q * H(I,EN)) / X H(I+1,EN) = (-SA - W * H(I,EN) - Q * H(I,NA)) / X GO TO 790 785 Z3 = CMPLX(-R-Y*H(I,NA),-S-Y*H(I,EN)) / CMPLX(ZZ,Q) H(I+1,NA) = T3(1) H(I+1,EN) = T3(2) 790 CONTINUE 800 CONTINUE DO 840 I = 1, N IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840 DO 820 J = I, N 820 Z(I,J) = H(I,J) 840 CONTINUE DO 880 JJ = LOW, N J = N + LOW - JJ M = MIN(J,IGH) DO 880 I = LOW, IGH ZZ = 0.0 DO 860 K = LOW, M 860 ZZ = ZZ + Z(I,K) * H(K,J) Z(I,J) = ZZ 880 CONTINUE GO TO 1001 1000 IERR = EN 1001 RETURN END