function [L,ierr] = Chol(A);
   %
   %
   %  Input:   A  an symmetric n by n matrix.
   %       
   %  Output:  L  a lower triangular matrix with A = LL' 
   %
   %        ierr  an error flag.  ierr = k >0 means failure at step k
   %                                   = 0    means normal exit 
   %
   %
   %
   %  D.C. Sorensen 
   %  25 Aug 03
   %-----------------------------------------------------------------
   %     
      [n,n] = size(A);
      ierr = 0;
   %
      for k = 1:n, 
   %
   %    exit if A is not positive definite
   %
        if (A(k,k) <= 0), ierr = k; return; end
   %
   %    Compute main diagonal elt. and then scale the k-th column 
   %

         A(k,k) = sqrt(A(k,k));
         A(k+1:n,k) = A(k+1:n,k)/A(k,k);
   %
   %     Update lower triangle of the trailing (n-k) by (n-k) block
   %
   
         for j = k+1:n,

                A(j:n,j) = A(j:n,j) - A(j:n,k)*A(j,k);
         end
      end
      L = tril(A);