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

[octave.git] / iterative / classical_iteration.m

diff --git a/iterative/classical_iteration.m b/iterative/classical_iteration.m
new file mode 100644 (file)
index 0000000..4e43ff1
--- /dev/null
+++ b/iterative/classical_iteration.m
@@ -0,0 +1,44 @@
+function [x, iterations, residual_norms] = ...
+        classical_iteration(A, b, x0, M_of_A, tolerance, max_iterations)
+
+  save_residual_norms = false;
+  if (nargout > 2)
+     save_residual_norms = true;
+     residual_norms = [];
+  end
+
+  % This creates a diagonal matrix from the diagonal entries of A.
+  M = M_of_A(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
+  end
+
+   iterations = k;
+   x = xk;
+end