  function [x,rVals] = cg1(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 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.
Tk = sparse(1,1);
beta0 = norm(b);
Qk = b/beta0;
beta = zeros(kMax,1);
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
    % Set up Tk Qk'*A*Qk using Lanczos recursions...
    alfa(k) = Qk(:,k)'*A*Qk(:,k);
    if k==1
        rk = A*Qk(:,k) - alfa(k)*Qk(:,k);
        Tk(k,k) = alfa(k);
    else
        rk = A*Qk(:,k) - alfa(k)*Qk(:,k) - beta(k-1)*Qk(:,k-1);
        Tk(k,k)   = alfa(k);
        Tk(k,k-1) = beta(k-1);
        Tk(k-1,k) = beta(k-1);
    end
    % To solve the tridiagonal system Tk*y = Qk'*b we note that the right
    % hand side has a special form since Qk(:,1) is a multiple of b and
    % orthogonal to Qk(:,2),...,Qk(:,k)....
    yk = Tk\[beta0;zeros(k-1,1)];
    xk = Qk*yk;
    rVals(k) = norm(A*xk-b);
    % Generate the next q-vector and augment it to Qk...
    if k<kMax
        beta(k) = norm(rk);
        Qk = [Qk rk/beta(k)];
    end
end
x = xk;



