function x=gauss(A,b)

% x = gauss0(A,b)
% 
% Risolve il sistema lineare
%
%  A*x = b
%
% applicando il metodo di eliminazione di Gauss
%
% Versione ottimizzata che sfrutta il prodotto
% righe per colonne tra matrici che Matlab rende
% disponibile.
                  
[n,m]=size(A);

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

b=b(:);

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

for j=1:n-1

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

  A(j+1:n,j+1:n) = A(j+1:n,j+1:n) - A(j+1:n,j)*A(j,j+1:n)/A(j,j);

  b(j+1:n) = b(j+1:n) - A(j+1:n,j)*b(j)/A(j,j);

end

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

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

x(n) = b(n)/A(n,n); x = x(:);

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

return
