  function [x,rVals] = cg0(A,b,kMax)
% x is the kMax iterate of the method of conjugate gradients applied to
% the nxn symmetric positive definite system Ax = b.
% This implementation is based on a ``literal'' interpretation of what
% the method does at each step. 
% rVals is a row vector with rVals(k) = norm(A*xk - b), k=1:kMax.

% The first q-vector...
beta0 = norm(b);
Qk = b/beta0;
for k=1:kMax
    % Let Qk be an nxk matrix with orthonormal columns whose range is the
    %    span of the vectors b, Ab, A^{2}*b,...,A^{k-1}b.
    % Compute the vector xk that minimizes 
    %
    %           phi(x) = (x'*A*x)/2 - x'*b
    %
    % over the range of Qk. This is equivalent to finding yk (kx1) so that 
    %
    %             phi(Qk*y)= y'(Qk'*A*Qk)y/2 - (Qk*y)'b 
    %
    % is minimized by setting y = yk. This vector is the solution to 
    % (Qk'*A*Qk)*y = Qk'*b. Thus
    Tk = Qk'*(A*Qk);
    yk = Tk\(Qk'*b);
    xk = Qk*yk;
    % Expand  span{Qk(:,1),...,Qk(:,k)} = ran(Qk) so that it includes the gradient 
    % of phi at xk....
    gk = A*xk - b;
    Qk = Extend(Qk,gk);   
    rVals(k) = norm(A*xk-b);
end
x = xk;