subroutine lnsrch(xold,fold,g,p,x,f,stpmax,check,func)
real(8), dimension(:), intent(in) :: xold,g
real(8), dimension(:), intent(inout) :: p
real(8), intent(in) :: fold,stpmax
real(8), dimension(:), intent(out) :: x
real(8), intent(out) :: f
logical, intent(out) :: check
interface
function func(x)
real(8) :: func
real(8), dimension(:), intent(in):: x
end function func
end interface
real(8), parameter :: alf=1.0e-4,tolx=epsilon(x)
integer :: ndum
real(8) :: a,alam,alam2,alamin,b,disc,f2,pabs,rhs1,rhs2,slope,tmplam
ndum=assert_eq4(size(g),size(p),size(x),size(xold),'lnsrch')
check=.false.
pabs=vabs(p(:))
if (pabs > stpmax) p(:)=p(:)*stpmax/pabs
slope=dot_product(g,p)
if (slope >= 0.0) then
print*,'roundoff problem in lnsrch'
stop
endif
alamin=tolx/maxval(abs(p(:))/max(abs(xold(:)),1.0d0))
alam=1.0
do
x(:)=xold(:)+alam*p(:)
f=func(x)
if (alam < alamin) then
x(:)=xold(:)
check=.true.
return
else if (f <= fold+alf*alam*slope) then
return
else
if (alam == 1.0) then
tmplam=-slope/(2.0d0*(f-fold-slope))
else
rhs1=f-fold-alam*slope
rhs2=f2-fold-alam2*slope
a=(rhs1/alam**2-rhs2/alam2**2)/(alam-alam2)
b=(-alam2*rhs1/alam**2+alam*rhs2/alam2**2)/&
(alam-alam2)
if (a == 0.0) then
tmplam=-slope/(2.0d0*b)
else
disc=b*b-3.0d0*a*slope
if (disc < 0.0) then
tmplam=0.5d0*alam
else if (b <= 0.0) then
tmplam=(-b+sqrt(disc))/(3.0d0*a)
else
tmplam=-slope/(b+sqrt(disc))
end if
end if
if (tmplam > 0.5d0*alam) tmplam=0.5d0*alam
end if
end if
alam2=alam
f2=f
alam=max(tmplam,0.1d0*alam)
end do
end subroutine lnsrch