function [L,d] = LDLT(A)
% LDLT implements the LDLT factorization for a symmetric matrix A, assumed
% to have an LU factorization.
%
% INPUT: A (Symmetric matrix to be factorized).
%
% OUTPUT: L (Unit lower triangular matrix).
%         d (Diagonal of D).
%
%
n = length(A);
L = eye(n);
d = zeros(1,n);
d(1) = A(1,1);

for j = 1:n-1    
    L(j+1:n,j) = A(j+1:n,j);
    for k = 1:j-1
        L(j+1:n,j) = L(j+1:n,j) - L(j,k)*d(k)*L(j+1:n,k);
    end
    L(j+1:n,j) = L(j+1:n,j)/d(j);
    
    d(j+1) = A(j+1,j+1) - sum(d(1:j).*L(j+1,1:j).^2);    
end


