dbrent_ Function

  function dbrent_(ax,bx,cx,func,fjac,tol,xmin,itmax) result(dbrent)
    real(8),intent(in)  :: ax,bx,cx,tol
    real(8),intent(out) :: xmin
    real(8)             :: dbrent
    integer,optional    :: itmax
    integer             :: itmax_=100
    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), intent(in) :: x
         real(8)             :: func
       end function func
       function fjac(x)
         real(8), intent(in) :: x
         real(8)             :: fjac
       end function fjac
    end interface
    itmax_=100;if(present(itmax))itmax_=itmax    
    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_