 function [L,d] = ldlt0 (A) 
%---------------------------------------------
%  function [L,U] = lufactor (A) 
%  performs LU factorization of a matrix A 
%---------------------------------------------
 n = size(A,1) ;
 L = eye(n);
 for k=1:n-1 
     for i=k+1:n
       piv =  A(k,i) / A(k,k) ;
       L(i,k)   = piv; 
       A(i,i:n) = A(i,i:n) - piv*A(k,i:n);
   end
 end 
 d = diag(A);
