function LD = myldl(A)
% LD = myldl(A) computes the LDL^T factorization of the symmetric matrix A.
% The diagonal entries of LD correspond to the entries of the matrix D, the
% entries below the diagonal to the (non-trivial) entries of L.

% Check whether the matrix is symmetric
if ~isequal(A,A')
    error('The matrix is not symmetric');
end

% Initialize the output as a matrix of correct size containing only zeros.
LD = zeros(size(A));

% Initialize the diagonal separately from LD. This is done in order to
% avoid loops during the main iteration.
d = zeros(1,length(A));

for j=1:length(A)
    % Compute the diagonal entries.
    d(j) = A(j,j)-sum(d(1:j-1).*LD(j,1:j-1).^2);
    % If the j-th diagonal entry becomes 0, stop the iteration. An LDL^T
    % factorization most probably does not exists; if it exists, it will be
    % useless for solving linear systems.
    % Note that it might also be intelligent to produce a warning whenever
    % the diagonal entry is very small compared to the other entries of LD,
    % as this will probably entail stability issues.
    if(d(j)==0)
        error('Computation stopped because of one diagonal entry being equal to 0.');
    end
    % Compute the off-diagonal entries in column j.
    LD(j+1:end,j) = (A(j+1:end,j) - ... 
        LD(j+1:end,1:j-1)*(d(1:j-1).*LD(j,1:j-1))')/d(j);
end

% Insert the diagonal entries into the final matrix LD.
LD = LD + diag(d);