function [P,L,U,X] = PLU1X(A,B) // resolution + decomposition LU ou PLU
                                // avec U unitaire
                                // car on divise la ligne du pivot
                                // on ne doit pas avoir de pivot nul
    // on peut appeler [L,U]=PLU1X(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)
perm=[1:n] 

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

for k = 1 : n ,   // etape k
    if A(k,k)==0 then // recherche d'un pivot non nul
        printf('pivot nul %d',k),
        if k<n then i=k+1,
                    while(A(i,k)==0 & i<n ), i=i+1,
                    end,
               else i=k,
        end,
        if A(i,k)==0 then  L='pas de pivot non nul',U='',
                           X='pas de solution', 
                           P='', return
                     else printf('on echange les lignes %d et %d',k,i),
                          A([k,i],:)=A([i,k],:);
                          perm([k,i])=perm([i,k]);
                          print(6,A),
        end,
    end,

  
   printf('pivot=%f',A(k,k)),
   de=de*A(k,k),

   A(k,k+1:n+1) = A(k,k+1:n+1)/A(k,k), // nouvelle ligne 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) 
   // 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)

   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: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)),
end,

print(6,perm)

printf('solutions'),

X=A(:,n+1)
for i=1:n, P(perm(i),i)=1, end
L=zeros(n,n)
for i=1:n,L(i,1:i)=A(i,1:i), end
U=diag(ones(1,n))
for i=1:n,U(i,i+1:n)=A(i,i+1: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