This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / f / bisect.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:02:33  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10       SUBROUTINE BISECT(N,EPS1,D,E,E2,LB,UB,MM,M,W,IND,IERR,RV4,RV5)
11       INTEGER I,J,K,L,M,N,P,Q,R,S,II,MM,M1,M2,TAG,IERR,ISTURM
12       REAL D(N),E(N),E2(N),W(MM),RV4(N),RV5(N)
13       REAL U,V,LB,T1,T2,UB,XU,X0,X1,EPS1,MACHEP
14       INTEGER IND(MM)
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       TAG = 0
23       T1 = LB
24       T2 = UB
25       DO 40 I = 1, N
26          IF (I .EQ. 1) GO TO 20
27          IF (ABS(E(I)) .GT. MACHEP * (ABS(D(I)) + ABS(D(I-1))))
28      X      GO TO 40
29    20    E2(I) = 0.0
30    40 CONTINUE
31       P = 1
32       Q = N
33       X1 = UB
34       ISTURM = 1
35       GO TO 320
36    60 M = S
37       X1 = LB
38       ISTURM = 2
39       GO TO 320
40    80 M = M - S
41       IF (M .GT. MM) GO TO 980
42       Q = 0
43       R = 0
44   100 IF (R .EQ. M) GO TO 1001
45       TAG = TAG + 1
46       P = Q + 1
47       XU = D(P)
48       X0 = D(P)
49       U = 0.0
50       DO 120 Q = P, N
51          X1 = U
52          U = 0.0
53          V = 0.0
54          IF (Q .EQ. N) GO TO 110
55          U = ABS(E(Q+1))
56          V = E2(Q+1)
57   110    XU = MIN(D(Q)-(X1+U),XU)
58          X0 = MAX(D(Q)+(X1+U),X0)
59          IF (V .EQ. 0.0) GO TO 140
60   120 CONTINUE
61   140 X1 = MAX(ABS(XU),ABS(X0)) * MACHEP
62       IF (EPS1 .LE. 0.0) EPS1 = -X1
63       IF (P .NE. Q) GO TO 180
64       IF (T1 .GT. D(P) .OR. D(P) .GE. T2) GO TO 940
65       M1 = P
66       M2 = P
67       RV5(P) = D(P)
68       GO TO 900
69   180 X1 = X1 * (Q-P+1)
70       LB = MAX(T1,XU-X1)
71       UB = MIN(T2,X0+X1)
72       X1 = LB
73       ISTURM = 3
74       GO TO 320
75   200 M1 = S + 1
76       X1 = UB
77       ISTURM = 4
78       GO TO 320
79   220 M2 = S
80       IF (M1 .GT. M2) GO TO 940
81       X0 = UB
82       ISTURM = 5
83       DO 240 I = M1, M2
84          RV5(I) = UB
85          RV4(I) = LB
86   240 CONTINUE
87       K = M2
88   250    XU = LB
89          DO 260 II = M1, K
90             I = M1 + K - II
91             IF (XU .GE. RV4(I)) GO TO 260
92             XU = RV4(I)
93             GO TO 280
94   260    CONTINUE
95   280    IF (X0 .GT. RV5(K)) X0 = RV5(K)
96   300    X1 = (XU + X0) * 0.5
97          IF ((X0 - XU) .LE. (2.0 * MACHEP *
98      X      (ABS(XU) + ABS(X0)) + ABS(EPS1))) GO TO 420
99   320    S = P - 1
100          U = 1.0
101          DO 340 I = P, Q
102             IF (U .NE. 0.0) GO TO 325
103             V = ABS(E(I)) / MACHEP
104             GO TO 330
105   325       V = E2(I) / U
106   330       U = D(I) - X1 - V
107             IF (U .LT. 0.0) S = S + 1
108   340    CONTINUE
109          GO TO (60,80,200,220,360), ISTURM
110   360    IF (S .GE. K) GO TO 400
111          XU = X1
112          IF (S .GE. M1) GO TO 380
113          RV4(M1) = X1
114          GO TO 300
115   380    RV4(S+1) = X1
116          IF (RV5(S) .GT. X1) RV5(S) = X1
117          GO TO 300
118   400    X0 = X1
119          GO TO 300
120   420    RV5(K) = X1
121       K = K - 1
122       IF (K .GE. M1) GO TO 250
123   900 S = R
124       R = R + M2 - M1 + 1
125       J = 1
126       K = M1
127       DO 920 L = 1, R
128          IF (J .GT. S) GO TO 910
129          IF (K .GT. M2) GO TO 940
130          IF (RV5(K) .GE. W(L)) GO TO 915
131          DO 905 II = J, S
132             I = L + S - II
133             W(I+1) = W(I)
134             IND(I+1) = IND(I)
135   905    CONTINUE
136   910    W(L) = RV5(K)
137          IND(L) = TAG
138          K = K + 1
139          GO TO 920
140   915    J = J + 1
141   920 CONTINUE
142   940 IF (Q .LT. N) GO TO 100
143       GO TO 1001
144   980 IERR = 3 * N + 1
145  1001 LB = T1
146       UB = T2
147       RETURN
148       END