From 2357d986b4d33f5d18681a97c246a141e94667c5 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Mon, 6 May 2013 12:54:07 -0400
Subject: [PATCH] Use infinity norms in jacobi().

---
 iterative/jacobi.m | 30 ++++++++++++++++++++----------
 1 file changed, 20 insertions(+), 10 deletions(-)

diff --git a/iterative/jacobi.m b/iterative/jacobi.m
index 3e4f946..347861d 100644
--- a/iterative/jacobi.m
+++ b/iterative/jacobi.m
@@ -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
 
-- 
2.44.2