Add classical_iteration() sans comments, refactor jacobi() to use it.

[octave.git] / iterative / jacobi.m

diff --git a/iterative/jacobi.m b/iterative/jacobi.m
index 3e4f946aa3f15cb2e728ec4befc752f8c7357382..ffee379a20ffddcdfb36b1db09a6ebd8ace83742 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.
-  %                    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,15 +45,11 @@ 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)
-     save_residual_norms = true;
-     residual_norms = [];
-  end
 
   if (nargin < 4)
     tolerance = 1e-10;
@@ -63,29 +59,12 @@ function [x, iterations, residual_norms] = ...
     max_iterations = intmax();
   end
 
-  % This creates a diagonal matrix from the diagonal entries of A.
-  M = diag(diag(A));
-  N = M - A;
-
-  % Avoid recomputing this each iteration.
-  b_norm = norm(b);
-
-  k = 0;
-  xk = x0;
-  rk = b - A*xk;
-
-  while (norm(rk) > tolerance*b_norm && k < max_iterations)
-    x_next = M \ (N*xk + b);
-    r_next = b - A*x_next;
-
-    k = k + 1;
-    xk = x_next;
-    rk = r_next;
-    if (save_residual_norms)
-      residual_norms = [ residual_norms; norm(rk) ];
-    end
+  M_of_A = @(A) diag(diag(A));
+  if (nargout > 2)
+    [x, iterations, residual_norms] = ...
+      classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);
+  else
+    [x, iterations] = ...
+      classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations);
   end
-
-   iterations = k;
-   x = xk;
 end