Add jacobi.m.

[octave.git] / iterative / jacobi.m

diff --git a/iterative/jacobi.m b/iterative/jacobi.m
new file mode 100644 (file)
index 0000000..3e4f946
--- /dev/null
+++ b/iterative/jacobi.m
@@ -0,0 +1,91 @@
+function [x, iterations, residual_norms] = ...
+        jacobi(A, b, x0, tolerance, max_iterations)
+  %
+  % Solve the system,
+  %
+  %   Ax = b
+  %
+  % iteratively using Jacobi iterations. That is, we let,
+  %
+  %   A = M - N = D - (L + U)
+  %
+  % where D is a diagonal matrix consisting of the diagonal entries of
+  % A (the rest zeros), and N = (L + U) are the remaining upper- and
+  % lower-triangular parts of A. Now,
+  %
+  %   Ax = (M - N)x = Mx - Nx = b
+  %
+  % has solution,
+  %
+  %   x = M^(-1)*(Nx + b)
+  %
+  % Thus, our iterations are of the form,
+  %
+  %   x_{k+1} = M^(-1)*(Nx_{k} + b)
+  %
+  % INPUT:
+  %
+  %   ``A`` -- The n-by-n coefficient matrix of the system.
+  %
+  %   ``b`` -- An n-by-1 vector; the right-hand side of the system.
+  %
+  %   ``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``.
+  %
+  %   ``max_iterations`` -- (optional; default: intmax) the maximum
+  %                         number of iterations we will perform.
+  %
+  % OUTPUT:
+  %
+  %   ``x`` -- An n-by-1 vector; the approximate solution to the system.
+  %
+  %   ``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.
+  %
+  save_residual_norms = false;
+  if (nargout > 2)
+     save_residual_norms = true;
+     residual_norms = [];
+  end
+
+  if (nargin < 4)
+    tolerance = 1e-10;
+  end
+
+  if (nargin < 5)
+    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
+  end
+
+   iterations = k;
+   x = xk;
+end