 module lnsrch

implicit none
contains
!-------------------------------------------------------------------------
SUBROUTINE LNESRCH(N,X0,F0,G0,FP0,D0,ALPHA,FCNG,TOL,STEPMX,ACC,EXTBND,MAXFCN,NFCALL,IERR,X,F,G) 
IMPLICIT REAL*8(A-H,O-Z) 
DOUBLE PRECISION X0(N),G0(N),D0(N),X(N),G(N) 
EXTERNAL FCNG 
! 
!***** 
! 
! The problem is of dimension N. The current iterate is X0 
! with associated function value F0 and gradient G0. The search 
! direction is D0 and the directional derivative along it is FP0. 
! ALPHA defines the initial step along D0, for example, 1.0. 
! All of these quantities are input to the subroutine. FCNG is 
! an externally supplied subroutine with calling sequence 
! FCNG(N,x,f,g), which returns the function value f and gradient 
! vector g at any given point x. 
! 
! TOLX, STEPMX, ACC, and EXTBND are tolerances that are used in 
! the line search subroutine (see the relevant statements below). 
! In particular, 0 < ACC < 1 determines the accuracy that must be 
! achieved on exit. The new iterate is returned in X0 with 
! associated information in F0 and G0. 
! 
! No more than MAXFCN calls to the subroutine FCNG are allowed. 
! The actual number of calls are returned in NFCALL, and IERR is 
! a flag that specifies error information (if nonzero on exit). 
! 
! X, F, and G are work space for the line search routine. 
! 
! The cubic interpolation formulae used below are given in 
! Chapter 4. See, in particular, the third item of Section 4.3 
! and expression (4.3). The LNESRCH Fortran code is derived from 
! Hillstrom [1976a], which is, in turn, derived from a routine 
! in the Harwell library. 
! 
!***** 
! 
IERR=0 
NFCALL=0 
! 
100 TOT=0.0D0 
CDERIV=FP0 
PDERIV=FP0 
FTEST=F0 
! 
105 IF (ALPHA.LE.TOL) THEN 
IERR=2 
GOTO 150 
ENDIF 
! 
! begin the line search by incrementing the solution vector and 
! calculating the function and gradient at the incremented vector. 
! 
X=X0
G=G0 
X0=X0 + ALPHA*D0 
F=FTEST 
! 
CALL FCNG(N,X0,FTEST,G0) 
! 
NFCALL=NFCALL+1 
IF (NFCALL.GT.MAXFCN) THEN 
IERR=1 
RETURN 
ENDIF 
! 
! compute the directional derivative 
! 
FP0=0.0D0 
DO 110 I=1,N 
FP0=FP0 + G0(I)*D0(I) 
110 CONTINUE 
! 
! test whether the function has decreased 
! 
IF (FTEST.GE.F) GOTO 120 
!
! sufficient progress and convergence criterion for exit 
! 
IF (DABS(FP0/PDERIV).LE.ACC) GOTO 140 
! 
116 CONTINUE 
! 
! test for extrapolation or interpolation 
! 
IF (FP0.GT. 0.0D0) GOTO 120 
! 
! revise ALPHA by extrapolation, limiting the length of 
! the multiplier by EXTBND 
! 
TOT=TOT+ALPHA 
IF (TOT.GT.STEPMX) GOTO 145 
TEMP=EXTBND 
IF (CDERIV.LT.FP0) TEMP=FP0/(CDERIV-FP0) 
IF (TEMP.GT.EXTBND) TEMP=EXTBND 
F=FTEST 
CDERIV=FP0 
ALPHA=ALPHA*TEMP 
GOTO 105 
! 
! reset solution vector and revise ALPHA by cubic interpolation 
! 
120 CONTINUE 
DO 130 I=1,N 
X0(I)=X(I) 
G0(I)=G(I) 
130 CONTINUE 
! 
TEMP=3.0D0*(F-FTEST)/ALPHA + FP0 + CDERIV 
WT=DMAX1(DABS(TEMP),DABS(FP0),DABS(CDERIV)) 
WW=TEMP/WT 
WW=WW*WW - (CDERIV/WT)*(FP0/WT) 
IF (WW.LT.0.0D0) WW=0.0D0 
WW=DSQRT(WW)*WT 
TEMP=1 - (FP0 + WW - TEMP)/(2*WW + FP0 - CDERIV) 
ALPHA=ALPHA*TEMP 
FTEST=F 
GOTO 105 
! 
! ALPHA is accepted 
! 
140 F0=FTEST 
145 ALPHA=TOT+ALPHA 
C
150 CONTINUE 
RETURN 
! 
END SUBROUTINE

END MODULE

