function retval = wolfelinesrch(fngrad,fn,x,p)
% usage: alpha = wolfelinesrch(fngrad,fn,x,p)
% description: returns the value alpha in the interval
% [a,b] at which strong Wolfe conditions are satisfied.
% Passed parameters:
% fngrad: gradient function name
% fn: function name
% x: starting point
% p: search direction

% shorthand: phi(alpha) = f(x + alpha*p) so that
% phip(alpha) = gradf(x+alpha*p)'*p.

% global variables:
global parms;

% local variables:
% rho: growth factor
% c1: Wolfe condition 1 parameter
% c2: Wolfe condition 2 parameter
% alpha1: initial step length
% alphamax: maximum allowable step size
% alphaold,alphanew: successive approxns to alpha
% phi0,phip0: phi(0), phip(0)
% phiold, phinew: phi(alpha*)
% phipold,phipnew: phip(alpha*)
% psi1,psi2: convenient constants
% doneflag: done flag

% no bulletproofing here
rho = parms.wolferho;
c1 = parms.c1;
c2 = parms.c2;
alpha1 = parms.alpha1;
alphamax = parms.alphamax;
alphaold = 0;
alphanew = alpha1;
phi0 = feval(fn,x);
phip0 = feval(fngrad,x)'*p;
psi1 = c1*phip0;
psi2 = -c2*phip0;
doneflag = 0;
phinew = phi0;
phipnew = phip0;

while (~doneflag)
  phiold = phinew;
  phipold = phipnew;
  phinew = feval(fn,x+alphanew*p);
  phipnew = feval(fngrad,x+alphanew*p)'*p;
  % test if new phi value is above Wolfe1 line
  if (phinew > phi0 + alphanew*psi1) | ((phinew >= phiold) & (alphanew ~= alpha1))
    retval = zoom(alphaold, alphanew, phi0, psi1, psi2, fn, fngrad, x, p);
    doneflag = 1;
  else
    if (abs(phipnew) <= psi2) % passes Wolfe2?
      retval = alphanew;
      doneflag = 1;
    else
      if (phipnew >= 0) % positive slope on right?
        retval = zoom(alphanew, alphaold, phi0, psi1, psi2, fn, fngrad, x, p);
        doneflag = 1;
      else
        alphaold = alphanew;
        alphanew = rho*alphanew;
        if (alphanew >= alphamax) | (alphanew <= 100*eps)
          retval = alphamax;
          doneflag = 1;  
        end
      end
    end  
  end
end
end % end this main function so that locals are hidden from function below
%
function alphaopt = zoom(alphalo, alphahi, phi0, psi1, psi2, fn, fngrad, x, p)
% Passed parameters:
% alphalo: alpha satisfying suff decrease and smallest fn value so far
% alphahi: chosen so that phip(alphalo)*(alphahi-alphalo) < 0
% phi0: as in main function
% psi1: as in main function
% psi2: as in main function
% fn: as in main function
% fngrad: as in main function
% x: as in main function
% p: as in main function
% local variables:
% doneflag: done flag
% philo,phihi,phitest: phi(*)
% phiplo,phiphi,phiptest: phip(*)
% alphatest,alphaopt: test and optimal values of alpha
% Initialize
doneflag = 0;
while (~doneflag)
  philo = feval(fn,x+alphalo*p);
  phiplo = feval(fngrad,x+alphalo*p)'*p;
  phihi = feval(fn,x+alphahi*p);
  phiphi = feval(fngrad,x+alphahi*p)'*p;
  if (alphalo < alphahi) % must do cubicmin in right order
    alphatest = cubicmin(alphalo,alphahi, philo, phihi, phiplo, phiphi);
  else
    alphatest = cubicmin(alphahi, alphalo, phihi, philo, phiphi, phiplo);
  end
  if (alphatest == alphalo) || (alphatest == alphahi) % bad news..min not interior
    disp('Warning: numeric failure in line search');
    alphaopt = alphahi;
    return;
  end
  phitest = feval(fn,x+alphatest*p);
  if ((phitest > phi0 + alphatest*psi1) | (phitest >= philo) )
    alphahi = alphatest;
  else
    phiptest = feval(fngrad,x+alphatest*p)'*p;
    if (abs(phiptest) <= psi2) % pass Wolfe2?
      alphaopt = alphatest;
      doneflag = 1;
    else
      if (phiptest*(alphahi-alphalo) >= 0) % keep entry condition
        alphahi = alphalo;
      end
      alphalo = alphatest;
    end
  end
end
end