 function [L,U] = gpd(A) 
%---------------------------------------------
% LU factorization of A using a *dense*  GP algorithm
%---------------------------------------------
 n = size(A,1) ;
 L = eye(n);
 U = [];
 for j=1:n
   w = L\A(:,j);
   U(:,j)  = [w(1:j); zeros(n-j,1)];
   L(j+1:n,j) = w(j+1:n)/w(j);
 end 
 %%  U(:,n)  = L\A(:,n); 
