  function P2()
% Checks out MaxDiag

% Small example checked against the naive method
n = 100; r = 3; d = rand(n,1); C = randn(n,r);
mu0 = max(diag(inv(diag(d)+C*C')));
mu = MaxDiag(d,C);
clc
fprintf('n = %1d, r = %1d, error = %10.3e\n\n',n,r,abs(mu-mu0))

% A large example to check efficiency

n = 100000; r = 10; d = rand(n,1); C = randn(n,r);
tic
muCVL = MaxDiagCVL(d,C);
tCVL = toc;
tic
muYOU = MaxDiag(d,C);
tYOU = toc;

fprintf('n = %1d, r = %1d, tYOU/tCVL = %6.2f,  error = %10.3e\n\n',n,r,tYOU/tCVL,abs(mu-mu0))

  function mu = MaxDiagCVL(d,C)
% d is column n vectors with positive entries and C is n-by-r with r<<n.
% mu is the largest entry along the diagonal of inv(diag(d) + C*C').

[n,r] = size(C);
% Sherman-Morrison-Woodbury says
%           inv(D + C*C') = inv(D) - inv(D)*C*inv(I+C'*inv(D)C)*C'*inv(D)
% 
% If

d1 = 1./sqrt(d);

% and we overwtite C with diag(d1)*C

C = C.*(d1*ones(1,r));

% then we want the max diagonal entry of
%               diag(d1)*(I - C*inv(M)*C')*diag(d1)
% where
M = eye(r,r)+C'*C;

% If
G = chol(M,'lower');

% and we overwrite C' with G\C'
C = (G\C')';

% then we need the maximum diagonal entry of
%               diag(d1)*(I - C*C')*diag(d1)

mu = max((1-sum(C'.^2))./d');


%  tYOU/tCVL  \approx 5     9/10
%  tYOU/tCVL  \approx 20    8/10
%  tYOU/tCVL  \approx 50    7/10

% Outof memory trying to do sherman morrison woodbury: 5/10
% Out of memory : 3/10