function [L,D,P] = ldl(A)
% A is symmetric and positive semidefinite
% L is unit lower triangular
% D is diagonal with D(1,1)>=D(2,2)>=...>=D(n,n)>= 0
% P is a permutation.
% P*A*P' = L*D*L';

[n,n] = size(A);
P = eye(n,n);
L = eye(n,n);
for k=1:n-1
    [m,j]  = max(diag(A(k:n,k:n)));
    j = j+k-1;
    P([k j],:) = P([j k],:);
    L([k j],1:k-1) = L([j k],1:k-1);
    A([k j],k:n) = A([j k],k:n);
    A(k:n,[k j]) = A(k:n,[j k]);
    L(k+1:n,k) = A(k+1:n,k)/A(k,k);
    D(k,k) = A(k,k);
    A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - L(k+1:n,k)*D(k,k)*  L(k+1:n,k)';
end
D(n,n) = A(n,n);