This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / 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