function [x,L,R] = fattLR(A,b)


% Risolve il sistema lineare
%
%        A*x = b
%
% applicando il metodo di eliminazione di Gauss
% e calcola la fattorizzazione 
%
%        A = L*R
%
% Nota: e' possibile ottimizzare questo codice
%       sfruttando il prodotto righe per colonne
%       tra matrici che Matlab rende disponibile.


[n,m]=size(A);

if(n~=m), 
   disp('La matrice deve essere quadrata'),
end


% Le seguenti istruzioni trasformano il sistema di
% equazioni in uno ad esso equivalente ma con matrice
% dei coefficienti triangolare superiore

b=b(:);

L=eye(n);

for j = 1:n-1

  if (A(j,j)==0), 
     error('Metodo di eliminazione di Gauss non applicabile'), 
  end

  for i = j+1:n 

     L(i,j) = A(i,j)/A(j,j);

     for k = j+1:n
        A(i,k) = A(i,k) - L(i,j)*A(j,k);
     end

     b(i) = b(i) - L(i,j)*b(j);
  end

end

if (A(n,n)==0), 
     error('Metodo di eliminazione di Gauss non applicabile'), 
end

% La seguente istruzione assegna ad R la parte triangolare 
% superiore della A ottenuta dopo la eliminazione di Gauss.
% Gli elementi della parte strettamente triangolare inferiore 
% di R vengono posti uguale a 0.

R = triu(A);

% Risoluzione del sistema triangolare superiore ottenuto
% mediante sostituzioni successive all'indietro

x(n) = b(n)/R(n,n); 

x= x(:); % Rende x un vettore colonna 

for i=n-1:-1:1

   x(i) = b(i);

   for j=i+1:n
     x(i) = x(i) - R(i,j)*x(j);
   end

   x(i) = x(i)/R(i,i);

end    

return
                  
