function [X,R,P,alpha,beta]=CG(A,b,x0,m);
% CG conjugate gradient method
%   [X,R,P,alpha,beta]=CG(A,b,x0,m) computes an approximation for the
%   solution of the linear system Ax=b performing m steps of the
%   conjugate gradient method starting with x0, and returns in the
%   matrix X the iterates, in the matrix R the residuals, in P the
%   conjugate directions and the coefficients alpha and beta.

x=x0; r=b-A*x; p=r;
R =r; P=p; X=x;
oldrho=r'*r;
for k=1:m
  Ap=A*p;
  alpha(k)=oldrho/(p'*Ap);
  x=x+alpha(k)*p; r=r-alpha(k)*Ap;
  rho=r'*r;
  beta(k)=rho/oldrho; oldrho=rho;
  p=r+beta(k)*p;
  X=[X,x]; R=[R,r]; P=[P,p];
end