dbrent_optimize Function

function dbrent_optimize(ax,bx,cx,func,fjac,tol,itmax,xmin) result(dbrent)
  real(8),intent(in)  :: ax,bx,cx,tol
  real(8),intent(out) :: xmin
  real(8)             :: dbrent
  integer             :: itmax
  real(8), parameter  :: zeps=1.d-3*epsilon(ax)
  integer             :: iter
  real(8)             :: a,b,d,d1,d2,du,dv,dw,dx,e,fu,fv,fw,fx,olde,tol1,tol2
  real(8)             :: u,u1,u2,v,w,x,xm
  logical             :: ok1,ok2
  interface
     function func(x)
       real(8) :: x
       real(8) :: func
     end function func
     function fjac(x)
       real(8) :: x
       real(8) :: fjac
     end function fjac
  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
  dx=fjac(x)
  dv=dx
  dw=dx
  do iter=1,ITMAX
     xm=0.5d0*(a+b)
     tol1=tol*abs(x)+ZEPS
     tol2=2.0d0*tol1
     if (abs(x-xm) <= (tol2-0.5d0*(b-a))) exit
     if (abs(e) > tol1) then
        d1=2.0d0*(b-a)
        d2=d1
        if (dw /= dx) d1=(w-x)*dx/(dx-dw)
        if (dv /= dx) d2=(v-x)*dx/(dx-dv)
        u1=x+d1
        u2=x+d2
        ok1=((a-u1)*(u1-b) > 0.d0) .and. (dx*d1 <= 0.d0)
        ok2=((a-u2)*(u2-b) > 0.d0) .and. (dx*d2 <= 0.d0)
        olde=e
        e=d
        if (ok1 .or. ok2) then
           if (ok1 .and. ok2) then
              d=merge(d1,d2, abs(d1) < abs(d2))
           else
              d=merge(d1,d2,ok1)
           end if
           if (abs(d) <= abs(0.5d0*olde)) then
              u=x+d
              if (u-a < tol2 .or. b-u < tol2) &
                   d=sign(tol1,xm-x)
           else
              e=merge(a,b, dx >= 0.d0)-x
              d=0.5d0*e
           end if
        else
           e=merge(a,b, dx >= 0.d0)-x
           d=0.5d0*e
        end if
     else
        e=merge(a,b, dx >= 0.d0)-x
        d=0.5d0*e
     end if
     if (abs(d) >= tol1) then
        u=x+d
        fu=func(u)
     else
        u=x+sign(tol1,d)
        fu=func(u)
        if (fu > fx) exit
     end if
     du=fjac(u)
     if (fu <= fx) then
        if (u >= x) then
           a=x
        else
           b=x
        end if
        call mov3(v,fv,dv,w,fw,dw)
        call mov3(w,fw,dw,x,fx,dx)
        call mov3(x,fx,dx,u,fu,du)
     else
        if (u < x) then
           a=u
        else
           b=u
        end if
        if (fu <= fw .or. w == x) then
           call mov3(v,fv,dv,w,fw,dw)
           call mov3(w,fw,dw,u,fu,du)
        else if (fu <= fv .or. v == x .or. v == w) then
           call mov3(v,fv,dv,u,fu,du)
        end if
     end if
  end do
  if (iter > ITMAX) stop 'dbrent: exceeded maximum iterations'
  xmin=x
  dbrent=fx
contains
  !bl
  subroutine mov3(a,b,c,d,e,f)
    real(8), intent(in) :: d,e,f
    real(8), intent(out) :: a,b,c
    a=d
    b=e
    c=f
  end subroutine mov3
end function dbrent_optimize