brent_optimize Function

function brent_optimize(ax,bx,cx,func,tol,itmax,xmin)
  real(8), intent(in)  :: ax,bx,cx,tol
  real(8), intent(out) :: xmin
  real(8)              :: brent_optimize
  integer              :: itmax
  real(8), parameter   :: cgold=0.3819660d0,zeps=1.0d-3*epsilon(ax)
  integer              :: iter
  real(8)              :: a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
  interface
     function func(x)
       real(8) :: x
       real(8) :: func
     end function func
  end interface
  a=min(ax,cx)
  b=max(ax,cx)
  v=bx
  w=v
  x=v
  e=0.d0
  fx=func(x)
  fv=fx
  fw=fx
  do iter=1,itmax
     xm=0.5d0*(a+b)
     tol1=tol*abs(x)+zeps
     tol2=2.0*tol1
     if (abs(x-xm) <= (tol2-0.5d0*(b-a))) then
        xmin=x
        brent_optimize=fx
        return
     end if
     if (abs(e) > tol1) then
        r=(x-w)*(fx-fv)
        q=(x-v)*(fx-fw)
        p=(x-v)*q-(x-w)*r
        q=2.d0*(q-r)
        if (q > 0.d0) p=-p
        q=abs(q)
        etemp=e
        e=d
        if (abs(p) >= abs(0.5d0*q*etemp) .or. &
             p <= q*(a-x) .or. p >= q*(b-x)) then
           e=merge(a-x,b-x, x >= xm )
           d=cgold*e
        else
           d=p/q
           u=x+d
           if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
        end if
     else
        e=merge(a-x,b-x, x >= xm )
        d=cgold*e
     end if
     u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
     fu=func(u)
     if (fu <= fx) then
        if (u >= x) then
           a=x
        else
           b=x
        end if
        call shft(v,w,x,u)
        call shft(fv,fw,fx,fu)
     else
        if (u < x) then
           a=u
        else
           b=u
        end if
        if (fu <= fw .or. w == x) then
           v=w
           fv=fw
           w=u
           fw=fu
        else if (fu <= fv .or. v == x .or. v == w) then
           v=u
           fv=fu
        end if
     end if
  end do
  !pause 'brent: exceed maximum iterations'
contains
  subroutine shft(a,b,c,d)
    real(8), intent(out) :: a
    real(8), intent(inout) :: b,c
    real(8), intent(in) :: d
    a=b
    b=c
    c=d
  end subroutine shft
end function brent_optimize