function [x,flag,relres,iter]=mypcg(A,b,tol,maxit,Mfun)
%Textbook page 133
%Mfun(r) solves M*z=r (i.e., compute z=M\r=M^{-1}r)
n=length(b);
x=zeros(n,1);
r=b-A*x; r0=norm(r);
d=Mfun(r); z=d;
flag=1;relres=[];
for iter=1:maxit
relres=[relres;norm(r)/r0]; %store relative residuals
if(norm(r)/r0)<=tol
flag=0;break;
end
Ad=A*d;
alf=r'*z/(d'*Ad);
x=x+alf*d;
r_old=r; r=r_old-alf*Ad;
z_old=z; z=Mfun(r);
beta=r'*z/(r_old'*z_old);
d=z+beta*d;
end