function R = cholesky(A)
% R = cholesky(A) if A is a symmetric positive definite matrix
%         returns a upper triangular matrix R and a diagonal matrix D
%         such that  A = R1'*D*R1'
%         R1 has all diagonal entries equals 1
%         the values of D are returned on the diagonal of R

TOL=1e-10; %% there certainly are better test than this!!

[n,m] = size(A);
if (n ~= m)
  error('cholesky: matrix has to be square')
end%if

for k = 1:n-1
  if ( abs(A(k,k)) <= TOL) 
    error ('cholesky:might be a singular matrix')
  end%if
  for j = k+1:n
     A(j,k) = A(k,j)/A(k,k);
     % row operations only
     A(j,j:n) = A(j,j:n) - A(j,k)*A(k,j:n);
  end%for
end%for
if ( abs(A(n,n)) <= TOL) 
  error ('cholesky:might be a singular matrix')
end%if

% return the lower triangular part of A. 
% Transpose it to obtain an upper triangular matrix
R = tril(A)';
