+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