function R=Cholesky(A)
% CHOLESKY computes the Cholesky decomposition of a matrix
%   R=Cholesky(A) computes the Cholesky decomposition A=R'R

n=length(A);
for j=1:n,
  v=A(j,j:n);
  if j>1,
    v=A(j,j:n)-R(1:j-1,j)'*R(1:j-1,j:n);
  end; 
  if v(1)<=0
    error('Matrix is not positive definite')
  else 
    h=1/sqrt(v(1));
  end
  R(j,j:n)=v*h;
end