function [C,k] = cholesky(A)
% [C,k] = cholesky(A)
% for symmetric positive definite A: 
%   find C upper triangular such that A=C'*C
%   if matrix is positive definite: k=0
%   if matrix is not positive definite: return column j where algorithm breaks down,
%      and current matrix C. Note that C'*C gives elements A(i,j) for i<k-1 or j<k-1
n = size(A,1);
C = zeros(n,n);
for k=1:n
  v = C(1:k-1,k);
  s = A(k,k) - v'*v;
  if s<=0 
    warning(sprintf('matrix is not positive definite\nalgorithm breaks down in column %g with s=%g ',k,s))
    return
  end
  C(k,k) = sqrt(s);
  C(k,k+1:n) = (A(k,k+1:n)-v'*C(1:k-1,k+1:n))/C(k,k);
end 
k=0;


