module oc_lnsrch_mod
  implicit none
contains


    subroutine oc_lnsrch(n,xold,fold,g,p,x,f,stpmax,check,ALF,TOLX,func)
      integer  :: n
      logical  :: check
      real     :: f,fold,stpmax,g(n),p(n),x(n),xold(n),func,ALF,TOLX
!      parameter (ALF=1.e-4,TOLX=1.e-7)
      external :: func
!     USES func
!        Given an n-dimensional point xold(1:n), the value of the function and gradient there,
!        fold and g(1:n), and a direction p(1:n), finds a new point x(1:n) along the direction
!        p from xold where the function func has decreased sufficiently. The new function value
!        is returned in f. stpmax is an input quantity that limits the length of the steps so that you
!        do not try to evaluate the function in regions where it is undefined or subject to overflow.
!        p is usually the Newton direction. The output quantity check is false on a normal exit.
!        It is true when x is too close to xold. In a minimization algorithm, this usually signals
!        convergence and can be ignored. However, in a zero-finding algorithm the calling program
!        should check whether the convergence is spurious.
!        Parameters: ALF ensures sufficient decrease in function valuse; TOLX is the convergence
!        criterion on dx 
      integer  :: i
      real     :: a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,&
                  test,tmplam

      check=.false.
      sum=0.
      do i=1,n
        sum=sum+p(i)*p(i)
      end do
      sum=sqrt(sum)
      if(sum.gt.stpmax)then
        do i=1,n
          p(i)=p(i)*stpmax/sum
        end do
      endif
      slope=0.
      do i=1,n
        slope=slope+g(i)*p(i)
      end do
      test=0.
      do i=1,n
        temp=abs(p(i))/max(abs(xold(i)),1.)
        if(temp.gt.test)test=temp
      end do
      alamin=TOLX/test
      alam=1.

      do while (.true.)
        do i=1,n
          x(i)=xold(i)+alam*p(i)
        end do
        f=func(x)
        if(alam.lt.alamin)then
          do i=1,n
            x(i)=xold(i)
          end do 
          check=.true.
          return
        else if(f.le.fold+ALF*alam*slope)then
          return
        else
          if(alam.eq.1.)then
            tmplam=-slope/(2.*(f-fold-slope))
          else
            rhs1=f-fold-alam*slope
            rhs2=f2-fold2-alam2*slope
            a=(rhs1/alam**2-rhs2/alam2**2)/(alam-alam2)
            b=(-alam2*rhs1/alam**2+alam*rhs2/alam2**2)/(alam-alam2)
            if(a.eq.0.)then
              tmplam=-slope/(2.*b)
            else
              disc=b*b-3.*a*slope
              if(disc.lt.0.) pause 'roundoff problem in lnsrch'
              tmplam=(-b+sqrt(disc))/(3.*a)
            endif
            if(tmplam.gt..5*alam)tmplam=.5*alam
          endif
        endif
        alam2=alam
        f2=f
        fold2=fold
        alam=max(tmplam,.1*alam)
      end do

    end subroutine
end module
