function [ sol, it_hist, ierr ] = nsola ( x, f, tol, parms )

%*****************************************************************************80
%
%% nsola() is a globally convergent Newton-Krylov solver.
%
%  Discussion:
%
%    This algorithm uses an inexact Newton-Armijo iteration,
%    an Eisenstat-Walker forcing term, and parabolic line search 
%    via three point interpolation.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    14 April 2021
%
%  Author:
%
%    Tim Kelley
%
%  Reference:
%
%    Tim Kelley,
%    Iterative Methods for Linear and Nonlinear Equations,
%    SIAM, 2004,
%    ISBN: 0898713528,
%    LC: QA297.8.K45.
%
%  Input:
%
%    real X(N): the initial iterate.
%
%    function fx = F ( x ): MATLAB code that evaluates the function.
%
%    real TOL(2): [ ATOL, RTOL ] error tolerances.
%    ATOL: absolute error tolerance.
%    RTOL: relative error tolerance.
%
%    real PARMS(5): [ MAXIT, MAXITL, ETAMAX, LMETH, RESTART_LIMIT]
%    MAXIT: maximum number of nonlinear iterations,default = 40.
%    MAXITL: maximum number of inner iterations before restart in GMRES(m),
%    m=maxitl.  default = 40.  For iterative methods other than GMRES(m),
%    maxitl is the upper bound on linear iterations.
%    ETAMAX: Maximum error tolerance for residual in inner iteration. 
%    The inner iteration terminates when the relative linear residual is
%    smaller than eta*| F(x_c) |. eta is determined by the modified 
%    Eisenstat-Walker formula if etamax > 0.  If etamax < 0, then
%    eta = |etamax| for the entire iteration.  default: etamax=.9
%    LMETH: choice of linear iterative method
%    1: GMRES, default.
%    2: GMRES(m) 
%    3: BICGSTAB
%    4: TFQMR
%    RESTART_LIMIT: max number of restarts for GMRES if lmeth = 2; default=20.
%
%  Output:
%
%    real SOL(N): the estimated solution.
%
%    real IT_HIST(maxit,3): [ FNRM, FNUM, ARM] for each iteration.
%    FNRM: scaled l2 norms of nonlinear residuals;
%    FNUM: number of function evaluations;
%    ARM: number of steplength reductions
%
%    integer IERR: termination flag.
%    0: successful termination
%    1: the termination criterion is not satsified after maxit iterations,
%       or the ratio of successive nonlinear residuals exceeded 1 and the
%       iteration was stopped early.
%    2: the line search failed.  Too many step reductions.
%
%  Local:
%
%    logical DEBUG: if TRUE, prints iteration statistics during execution.
%
%    real ALPHA: parameter to measure sufficient decrease.
%
%    real SIGMA0, SIGMA1: safeguarding bounds for the linesearch.
%
%    integer MAXARM: maximum number of steplength reductions.
%
  debug = false;

  alpha = 1.0e-4; 
  sigma0 = 0.1;
  sigma1 = 0.5;
  maxarm = 50; 

  gamma = 0.9;
  ierr = 0;
  maxit = 40;
  lmaxit = 40;
  etamax = 0.9;
  it_histx = zeros ( maxit, 3 );
  lmeth = 1;
  restart_limit = 20;
%
%  initialize parameters for the iterative methods
%
  gmparms = [ abs(etamax), lmaxit ];
  if ( nargin == 4 )
    maxit = parms(1);
    lmaxit = parms(2);
    etamax = parms(3);
    gmparms = [ abs(etamax), lmaxit ];
    if ( 4 <= length ( parms ) )
      lmeth = parms(4);
    end
    if ( length ( parms ) == 5 )
      gmparms = [ abs(etamax), lmaxit, parms(5), 1 ];
    end
  end

  atol = tol(1);
  rtol = tol(2);

  n = length ( x );
  fnrm = 1.0;
  itc = 0;
%
%  Evaluate f at the initial iterate
%
  f0 = feval ( f, x );
  fnrm = norm ( f0 ) / sqrt ( n );
  it_histx(itc+1,1) = fnrm;
  it_histx(itc+1,2) = 0;
  it_histx(itc+1,3) = 0;
  fnrmo = 1.0;
  outstat(itc+1, :) = [ itc, fnrm, 0, 0, 0 ];
%
%  Compute the stop tolerance
%
  stop_tol = atol + rtol * fnrm;
%
%  Iteration.
%
  while ( stop_tol < fnrm & itc < maxit )
%
%  Track the ratio of successive residual norms.
%
    ratio = fnrm / fnrmo;
    fnrmo = fnrm; 
    itc = itc + 1;
    [ step, errstep, inner_it_count, inner_f_evals ] = ...
      fdkrylov ( f0, f, x, gmparms, lmeth );
%
%  The line search starts here.
%
    xold = x;
    lambda = 1.0;
    lamm = 1;
    lamc = lambda;
    iarm = 0;
    xt = x + lambda * step;
    ft = feval ( f, xt );
    nft = norm ( ft );
    nf0 = norm ( f0 );
    ff0 = nf0 * nf0;
    ffc = nft * nft;
    ffm = nft * nft;

    while ( ( 1.0 - alpha * lambda) * nf0 <= nft )
%
%  apply the three point parabolic model
%
      if ( iarm == 0 )
        lambda = sigma1 * lambda; 
      else
        lambda = parab3p ( lamc, lamm, ff0, ffc, ffm ); 
      end
%
% update x; keep the books on lambda
%
      xt = x + lambda * step;
      lamm = lamc;
      lamc = lambda;
%
%  keep the books on the function norms
%
      ft = feval ( f, xt );
      nft = norm ( ft );
      ffm = ffc;
      ffc = nft * nft;
      iarm = iarm + 1;
      if ( maxarm < iarm )
        disp ( 'nsola: Armijo failure, too many reductions ');
        ierr = 2;
        disp ( outstat )
        it_hist = it_histx(1:itc+1,:);
        sol = xold;
        return;
      end

    end

    x = xt;
    f0 = ft;
%
%   end of line search
%
    fnrm = norm ( f0 ) / sqrt ( n );
    it_histx(itc+1,1) = fnrm; 
%
%   How many function evaluations did this iteration require?
%
    it_histx(itc+1,2) = it_histx(itc,2) + inner_f_evals + iarm + 1;
    if ( itc == 1 )
      it_histx(itc+1,2) = it_histx(itc+1,2) + 1;
    end;
    it_histx(itc+1,3) = iarm;

    ratio = fnrm / fnrmo;
%
%  adjust eta as per Eisenstat-Walker
%
    if ( 0 < etamax )
      etaold = gmparms(1);
      etanew = gamma * ratio * ratio;
      if ( 0.1 <  gamma * etaold * etaold )
        etanew = max ( etanew, gamma * etaold * etaold );
      end
      gmparms(1) = min ( [ etanew, etamax ] );
      gmparms(1) = max ( gmparms(1), 0.5 * stop_tol / fnrm );
    end

    outstat(itc+1, :) = [ itc, fnrm, inner_it_count, ratio, iarm ];

    if ( debug )
      disp ( outstat(itc+1,:) )
    end

  end

  sol = x;
  it_hist = it_histx(1:itc+1,:);
  if ( debug )
    disp ( outstat )
    it_hist = it_histx(1:itc+1,:);
  end
%
%  on failure, set the error flag
%
  if ( stop_tol < fnrm )
    ierr = 1;
  end

  return
end

