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 infinity norm
  %                     of the residual divided by the infinity 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
  %                         (infinity)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 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
