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

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