% This program performs the preconditioned conjugate gradient method,
% by calling pcg(A,b,x,M), where x is the initial guess to the solution
% and M is the M-matrix of the splitting of A, the preconditioning
% matrix.
% 
% Written for MA221, by Ming Gu, Fall 2007
%

function [conv, x] = pcg(A,b,x,M);

r = b - A*x;
z = M\r;
p = z;
j = 1;
conv(1) = norm(r);
if b == 0
   err(1) = norm(x,inf);
end

while (norm(r)/norm(b) > 1.e-14) & (j < 500)
    alpha = (z'*r)/(p'*A*p);
    x = x + alpha*p;
    rnew = r - alpha*A*p;
    znew = M\rnew;
    beta = (znew'*rnew)/(z'*r);
    p = znew + beta*p;
    r = rnew;
    z = znew;
    j = j + 1;
    conv(j) = norm(r)/norm(b);
    if b == 0
       err(j) = norm(x,inf);
    end
    disp(sprintf('%5.0f %20.16e', j,conv(j)));
%    disp(sprintf('%5.0f %20.16e %20.16e', j,conv(j),err(j)));
end
subplot(2,1,1); plot((1:j),conv,'o',(1:j),conv,'--');
subplot(2,1,2); semilogy((1:j),conv,'o',(1:j),conv,'--');
%figure(2)
%plot((1:j),err,'o',(1:j),err,'--');
%semilogy((1:j),err,'o',(1:j),err,'--');
return