function [G11,G21,P] = PosSemiDefFact(A,tol)
% A is an n-by-n symmetric positive semidefinite matrix.
% tol is a small positive number with the property that any number
% less than tol*norm(A,1) can be declared to equal zero
% G11 is r-by-r and lower triangular, G21 is (n-r)-by-r, and P is
% a permutation matrix such that
%
%     norm( PAP' -  [G11 0;G21 0]*[G11 0 ; G21 0]' ) = O(tol*norm(A,1))

[n,n] = size(A);
[L,D,P] = ldl(A);
r = sum(diag(D)>tol*norm(A,1));
for k=1:n
    L(:,k) = L(:,k)*sqrt(D(k,k));
end
d = diag(D);
m = max(max(abs(L)));
G11 = L(1:r,1:r);
G21 = L(r+1:n,1:r);

errL = norm(P*A*P' - L*L',1);
errG = norm(P*A*P' - L(:,1:r)*L(:,1:r)',1);