real function fmin(ax,bx,f,tol)
real ax,bx,f,tol
c
c an approximation x to the point where f attains a minimum on
c the interval (ax,bx) is determined.
c
c
c input..
c
c ax left endpoint of initial interval
c bx right endpoint of initial interval
c f function subprogram which evaluates f(x) for any x
c in the interval (ax,bx)
c tol desired length of the interval of uncertainty of the final
c result ( .ge. 0.0)
c
c
c output..
c
c fmin abcissa approximating the point where f attains a minimum
c
c
c the method used is a combination of golden section search and
c successive parabolic interpolation. convergence is never much slower
c than that for a fibonacci search. if f has a continuous second
c derivative which is positive at the minimum (which is not at ax or
c bx), then convergence is superlinear, and usually of the order of
c about 1.324....
c the function f is never evaluated at two points closer together
c than eps*abs(fmin) + (tol/3), where eps is approximately the square
c 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(x) + (tol/3), then fmin approximates
c the abcissa of the global minimum of f on the interval ax,bx with
c an error less than 3*eps*abs(fmin) + tol. if f is not unimodal,
c then fmin may approximate a local, but perhaps non-global, minimum to
c the same accuracy.
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
real a,b,c,d,e,eps,xm,p,q,r,tol1,tol2,u,v,w
real fu,fv,fw,fx,x
c
c c is the squared inverse of the golden ratio
c
c = 0.5*(3. - sqrt(5.0))
c
c eps is approximately the square root of the relative machine
c precision.
c
eps = 1.00
10 eps = eps/2.00
tol1 = 1.0 + eps
if (tol1 .gt. 1.00) 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.00*(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