function lambdap = parab3p ( lambdac, lambdam, ff0, ffc, ffm )

% Apply three-point safeguarded parabolic model for a line search.
%
% C. T. Kelley, June 29, 1994
%
% function lambdap = parab3p(lambdac, lambdam, ff0, ffc, ffm)
%
%  Disclaimer:
%
%    This code comes with no guarantee or warranty of any kind.
%
%  Author:
%
%    Tim Kelley
%
%  Reference:
%
%    Tim Kelley,
%    Iterative Methods for Linear and Nonlinear Equations,
%    SIAM, 2004,
%    ISBN: 0898713528,
%    LC: QA297.8.K45.
%
% input:
%       lambdac = current steplength
%       lambdam = previous steplength
%       ff0 = value of \| F(x_c) \|^2
%       ffc = value of \| F(x_c + \lambdac d) \|^2
%       ffm = value of \| F(x_c + \lambdam d) \|^2
%
% output:
%       lambdap = new value of lambda given parabolic model
%
% internal parameters:
%       sigma0 = .1, sigma1=.5, safeguarding bounds for the linesearch
%

%
% set internal parameters
%
sigma0=.1; sigma1=.5;
%
% compute coefficients of interpolation polynomial
%
% p(lambda) = ff0 + (c1 lambda + c2 lambda^2)/d1
%
% d1 = (lambdac - lambdam)*lambdac*lambdam < 0
%      so if c2 > 0 we have negative curvature and default to
%      lambdap = sigam1 * lambda
%
c2 = lambdam*(ffc-ff0)-lambdac*(ffm-ff0);
if c2 >= 0
    lambdap = sigma1*lambdac; return
end
c1=lambdac*lambdac*(ffm-ff0)-lambdam*lambdam*(ffc-ff0);
lambdap=-c1*.5/c2;
if (lambdap < sigma0*lambdac) lambdap=sigma0*lambdac; end
if (lambdap > sigma1*lambdac) lambdap=sigma1*lambdac; end
