function [ x, err, total_iters ] = cg ( x0, b, atv, params )

%*****************************************************************************80
%
%% cg() is a conjugate gradient solver for linear systems.
%
%  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 x0(n): the initial solution estimate.
%
%    real b(n): the right hand side.
%
%    function ax = atv(x):, a matrix-vector product routine computing A*x.
%
%    real params(2): iteration parameters.
%    params(1): relative residual reduction factor.
%    params(2): kmax = maximum number of iterations.
%
%  Output:
%
%    real x(n): the estimated solution.
%
%    real err(*): vector of iteration residual norms.
%
%    integer total_iters: number of iterations.
%
  x = x0;
  n = length ( b );
  errtol = params(1) * norm ( b );
  kmax = params(2);
  rho = zeros(kmax+1,1);

  r = b - feval ( atv, x );

  k = 1;
  rho(k) = r' * r;
  zeta = norm ( r );
  total_iters = 1;
  err = [ zeta ];
%
%  Iteration
%
  while ( errtol < sqrt ( rho(k) ) && k < kmax + 1 )

    if ( k == 1 )
      p = r;
    else
      beta = rho(k) / rho(k-1);
      p = r + beta * p;
    end
    w = feval ( atv, p );
    alpha = rho(k) / ( p' * w );
    x = x + alpha * p;
    r = r - alpha * w;

    k = k + 1;
    rho(k) = r' * r;
    zeta = norm ( r );
    total_iters = total_iters + 1;
    err = [ err, zeta ];

  end

  return
end

