function [ step, errstep, total_iters, f_evals ] = fdkrylov ( f0, f, x, ...
  params, lmeth )

%*****************************************************************************80
%
%% fdkrylov() implements a Krylov solver for use in nsola().
%
%  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 F0(N): function value at current point.
%
%    function fx = F ( x ): MATLAB code that evaluates the function.
%
%    real X(N): the current point.
%
%    real PARAMS(4): [ERRTOL,LMAXIT,RESTART,REORTH], iteration parameters.
%    ERRTOL: relative residual reduction factor
%    LMAXIT: maximum number of iterations
%    RESTART: maximum number of restarts for GMRES(m)
%    REORTH: reorthogonalization method in GMRES
%    1: Brown/Hindmarsh condition (default)
%    2: Never reorthogonalize (not recommended)
%    3: Always reorthogonalize (not cheap!)
%
%    integer LMETH: method choice
%    1: GMRES without restarts (default)
%    2: GMRES(m), m=params(2) and the maximum number of restarts is params(3) 
%    3: Bi-CGSTAB
%    4: TFQMR
%
%  Output:
%
%    real X(N): the estimated solution.
%
%    real ERRSTEP(total_iters): residual norms for the iteration.
%
%    integer TOTAL_ITERS: the number of iterations taken.
%
  lmaxit = params(2);

  if ( lmeth == 1 )
    restart_limit = 0;
  elseif ( 3 <= length ( params ) )
    restart_limit = params(3);
  else
    restart_limit = 20;
  end
%
%  default reorthogonalization
%
  if ( length ( params ) == 3 )
    gmparms = [ params(1), params(2), 1 ];
%
%  reorthogonalization method is params(4)
%
  elseif ( length ( params ) == 4 )
    gmparms = [ params(1), params(2), params(4) ];
  else
    gmparms = [ params(1), params(2) ];
  end
%
%  Linear iterative methods
%
  if ( lmeth == 1 || lmeth == 2 ) 
%
%  Compute the step using a special GMRES routine.
%
    [ step, errstep, total_iters ] = fdgmres ( f0, f, x, gmparms );
    kinn = 0;
%
%  Restart at most restart_limit times
%
    while ( total_iters == lmaxit && ...
          gmparms(1) * norm(f0) < errstep(total_iters) && ...
          kinn < restart_limit )
      kinn = kinn + 1;
      [ step, errstep, total_iters ] = fdgmres ( f0, f, x, gmparms, step );
    end
    total_iters = total_iters + kinn * lmaxit;
    f_evals = total_iters + kinn;
%
% Bi-CGSTAB
%
  elseif ( lmeth == 3 )
    [ step, errstep, total_iters ] = fdcgstab ( f0, f, x, gmparms );
    f_evals = 2 * total_iters;
%
% TFQMR
%
  elseif ( lmeth == 4 ) 
    [ step, errstep, total_iters ] = fdtfqmr ( f0, f, x, gmparms );
    f_evals = 2 * total_iters;
  else
    error ( 'fdkrylov: Illegal value of LMETH.' )
  end

  return
end

