SUBROUTINE DFZERO(F,B,C,R,RE,AE,IFLAG)
C***BEGIN PROLOGUE DFZERO
C***DATE WRITTEN 700901 (YYMMDD)
C***REVISION DATE 861211 (YYMMDD)
C***CATEGORY NO. F1B
C***KEYWORDS LIBRARY=SLATEC,TYPE=DOUBLE PRECISION(FZERO-S DFZERO-D),
C BISECTION,NONLINEAR,NONLINEAR EQUATIONS,ROOTS,ZEROES,
C ZEROS
C***AUTHOR SHAMPINE,L.F.,SNLA
C WATTS,H.A.,SNLA
C***PURPOSE Search for a zero of a function F(X) in a given
C interval (B,C). It is designed primarily for problems
C where F(B) and F(C) have opposite signs.
C***DESCRIPTION
C
C **** Double Precision version of FZERO ****
C
C Based on a method by T J Dekker
C written by L F Shampine and H A Watts
C
C DFZERO searches for a zero of a function F(X) between
C the given values B and C until the width of the interval
C (B,C) has collapsed to within a tolerance specified by
C the stopping criterion, DABS(B-C) .LE. 2.*(RW*DABS(B)+AE).
C The method used is an efficient combination of bisection
C and the secant rule.
C
C Description Of Arguments
C
C F,B,C,R,RE and AE are DOUBLE PRECISION input parameters.
C B and C are DOUBLE PRECISION output parameters and IFLAG (flagged
C by an * below).
C
C F - Name of the DOUBLE PRECISION valued external function.
C This name must be in an EXTERNAL statement in the
C calling program. F must be a function of one double
C precision argument.
C
C *B - One end of the interval (B,C). The value returned for
C B usually is the better approximation to a zero of F.
C
C *C - The other end of the interval (B,C)
C
C R - A (better) guess of a zero of F which could help in
C speeding up convergence. If F(B) and F(R) have
C opposite signs, a root will be found in the interval
C (B,R); if not, but F(R) and F(C) have opposite
C signs, a root will be found in the interval (R,C);
C otherwise, the interval (B,C) will be searched for a
C possible root. When no better guess is known, it is
C recommended that r be set to B or C; because if R is
C not interior to the interval (B,C), it will be ignored.
C
C RE - Relative error used for RW in the stopping criterion.
C If the requested RE is less than machine precision,
C then RW is set to approximately machine precision.
C
C AE - Absolute error used in the stopping criterion. If the
C given interval (B,C) contains the origin, then a
C nonzero value should be chosen for AE.
C
C *IFLAG - A status code. User must check IFLAG after each call.
C Control returns to the user from FZERO in all cases.
C XERROR does not process diagnostics in these cases.
C
C 1 B is within the requested tolerance of a zero.
C The interval (B,C) collapsed to the requested
C tolerance, the function changes sign in (B,C), and
C F(X) decreased in magnitude as (B,C) collapsed.
C
C 2 F(B) = 0. However, the interval (B,C) may not have
C collapsed to the requested tolerance.
C
C 3 B may be near a singular point of F(X).
C The interval (B,C) collapsed to the requested tol-
C erance and the function changes sign in (B,C), but
C F(X) increased in magnitude as (B,C) collapsed,i.e.
C abs(F(B out)) .GT. max(abs(F(B in)),abs(F(C in)))
C
C 4 No change in sign of F(X) was found although the
C interval (B,C) collapsed to the requested tolerance.
C The user must examine this case and decide whether
C B is near a local minimum of F(X), or B is near a
C zero of even multiplicity, or neither of these.
C
C 5 Too many (.GT. 500) function evaluations used.
C***REFERENCES L. F. SHAMPINE AND H. A. WATTS, *FZERO, A ROOT-SOLVING
C CODE*, SC-TM-70-631, SEPTEMBER 1970.
C T. J. DEKKER, *FINDING A ZERO BY MEANS OF SUCCESSIVE
C LINEAR INTERPOLATION*, 'CONSTRUCTIVE ASPECTS OF THE
C FUNDAMENTAL THEOREM OF ALGEBRA', EDITED BY B. DEJON
C P. HENRICI, 1969.
C***ROUTINES CALLED D1MACH
C***END PROLOGUE DFZERO
INTEGER I, IC, ICNT, IERR, IFLAG, IPASS, IPSS, ITEST(36),
* ITMP(15), J, KLUS, KOUNT, KPRINT, LUN, NDEG
DOUBLE PRECISION A, ACBS, ACMB, AE, AW, B, C,
* CMB, D1MACH, DABS,
* DMAX1, DMIN1, DSIGN1, DSQRT, ER, F, FA, FB, FC,
* FX, FZ, P, Q, R, RE, REL, RW, T, TOL, WI, WORK, WR, Z
C BEGIN BLOCK PERMITTING ...EXITS TO 200
C ER IS TWO TIMES THE COMPUTER UNIT ROUNDOFF VALUE WHICH IS
C DEFINED HERE BY THE FUNCTION D1MACH.
C***FIRST EXECUTABLE STATEMENT DFZERO
ccc ER = 2.0D0*D1MACH(4)
ccc For IEEE 754 double precision arithmetic
ccc d1mach(4) should be 2^{-53} x 2 = 2.22e-16.
ccc Added by Steve Verrill on 5/25/02.
er = 2.0d0*2.22d-16
C
C INITIALIZE
C
Z = R
IF (R .LE. DMIN1(B,C) .OR. R .GE. DMAX1(B,C)) Z = C
RW = DMAX1(RE,ER)
AW = DMAX1(AE,0.0D0)
IC = 0
T = Z
FZ = F(T)
FC = FZ
T = B
FB = F(T)
KOUNT = 2
IF (DSIGN(1.0D0,FZ) .EQ. DSIGN(1.0D0,FB)) GO TO 10
C = Z
GO TO 30
10 CONTINUE
C BEGIN BLOCK PERMITTING ...EXITS TO 20
C ...EXIT
IF (Z .EQ. C) GO TO 20
T = C
FC = F(T)
KOUNT = 3
C ...EXIT
IF (DSIGN(1.0D0,FZ) .EQ. DSIGN(1.0D0,FC)) GO TO 20
B = Z
FB = FZ
20 CONTINUE
30 CONTINUE
A = C
FA = FC
ACBS = DABS(B-C)
FX = DMAX1(DABS(FB),DABS(FC))
C
40 CONTINUE
C BEGIN BLOCK PERMITTING ...EXITS TO 180
IF (DABS(FC) .GE. DABS(FB)) GO TO 50
C PERFORM INTERCHANGE
A = B
FA = FB
B = C
FB = FC
C = A
FC = FA
50 CONTINUE
C
CMB = 0.5D0*(C - B)
ACMB = DABS(CMB)
TOL = RW*DABS(B) + AW
C
C TEST STOPPING CRITERION AND FUNCTION COUNT
C
IF (ACMB .GT. TOL) GO TO 90
C
C
C FINISHED. PROCESS RESULTS FOR PROPER SETTING OF IFLAG
C
IF (DSIGN(1.0D0,FB) .NE. DSIGN(1.0D0,FC)) GO TO 60
IFLAG = 4
GO TO 80
60 CONTINUE
IF (DABS(FB) .LE. FX) GO TO 70
IFLAG = 3
GO TO 80
70 CONTINUE
IFLAG = 1
80 CONTINUE
C ............EXIT
GO TO 200
90 CONTINUE
IF (FB .NE. 0.0D0) GO TO 100
IFLAG = 2
C ............EXIT
GO TO 200
100 CONTINUE
IF (KOUNT .LT. 500) GO TO 110
IFLAG = 5
C ............EXIT
GO TO 200
110 CONTINUE
C
C CALCULATE NEW ITERATE IMPLICITLY AS B+P/Q
C WHERE WE ARRANGE P .GE. 0.
C THE IMPLICIT FORM IS USED TO PREVENT OVERFLOW.
C
P = (B - A)*FB
Q = FA - FB
IF (P .GE. 0.0D0) GO TO 120
P = -P
Q = -Q
120 CONTINUE
C
C UPDATE A AND CHECK FOR SATISFACTORY REDUCTION
C IN THE SIZE OF THE BRACKETING INTERVAL.
C IF NOT, PERFORM BISECTION.
C
A = B
FA = FB
IC = IC + 1
IF (IC .LT. 4) GO TO 140
IF (8.0D0*ACMB .LT. ACBS) GO TO 130
C
C USE BISECTION
C
B = 0.5D0*(C + B)
C .........EXIT
GO TO 180
130 CONTINUE
IC = 0
ACBS = ACMB
140 CONTINUE
C
C TEST FOR TOO SMALL A CHANGE
C
IF (P .GT. DABS(Q)*TOL) GO TO 150
C
C INCREMENT BY TOLERANCE
C
B = B + DSIGN(TOL,CMB)
GO TO 170
150 CONTINUE
C
C ROOT OUGHT TO BE BETWEEN B AND (C+B)/2.
C
IF (P .LT. CMB*Q) GO TO 160
C
C USE BISECTION
C
B = 0.5D0*(C + B)
GO TO 170
160 CONTINUE
C
C USE SECANT RULE
C
B = B + P/Q
170 CONTINUE
180 CONTINUE
C
C HAVE COMPLETED COMPUTATION FOR NEW ITERATE B
C
T = B
FB = F(T)
KOUNT = KOUNT + 1
C
C DECIDE WHETHER NEXT STEP IS INTERPOLATION OR EXTRAPOLATION
C
IF (DSIGN(1.0D0,FB) .NE. DSIGN(1.0D0,FC)) GO TO 190
C = A
FC = FA
190 CONTINUE
GO TO 40
200 CONTINUE
RETURN
END