%LDL^T factorization of a symmetric matrix A
%uses naive Gaussian elimination to build L

A = [4,2,-2;2,10,14;-2,14,42]

[n,n] = size(A);      %A is nxn (A must be square)
L = eye(n);           %start with L = nxn identity matrix

%forward elimination on U = A

for k=1:n-1
    for i=k+1:n
        c = A(i,k)/A(k,k);   %c = multiple of row k to subtract from row i
        L(i,k) = c;
        for j=k:n
            A(i,j) = A(i,j)-c*A(k,j);
        end
    end
end

D = zeros(n,n);
for i=1:n
    D(i,i) = A(i,i);
end

disp('the L and D matrices are...')
L
D
