% ``x0`` -- An n-by-1 vector; an initial guess to the solution.
%
% ``tolerance`` -- (optional; default: 1e-10) the stopping tolerance.
- % we stop when the relative error (the norm of the
- % residual divided by the norm of ``b``) is less
- % than ``tolerance``.
+ % we stop when the relative error (the infinity norm
+ % of the residual divided by the infinity norm of
+ % ``b``) is less than ``tolerance``.
%
% ``max_iterations`` -- (optional; default: intmax) the maximum
% number of iterations we will perform.
%
% ``iterations`` -- The number of iterations taken.
%
- % ``residual_norms`` -- An n-by-iterations vector of the residual norms
- % at each iteration. If not requested, they will
- % not be computed to save space.
+ % ``residual_norms`` -- An n-by-iterations vector of the residual
+ % (infinity)norms at each iteration. If not
+ % requested, they will not be computed to
+ % save space.
%
save_residual_norms = false;
if (nargout > 2)
M = diag(diag(A));
N = M - A;
- % Avoid recomputing this each iteration.
- b_norm = norm(b);
+ % Avoid recomputing this at the beginning of each iteration.
+ b_norm = norm(b, 'inf');
+ relative_tolerance = tolerance*b_norm;
k = 0;
xk = x0;
rk = b - A*xk;
+ rk_norm = norm(rk, 'inf');
- while (norm(rk) > tolerance*b_norm && k < max_iterations)
+ while (rk_norm > relative_tolerance && k < max_iterations)
+ % This is almost certainly slower than updating the individual
+ % components of xk, but it's "fast enough."
x_next = M \ (N*xk + b);
r_next = b - A*x_next;
k = k + 1;
xk = x_next;
rk = r_next;
+
+ % We store the current residual norm so that, if we're saving
+ % them, we don't need to recompute it at the beginning of the next
+ % iteration.
+ rk_norm = norm(rk, 'inf');
if (save_residual_norms)
- residual_norms = [ residual_norms; norm(rk) ];
+ residual_norms = [ residual_norms; rk_norm ];
end
end