function [L,U] = LUdiag(A,D)     // decomposition LU
                                // avec D la diagonale de U imposée
                                // on ne doit pas avoir de pivot nul

n=size(A,'c') // calcul de la taille de A et vérif carrée
if ~ size(A,'r')==n  then printf('matrice non carree'), return,end

print(6,A)

for k = 1 : n ,   // etape k
    if A(k,k)==0 then 
        printf('pivot nul à l''étape %d',k),
        L='pas de calcul',
        U='',
        return
    end,
  
   printf('pivot=%f',A(k,k)),

   racpivot = sqrt(A(k,k))
   A(k,k+1:n) = A(k,k+1:n)*D(k)/A(k,k) // nouvelle ligne k
   A(k:n,k) = A(k:n,k)/D(k), // nouvelle colonne k
   
   // for i = k+1 : n ,   // nouvelle ligne i
   //     A(i,k+1:n) = A(i,k+1:n) - A(i,k)*A(k,k+1:n)) 
   // end,
   // ou
   A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n)
   
   printf('etape %d',k),
   print(6,A),
end,

L=zeros(n,n)
for i=1:n,L(i,1:i)=A(i,1:i), end
U=diag(D)
for i=1:n,U(i,i+1:n)=A(i,i+1:n), end

// A=[4 5 3;3 4 2;2 3 -1] 
// D= [10 100 1000]