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

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:02:38  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      SUBROUTINE TSTURM(NM,N,EPS1,D,E,E2,LB,UB,MM,M,W,Z,
     X                  IERR,RV1,RV2,RV3,RV4,RV5,RV6)
      INTEGER I,J,K,M,N,P,Q,R,S,II,IP,JJ,MM,M1,M2,NM,ITS,
     X        IERR,GROUP,ISTURM
      REAL D(N),E(N),E2(N),W(MM),Z(NM,MM),
     X       RV1(N),RV2(N),RV3(N),RV4(N),RV5(N),RV6(N)
      REAL U,V,LB,T1,T2,UB,UK,XU,X0,X1,EPS1,EPS2,EPS3,EPS4,
     X       NORM,MACHEP
#if defined(CERNLIB_CDC)
      MACHEP=2.**(-47)
#endif
#if !defined(CERNLIB_CDC)
      MACHEP=2.**(-23)
#endif
      IERR = 0
      T1 = LB
      T2 = UB
      DO 40 I = 1, N
         IF (I .EQ. 1) GO TO 20
         IF (ABS(E(I)) .GT. MACHEP * (ABS(D(I)) + ABS(D(I-1))))
     X      GO TO 40
   20    E2(I) = 0.0
   40 CONTINUE
      P = 1
      Q = N
      X1 = UB
      ISTURM = 1
      GO TO 320
   60 M = S
      X1 = LB
      ISTURM = 2
      GO TO 320
   80 M = M - S
      IF (M .GT. MM) GO TO 980
      Q = 0
      R = 0
  100 IF (R .EQ. M) GO TO 1001
      P = Q + 1
      XU = D(P)
      X0 = D(P)
      U = 0.0
      DO 120 Q = P, N
         X1 = U
         U = 0.0
         V = 0.0
         IF (Q .EQ. N) GO TO 110
         U = ABS(E(Q+1))
         V = E2(Q+1)
  110    XU = MIN(D(Q)-(X1+U),XU)
         X0 = MAX(D(Q)+(X1+U),X0)
         IF (V .EQ. 0.0) GO TO 140
  120 CONTINUE
  140 X1 = MAX(ABS(XU),ABS(X0)) * MACHEP
      IF (EPS1 .LE. 0.0) EPS1 = -X1
      IF (P .NE. Q) GO TO 180
      IF (T1 .GT. D(P) .OR. D(P) .GE. T2) GO TO 940
      R = R + 1
      DO 160 I = 1, N
  160 Z(I,R) = 0.0
      W(R) = D(P)
      Z(P,R) = 1.0
      GO TO 940
  180 X1 = X1 * (Q-P+1)
      LB = MAX(T1,XU-X1)
      UB = MIN(T2,X0+X1)
      X1 = LB
      ISTURM = 3
      GO TO 320
  200 M1 = S + 1
      X1 = UB
      ISTURM = 4
      GO TO 320
  220 M2 = S
      IF (M1 .GT. M2) GO TO 940
      X0 = UB
      ISTURM = 5
      DO 240 I = M1, M2
         RV5(I) = UB
         RV4(I) = LB
  240 CONTINUE
      K = M2
  250    XU = LB
         DO 260 II = M1, K
            I = M1 + K - II
            IF (XU .GE. RV4(I)) GO TO 260
            XU = RV4(I)
            GO TO 280
  260    CONTINUE
  280    IF (X0 .GT. RV5(K)) X0 = RV5(K)
  300    X1 = (XU + X0) * 0.5
         IF ((X0 - XU) .LE. (2.0 * MACHEP *
     X      (ABS(XU) + ABS(X0)) + ABS(EPS1))) GO TO 420
  320    S = P - 1
         U = 1.0
         DO 340 I = P, Q
            IF (U .NE. 0.0) GO TO 325
            V = ABS(E(I)) / MACHEP
            GO TO 330
  325       V = E2(I) / U
  330       U = D(I) - X1 - V
            IF (U .LT. 0.0) S = S + 1
  340    CONTINUE
         GO TO (60,80,200,220,360), ISTURM
  360    IF (S .GE. K) GO TO 400
         XU = X1
         IF (S .GE. M1) GO TO 380
         RV4(M1) = X1
         GO TO 300
  380    RV4(S+1) = X1
         IF (RV5(S) .GT. X1) RV5(S) = X1
         GO TO 300
  400    X0 = X1
         GO TO 300
  420    RV5(K) = X1
      K = K - 1
      IF (K .GE. M1) GO TO 250
      NORM = ABS(D(P))
      IP = P + 1
      DO 500 I = IP, Q
  500 NORM = NORM + ABS(D(I)) + ABS(E(I))
      EPS2 = 1.0E-3 * NORM
      EPS3 = MACHEP * NORM
      UK = Q-P+1
      EPS4 = UK * EPS3
      UK = EPS4 / SQRT(UK)
      GROUP = 0
      S = P
      DO 920 K = M1, M2
         R = R + 1
         ITS = 1
         W(R) = RV5(K)
         X1 = RV5(K)
         IF (K .EQ. M1) GO TO 520
         IF (X1 - X0 .GE. EPS2) GROUP = -1
         GROUP = GROUP + 1
         IF (X1 .LE. X0) X1 = X0 + EPS3
  520    V = 0.0
         DO 580 I = P, Q
            RV6(I) = UK
            IF (I .EQ. P) GO TO 560
            IF (ABS(E(I)) .LT. ABS(U)) GO TO 540
            XU = U / E(I)
            RV4(I) = XU
            RV1(I-1) = E(I)
            RV2(I-1) = D(I) - X1
            RV3(I-1) = 0.0
            IF (I .NE. Q) RV3(I-1) = E(I+1)
            U = V - XU * RV2(I-1)
            V = -XU * RV3(I-1)
            GO TO 580
  540       XU = E(I) / U
            RV4(I) = XU
            RV1(I-1) = U
            RV2(I-1) = V
            RV3(I-1) = 0.0
  560       U = D(I) - X1 - XU * V
            IF (I .NE. Q) V = E(I+1)
  580    CONTINUE
         IF (U .EQ. 0.0) U = EPS3
         RV1(Q) = U
         RV2(Q) = 0.0
         RV3(Q) = 0.0
  600    DO 620 II = P, Q
            I = P + Q - II
            RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
            V = U
            U = RV6(I)
  620    CONTINUE
         IF (GROUP .EQ. 0) GO TO 700
         DO 680 JJ = 1, GROUP
            J = R - GROUP - 1 + JJ
            XU = 0.0
            DO 640 I = P, Q
  640       XU = XU + RV6(I) * Z(I,J)
            DO 660 I = P, Q
  660       RV6(I) = RV6(I) - XU * Z(I,J)
  680    CONTINUE
  700    NORM = 0.0
         DO 720 I = P, Q
  720    NORM = NORM + ABS(RV6(I))
         IF (NORM .GE. 1.0) GO TO 840
         IF (ITS .EQ. 5) GO TO 960
         IF (NORM .NE. 0.0) GO TO 740
         RV6(S) = EPS4
         S = S + 1
         IF (S .GT. Q) S = P
         GO TO 780
  740    XU = EPS4 / NORM
         DO 760 I = P, Q
  760    RV6(I) = RV6(I) * XU
  780    DO 820 I = IP, Q
            U = RV6(I)
            IF (RV1(I-1) .NE. E(I)) GO TO 800
            U = RV6(I-1)
            RV6(I-1) = RV6(I)
  800       RV6(I) = U - RV4(I) * RV6(I-1)
  820    CONTINUE
         ITS = ITS + 1
         GO TO 600
  840    U = 0.0
         DO 860 I = P, Q
  860    U = U + RV6(I)**2
         XU = 1.0 / SQRT(U)
         DO 880 I = 1, N
  880    Z(I,R) = 0.0
         DO 900 I = P, Q
  900    Z(I,R) = RV6(I) * XU
         X0 = X1
  920 CONTINUE
  940 IF (Q .LT. N) GO TO 100
      GO TO 1001
  960 IERR = 4 * N + R
      GO TO 1001
  980 IERR = 3 * N + 1
 1001 LB = T1
      UB = T2
      RETURN
      END
Commit	Line	Data
fe4da5cc	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