function L = Cholesky(A)
% Cholesky implements the Cholesky factorization for a symmetrc
% positive definite matrix A.
%
% INPUT: A (SPD matrix to be factorized).
%
% OUTPUT: L (Lower triangular matrix).
%
%
n = length(A);
L = zeros(n);

L(1,1) = sqrt(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)*L(j+1:n,k);
    end
    L(j+1:n,j) = L(j+1:n,j)/L(j,j);
    
    L(j+1,j+1) = sqrt(A(j+1,j+1) - sum(L(j+1,1:j).^2));    
end


