function [x,iter,err] = gs(A,x,b,iter_max,tol)
%
% Cerca di risolve un sistema lineare Ax = b
% usando il metodo di Gauss-Seidel.
%
% USO: 
%             [x,iter,err] = gs(A,x,b,iter_max,tol);
%
% In ingresso:
% A           matrice dei coefficienti
% x           x_0
% b           vettore dei termini noti
% iter_max    numero massimo di iterazioni
% tol         tolleranza
%
% In uscita:
% x           soluzione 
% iter        numero di iterazioni
% err         errore calcolato ad ogni iterazione 

iter = 0;

bn2 = norm(b);
if bn2==0.0
  bn2 = 1.0;
end

r = b-A*x;
err(1) = norm(r)/bn2;
if err(1)<=tol
  return;
end

E = tril(A,-1);    % estrae la parte triangolare inferiore da A
D = diag(diag(A)); % costruisce la matrice diagonale prendendo
                   % la diagonale principale di A
F = triu(A,1);     % estrae la parte triangolare superiore da A

M = E+D;
N = -F;

for iter = 1:iter_max
   x_old = x;
   x = M\(N*x+b); % matlab riconosce automaticamente che la matrice M
                  % e` traingolare inferiore ed utilizza il metodo
                  % piu` opportuno per "risolvere" il sistema lineare
   err(iter) = norm(x-x_old)/norm(x);
   if err(iter)<=tol
     break;
   end
end

if iter>=iter_max
  disp('ERRORE: numero eccessivo di iterazioni in Gauss-Seidel\n')
end