function L = cholesky(A)

n=size(A,1);
L=zeros(size(A));

%*** test if A is symmetric
if max(max(abs(A-A')))>0
  disp('Matrix is not symmetric')
  return
end
%*** test 2nd and 3rd criterion
if max(diag(A))<max(max(abs(A))) || sum(diag(A)>0)<n
  disp('2nd and 3rd criterion not valid')
  return
end

for k=1:n
  tmp = A(k,k)- sum(L(k,1:k-1).^2);
  if tmp<=0
    disp('Matrix is not p.d.')
    L=zeros(size(A));
    return;
  end
  L(k,k) = sqrt(tmp);
  for i=k+1:n
    L(i,k) = 1/L(k,k)*(A(i,k)-sum(L(i,1:k-1).*L(k,1:k-1)));
  end
end
