clear

// DONNEES 

// EXEMPLE 1 
n = 3;  // dimension de la matrice carrée
b = [1;2;1]; // second membre 
A = [2,3,1;-1,2,1;2,3,2]; // matrice 
 

// EXEMPLE 2 AVEC PIVOT PETIT 
//delta = 1.0E-16;
//n = 2;  // dimension de la matrice carrée
//b = [delta+1;2]; // second membre 
//A = [delta,1;1,1]; // matrice 


// solution avec la fonction backslash de scilab 
y = A\b



// FACTORISATION A=LU de A sans pivotage
function A=LUSansPivotage(A,n)
  for k=1:n-1 
        
    if (abs(A(k,k))<1.0E-17) then
        mprintf( 'k=%d, pivot Akk trop petit %e\n', k, A(k,k) )       
       abort    
    end     
    A(k+1:n,k)= A(k+1:n,k)/A(k,k); // mise a jour de L
    A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); // mise a jour de U
  end  
endfunction



function x=SolLUSansPivotage(A,n,b)
    
x = zeros(n,1)
z = zeros(n,1)
    
// Lz = b descente 

z(1) = b(1);
for i=2:n
    z(i) = b(i) - A(i,1:i-1)*z(1:i-1);
end

// Ux = z remontée

x(n) = z(n)/A(n,n);
for i=n-1:-1:1
    x(i) = (z(i) - A(i,i+1:n)*x(i+1:n))/A(i,i);
end        
endfunction



A=LUSansPivotage(A,n)
x=SolLUSansPivotage(A,n,b)

//erreur 
e=y-x

mprintf( 'Norme erreur : %e\n', norm(e) ) ;

disp('x',x)

disp('y',y)

disp('x-y',x-y)

  
  




