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

%*****************************************************************************80
%
%% nsolgm() is a locally convergent Newton-GMRES solver.
%
%  Discussion:
%
%    This algorithm uses an Eisenstat-Walker forcing term.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 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(3): [maxit, maxitl, etamax], iteration control parameters:
%    maxit: maximum number of nonlinear iterations, default = 40
%    maxitl: maximum number of inner iterations, default = 40
%    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
%
%  Output:
%
%    real SOL(N): the estimated solution.
%
%    real IT_HIST(maxit,3): [ fnorm, fnum, arm ] for each iteration.
%    FNORM:  scaled l2 norms of nonlinear residuals
%    FNUM: function evaluations,
%    ARM: number of steplength reductions
%
%    integer IERR: termination code.
%    0: successful termination
%    1: the error tolerance was not acheived after maxit iterations,
%       or the iteration was terminated early because the ratio of 
%       successive nonlinear residuals exceeded 1.
%
%  Local:
%
%    logical DEBUG: if TRUE, prints iteration statistics during execution.
%
  debug = true;

  gamma = 0.9;
  ierr = 0;
  maxit = 40;
  lmaxit = 40;
  etamax = 0.9;
  it_histx = zeros ( maxit, 2 );

  if ( nargin == 4 )
    maxit = parms(1);
    lmaxit = parms(2);
    etamax = parms(3);
  end
  gmparms = [ abs(etamax), lmaxit ];
  rtol = tol(2);
  atol = tol(1);
  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; 
  fnrmo = 1.0;
  outstat(itc+1, :) = [ itc, fnrm, 0.0, 0.0 ];
%
%  Compute the stop tolerance
%
  stop_tol = atol + rtol * fnrm;
%
%  Iteration.
%
  while ( stop_tol < fnrm && itc < maxit )
%
% keep track of the ratio (rat = fnrm/frnmo)
% of successive residual norms and 
% the iteration counter (itc)
%
    rat = fnrm / fnrmo;
    fnrmo = fnrm; 
    itc = itc + 1;
%
%  Compute the step using a special GMRES routine.
%
    [ step, errstep, inner_it_count ] = fdgmres ( f0, f, x, gmparms ); 
    xold = x;
    x = x + step;
    f0 = feval ( f, x );
    fnrm = norm ( f0 ) / sqrt ( n );
    it_histx(itc+1,1) = fnrm;
    rat = fnrm / fnrmo;
    outstat(itc+1, :) = [ itc, fnrm, inner_it_count, rat ];

    if ( debug )
      disp ( outstat(itc+1,:) )
    end
%
%  How many function evaluations did this iteration require?
%
    it_histx(itc+1,2) = it_histx(itc,2) + inner_it_count + 1;
    if ( itc == 1 ) 
      it_histx(itc+1,2) = it_histx(itc+1,2) + 1;
    end;
%
%  adjust eta 
%
    if ( 0.0 < etamax )
      etaold = gmparms(1);
      etanew = gamma * rat * rat;
      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

  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

