function G = CholBlockDiag(U)
% U is n-by-r and n is an integral multiple of r
% G is an m-by-1 cell array and G{k} = G(idx,idx) where
% GG' = I + UU' is the Cholesky factorization and idx = 1+(k-1)*r:k*r

[n,r] = size(U);
m = n/r;
Ir = eye(r,r);
G = cell(m,1);
% G{1}G{1}' = Ir + U(1:r,:)U(1:r,:)' = Ir + F where
F = zeros(r,r);
G{1} = chol(Ir+U(1:r,:)*U(1:r,:)','lower');
for k=2:m
idx = 1+(k-1)*r:k*r;
% Consider
% Iq+U1*U1' U1*U2' G0 0 G0' G21'
% = *
% U2*U1' Ir+U2*U2' G21 Gk 0 Gk'
% where
% U1 = U(1:q,:)
% U2 = U(idx,:)
% idx = 1+(k-1)*r:k*r;
% q = (k-1)*r
% Note that
% Gk*Gk' = Ir + U2*U2' - G21*G21'
% = Ir + U2*U2' - (G0\(U1*U2'))' *(G0\(U1*U2'))
% = Ir + U2*U2' - (U1*U2')'*inv(Iq+U1*U1')*U1*U2'
% = Ir + U2*U2' - U2*(U1'*inv(Iq+U1*U1')*U1)*U2'
% = Ir + U2*M*U2'
% where
% M = Ir - U1'*inv(Iq+U1*U1')*U1
%
% From Sherman-Morrison,
% M = Ir - U1'*(Iq - U1*inv(Ir + U1'*U1)*U1')U1
% = Ir - U1'*U1 + (U1'*U1)*inv(Ir+U1'*U1)*(U1'*U1)
% Note that if
F = F + U(idx-r,:)'*U(idx-r,:); %= U(1:(k-1)*r,:)'*U(1:(k-1)*r,:)= U1'*U1

% then
% M = Ir - F + F*inv(Ir + F)*F
% = Ir - F + inv(Ir+F)*F*F
% = inv(Ir+F)*( (Ir - F*F + F*F) = inv(Ir+F)
%
% So if
LM = chol(Ir+F,'lower');
% then
% Gk*Gk' = Ir + U2*inv(LM*LM')*U2'
% = Ir + U2*inv(LM')*inv(LM)*U2' = Ir + H'*H where H = LM\U2

H = LM\U(idx,:)'; % = LM\U2'
G{k} = chol(Ir + H'*H,'lower');

% Each pass there is are two r-by-r choleskies and one r-by-r mat mult
end

% Work = O(qr^3) = O(nr^2)