      DOUBLE PRECISION FUNCTION DFMIN(AX,BX,F,TOL)
      DOUBLE PRECISION AX,BX,F,TOL
C***BEGIN PROLOGUE  DFMIN
C***CATEGORY NO.    G1A2
C***KEYWORD(S)  ONE-DIMENSIONAL MINIMIZATION, UNIMODAL FUNCTION
C***AUTHOR  R. BRENT
C***DATE WRITTEN    730101  (YYMMDD)
C***PURPOSE
C     An approximation to the point where  F  attains a minimum  on
C     the interval (AX,BX) is determined as the value of the function
C     DFMIN.
C
C PARAMETERS
C
C INPUT
C
C  AX    (double precision)  left endpoint of initial interval
C  BX    (double precision) right endpoint of initial interval
C  F     function subprogram which evaluates F(X)  for any  X
C        in the interval  (AX,BX)
C  TOL   (double precision) desired length of the interval of uncertainty 
C        of the final result ( .ge. 0.0)
C
C
C OUTPUT
C
C DFMIN   abcissa approximating the minimizer of F
C AX     lower bound for minimizer
C BX     upper bound for minimizer
C
C***DESCRIPTION
C
C     The method used is a combination of golden section search and
C     successive parabolic interpolation.  Convergence is never much 
C     slower than that for a Fibonacci search.  If F has a continuous 
C     second derivative which is positive at the minimum (which is not
C     at AX or BX), then convergence is superlinear, and usually of the 
C     order of about 1.324....
C
C     The function F is never evaluated at two points closer together
C     than EPS*ABS(DFMIN) + (TOL/3), where EPS is approximately the 
C     square root of the relative machine precision.  If F is a unimodal
C     function and the computed values of F are always unimodal when
C     separated by at least EPS*ABS(XSTAR) + (TOL/3), then DFMIN 
C     approximates the abcissa of the global minimum of F on the 
C     interval AX,BX with an error less than 3*EPS*ABS(DFMIN) + TOL.  
C     If F is not unimodal, then DFMIN may approximate a local, but 
C     perhaps non-global, minimum to the same accuracy.
C
C     This function subprogram is a slightly modified version of the
C     ALGOL 60 procedure LOCALMIN given in Richard Brent, Algorithms for
C     Minimization Without Derivatives, Prentice-Hall, Inc. (1973).
C
C***REFERENCE(S)
C     Richard Brent, Algorithms for Minimization Without Derivatives,
C     Prentice-Hall, Inc. (1973).
C***ROUTINES CALLED   NONE
C***END PROLOGUE
      DOUBLE PRECISION  A,B,C,D,E,EPS,XM,P,Q,R,TOL1,TOL2,U,V,W
      DOUBLE PRECISION  FU,FV,FW,FX,X
      DOUBLE PRECISION  DABS,DSQRT,DSIGN
C
C  C is the squared inverse of the golden ratio
C
      C = 0.5D0*(3. - DSQRT(5.0D0))
C
C  EPS is approximately the square root of the relative machine
C  precision.
C
      EPS = 1.0D00
   10 EPS = EPS/2.0D00
      TOL1 = 1.0D0 + EPS
      IF (TOL1 .GT. 1.0D00) GO TO 10
      EPS = DSQRT(EPS)
C
C  initialization
C
      A = AX
      B = BX
      V = A + C*(B - A)
      W = V
      X = V
      E = 0.0D0
      FX = F(X)
      FV = FX
      FW = FX
C
C  main loop starts here
C
   20 XM = 0.5D0*(A + B)
      TOL1 = EPS*DABS(X) + TOL/3.0D0
      TOL2 = 2.0D0*TOL1
C
C  check stopping criterion
C
      IF (DABS(X - XM) .LE. (TOL2 - 0.5D0*(B - A))) GO TO 90
C
C is golden-section necessary
C
      IF (DABS(E) .LE. TOL1) GO TO 40
C
C  fit parabola
C
      R = (X - W)*(FX - FV)
      Q = (X - V)*(FX - FW)
      P = (X - V)*Q - (X - W)*R
      Q = 2.0D00*(Q - R)
      IF (Q .GT. 0.0D0) P = -P
      Q =  DABS(Q)
      R = E
      E = D
C
C  is parabola acceptable
C
   30 IF (DABS(P) .GE. DABS(0.5D0*Q*R)) GO TO 40
      IF (P .LE. Q*(A - X)) GO TO 40
      IF (P .GE. Q*(B - X)) GO TO 40
C
C  a parabolic interpolation step
C
      D = P/Q
      U = X + D
C
C  F must not be evaluated too close to AX or BX
C
      IF ((U - A) .LT. TOL2) D = DSIGN(TOL1, XM - X)
      IF ((B - U) .LT. TOL2) D = DSIGN(TOL1, XM - X)
      GO TO 50
C
C  a golden-section step
C
   40 IF (X .GE. XM) E = A - X
      IF (X .LT. XM) E = B - X
      D = C*E
C
C  F must not be evaluated too close to X
C
   50 IF (DABS(D) .GE. TOL1) U = X + D
      IF (DABS(D) .LT. TOL1) U = X + DSIGN(TOL1, D)
      FU = F(U)
C
C  update  A, B, V, W, and X
C
      IF (FU .GT. FX) GO TO 60
      IF (U .GE. X) A = X
      IF (U .LT. X) B = X
      V = W
      FV = FW
      W = X
      FW = FX
      X = U
      FX = FU
      GO TO 20
   60 IF (U .LT. X) A = U
      IF (U .GE. X) B = U
      IF (FU .LE. FW) GO TO 70
      IF (W .EQ. X) GO TO 70
      IF (FU .LE. FV) GO TO 80
      IF (V .EQ. X) GO TO 80
      IF (V .EQ. W) GO TO 80
      GO TO 20
   70 V = W
      FV = FW
      W = U
      FW = FU
      GO TO 20
   80 V = U
      FV = FU
      GO TO 20
C
C  end of main loop
C
   90 DFMIN = X
      RETURN
      END
