function R = cholesky(A)
%cholesky Cholesky factoriztion
%
%   R = cholesky(A) computes the Cholesky factorization
%   of an n x n positive definite matrix A, where R is an
%   upper triangular matrix such that A = R'R.
%   If A appears not to be positive definite, the
%   function prints an error message.

[m, n] = size(A);
if m ~= n
	error('The matrix is not square.');
end

R = zeros(n,n);
for i = 1:n
	% we need to take the square root of this expression
	tmp = A(i,i) - sum(R(1:i-1,i).^2);
	if (tmp <= 0)
		disp('The matrix is not positive definite.');
      R = [];
      return;
	end
	R(i,i) = sqrt(tmp);
	% now compute the off-diagonal entries
	for j = i+1:n
		R(i,j) = (A(i,j) - sum(R(1:i-1,i).*R(1:i-1,j)))/R(i,i);
	end
end
