//  Factorisation LU

  a=[1 1+0.5*10^(-15) 3; 2 2 20; 3 6 4];
// Comparer avec   a=[1 1+0.5*10^(-2) 3; 2 2 20; 3 6 4]
  
  b = a;
  
  n = size(a,1);
  for j = 1:n-1
    for i = j+1:n
      a(i,j) = a(i,j)/a(j,j);
      for k = j+1:n
        a(i,k) = a(i,k)-a(i,j)*a(j,k);
      end;
    end;
  end;
  
  l = eye(n,n);
  u = zeros(n,n);
  for i =2:n
    for j = 1:i-1
      l(i,j) = a(i,j);
    end
  end;
  for i = 1:n
    for j = i:n
      u(i,j) = a(i,j);
     end;
  end;
  
  b - l*u
  
