function [x, iterations, residual_norms] = ...
	 successive_over_relaxation(omega, A, b, x0, ...
				    tolerance, max_iterations)
  %
  % Solve the system,
  %
  %   Ax = b
  %
  % iteratively using SOR. That is, we let,
  %
  %   A = M - N = ((1/omega)*D + L) - U
  %
  % where D is a diagonal matrix consisting of the diagonal entries of
  % A (the rest zeros), and L,U are the upper- and lower-triangular
  % parts of A respectively. 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:
  %
  %   ``omega`` -- The relaxation factor.
  %
  %   ``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.
  %

  if (nargin < 5)
    tolerance = 1e-10;
  end

  if (nargin < 6)
    max_iterations = intmax();
  end

  M_of_A = @(A) (1/omega)*diag(diag(A)) + tril(A,-1);

  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
end
