+function [x, iterations, residual_norms] = ...
+ successive_over_relaxation(A, b, omega, 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:
+ %
+ % ``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.
+ %
+ % ``omega`` -- The relaxation factor.
+ %
+ % ``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