1 REAL FUNCTION ZEROIN (AX,BX,F,TOL)
2 REAL AX,BX,F,TOL
3 C
4 C
5 REAL A, B, C, D, E, EPS, FA, FB, FC, TOL1, XM, P, Q, R, S
6 C
7 C COMPUTER EPS, THE RELATIVE MACHINE PRECISION
8 C
9 EPS = 1.0
10 10 EPS = EPS/2.0
11 TOL1 = 1.0 + EPS
12 IF (TOL1 .GT. 1.0) GO TO 10
13 C
14 C INITIALIZATION
15 C
18 A = AX
19 B = BX
20 FA = F(A)
21 FB = F(B)
22 C
23 C BEGIN STEP
24 C
25 20 C = A
26 FC = FA
27 D = B - A
28 E = D
31 IF (ABS(FC) .GE. ABS(FB)) GO TO 40
32 A = B
33 B = C
34 C = A
35 FA = FB
36 FB = FC
37 FC = FA
38 C
39 C CONVERGENCE TEST
40 C
41 40 TOL1 = 2.0*EPS*ABS(B) + 0.5*TOL
42 XM = .5*(C-B)
43 IF (ABS(XM) .LE. TOL1) GO TO 90
44 IF (FB .EQ. 0.0) GO TO 90
45 C
46 C IS BISECTION NECESSARY
47 C
48 IF (ABS(E) .LT. TOL1) GO TO 70
49 IF (ABS(FA) .LE. ABS(FB)) GO TO 70
50 C
51 C IS QUADRATIC INTERPOLATION POSSIBLE
52 C
53 IF (A .NE. C) GO TO 50
54 C
55 C LINEAR INTERPOLATION
56 C
57 S = FB/FA
58 P = 2.0*XM*S
59 Q = 1.0 - S
60 GO TO 60
61 C
62 C INVERSE QUADRATIC INTERPOLATION
63 C
64 50 Q = FA/FC
65 R = FB/FC
66 S = FB/FA
67 P = S*(2.0*XM*Q*(Q-R) - (B-A) * (R-1.0))
68 Q = (Q-1.0)*(R-1.0)*(S-1.0)
69 C
70 C ADJUST SIGNS
71 C
72 60 IF (P .GT. 0.0) Q = -Q
73 P = ABS(P)
74 C
75 C IS INTERPOLATION ACCEPTABLE
76 C
77 IF ((2.0*P) .GE. (3.0*XM*Q - ABS(TOL1*Q))) GO TO 70
78 IF (P .GE. ABS(0.5*E*Q)) GO TO 70
79 E = D
80 D = P/Q
81 GO TO 80
82 C
83 C BISECTION
84 C
85 70 D = XM
86 E = D
87 C
88 C COMPLETE STEP
89 C
90 80 A = B
91 FA = FB
92 IF (ABS(D) .GT. TOL1) B = B + D
93 IF (ABS(D) .LE. TOL1) B = B + SIGN(TOL1,XM)
94 FB = F(B)
95 IF ((FB*(FC/ABS (FC))) .GT. 0.0) GO TO 20
96 GO TO 30
97 C
98 C DONE
99 C
100 90 ZEROIN = B
101 RETURN
102 END