function R = icholesky(A)
%ICHOLESKY computes the incomplete Cholesky factorization of
%the symmetric positive definite matrix A.
%
%   R = icholesky(A) returns an n x n upper triangular matrix R
%   that has the same distribution of nonzeros as A above the diagonal.

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

R = spalloc(n,n,floor(nnz(A)/2));
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)
		error('The matrix is not positive definite');
	end
	R(i,i) = sqrt(tmp);
	%now compute the off-diagonal entries
	for j = i+1:n
		if A(i,j) == 0
			R(i,j) = 0;
		else
			R(i,j) = (A(i,j) - sum(R(1:i-1,i).*R(1:i-1,j)))/R(i,i);
		end
	end
end


