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 347861d5aee459bb10b74c6457ad215c8e32b567..ffee379a20ffddcdfb36b1db09a6ebd8ace83742 100644 (file)
--- a/iterative/jacobi.m
+++ b/iterative/jacobi.m
@@ -50,11 +50,6 @@ function [x, iterations, residual_norms] = ...
   %                         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;
@@ -64,38 +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 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 (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; rk_norm ];
-    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