[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / f / hqr2.F

*
* $Id$
*
* $Log$
* 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
Commit	Line	Data
fe4da5cc	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 HQR2(NM,N,LOW,IGH,H,WR,WI,Z,IERR)
	11	INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NA,NM,NN,
	12	X IGH,ITS,LOW,MP2,ENM2,IERR
	13	REAL H(NM,N),WR(N),WI(N),Z(NM,N)
	14	REAL P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,MACHEP
	15	LOGICAL NOTLAS
	16	COMPLEX Z3
	17	REAL T3(2)
	18	EQUIVALENCE (Z3,T3(1))
	19	#if defined(CERNLIB_CDC)
	20	MACHEP=2.**(-47)
	21	#endif
	22	#if !defined(CERNLIB_CDC)
	23	MACHEP=2.**(-23)
	24	#endif
	25	IERR = 0
	26	DO 50 I = 1, N
	27	IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50
	28	WR(I) = H(I,I)
	29	WI(I) = 0.0
	30	50 CONTINUE
	31	EN = IGH
	32	T = 0.0
	33	60 IF (EN .LT. LOW) GO TO 340
	34	ITS = 0
	35	NA = EN - 1
	36	ENM2 = NA - 1
	37	70 DO 80 LL = LOW, EN
	38	L = EN + LOW - LL
	39	IF (L .EQ. LOW) GO TO 100
	40	IF (ABS(H(L,L-1)) .LE. MACHEP * (ABS(H(L-1,L-1))
	41	X + ABS(H(L,L)))) GO TO 100
	42	80 CONTINUE
	43	100 X = H(EN,EN)
	44	IF (L .EQ. EN) GO TO 270
	45	Y = H(NA,NA)
	46	W = H(EN,NA) * H(NA,EN)
	47	IF (L .EQ. NA) GO TO 280
	48	IF (ITS .EQ. 30) GO TO 1000
	49	IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130
	50	T = T + X
	51	DO 120 I = LOW, EN
	52	120 H(I,I) = H(I,I) - X
	53	S = ABS(H(EN,NA)) + ABS(H(NA,ENM2))
	54	X = 0.75 * S
	55	Y = X
	56	W = -0.4375 * S * S
	57	130 ITS = ITS + 1
	58	DO 140 MM = L, ENM2
	59	M = ENM2 + L - MM
	60	ZZ = H(M,M)
	61	R = X - ZZ
	62	S = Y - ZZ
	63	P = (R * S - W) / H(M+1,M) + H(M,M+1)
	64	Q = H(M+1,M+1) - ZZ - R - S
65	R = H(M+2,M+1)
66	S = ABS(P) + ABS(Q) + ABS(R)
67	P = P / S
68	Q = Q / S
69	R = R / S
70	IF (M .EQ. L) GO TO 150
71	IF (ABS(H(M,M-1)) * (ABS(Q) + ABS(R)) .LE. MACHEP * ABS(P)
72	X * (ABS(H(M-1,M-1)) + ABS(ZZ) + ABS(H(M+1,M+1)))) GO TO 150
73	140 CONTINUE
74	150 MP2 = M + 2
75	DO 160 I = MP2, EN
76	H(I,I-2) = 0.0
77	IF (I .EQ. MP2) GO TO 160
78	H(I,I-3) = 0.0
79	160 CONTINUE
80	DO 260 K = M, NA
81	NOTLAS = K .NE. NA
82	IF (K .EQ. M) GO TO 170
83	P = H(K,K-1)
84	Q = H(K+1,K-1)
85	R = 0.0
86	IF (NOTLAS) R = H(K+2,K-1)
87	X = ABS(P) + ABS(Q) + ABS(R)
88	IF (X .EQ. 0.0) GO TO 260
89	P = P / X
90	Q = Q / X
91	R = R / X
92	170 S = SIGN(SQRT(PP+QQ+R*R),P)
93	IF (K .EQ. M) GO TO 180
94	H(K,K-1) = -S * X
95	GO TO 190
96	180 IF (L .NE. M) H(K,K-1) = -H(K,K-1)
97	190 P = P + S
98	X = P / S
99	Y = Q / S
100	ZZ = R / S
101	Q = Q / P
102	R = R / P
103	DO 210 J = K, N
104	P = H(K,J) + Q * H(K+1,J)
105	IF (.NOT. NOTLAS) GO TO 200
106	P = P + R * H(K+2,J)
107	H(K+2,J) = H(K+2,J) - P * ZZ
108	200 H(K+1,J) = H(K+1,J) - P * Y
109	H(K,J) = H(K,J) - P * X
110	210 CONTINUE
111	J = MIN(EN,K+3)
112	DO 230 I = 1, J
113	P = X * H(I,K) + Y * H(I,K+1)
114	IF (.NOT. NOTLAS) GO TO 220
115	P = P + ZZ * H(I,K+2)
116	H(I,K+2) = H(I,K+2) - P * R
117	220 H(I,K+1) = H(I,K+1) - P * Q
118	H(I,K) = H(I,K) - P
119	230 CONTINUE
120	DO 250 I = LOW, IGH
121	P = X * Z(I,K) + Y * Z(I,K+1)
122	IF (.NOT. NOTLAS) GO TO 240
123	P = P + ZZ * Z(I,K+2)
124	Z(I,K+2) = Z(I,K+2) - P * R
125	240 Z(I,K+1) = Z(I,K+1) - P * Q
126	Z(I,K) = Z(I,K) - P
127	250 CONTINUE
128	260 CONTINUE
129	GO TO 70
130	270 H(EN,EN) = X + T
131	WR(EN) = H(EN,EN)
132	WI(EN) = 0.0
133	EN = NA
134	GO TO 60
135	280 P = (Y - X) / 2.0
136	Q = P * P + W
137	ZZ = SQRT(ABS(Q))
138	H(EN,EN) = X + T
139	X = H(EN,EN)
140	H(NA,NA) = Y + T
141	IF (Q .LT. 0.0) GO TO 320
142	ZZ = P + SIGN(ZZ,P)
143	WR(NA) = X + ZZ
144	WR(EN) = WR(NA)
145	IF (ZZ .NE. 0.0) WR(EN) = X - W / ZZ
146	WI(NA) = 0.0
147	WI(EN) = 0.0
148	X = H(EN,NA)
149	R = SQRT(XX+ZZZZ)
150	P = X / R
151	Q = ZZ / R
152	DO 290 J = NA, N
153	ZZ = H(NA,J)
154	H(NA,J) = Q * ZZ + P * H(EN,J)
155	H(EN,J) = Q * H(EN,J) - P * ZZ
156	290 CONTINUE
157	DO 300 I = 1, EN
158	ZZ = H(I,NA)
159	H(I,NA) = Q * ZZ + P * H(I,EN)
160	H(I,EN) = Q * H(I,EN) - P * ZZ
161	300 CONTINUE
162	DO 310 I = LOW, IGH
163	ZZ = Z(I,NA)
164	Z(I,NA) = Q * ZZ + P * Z(I,EN)
165	Z(I,EN) = Q * Z(I,EN) - P * ZZ
166	310 CONTINUE
167	GO TO 330
168	320 WR(NA) = X + P
169	WR(EN) = X + P
170	WI(NA) = ZZ
171	WI(EN) = -ZZ
172	330 EN = ENM2
173	GO TO 60
174	340 NORM = 0.0
175	K = 1
176	DO 360 I = 1, N
177	DO 350 J = K, N
178	350 NORM = NORM + ABS(H(I,J))
179	K = I
180	360 CONTINUE
181	IF (NORM .EQ. 0.0) GO TO 1001
182	DO 800 NN = 1, N
183	EN = N + 1 - NN
184	P = WR(EN)
185	Q = WI(EN)
186	NA = EN - 1
187	IF (Q) 710, 600, 800
188	600 M = EN
189	H(EN,EN) = 1.0
190	IF (NA .EQ. 0) GO TO 800
191	DO 700 II = 1, NA
192	I = EN - II
193	W = H(I,I) - P
194	R = H(I,EN)
195	IF (M .GT. NA) GO TO 620
196	DO 610 J = M, NA
197	610 R = R + H(I,J) * H(J,EN)
198	620 IF (WI(I) .GE. 0.0) GO TO 630
199	ZZ = W
200	S = R
201	GO TO 700
202	630 M = I
203	IF (WI(I) .NE. 0.0) GO TO 640
204	T = W
205	IF (W .EQ. 0.0) T = MACHEP * NORM
206	H(I,EN) = -R / T
207	GO TO 700
208	640 X = H(I,I+1)
209	Y = H(I+1,I)
210	Q = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I)
211	T = (X * S - ZZ * R) / Q
212	H(I,EN) = T
213	IF (ABS(X) .LE. ABS(ZZ)) GO TO 650
214	H(I+1,EN) = (-R - W * T) / X
215	GO TO 700
216	650 H(I+1,EN) = (-S - Y * T) / ZZ
217	700 CONTINUE
218	GO TO 800
219	710 M = NA
220	IF (ABS(H(EN,NA)) .LE. ABS(H(NA,EN))) GO TO 720
221	H(NA,NA) = Q / H(EN,NA)
222	H(NA,EN) = -(H(EN,EN) - P) / H(EN,NA)
223	GO TO 730
224	720 Z3 = CMPLX(0.0,-H(NA,EN)) / CMPLX(H(NA,NA)-P,Q)
225	H(NA,NA) = T3(1)
226	H(NA,EN) = T3(2)
227	730 H(EN,NA) = 0.0
228	H(EN,EN) = 1.0
229	ENM2 = NA - 1
230	IF (ENM2 .EQ. 0) GO TO 800
231	DO 790 II = 1, ENM2
232	I = NA - II
233	W = H(I,I) - P
234	RA = 0.0
235	SA = H(I,EN)
236	DO 760 J = M, NA
237	RA = RA + H(I,J) * H(J,NA)
238	SA = SA + H(I,J) * H(J,EN)
239	760 CONTINUE
240	IF (WI(I) .GE. 0.0) GO TO 770
241	ZZ = W
242	R = RA
243	S = SA
244	GO TO 790
245	770 M = I
246	IF (WI(I) .NE. 0.0) GO TO 780
247	Z3 = CMPLX(-RA,-SA) / CMPLX(W,Q)
248	H(I,NA) = T3(1)
249	H(I,EN) = T3(2)
250	GO TO 790
251	780 X = H(I,I+1)
252	Y = H(I+1,I)
253	VR = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) - Q * Q
254	VI = (WR(I) - P) * 2.0 * Q
255	IF (VR .EQ. 0.0 .AND. VI .EQ. 0.0) VR = MACHEP * NORM
256	X * (ABS(W) + ABS(Q) + ABS(X) + ABS(Y) + ABS(ZZ))
257	Z3 = CMPLX(XR-ZZRA+QSA,XS-ZZSA-QRA) / CMPLX(VR,VI)
258	H(I,NA) = T3(1)
259	H(I,EN) = T3(2)
260	IF (ABS(X) .LE. ABS(ZZ) + ABS(Q)) GO TO 785
261	H(I+1,NA) = (-RA - W * H(I,NA) + Q * H(I,EN)) / X
262	H(I+1,EN) = (-SA - W * H(I,EN) - Q * H(I,NA)) / X
263	GO TO 790
264	785 Z3 = CMPLX(-R-YH(I,NA),-S-YH(I,EN)) / CMPLX(ZZ,Q)
265	H(I+1,NA) = T3(1)
266	H(I+1,EN) = T3(2)
267	790 CONTINUE
268	800 CONTINUE
269	DO 840 I = 1, N
270	IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840
271	DO 820 J = I, N
272	820 Z(I,J) = H(I,J)
273	840 CONTINUE
274	DO 880 JJ = LOW, N
275	J = N + LOW - JJ
276	M = MIN(J,IGH)
277	DO 880 I = LOW, IGH
278	ZZ = 0.0
279	DO 860 K = LOW, M
280	860 ZZ = ZZ + Z(I,K) * H(K,J)
281	Z(I,J) = ZZ
282	880 CONTINUE
283	GO TO 1001
284	1000 IERR = EN
285	1001 RETURN
286	END