ctmc_gmres.h File Reference

Restarted GMRES with an ILUT preconditioner, for the linear systems a generator produces. More...

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

Include dependency graph for ctmc_gmres.h:

Go to the source code of this file.

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

Namespaces
namespace	line
namespace	line::mc

Functions
template<class T>
GmresResult< T >	line::mc::ctmc_gmres (const Matrix< T > &A, const std::vector< T > &b, double tol=1e-12, long restart=0, long maxit=0, const std::vector< T > &x0=std::vector< T >())
	Restarted GMRES with an ILUT preconditioner, for the linear systems a generator produces.

Variables
constexpr std::size_t	line::mc::GMRES_MIN_STATES = 6000
	Order above which the direct sparse factorization is abandoned in favour of the Krylov path.

Detailed Description

Restarted GMRES with an ILUT preconditioner, for the linear systems a generator produces.

Templated port of matlab/src/api/mc/ctmc_gmres.m and jar/src/main/java/jline/api/mc/Ctmc_gmres.java. This is the iterative counterpart of the direct solve in ctmc_solve, meant for generators whose LU fill-in exceeds the memory available; no CTMC-specific processing happens here, so the same routine serves the stochastic complement and the aggregation methods.

Two preparation steps are not optional on a generator, and both are ported. Rows are equilibrated to unit max norm, so the O(1) normalization row does not mix with rows carrying rates of a wholly different magnitude. The states are then reordered by reverse Cuthill-McKee: in the natural ordering of a birth-death chain the unpivoted incomplete elimination has growth factor (mu/lambda)^n, which overflows within a few thousand states, and a fill-reducing ordering rather than pivoting is what removes it.

The preconditioner is ILUT(p, tau) of Saad (1994): the row is expanded into a dense workspace, entries below tau times the mean magnitude of the original row are dropped as they are produced, and the p largest survivors are kept in each of the L and U parts. A vanishing pivot is replaced by a value of the order of the row threshold, keeping its sign, which is Saad's remedy and avoids abandoning the factorization for one rate-free state. When the factorization breaks down entirely the preconditioner degrades to Jacobi.

MATLAB-VS-JAVA DISAGREEMENT, resolved in favour of Java. MATLAB calls the built-in gmres with (L,U), which applies the preconditioner on the LEFT, so its RELRES is the preconditioned residual norm(M\(b-Ax))/norm(M\b) and depends on the preconditioner. The JAR expands the Krylov space of A M^-1 instead, and reports the true residual norm(b-Ax)/norm(b). The solution both converge to is the same; the reported relres is not, and a tolerance on the true residual is the one a caller can act on, so this port follows the JAR. FLAG keeps the MATLAB convention: 0 converged, 1 iteration limit, 3 stagnation or divergence or a non-finite iterate.

GATED ON TRANSCENDENTAL ARITHMETIC. GMRES stops on a residual tolerance and its Arnoldi step normalizes by a Euclidean norm, so square roots appear in every iteration and there is no exact answer to converge to: at Rational the iteration would run to the iteration limit with exploding denominators. The exact solve of the same system is ctmc_solve, which is what a Rational caller should use.

Definition in file ctmc_gmres.h.