function [ x, error, total_iters ] = fdgmres ( f0, f, xc, params, xinit )

%*****************************************************************************80
%
%% fdgmres() implements GMRES for a Newton-GMRES solver.
%
%  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 F0(N): the function value at the current point.
%
%    function fx = F ( x ): MATLAB code that evaluates the function.
%
%    real XC(N): the current point.
%
%    real PARAMS(3):
%    params(1): relative residual reduction factor
%    params(2): max number of iterations
%    params(3): reorthogonalization method
%    1: Brown/Hindmarsh condition (default)
%    2: Never reorthogonalize (not recommended)
%    3: Always reorthogonalize (not cheap!)
%
%    real xinit(N): the initial iterate.
%
%  Output:
%
%    real X(N): the estimated solution.
%
%    real ERROR(total_iters): residual norms for the iteration.
%
%    integer TOTAL_ITERS: the number of iterations taken.
%
  errtol = params(1);
  kmax = params(2);
  reorth = 1;
  if ( length ( params ) == 3 )
    reorth = params(3);
  end
%
%  right side of linear equation for the step is -f0 if the
%  default initial iterate is used
%
  b = - f0;
  n = length ( b );
%
%  Use zero vector as initial iterate for Newton step unless
%  the calling routine has a better idea (useful for GMRES(m)).
%
  x = zeros(n,1);
  r = b;
  if ( nargin == 5 )
    x = xinit;
    r = - dirder ( xc, x, f, f0 ) - f0;
  end

  h = zeros(kmax,kmax);
  v = zeros(n,kmax);
  c = zeros(kmax+1,1);
  s = zeros(kmax+1,1);
  rho = norm ( r );
  g = rho * eye ( kmax+1, 1 );
  errtol = errtol * norm ( b );
  error = [];
%
%  test for termination on entry
%
  error = [ error, rho ];
  total_iters = 0;
  if ( rho < errtol ) 
%   disp(' early termination ')
    return
  end

  v(:,1) = r / rho;
  beta = rho;
  k = 0;
%
%  Iteration.
%
  while ( ( errtol < rho ) && ( k < kmax ) )

    k = k + 1;
%
%  Call directional derivative function
%
    v(:,k+1) = dirder ( xc, v(:,k), f, f0 );
    normav = norm ( v(:,k+1) );
%
%  Modified Gram-Schmidt
%
    for j = 1 : k
      h(j,k) = v(:,j)' * v(:,k+1);
      v(:,k+1) = v(:,k+1) - h(j,k) * v(:,j);
    end
    h(k+1,k) = norm ( v(:,k+1) );
    normav2 = h(k+1,k);
%
%  Reorthogonalize?
%
    if ( ( reorth == 1 && normav + 0.001 * normav2 == normav) || reorth ==  3 )
      for j = 1 : k
        hr = v(:,j)' * v(:,k+1);
	      h(j,k) = h(j,k) + hr;
        v(:,k+1) = v(:,k+1) - hr * v(:,j);
      end
      h(k+1,k) = norm ( v(:,k+1) );
    end
%
%  Watch out for happy breakdown
%
    if ( h(k+1,k) ~= 0.0 )
      v(:,k+1) = v(:,k+1) / h(k+1,k);
    end
%
%  Form and store the information for the new Givens rotation
%
    if ( 1 < k )
      h(1:k,k) = givapp ( c(1:k-1), s(1:k-1), h(1:k,k), k - 1 );
    end
%
%  Don't divide by zero if solution has  been found
%
    nu = norm ( h(k:k+1,k) );

    if ( nu ~= 0.0 )
      c(k) = conj ( h(k,k) / nu );
      s(k) = - h(k+1,k) / nu;
      h(k,k) = c(k) * h(k,k) - s(k) * h(k+1,k);
      h(k+1,k) = 0.0;
      g(k:k+1) = givapp ( c(k), s(k), g(k:k+1), 1 );
    end
%
%  Update the residual norm
%
    rho = abs ( g(k+1) );
    error = [ error, rho ];

  end

  return
end

