ctmc_classical.h File Reference

Classical stationary iterative solvers for pi Q = 0, sum(pi) = 1: the power method on the uniformised DTMC, Jacobi, Gauss-Seidel and SOR (Stewart 1994, ch. More...

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <string>
#include <vector>
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"

Include dependency graph for ctmc_classical.h:

Go to the source code of this file.

Classes
struct	line::mc::ClassicalResult< T >

Namespaces
namespace	line
	Conservation laws of a layered queueing network, enumerated from its structure.
namespace	line::mc

Functions
template<class T>
ClassicalResult< T >	line::mc::ctmc_power (const Matrix< T > &Q, double tol=CLASSICAL_DEFAULT_TOL, long maxiter=CLASSICAL_DEFAULT_MAXITER, const std::vector< T > &p0={})
	Power method on the uniformised DTMC P = I + Q/q, q = 1.1 max|q_ii|.
template<class T>
ClassicalResult< T >	line::mc::ctmc_jacobi (const Matrix< T > &Q, double tol=CLASSICAL_DEFAULT_TOL, long maxiter=CLASSICAL_DEFAULT_MAXITER, const std::vector< T > &p0={})
	Jacobi iteration on Q' x = 0 with renormalisation (see file comment).
template<class T>
ClassicalResult< T >	line::mc::ctmc_gaussseidel (const Matrix< T > &Q, double tol=CLASSICAL_DEFAULT_TOL, long maxiter=CLASSICAL_DEFAULT_MAXITER, const std::vector< T > &p0={})
	Forward Gauss-Seidel iteration on Q' x = 0 with renormalisation.
template<class T>
ClassicalResult< T >	line::mc::ctmc_sor (const Matrix< T > &Q, double tol=CLASSICAL_DEFAULT_TOL, long maxiter=CLASSICAL_DEFAULT_MAXITER, double omega=0.0, const std::vector< T > &p0={})
	Successive over-relaxation on Q' x = 0 with renormalisation.
bool	line::mc::is_classical_linsolver (const std::string &name)
	Whether a linsolver name selects one of the classical solvers, as accepted by ctmc_solve in all four codebases: "power", "jacobi", "gaussseidel", "sor".
template<class T>
ClassicalResult< T >	line::mc::ctmc_classical_by_name (const std::string &name, const Matrix< T > &Q)
	Run the classical solver named by is_classical_linsolver with defaults.

Variables
constexpr double	line::mc::CLASSICAL_DEFAULT_TOL = 1e-12
constexpr long	line::mc::CLASSICAL_DEFAULT_MAXITER = 10000

Detailed Description

Classical stationary iterative solvers for pi Q = 0, sum(pi) = 1: the power method on the uniformised DTMC, Jacobi, Gauss-Seidel and SOR (Stewart 1994, ch.

3; Haverkort 1998, section 15.1.3).

Templated port of matlab/src/api/mc/ctmc_power.m, ctmc_jacobi.m, ctmc_gaussseidel.m and ctmc_sor.m, which share the loop of matlab/src/api/mc/private/stationary_sweep_iterate.m. The JAR twin is jline.api.mc.StationarySweep and the Python twin line_solver.api.mc.classical. All four codebases run the same loop:

x = p0 / sum(p0) (uniform by default); for k = 0..maxiter: relres = ||x Q||_1 / max|q_ii|; stop with flag 0 if relres <= tol, flag 1 if k == maxiter; otherwise sweep, renormalise to sum one, and stop with flag 4 if the sweep sum is zero or not finite.

Sweeps, with A = Q' split as D + L + U: power x <- x + (x Q)'/q, q = 1.1 max|q_ii| (aperiodic by the self-loop of P = I + Q/q, so periodic chains converge) jacobi x <- -D^{-1} (L + U) x gaussseidel x <- (D + L)^{-1} (-U x) sor x <- (D + w L)^{-1} ((1 - w) D - w U) x

Flag 2 reports a zero diagonal (an absorbing state) for every method but the power method, which needs no diagonal inverse. Jacobi does not converge when the embedded jump chain is periodic, for instance on any birth-death chain; it then reports flag 1. On a reducible generator the limit depends on p0.

SOR with omega <= 0 selects the adaptive rule: start at omega = 1 and every 10 sweeps estimate the Gauss-Seidel spectral radius from the ratio of successive update norms, lambda, then set mu^2 = (lambda + omega - 1)^2 / (lambda omega^2), omega = min(2 / (1 + sqrt(1 - mu^2)), 1.95) when that increases omega (Young's formula for consistently ordered matrices, Hageman and Young 1981). If the residual later grows above 100 times the reference residual recorded at the last increase, omega reverts to 1 and adaptation stops.

Like GMRES these stop on a residual tolerance, so the exact instantiation is refused at compile time: there is nothing exact for the iteration to reach.

Definition in file ctmc_classical.h.