function R = cholesky(S)

% function R = cholesky(S)
% Cholesky factorization S = R'*R 
% S symmetric positive definite, R upper triangular 

n = size(S,1); 
R = zeros(size(S));
R(1,1) = sqrt(S(1,1));

for j = 1:n-1;
   I = [1:j-1]; K = [j+1:n];
   R(j,K) = (S(j,K) - R(I,j)'*R(I,K))/R(j,j); 
   R(j+1,j+1) = sqrt(S(j+1,j+1)-sum(R(1:j,j+1).^2));
end;

