 function [L, d]  = ldlt (A) 
%---------------------------------------------
%  performs LDLT factorization of a matrix  A 
%---------------------------------------------
 n = size(A,1) ;
 for k=1:n-1 
   if (A(k,k)==0), error(' zero pivot'), end
   for i=k+1:n
      piv =  A(k,i)/A(k,k);
      A(i,i:n)=A(i,i:n)-piv*A(k,i:n);
  end 
 end
 d = diag(A);
 A  = (1 ./d) .* triu(A);
 L = A';
