% cholesky.m
% Calculates the choleksy decompostion of a symmetric positive definite
% matrix A, s.t. A = L * L'
%
% input:  A - s.p.d. matrix
% output: L - lower triangular matrix, s.t. A = L*L'
%
% author: Pascal F. Heiter
% date:   2013/12/03

function L = cholesky(A)

%get problem dimension and allocate memory
n = size(A,1);
L = zeros(n,n);

%k = 1
L(1,1) = sqrt(A(1,1));
L(2:n,1) = 1/L(1,1) * A(2:n,1);

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