function x = CG(A,b)
% Solves Linear Equations by Conjugate Gradient Method for Symmetric matrix
%               H. Hasegawa; July 2, 2006
s=size(A); x = rand(s(1),1); k= 0; res = 1; r = b-A*x; p = r;
while ( res > 1.0e-8)
   k = k+1;
   rtr = r'*r; Ap = A*p;
   alpha = rtr/(p'*Ap);
   x = x+alpha*p;
   r = r-alpha*Ap; R(:,k) = r;
   res = norm(b-A*x,2);
   beta = (r'*r)/rtr;
   p = r+beta*p;   P(:,k) = p;
   out = [k res]; 
   disp(sprintf('%5d %15.3e',out))
end
R'*R
P'*(A*P)