function [L,U] = L1U(A,B)     // decomposition LU
                                  // avec L unitaire
                                  // on ne divise pas la ligne du pivot
                                  // mais on doit diviser la colonne du pivot
                                  // 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'), end

print(6,A),

de=1, // pour calculer le determinant en meme temps

for k = 1 : n ,   // etape k
    if A(k,k)==0 then 
        printf('pivot nul %d',k),
        L='pas de decomposition LU',
        U=L,
        return
    end,
  
   printf('pivot=%f',A(k,k)),
   de=de*A(k,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)/A(k,k) 
   // 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)/A(k,k)
  
   A(k+1:n,k) = A(k+1:n,k)/A(k,k), // nouvelle colonne k
   
   printf('etape %d',k),
   print(6,A),
end,

printf('determinant=%f',de),

L=diag(ones(1,n))
for i=1:n,L(i,1:i-1)=A(i,1:i-1), end
U=zeros(n,n)
for i=1:n,U(i,i:n)=A(i,i:n), end


// A=[2,1,-4;4,2,1;2,1,1]         plusieurs decompositions LU
// A=[2,1,-4;4,2,1;2,2,1]         pas de decomposition LU

// A= [4 5 3;3 4 2;2 3 -1]
// A=[9 9 6;9 13 8;6 8 6]