X-Git-Url: http://gitweb.michael.orlitzky.com/?p=octave.git;a=blobdiff_plain;f=iterative%2Fjacobi.m;h=ffee379a20ffddcdfb36b1db09a6ebd8ace83742;hp=347861d5aee459bb10b74c6457ad215c8e32b567;hb=40a26bb2dfcf53296d9996b88008a605ad59573c;hpb=2357d986b4d33f5d18681a97c246a141e94667c5

diff --git a/iterative/jacobi.m b/iterative/jacobi.m
index 347861d..ffee379 100644
--- 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