% cg.m - Simple CG iteration for 10-digit accuracy.
%        Remember, A must be symmetric positive definite!
 
  function x = cg(A,b)
  N = length(A);
  x = zeros(N,1); r = b; p = r;
  while norm(r) > 1e-10

     Ap = A*p;
     rr = r'*r;
     a  = rr/(p'*Ap);    % step length
     x  = x + a*p;       % approximate solution
     r  = r - a*Ap;      % residual
     c  = r'*r/rr;       % improvement this step
     p  = r + c*p;       % search direction

     fprintf('%20.15f\n',norm(r))
  end