c To get d1mach, mail netlib c send d1mach from core double precision function fmin(ax,bx,f,tol) double precision 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 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.) c c output.. c c fmin abcissa approximating the point where f attains a c minimum 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 double precision a,b,c,d,e,eps,xm,p,q,r,tol1,t2,u,v,w,fu,fv,fw, 2 fx,x,tol3 double precision dabs,dsqrt,d1mach c c c is the squared inverse of the golden ratio c=0.5d0*(3.0d0-dsqrt(5.0d0)) c c eps is approximately the square root of the relative machine c precision. c 10 eps=d1mach(4) tol1=eps+1.0d0 eps=dsqrt(eps) c a=ax b=bx v=a+c*(b-a) w=v x=v e=0.0d0 fx=f(x) fv=fx fw=fx tol3=tol/3.0d0 c c main loop starts here c 20 xm=0.5d0*(a+b) tol1=eps*dabs(x)+tol3 t2=2.0d0*tol1 c c check stopping criterion c if (dabs(x-xm).le.(t2-0.5d0*(b-a))) go to 190 p=0.0d0 q=0.0d0 r=0.0d0 if (dabs(e).le.tol1) go to 50 c c fit parabola c r=(x-w)*(fx-fv) q=(x-v)*(fx-fw) p=(x-v)*q-(x-w)*r q=2.0d0*(q-r) if (q.le.0.0d0) go to 30 p=-p go to 40 30 q=-q 40 r=e e=d 50 if ((dabs(p).ge.dabs(0.5d0*q*r)).or.(p.le.q*(a-x)) 2 .or.(p.ge.q*(b-x))) go to 60 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).ge.t2).and.((b-u).ge.t2)) go to 90 d=tol1 if (x.ge.xm) d=-d go to 90 c c a golden-section step c 60 if (x.ge.xm) go to 70 e=b-x go to 80 70 e=a-x 80 d=c*e c c f must not be evaluated too close to x c 90 if (dabs(d).lt.tol1) go to 100 u=x+d go to 120 100 if (d.le.0.0d0) go to 110 u=x+tol1 go to 120 110 u=x-tol1 120 fu=f(u) c c update a, b, v, w, and x c if (fx.gt.fu) go to 140 if (u.ge.x) go to 130 a=u go to 140 130 b=u 140 if (fu.gt.fx) go to 170 if (u.ge.x) go to 150 b=x go to 160 150 a=x 160 v=w fv=fw w=x fw=fx x=u fx=fu go to 20 170 if ((fu.gt.fw).and.(w.ne.x)) go to 180 v=w fv=fw w=u fw=fu go to 20 180 if ((fu.gt.fv).and.(v.ne.x).and.(v.ne.w)) go to 20 v=u fv=fu go to 20 c c end of main loop c 190 fmin=x return end