function P2Grade()
% Checks out DiagChol

% Generate an example and solve it using the O(n^3) method...
n = 600;
u = randn(n,1);
tau = rand(1);
A = eye(n,n) + tau*u*u';
G = chol(A,'lower');
dTrue = diag(G);
% Now compare it to the quick method for computing d..
d = DiagChol(u,tau);
clc
fprintf('|| dTrue - DiagChol(u,tau) || = %10.2e\n',norm(dTrue - d))



  function d = DiagCholCVL(u,tau)
% Suppose u is a  column n-vector and tau >= 0
% I + tau*u*u' = G*G' is the Cholesky factorization of I + tau*u*u'.
% d = diag(G), i.e., d(i) = G(i,i), i=1:n

n = length(u);
d = zeros(n,1);
%
% Comparing blocks in
%
%    d(1)  0    d(1) v'                      1+tau*u(1)*u(1)    tau*u(1)*u(2:n)'
%             *           = I+ tau*u*u'  = 
%    v    G1     0   G1'                     tau*u(1)*u(2:n)    I + tau*u(2:n)*u(2:n)'  
% we see that
%            
%          d(1) = sqrt(1+tau*u(1)^2)
%             v = (tau*u(1)/d(1))*u(2:n) 
%        G1*G1' =  I + tau*u(2:n)*u(2:n)' - v*v'
%               =  I + (tau/(1+tau*u(1)^2)*u(2:n)*u(2:n)
%               =  I + tau1*u(2:n)*u(2:n)'
%
%  where tau1 = (tau/(1+tau*u(1)^2).
%
%  Since we want d = [d(1);diag(G1)] we repeat this process on the reduced
%  problem G1*G1' = I + tau1*u(2:n)*u(2:n)' 

for k=1:n
    d(k) = sqrt(1 + tau*u(k)^2);
    tau = tau/(1+tau*u(k)^2);      % The "next" tau
end

% -2 recursive implementations involve O(n^2) storage and dont work
% if n is bigger than the recursion limit
    