function [x, its, norm_r] = gradott (A, b, tol, max_its, x)
%
% [x, its, norm_r] = gradott (A, b, tol, max_its, x0)
%
% Compute the solution of A*x=b through the gradient 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. A has to be symmetric and positive definite.
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
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;
  Ar = A * r;
  alpha = r' * r / (r' * Ar);
  if (alpha < 0)
     warning ('La matrice non sembra definita positiva')
  end
  x = x + alpha * r;
  r = r - alpha * Ar;
  norm_r = norm (r);
end
if (norm_r > tol(1) * norm_b + tol(2))
  warning ('Impossibile raggiungere la tolleranza richiesta')
end
%!test
%! A = sparse ([2, 1, 0;1, 2, 1;0, 1, 2]);
%! b = A * (1:3)';
%! x = gradott (A, b, 1e-8);
%! assert (x, (1:3)', 1e-7)
%!test
%! A = diag ([1, 100]);
%! b = A * (1:2)';
%! x0 = [100; 100];
%! [x, its] = gradott (A, b, [], [], x0);
%! assert (its, 7)
%!test % matrice non definita positiva, non se ne accorge, non convergenza
%! A = [1, -1;1, 1];
%! b = A * (1:2)';
%! x0 = [100; 100];
%! [x, its] = gradott (A, b, [], [], x0);
%! assert (its, 100)
%!test % matrice non definita positiva, se ne accorge, convergenza
%! A = [1, 0;1, -1];
%! b = A * (1:2)';
%! x0 = [100; 100];
%! [x, its] = gradott (A, b, [], [], x0);
%! assert (x, (1:2)', 1e-6)
