REAL FUNCTION FMIN(AX,BX,F,TOL)
C***BEGIN PROLOGUE FMIN
C***DATE WRITTEN 730101 (YYMMDD)
C***REVISION DATE 730101 (YYMMDD)
C***CATEGORY NO. G1A2
C***KEYWORDS ONE-DIMENSIONAL MINIMIZATION, UNIMODAL FUNCTION
C***AUTHOR BRENT, R.
C***PURPOSE An approximation to the point where F attains a minimum on
C the interval (AX,BX) is determined as the value of the
C function FMIN.
C***DESCRIPTION
C
C From the book, "Numerical Methods and Software" by
C D. Kahaner, C. Moler, S. Nash
C Prentice Hall, 1988
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(FMIN) + (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 FMIN
C approximates the abcissa of the global minimum of F on the
C interval AX,BX with an error less than 3*EPS*ABS(FMIN) + TOL.
C If F is not unimodal, then FMIN 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 INPUT PARAMETERS
C
C AX (real) left endpoint of initial interval
C BX (real) right endpoint of initial interval
C F Real function of the form REAL FUNCTION F(X) which evaluates
C F(X) for any X in the interval (AX,BX)
C Must be declared EXTERNAL in calling routine.
C TOL (real) desired length of the interval of uncertainty of the
C final result ( .ge. 0.0)
C
C
C OUTPUT PARAMETERS
C
C FMIN abcissa approximating the minimizer of F
C AX lower bound for minimizer
C BX upper bound for minimizer
C
C***REFERENCES RICHARD BRENT, ALGORITHMS FOR MINIMIZATION WITHOUT
C DERIVATIVES, PRENTICE-HALL, INC. (1973).
C***ROUTINES CALLED (NONE)
C***END PROLOGUE FMIN
REAL AX,BX,F,TOL
REAL A,B,C,D,E,EPS,XM,P,Q,R,TOL1,TOL2,U,V,W
REAL FU,FV,FW,FX,X
REAL ABS,SQRT,SIGN
C***FIRST EXECUTABLE STATEMENT FMIN
C = 0.5*(3. - SQRT(5.0))
C
C C is the squared inverse of the golden ratio
C
C EPS is approximately the square root of the relative machine
C precision.
C
EPS = 1.0
10 EPS = EPS/2.0
TOL1 = 1.0 + EPS
IF (TOL1 .GT. 1.0) GO TO 10
EPS = SQRT(EPS)
C
C initialization
C
A = AX
B = BX
V = A + C*(B - A)
W = V
X = V
E = 0.0
FX = F(X)
FV = FX
FW = FX
C
C main loop starts here
C
20 XM = 0.5*(A + B)
TOL1 = EPS*ABS(X) + TOL/3.0
TOL2 = 2.0*TOL1
C
C check stopping criterion
C
IF (ABS(X - XM) .LE. (TOL2 - 0.5*(B - A))) GO TO 90
C
C is golden-section necessary
C
IF (ABS(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.0*(Q - R)
IF (Q .GT. 0.0) P = -P
Q = ABS(Q)
R = E
E = D
C
C is parabola acceptable
C
30 IF (ABS(P) .GE. ABS(0.5*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 = SIGN(TOL1, XM - X)
IF ((B - U) .LT. TOL2) D = SIGN(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 (ABS(D) .GE. TOL1) U = X + D
IF (ABS(D) .LT. TOL1) U = X + SIGN(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 FMIN = X
RETURN
END