MINICERN/mathlib/gen/f/hqr.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:02:35  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE HQR(NM,N,LOW,IGH,H,WR,WI,IERR)
  11       INTEGER I,J,K,L,M,N,EN,LL,MM,NA,NM,IGH,ITS,LOW,MP2,ENM2,IERR
  12       REAL H(NM,N),WR(N),WI(N)
  13       REAL P,Q,R,S,T,W,X,Y,ZZ,MACHEP
  14       LOGICAL NOTLAS
  15 #if defined(CERNLIB_CDC)
  16       MACHEP=2.**(-47)
  17 #endif
  18 #if !defined(CERNLIB_CDC)
  19       MACHEP=2.**(-23)
  20 #endif
  21       IERR = 0
  22       DO 50 I = 1, N
  23          IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50
  24          WR(I) = H(I,I)
  25          WI(I) = 0.0
  26    50 CONTINUE
  27       EN = IGH
  28       T = 0.0
  29    60 IF (EN .LT. LOW) GO TO 1001
  30       ITS = 0
  31       NA = EN - 1
  32       ENM2 = NA - 1
  33    70 DO 80 LL = LOW, EN
  34          L = EN + LOW - LL
  35          IF (L .EQ. LOW) GO TO 100
  36          IF (ABS(H(L,L-1)) .LE. MACHEP * (ABS(H(L-1,L-1))
  37      X      + ABS(H(L,L)))) GO TO 100
  38    80 CONTINUE
  39   100 X = H(EN,EN)
  40       IF (L .EQ. EN) GO TO 270
  41       Y = H(NA,NA)
  42       W = H(EN,NA) * H(NA,EN)
  43       IF (L .EQ. NA) GO TO 280
  44       IF (ITS .EQ. 30) GO TO 1000
  45       IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130
  46       T = T + X
  47       DO 120 I = LOW, EN
  48   120 H(I,I) = H(I,I) - X
  49       S = ABS(H(EN,NA)) + ABS(H(NA,ENM2))
  50       X = 0.75 * S
  51       Y = X
  52       W = -0.4375 * S * S
  53   130 ITS = ITS + 1
  54       DO 140 MM = L, ENM2
  55          M = ENM2 + L - MM
  56          ZZ = H(M,M)
  57          R = X - ZZ
  58          S = Y - ZZ
  59          P = (R * S - W) / H(M+1,M) + H(M,M+1)
  60          Q = H(M+1,M+1) - ZZ - R - S
  61          R = H(M+2,M+1)
  62          S = ABS(P) + ABS(Q) + ABS(R)
  63          P = P / S
  64          Q = Q / S
  65          R = R / S
  66          IF (M .EQ. L) GO TO 150
  67          IF (ABS(H(M,M-1)) * (ABS(Q) + ABS(R)) .LE. MACHEP * ABS(P)
  68      X    * (ABS(H(M-1,M-1)) + ABS(ZZ) + ABS(H(M+1,M+1)))) GO TO 150
  69   140 CONTINUE
  70   150 MP2 = M + 2
  71       DO 160 I = MP2, EN
  72          H(I,I-2) = 0.0
  73          IF (I .EQ. MP2) GO TO 160
  74          H(I,I-3) = 0.0
  75   160 CONTINUE
  76       DO 260 K = M, NA
  77          NOTLAS = K .NE. NA
  78          IF (K .EQ. M) GO TO 170
  79          P = H(K,K-1)
  80          Q = H(K+1,K-1)
  81          R = 0.0
  82          IF (NOTLAS) R = H(K+2,K-1)
  83          X = ABS(P) + ABS(Q) + ABS(R)
  84          IF (X .EQ. 0.0) GO TO 260
  85          P = P / X
  86          Q = Q / X
  87          R = R / X
  88   170    S = SIGN(SQRT(P*P+Q*Q+R*R),P)
  89          IF (K .EQ. M) GO TO 180
  90          H(K,K-1) = -S * X
  91          GO TO 190
  92   180    IF (L .NE. M) H(K,K-1) = -H(K,K-1)
  93   190    P = P + S
  94          X = P / S
  95          Y = Q / S
  96          ZZ = R / S
  97          Q = Q / P
  98          R = R / P
  99          DO 210 J = K, EN
 100             P = H(K,J) + Q * H(K+1,J)
 101             IF (.NOT. NOTLAS) GO TO 200
 102             P = P + R * H(K+2,J)
 103             H(K+2,J) = H(K+2,J) - P * ZZ
 104   200       H(K+1,J) = H(K+1,J) - P * Y
 105             H(K,J) = H(K,J) - P * X
 106   210    CONTINUE
 107          J = MIN(EN,K+3)
 108          DO 230 I = L, J
 109             P = X * H(I,K) + Y * H(I,K+1)
 110             IF (.NOT. NOTLAS) GO TO 220
 111             P = P + ZZ * H(I,K+2)
 112             H(I,K+2) = H(I,K+2) - P * R
 113   220       H(I,K+1) = H(I,K+1) - P * Q
 114             H(I,K) = H(I,K) - P
 115   230    CONTINUE
 116   260 CONTINUE
 117       GO TO 70
 118   270 WR(EN) = X + T
 119       WI(EN) = 0.0
 120       EN = NA
 121       GO TO 60
 122   280 P = (Y - X) / 2.0
 123       Q = P * P + W
 124       ZZ = SQRT(ABS(Q))
 125       X = X + T
 126       IF (Q .LT. 0.0) GO TO 320
 127       ZZ = P + SIGN(ZZ,P)
 128       WR(NA) = X + ZZ
 129       WR(EN) = WR(NA)
 130       IF (ZZ .NE. 0.0) WR(EN) = X - W / ZZ
 131       WI(NA) = 0.0
 132       WI(EN) = 0.0
 133       GO TO 330
 134   320 WR(NA) = X + P
 135       WR(EN) = X + P
 136       WI(NA) = ZZ
 137       WI(EN) = -ZZ
 138   330 EN = ENM2
 139       GO TO 60
 140  1000 IERR = EN
 141  1001 RETURN
 142       END