*DECK DFZERO
SUBROUTINE DFZERO (F, B, C, R, RE, AE, IFLAG)
C***BEGIN PROLOGUE DFZERO
C***PURPOSE Search for a zero of a function F(X) in a given interval
C (B,C). It is designed primarily for problems where F(B)
C and F(C) have opposite signs.
C***LIBRARY SLATEC
C***CATEGORY F1B
C***TYPE DOUBLE PRECISION (FZERO-S, DFZERO-D)
C***KEYWORDS BISECTION, NONLINEAR, ROOTS, ZEROS
C***AUTHOR Shampine, L. F., (SNLA)
C Watts, H. A., (SNLA)
C***DESCRIPTION
C
C DFZERO searches for a zero of a DOUBLE PRECISION function F(X)
C between the given DOUBLE PRECISION values B and C until the width
C of the interval (B,C) has collapsed to within a tolerance
C specified by the stopping criterion,
C ABS(B-C) .LE. 2.*(RW*ABS(B)+AE).
C The method used is an efficient combination of bisection and the
C secant rule and is due to T. J. Dekker.
C
C Description Of Arguments
C
C F :EXT - Name of the DOUBLE PRECISION external function. This
C name must be in an EXTERNAL statement in the calling
C program. F must be a function of one DOUBLE
C PRECISION argument.
C
C B :INOUT - One end of the DOUBLE PRECISION interval (B,C). The
C value returned for B usually is the better
C approximation to a zero of F.
C
C C :INOUT - The other end of the DOUBLE PRECISION interval (B,C)
C
C R :IN - A (better) DOUBLE PRECISION guess of a zero of F
C which could help in speeding up convergence. If F(B)
C and F(R) have opposite signs, a root will be found in
C the interval (B,R); if not, but F(R) and F(C) have
C opposite signs, a root will be found in the interval
C (R,C); otherwise, the interval (B,C) will be
C searched for a possible root. When no better guess
C is known, it is recommended that R be set to B or C,
C since if R is not interior to the interval (B,C), it
C will be ignored.
C
C RE :IN - 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 :IN - Absolute error used in the stopping criterion. If
C the given interval (B,C) contains the origin, then a
C nonzero value should be chosen for AE.
C
C IFLAG :OUT - A status code. User must check IFLAG after each
C call. Control returns to the user from DFZERO in all
C 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
C***REFERENCES L. F. Shampine and H. A. Watts, FZERO, a root-solving
C code, Report SC-TM-70-631, Sandia Laboratories,
C 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 and P. Henrici, Wiley-Interscience, 1969.
C***ROUTINES CALLED D1MACH
C***REVISION HISTORY (YYMMDD)
C 700901 DATE WRITTEN
C 890531 Changed all specific intrinsics to generic. (WRB)
C 890531 REVISION DATE from Version 3.2
C 891214 Prologue converted to Version 4.0 format. (BAB)
C 920501 Reformatted the REFERENCES section. (WRB)
C***END PROLOGUE DFZERO
DOUBLE PRECISION A,ACBS,ACMB,AE,AW,B,C,CMB,D1MACH,ER,
+ F,FA,FB,FC,FX,FZ,P,Q,R,RE,RW,T,TOL,Z
INTEGER IC,IFLAG,KOUNT
C
C***FIRST EXECUTABLE STATEMENT DFZERO
C
C ER is two times the computer unit roundoff value which is defined
C here by the function D1MACH.
C
ER = 2.0D0 * D1MACH(4)
C
C Initialize.
C
Z = R
IF (R .LE. MIN(B,C) .OR. R .GE. MAX(B,C)) Z = C
RW = MAX(RE,ER)
AW = MAX(AE,0.D0)
IC = 0
T = Z
FZ = F(T)
FC = FZ
T = B
FB = F(T)
KOUNT = 2
IF (SIGN(1.0D0,FZ) .EQ. SIGN(1.0D0,FB)) GO TO 1
C = Z
GO TO 2
1 IF (Z .EQ. C) GO TO 2
T = C
FC = F(T)
KOUNT = 3
IF (SIGN(1.0D0,FZ) .EQ. SIGN(1.0D0,FC)) GO TO 2
B = Z
FB = FZ
2 A = C
FA = FC
ACBS = ABS(B-C)
FX = MAX(ABS(FB),ABS(FC))
C
3 IF (ABS(FC) .GE. ABS(FB)) GO TO 4
C
C Perform interchange.
C
A = B
FA = FB
B = C
FB = FC
C = A
FC = FA
C
4 CMB = 0.5D0*(C-B)
ACMB = ABS(CMB)
TOL = RW*ABS(B) + AW
C
C Test stopping criterion and function count.
C
IF (ACMB .LE. TOL) GO TO 10
IF (FB .EQ. 0.D0) GO TO 11
IF (KOUNT .GE. 500) GO TO 14
C
C Calculate new iterate implicitly as B+P/Q, where we arrange
C P .GE. 0. The implicit form is used to prevent overflow.
C
P = (B-A)*FB
Q = FA - FB
IF (P .GE. 0.D0) GO TO 5
P = -P
Q = -Q
C
C Update A and check for satisfactory reduction in the size of the
C bracketing interval. If not, perform bisection.
C
5 A = B
FA = FB
IC = IC + 1
IF (IC .LT. 4) GO TO 6
IF (8.0D0*ACMB .GE. ACBS) GO TO 8
IC = 0
ACBS = ACMB
C
C Test for too small a change.
C
6 IF (P .GT. ABS(Q)*TOL) GO TO 7
C
C Increment by TOLerance.
C
B = B + SIGN(TOL,CMB)
GO TO 9
C
C Root ought to be between B and (C+B)/2.
C
7 IF (P .GE. CMB*Q) GO TO 8
C
C Use secant rule.
C
B = B + P/Q
GO TO 9
C
C Use bisection (C+B)/2.
C
8 B = B + CMB
C
C Have completed computation for new iterate B.
C
9 T = B
FB = F(T)
KOUNT = KOUNT + 1
C
C Decide whether next step is interpolation or extrapolation.
C
IF (SIGN(1.0D0,FB) .NE. SIGN(1.0D0,FC)) GO TO 3
C = A
FC = FA
GO TO 3
C
C Finished. Process results for proper setting of IFLAG.
C
10 IF (SIGN(1.0D0,FB) .EQ. SIGN(1.0D0,FC)) GO TO 13
IF (ABS(FB) .GT. FX) GO TO 12
IFLAG = 1
RETURN
11 IFLAG = 2
RETURN
12 IFLAG = 3
RETURN
13 IFLAG = 4
RETURN
14 IFLAG = 5
RETURN
END