function X = pcg(M,A,b,n)
% usage: X = pcg(M,A,b,n)
% description: given SPD mxm A and preconditioner M,
% and mx1 b != 0, this function calculates n steps
% of conjugate gradients, and returns mxn X containing
% the first n approximates to solution of A*x = b.
% We do not count the 0th, which is x^(0) = 0.

% local variables:
% ii: index variables

m = length(b);
X = zeros(m,n);
x = zeros(m,1);
r = b;
z = M\r;% in practice, we would *never* do this
p = z;
for jj = 1:n
    rold = r;
    zold = z;
    alpha = rold'*zold/(p'*A*p);
    x = x + alpha*p;
    r = rold - alpha*A*p;
    z = M\r; % see above
    bta = r'*z/(rold'*zold);
    p = z + bta*p;
    X(:,jj) = x;
end