module wolfe

implicit none
contains

!==================================================================================================
subroutine wolfe_search2(x,p,alpha0,c1,c2,f0,df0,alpha,fval,dfval,nfe,nge)

real,intent(in)            :: f,df(:),x(:),p(:),alpha0,c1,c2,f0,df0(:)
real, intent(out)          :: alpha,fval,dfval
integer, intent(out)       :: nfe, nge

! internal variables
real                       :: slope0,slope,xplus(size(x)),old_alpha,old_fval,old_dfval,old_slope
real                       :: temp,alpha_lo,alpha_hi,f_lo,f_hi,df_lo,slope_lo,a,b,dalpha,c
integer                    :: lsiter,iter_cnt

! [ALPHA,FVAL,DFVAL,NFE,NGE] = wolfe_search2(F,DF,PARAMS,X,P,ALPHA0,C1,C2,F0,DF0,TRACE)
! Performs search for point satisfying the two strong Wolfe conditions:
!
!	f(x+alpha.p) <= f(x) + c1.alpha p'.grad f(0)		(WC1)
!	|p'.grad f(x+alpha.p)| <= c2|p'.grad f'(x)|		(WC2b)
!
! ALPHA0 is the initial value of alpha in the search.  If ALPHA0
!	satisfies the Wolfe conditions, then the search stops there. 
! F0 is f(x) at the initial x.
! DF0 is grad f(x) at the initial x.
! DF(PARAMS,X) is the gradient of F(PARAMS,X) with
!	respect to X
! Note that PARAMS is the parameter for F: F(PARAMS,X)
! If TRACE is non-zero then print out information about the process
!
! The algorithm follows that of Wright & Nocedal "Numerical Optimization",
! pp. 60-62, section 3.4
!
! Returns
!	ALPHA	-- step length parameter
!	FVAL	-- function value at end
!	DFVAL	-- gradient value at end
!	NFE	-- # function evaluations
!	NGE	-- # gradient evaluations

!
! Check parameters
!
alpha = 0;
fval = f0;
dfval = df0;
nfe = 0;
nge = 0;

if (not(0<c1.and.c1<c2.and.c2<1)) then
  write(0,*) 'wolfe_search2: Error: Need 0 < c1 < c2 < 1'
  write(0,*) 'wolfe_search2: c1 = ',c1,' c2 = ', c2
  return
end if

! Initialization

slope0 = dot_product(p,df0)
if (slope0 >= 0) then
  write(0,*) 'wolfe_search2: Error: Need a descent direction'
  write(0,*) 'wolfe_search2: p*grad f(x) = ',slope0
  return
end

!
! Bracketing phase
!
write(0,*) 'wolfe_search2: Bracketing phase: alpha = ', alpha0
alpha = alpha0
old_alpha = 0
!!! f0  = feval(f, params,x);
!!! df0 = feval(df,params,x);
old_fval = f0
old_dfval = df0
old_slope = slope0
write(0,*) 'wolfe_search2: f(x) = ',f0,' p*grad f(x) = ',slope0

! Main loop
!!firsttrip = 1;
lsiter=1
do while (lsiter < ls_iter ) !!! forever do...
  xplus = x+alpha*p

  fval = feval(f,params,xplus)

  nfe = nfe + 1
  write(0,*) 'wolfe_search2: alpha = ',alpha,' f(x+alpha*p) = ',fval
  if (fval > f0 + c1*alpha*slope0 .or. fval >= old_fval ) then
     write(0,*) 'wolfe_search2: (WC1) failed or f increased'
     exit;
  end if
  write(0,*) 'wolfe_search2: (WC1) holds & f decreased'

  dfval = feval(df,params,xplus);

  nge = nge + 1;
  slope = p*dfval;
  write(0,*) 'wolfe_search2: p''.grad f(x+alpha*p) = ', slope
  if ( abs(slope) <= c2*abs(slope0) ) then
     write(0,*)  'wolfe_search2: (WC2) holds'
    return
  end if
  if ( slope >= 0 ) then
    write(0,*) 'wolfe_search2: f''(alpha) >= 0'
    exit
  end

  ! Update variables -- note no upper limit on alpha
  temp = alpha
  ! alpha = alpha + old_alpha
  alpha = 2*alpha
  old_alpha = temp
  old_fval = fval
  old_slope = slope
  lsiter=lsiter+1
end do

write(0,*) 'wolfe_search2: Entering ''zoom'' phase.'

!
! "Zoom" phase
!
alpha_lo = old_alpha;
alpha_hi =     alpha;
f_lo     = old_fval;
f_hi     =     fval;
df_lo    = old_dfval;
slope_lo = old_slope;
!!df_hi   =     dfval;
!!slope_hi=     slope;

write(0,*) 'wolfe_search2: zoom phase: alpha_lo = ',alpha_lo,' alpha_hi = ',alpha_hi
iter_cnt = 0
lsiter=1
do while (lsiter < ls_iter ) !!! forever do...

  ! form quadratic interpolant of function values on alpha_lo, alpha_hi
  ! and the derivative at alpha_lo...
  ! and find min of interpolant within interval [alpha_lo,alpha_hi]
  a = f_lo
  b = slope_lo
  dalpha = alpha_hi-alpha_lo;
  c = (f_hi - f_lo - dalpha*slope_lo)/dalpha**2
  ! iter_cnt
  if ( c <= 0 .or. mod(iter_cnt,3) == 2 ) then
    ! Use bisection
    alpha = alpha_lo + 0.5*dalpha
    write(0,*) 'wolfe_search2: using bisection: c = ', c
  else
    ! Use min of quadratic
    alpha = alpha_lo - 0.5*b/c
    write(0,*) 'wolfe_search2: using quadratic: a = ', a,' b = ',b, ' c = ', c
  end if
  write(0,*) 'wolfe_search2: alpha = ', alpha
  end

  ! main part of loop
  xplus = x + alpha*p

  fval = feval(f,params,xplus);

  nfe = nfe + 1;
  write(0,*) 'wolfe_search2: f(x+alpha*p) = ', fval
  if (fval > f0 + c1*alpha*slope0 .or. fval >= f_lo ) then
     write(0,*) ('wolfe_search2: zoom: (WC1) fails or f increased')
     write(0,*) ('wolfe_search2: zoom: update alpha_hi')
     alpha_hi = alpha;
     f_hi   = fval;
    ! df_hi = feval(df,params,xplus);
  else
    write(0,*) 'wolfe_search2: zoom: (WC1) holds and f decreased'

    dfval = feval(df,params,xplus);

    nge = nge + 1
    slope = p*dfval
    write(0,*) 'wolfe_search2: p''.grad f(x+alpha*p) = ', slope
    if (abs(slope) <= c2*abs(slope0)) then
       write(0,*) 'wolfe_search2: zoom: (WC2b) holds'
       return
    end if
    if ( slope*dalpha >= 0 ) then
       write(0,*) 'wolfe_search2: zoom: alpha_hi <- alpha_lo & update alpha_lo'
       alpha_hi = alpha_lo;
       f_hi     = f_lo;
       !!! df_hi  = df_lo;
       !!! slope_hi= slope_lo;
       alpha_lo = alpha;
       f_lo     = fval;
       df_lo    = dfval;
       slope_lo = slope;
    else
       write(0,*) 'wolfe_search2: zoom: update alpha_lo'
       alpha_lo = alpha;
       f_lo     = fval;
       df_lo    = dfval;
       slope_lo = slope;
    end if
  end if
  ! Update iteration count
  iter_cnt = iter_cnt + 1;
  lsiter=lsiter+1
end do
end subroutine
end module
