% conjgrad.m  solve Ax=b for x by conjugate gradient method (initial x 0 OK)
%                   A must be symmetric positive definite for convergence
% added file output of debug
function [x] = conjgrad(A, b, x, tol, maxit)
  debug = 2;
  r = b - A * x;
  p = r;
  rsold = r' * r;

  for i = 1:maxit
    Ap = A * p;
    alpha = rsold / (p' * Ap);
    x = x + alpha * p;
    r = r - alpha * Ap;
    rsnew = r' * r;
    if(debug > 1)
      disp('conjgrad.m iter, rsnew =')
      disp(i)
      disp(rsnew)
    end % if
    if sqrt(rsnew) < tol
      break;
    end % if
    p = r + (rsnew / rsold) * p;
    rsold = rsnew;
  end % end for
  if(debug > 0)
      disp('conjgrad final iter, rsnew =')
      disp(i)
      disp(rsnew)
  end % if
end % conjgrad
