matlab code

function [x,flag,relres,iter]=mypcgfun(Afun,b,tol,maxit,Mfun)

%Textbook page 133

%Mfun(r) solves M*z=r (i.e., compute z=M\r=M^{-1}r)

if(nargin<5)

Mfun=@(r) r; %no preconditioner

end

n=length(b);

x=zeros(n,1); %initial guess

r=b-Afun(x); r0=norm(r);

d=Mfun(r); z=d;

flag=1;relres=[];

for iter=1:maxit

Ad=Afun(d);

alf=r'*z/(d'*Ad);

x=x+alf*d;

r_old=r; r=r_old-alf*Ad;

relres=[relres;norm(r)/r0]; %store relative residuals

if(norm(r)/r0)<=tol

flag=0;break;

end

z_old=z; z=Mfun(r);

beta=r'*z/(r_old'*z_old);

d=z+beta*d;

end