function [ x, error, total_iters, it_hist ] = pcg ( x0, b, atv, params, pcv )

%*****************************************************************************80
%
%% pcg() is a preconditioned Conjugate-Gradient 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 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
%
%    function px = PCV ( x ): MATLAB code that applies the preconditioner
%    If omitted, the identity is used.
%
%  Output:
%
%    real X(N): the estimated solution.
%
%    real ERROR(total_iters): iteration residual norms.
%
%    integer TOTAL_ITERS: the number of iterations taken.
%
%    real IT_HIST(N,TOTAL_ITERS): matrix of the sequence of solutions X.
%

%
%  Ensure input vectors are columnwise.
%
  x0 = x0(:);

  if ( nargout == 4 )
    it_hist = [];
  end

  n = length ( b );
  errtol = params(1);
  maxiters = params(2);
  error = [];
  x = x0;

  if ( nargout == 4 )
    it_hist = [ it_hist, x ];
  end

  r = b - feval ( atv, x );

  if ( nargin == 4 )
    z = r;
  else
    z = feval ( pcv, r );
  end

  rho = z' * r;
  tst = norm ( r );
  terminate = errtol * norm ( b );
  error = [ error, tst ];
  it = 1;

  while ( ( terminate < tst ) && ( it <= maxiters ) )

    if ( it == 1 ) 
      p = z;
    else
      beta = rho / rhoold;
      p = z + beta * p;
    end

    w = feval ( atv, p );
    alpha = p' * w;
%
%  Test here to make sure the linear transformation is positive definite.
%  A non-positive value of alpha is a very bad sign.
%
    if ( alpha <= 0.0 )
      [ alpha, rho, it ]
      error ( 'pcg: negative curvature' ) 
    end

    alpha = rho / alpha;
    x = x + alpha * p;

    if ( nargout == 4 )
      it_hist = [ it_hist, x ]; 
    end

    r = r - alpha * w;
    tst = norm ( r );
    rhoold = rho;

    if ( nargin == 4 )
      z = r;
    else
      z = feval ( pcv, r );
    end

    rho = z' * r;
    it = it + 1;
    error = [ error, tst ];
    total_iters = it - 1;

  end

  return
end

