function [R, f] = mychol(A)
%MYCHOL My own version of Cholesky factorization.
%   R = mychol(A), for a positive definite matrix A,
%   is an lower triangular matrix R so that R*R' = A.
%   f, flag; f=1 if is s.p.d, f=0 if not

if ~isequal(A',A)
   error('Matrix must be symmetric or Hermitian')
end
f=1;
[n,n] = size(A);
R = zeros(n,n);
for k = 1:n
   i = (1:k-1)';
   j = k+1:n;
   s = A(k,k) - R(i,k)'*R(i,k);
   if s > 0
      R(k,k) = sqrt(s);
   else
       f=0;
      error('Matrix must be positive definite')
   end
   R(k,j) = (A(k,j) - R(i,k)'*R(i,j))/R(k,k);
end
R=R';