function [ x, error, total_iters ] = gmresb ( x0, b, atv, params )

%*****************************************************************************80
%
%% gmresb() implements brute-force GMRES solver.
%
%  Discussion:
%
%    This version uses the modified Gram-Schmidt orthogonalization.
%
%  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 X0(N): the initial iterate.
%
%    real B(N): the right hand side.
%
%    function ax = ATV ( X ): returns A * X.
%
%    real PARAMS(2): iteration control parameters.
%    params(1) = relative residual reduction factor
%    params(2) = max number of iterations
%
%  Output:
%
%    real X(N): the estimated solution.
%
%    real ERROR(total_iters): iteration residual norms.
%
%    integer TOTAL_ITERS: the number of iterations taken.
%
  n = length ( b );
  errtol = params(1);
  kmax = params(2);
  x = x0;
  h = zeros(kmax,kmax);
  v = zeros(n,kmax);
  if ( norm ( x ) ~= 0.0 )
    r = b - feval ( atv, x );
  else
    r = b;
  end
  rho = norm ( r );
  errtol = errtol * norm ( b );
  error = [];
%
%  test for termination on entry
%
  error = [ error, rho ];
  total_iters = 0;
  if ( rho < errtol ) 
    return
  end
  v(:,1) = r / rho;
  beta = rho;
  k = 0;
%
%  Iteration
%
  while ( ( errtol < rho ) && ( k < kmax ) )

    k = k + 1;
    v(:,k+1) = feval ( atv, v(:,k) );
%
%  Modified Gram-Schmidt
%
    for j = 1 : k
      h(j,k) = v(:,k+1)' * v(:,j);
      v(:,k+1) = v(:,k+1) - h(j,k) * v(:,j);
    end
    h(k+1,k) = norm ( v(:,k+1) );
%
%   Watch out for happy breakdown.
%
    if ( h(k+1,k) ~= 0.0 )
      v(:,k+1) = v(:,k+1) / h(k+1,k);
    end
%
%  Formulate and solve the least squares problem.
%  Update the residual approximation.
%
    y = h(1:k+1,1:k) \ ( beta * eye(k+1,1) );
    rho = norm ( beta * eye(k+1,1) - h(1:k+1,1:k) * y );
    error = [ error, rho ];
  end
%
%  At this point either k > kmax or rho < errtol.
%  It's time to compute x and leave.
%
  total_iters = k;
  x = x + v(1:n,1:k) * y;

  return
end

