This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / f / tsturm.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:02:38  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10       SUBROUTINE TSTURM(NM,N,EPS1,D,E,E2,LB,UB,MM,M,W,Z,
11      X                  IERR,RV1,RV2,RV3,RV4,RV5,RV6)
12       INTEGER I,J,K,M,N,P,Q,R,S,II,IP,JJ,MM,M1,M2,NM,ITS,
13      X        IERR,GROUP,ISTURM
14       REAL D(N),E(N),E2(N),W(MM),Z(NM,MM),
15      X       RV1(N),RV2(N),RV3(N),RV4(N),RV5(N),RV6(N)
16       REAL U,V,LB,T1,T2,UB,UK,XU,X0,X1,EPS1,EPS2,EPS3,EPS4,
17      X       NORM,MACHEP
18 #if defined(CERNLIB_CDC)
19       MACHEP=2.**(-47)
20 #endif
21 #if !defined(CERNLIB_CDC)
22       MACHEP=2.**(-23)
23 #endif
24       IERR = 0
25       T1 = LB
26       T2 = UB
27       DO 40 I = 1, N
28          IF (I .EQ. 1) GO TO 20
29          IF (ABS(E(I)) .GT. MACHEP * (ABS(D(I)) + ABS(D(I-1))))
30      X      GO TO 40
31    20    E2(I) = 0.0
32    40 CONTINUE
33       P = 1
34       Q = N
35       X1 = UB
36       ISTURM = 1
37       GO TO 320
38    60 M = S
39       X1 = LB
40       ISTURM = 2
41       GO TO 320
42    80 M = M - S
43       IF (M .GT. MM) GO TO 980
44       Q = 0
45       R = 0
46   100 IF (R .EQ. M) GO TO 1001
47       P = Q + 1
48       XU = D(P)
49       X0 = D(P)
50       U = 0.0
51       DO 120 Q = P, N
52          X1 = U
53          U = 0.0
54          V = 0.0
55          IF (Q .EQ. N) GO TO 110
56          U = ABS(E(Q+1))
57          V = E2(Q+1)
58   110    XU = MIN(D(Q)-(X1+U),XU)
59          X0 = MAX(D(Q)+(X1+U),X0)
60          IF (V .EQ. 0.0) GO TO 140
61   120 CONTINUE
62   140 X1 = MAX(ABS(XU),ABS(X0)) * MACHEP
63       IF (EPS1 .LE. 0.0) EPS1 = -X1
64       IF (P .NE. Q) GO TO 180
65       IF (T1 .GT. D(P) .OR. D(P) .GE. T2) GO TO 940
66       R = R + 1
67       DO 160 I = 1, N
68   160 Z(I,R) = 0.0
69       W(R) = D(P)
70       Z(P,R) = 1.0
71       GO TO 940
72   180 X1 = X1 * (Q-P+1)
73       LB = MAX(T1,XU-X1)
74       UB = MIN(T2,X0+X1)
75       X1 = LB
76       ISTURM = 3
77       GO TO 320
78   200 M1 = S + 1
79       X1 = UB
80       ISTURM = 4
81       GO TO 320
82   220 M2 = S
83       IF (M1 .GT. M2) GO TO 940
84       X0 = UB
85       ISTURM = 5
86       DO 240 I = M1, M2
87          RV5(I) = UB
88          RV4(I) = LB
89   240 CONTINUE
90       K = M2
91   250    XU = LB
92          DO 260 II = M1, K
93             I = M1 + K - II
94             IF (XU .GE. RV4(I)) GO TO 260
95             XU = RV4(I)
96             GO TO 280
97   260    CONTINUE
98   280    IF (X0 .GT. RV5(K)) X0 = RV5(K)
99   300    X1 = (XU + X0) * 0.5
100          IF ((X0 - XU) .LE. (2.0 * MACHEP *
101      X      (ABS(XU) + ABS(X0)) + ABS(EPS1))) GO TO 420
102   320    S = P - 1
103          U = 1.0
104          DO 340 I = P, Q
105             IF (U .NE. 0.0) GO TO 325
106             V = ABS(E(I)) / MACHEP
107             GO TO 330
108   325       V = E2(I) / U
109   330       U = D(I) - X1 - V
110             IF (U .LT. 0.0) S = S + 1
111   340    CONTINUE
112          GO TO (60,80,200,220,360), ISTURM
113   360    IF (S .GE. K) GO TO 400
114          XU = X1
115          IF (S .GE. M1) GO TO 380
116          RV4(M1) = X1
117          GO TO 300
118   380    RV4(S+1) = X1
119          IF (RV5(S) .GT. X1) RV5(S) = X1
120          GO TO 300
121   400    X0 = X1
122          GO TO 300
123   420    RV5(K) = X1
124       K = K - 1
125       IF (K .GE. M1) GO TO 250
126       NORM = ABS(D(P))
127       IP = P + 1
128       DO 500 I = IP, Q
129   500 NORM = NORM + ABS(D(I)) + ABS(E(I))
130       EPS2 = 1.0E-3 * NORM
131       EPS3 = MACHEP * NORM
132       UK = Q-P+1
133       EPS4 = UK * EPS3
134       UK = EPS4 / SQRT(UK)
135       GROUP = 0
136       S = P
137       DO 920 K = M1, M2
138          R = R + 1
139          ITS = 1
140          W(R) = RV5(K)
141          X1 = RV5(K)
142          IF (K .EQ. M1) GO TO 520
143          IF (X1 - X0 .GE. EPS2) GROUP = -1
144          GROUP = GROUP + 1
145          IF (X1 .LE. X0) X1 = X0 + EPS3
146   520    V = 0.0
147          DO 580 I = P, Q
148             RV6(I) = UK
149             IF (I .EQ. P) GO TO 560
150             IF (ABS(E(I)) .LT. ABS(U)) GO TO 540
151             XU = U / E(I)
152             RV4(I) = XU
153             RV1(I-1) = E(I)
154             RV2(I-1) = D(I) - X1
155             RV3(I-1) = 0.0
156             IF (I .NE. Q) RV3(I-1) = E(I+1)
157             U = V - XU * RV2(I-1)
158             V = -XU * RV3(I-1)
159             GO TO 580
160   540       XU = E(I) / U
161             RV4(I) = XU
162             RV1(I-1) = U
163             RV2(I-1) = V
164             RV3(I-1) = 0.0
165   560       U = D(I) - X1 - XU * V
166             IF (I .NE. Q) V = E(I+1)
167   580    CONTINUE
168          IF (U .EQ. 0.0) U = EPS3
169          RV1(Q) = U
170          RV2(Q) = 0.0
171          RV3(Q) = 0.0
172   600    DO 620 II = P, Q
173             I = P + Q - II
174             RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
175             V = U
176             U = RV6(I)
177   620    CONTINUE
178          IF (GROUP .EQ. 0) GO TO 700
179          DO 680 JJ = 1, GROUP
180             J = R - GROUP - 1 + JJ
181             XU = 0.0
182             DO 640 I = P, Q
183   640       XU = XU + RV6(I) * Z(I,J)
184             DO 660 I = P, Q
185   660       RV6(I) = RV6(I) - XU * Z(I,J)
186   680    CONTINUE
187   700    NORM = 0.0
188          DO 720 I = P, Q
189   720    NORM = NORM + ABS(RV6(I))
190          IF (NORM .GE. 1.0) GO TO 840
191          IF (ITS .EQ. 5) GO TO 960
192          IF (NORM .NE. 0.0) GO TO 740
193          RV6(S) = EPS4
194          S = S + 1
195          IF (S .GT. Q) S = P
196          GO TO 780
197   740    XU = EPS4 / NORM
198          DO 760 I = P, Q
199   760    RV6(I) = RV6(I) * XU
200   780    DO 820 I = IP, Q
201             U = RV6(I)
202             IF (RV1(I-1) .NE. E(I)) GO TO 800
203             U = RV6(I-1)
204             RV6(I-1) = RV6(I)
205   800       RV6(I) = U - RV4(I) * RV6(I-1)
206   820    CONTINUE
207          ITS = ITS + 1
208          GO TO 600
209   840    U = 0.0
210          DO 860 I = P, Q
211   860    U = U + RV6(I)**2
212          XU = 1.0 / SQRT(U)
213          DO 880 I = 1, N
214   880    Z(I,R) = 0.0
215          DO 900 I = P, Q
216   900    Z(I,R) = RV6(I) * XU
217          X0 = X1
218   920 CONTINUE
219   940 IF (Q .LT. N) GO TO 100
220       GO TO 1001
221   960 IERR = 4 * N + R
222       GO TO 1001
223   980 IERR = 3 * N + 1
224  1001 LB = T1
225       UB = T2
226       RETURN
227       END