Use infinity norms in jacobi().

diff --git a/iterative/jacobi.m b/iterative/jacobi.m
index 3e4f946aa3f15cb2e728ec4befc752f8c7357382..347861d5aee459bb10b74c6457ad215c8e32b567 100644 (file)
--- 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.
   %   ``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.
   %
   %   ``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.
   %
   %
   %   ``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)
   %
   save_residual_norms = false;
   if (nargout > 2)


@@ -67,22 +68,31 @@ function [x, iterations, residual_norms] = ...
   M = diag(diag(A));
   N = M - A;
 
   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;
 
   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;
     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)
     if (save_residual_norms)

-      residual_norms = [ residual_norms; norm(rk) ];
+      residual_norms = [ residual_norms; rk_norm ];

     end
   end
 
     end
   end