  function [x,rVals] = pcg2(A,b,kMax,G_M)
% x is the kMax iterate of the method of conjugate gradients applied to
% the nxn symmetric positive definite system Ax = b.
% This implementation is identical to cg0 except that it uses the Lanczos 
% iteration to generate the q-vectors and exploits the fact that the
% matrix Tk = Qk'*A*Qk is tridiagonal.
% rVals is a row vector with rVals(k) = norm(A*xk - b), k=1:kMax.
r = b;
z = G_M'\(G_M\r);            
rho_c = z'*r;
x = zeros(length(b),1);
for k=1:kMax
    if k==1
        p = z;
    else
        tau = rho_c/rho_minus;
        p = z + tau*p;
    end
    w = A*p;
    mu = rho_c/(p'*w);
    x = x + mu*p;
    r = r-mu*w;
    z = G_M'\(G_M\r); 
    rho_minus = rho_c;
    rho_c = z'*r; 
    rVals(k) = norm(A*x-b);
end
