function [L,U,X] = L1UX(A,B)     // resolution + 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
    // on peut appeler [L,U]= L1UX(A) si on ne veut que la decomposition

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

// B est recopié dans la colonne n+1 de A
A(:,n+1) = B
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,X='pas de calcul'
        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+1) = A(i,k+1:n+1) - A(i,k)*A(k,k+1:n+1)/A(k,k) 
   // end,
   // ou
   A(k+1:n,k+1:n+1) = A(k+1:n,k+1:n+1) - A(k+1:n,k)*A(k,k+1:n+1)/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),

// remontée, calcul des xi  dans A(:,n+1)
for i=[n:-1:1],   // remontee : calcul des xi dans A(i,n+1)
    A(i,n+1) = (A(i,n+1)  - A(i,i+1:n)*A(i+1:n,n+1))/A(i,i),
end,
printf('solutions'),
// print(6,A(:,n+1)'),
X=A(:,n+1)

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;3,3,-5;4,5,-2], B=[8;14;16]       premier exemple du cours
// A=[2,1,-4;4,2,-7;2,1,-1], B=[8;15;9]        resol impossible
// A=[2,1,-4;4,2,-7;2,1,-1], B=[8;15;5]        resol indeterminee
// A=[2,1,-4;4,2,-7;2,2,-1], B=[8;15;5]        echange de lignes obligatoire
//                                             pas de LU mais PLU

// 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]