author Michael Orlitzky Mon, 6 May 2013 16:54:07 +0000 (12:54 -0400) committer Michael Orlitzky Mon, 6 May 2013 16:54:07 +0000 (12:54 -0400)

@@ -32,9 +32,9 @@ function [x, iterations, residual_norms] = ...
%   ``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.
@@ -45,9 +45,10 @@ function [x, iterations, residual_norms] = ...
%
%   ``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)
@@ -67,22 +68,31 @@ function [x, iterations, residual_norms] = ...
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