function [x, its, norm_r] = gaussseidel (A, b, tol, max_its, x)
%
% [x, its, norm_r] = gaussseidel (A, b, tol, max_its, x0)
%
% Compute the solution of A*x=b through the Gauss-Seidel method. tol is a
% vector of two components with the relative and the absolute
% tolerances (defaults [1e-6, 1e-6]). max_its is the maximum number of
% iterations (default 100) and x0 the initial solution (default zero).
% x is the computed solution, its the number of iterations and norm_r
% the norm of the final residual.
n = length (A);
if ((nargin <= 2) || isempty (tol))
  tol = [1e-6, 1e-6];
end
if ((nargin <= 3) || isempty (max_its))
  max_its = 100;
end
if ((nargin <= 4) || isempty (x))
  x = zeros (n, 1);
end
if (isscalar (tol))
  tol = [tol, tol];
end
n = length (A);
norm_b = norm (b);
r = b - A * x;
norm_r = norm (r);
its = 0;
while ((norm_r > tol(1) * norm_b + tol(2)) && (its < max_its))
  its = its + 1;
  for i = 1:n
    x(i) = (b(i) - ...
         A(i, 1:i - 1) * x(1:i - 1) - ...
         A(i, i + 1:n) * x(i + 1:n)) / A(i, i);
  end
  r = b - A * x;
  norm_r = norm (r);
end
if (norm_r > tol(1) * norm_b + tol(2))
  warning ('Impossibile raggiungere la tolleranza richiesta')
end
%!test
%! A = [2, 1, 0; 1, 2, 1; 0, 1, 2];
%! b = A * ones (3, 1);
%! x = gaussseidel (A, b, 1e-8);
%! assert (x, ones (3, 1), 1e-7)
%!test
%! A = [1, 1, 1; 1, 2, 1; 0, 1, 1];
%! b = A * ones (3, 1);
%! x = gaussseidel (A, b, 1e-8);
%! assert (x, ones (3, 1), 1e-7)
