% This program performs the conjugate gradient method,
% by calling pcg(A,b,x), where x is the initial guess to the solution
% 
% Written for MA128B, by Ming Gu, Spring 2010
%

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

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

while (norm(r)/norm(b) > 1.e-14) & (j < 500)
    alpha = (r'*r)/(p'*A*p);
    x = x + alpha*p;
    rnew = r - alpha*A*p;
    beta = (rnew'*rnew)/(r'*r);
    p = rnew + beta*p;
    r = rnew;
    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
