%Cholesky 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

%first we find the LU factorization via 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

%now we modify L so that A = LL^T
%columns of L are multiplied by square roots of diagonal entries of U

for j=1:n
    c = sqrt(A(j,j));
    for i=j:n
        L(i,j) = L(i,j)*c;
    end
end

disp('the L matrix is...')
L
