 function [U, d]  = ldltU (A)
% function [U, d]  =  ldlt (A)
% LDLT factorization of A . basic version. Gives
% the U matrix and the diagonal
%--------------------------------------------
 n = size(A,1) ;
 for k=1:n-1 
    for i=k+1:n
      piv = A(k,i) / A(k,k) ;
%%-------------------- operate only on U-part
      A(i,i:n)=A(i,i:n)-piv*A(k,i:n);
    end
 end 
%%-------------------- return U and diagonal
 d = diag(A);
 U = triu(A);
%% M = diag(d) \ triu(A);    L = M';

