function [LU,d]=FacLU(A)
  n=size(A,"r");
 
// Transformation de la matrice    
for(k=1:1:n-1)
  for (i=k+1:1:n)
    A(i,k)=A(i,k)/A(k,k);
    for (j=k+1:1:n)
     A(i,j)=A(i,j)-A(i,k)*A(k,j);     
  end 
 end
end
LU=A;
//Calcul du d�terminant 
d=1;
for (k=1:1:n)
 d=d*A(k,k);
end
endfunction

function y=DescenteL(L,b)
n=size(L,"r");  
//Resolution de Ly=b
for i=1:1:n
  S=0;
  for j=1:1:i-1
    S=S+L(i,j)*y(j);
    end  
y(i)=(b(i)-S);
end
endfunction

function x=RemonteeU(U,y)
  n=size(U,"r");
//Resolution de Ux=y
for i=n:-1:1
  S=0;
  for j=i+1:1:n
    S=S+U(i,j)*x(j);
    end  
x(i)=(y(i)-S)/U(i,i);
end
endfunction

function [x]=ResoutLU(A,b) 
n=size(A,"r"); 
LU=FacLU(A);
y=DescenteL(LU,b);
x=RemonteeU(LU,y);
endfunction

function X=ResoutLU2(A,B)
LU=FacLU(A);
p=size(B,"c"); 
for(j=1:1:p)
y(:,j)=DescenteL(LU,B(:,j));
X(:,j)=RemonteeU(LU,y(:,j));    
  end    
endfunction

function X=Inverse(A)
  n=size(A,"r");
  B=eye(n,n);   
  X=ResoutLU2(A,B);
  endfunction

       



   
