 |
LINE Solver (C++)
Templated C++ port of the LINE queueing solver
|
|
LINE Solver (C++)
Templated C++ port of the LINE queueing solver
|
The M/G/1-type and GI/M/1-type fundamental-matrix solvers of MAMSolver / SMCSolver, ported from matlab/lib/thirdparty/MG1files: stat.m, MG1_EG.m, MG1_Decay.m, GIM1_Caudal.m, MG1_Shifts.m, MG1_CR.m, MG1_FI.m, MG1_NI.m (with solveSylvPowersDirectSum.m and solveSylvPowersRealSchur_FW.m), MG1_RR.m (with MG1_RR_Btemp.m and MG1_RR_tempB.m), MG1_IS.m and GIM1_R.m. More...
#include <algorithm>#include <cmath>#include <complex>#include <cstddef>#include <limits>#include <string>#include <vector>#include "line/num/complex_number.h"#include "line/num/number.h"#include "line/util/eig.h"#include "line/util/error.h"#include "line/util/fft.h"#include "line/util/linalg.h"#include "line/util/lstsq.h"#include "line/util/lu.h"#include "line/util/matrix.h"#include "line/util/svd.h"Go to the source code of this file.
Classes | |
| struct | line::smc::Drift |
| The drift of an M/G/1-type sequence, and the invariant vector it uses. More... | |
| struct | line::smc::ShiftResult |
| What MG1_Shifts returns: the shifted sequence and the drift it measured. More... | |
| struct | line::smc::Mg1CrOptions |
| Options of MG1_CR, with the reference's defaults. More... | |
| struct | line::smc::Mg1FiOptions |
| Options of MG1_FI, with the reference's defaults. More... | |
| struct | line::smc::Mg1NiOptions |
| Options of MG1_NI, with the reference's defaults. More... | |
| struct | line::smc::Mg1RrOptions |
| Options of MG1_RR, with the reference's defaults. More... | |
| struct | line::smc::rr_detail::Disp |
| The displacement representation B = L(b) + L(c1) L(Z r1)' + L(c2) L(Z r2)' of the Ramaswami-reduction matrix, with L(x) the block lower-triangular Toeplitz matrix whose first block column is x (N*m x m) and Z the block down-shift. More... | |
| struct | line::smc::Mg1IsOptions |
| Options of MG1_IS, with the reference's defaults. More... | |
Namespaces | |
| namespace | line |
| Conservation laws of a layered queueing network, enumerated from its structure. | |
| namespace | line::smc |
| namespace | line::smc::mg1x_detail |
| namespace | line::smc::rr_detail |
Typedefs | |
| using | line::smc::Blocks = std::vector<Matrix<double>> |
Functions | |
| template<class T> | |
| Matrix< T > | line::smc::madd (const Matrix< T > &A, const Matrix< T > &B) |
| A + B. | |
| template<class T> | |
| Matrix< T > | line::smc::msub (const Matrix< T > &A, const Matrix< T > &B) |
| A - B. | |
| Matrix< double > | line::smc::mscale (const Matrix< double > &A, double c) |
| c * A. | |
| std::vector< double > | line::smc::rowsums (const Matrix< double > &A) |
| sum(A,2), the row sums, as a column held in a vector. | |
| double | line::smc::inf_norm (const Matrix< double > &A) |
| norm(A,inf), the largest absolute row sum. | |
| double | line::smc::inf_norm (const Blocks &blk, std::size_t from) |
| norm(A,inf) over a whole block sequence stacked vertically. | |
| double | line::smc::max_abs_diff (const Matrix< double > &A, const Matrix< double > &B) |
| max(max(abs(A-B))). | |
| double | line::smc::max_col_sum (const Matrix< double > &A) |
| max(sum(A)), the largest column sum WITHOUT absolute values, as in MATLAB. | |
| Blocks | line::smc::blocks_of (const Matrix< double > &A, std::size_t m) |
| Splits the wide [A0 A1 ... Amax] into its m x m blocks. | |
| Matrix< double > | line::smc::hcat (const Blocks &blk) |
| Re-assembles a block sequence into the wide [A0 A1 ... Amax]. | |
| Matrix< double > | line::smc::vcat (const Blocks &blk) |
| Stacks a block sequence vertically, [A0; A1; ...; Amax]. | |
| Blocks | line::smc::vblocks_of (const Matrix< double > &A, std::size_t r) |
| Splits a vertical stack into its blocks of r rows. | |
| std::vector< double > | line::smc::stat (const Matrix< double > &A) |
| Stationary distribution of a stochastic matrix: the left eigenvector for eigenvalue 1, nonnegative and summing to one. | |
| std::vector< double > | line::smc::rowvec_times (const std::vector< double > &v, const Matrix< double > &A) |
| theta A, the row vector times matrix product used throughout. | |
| double | line::smc::dot (const std::vector< double > &a, const std::vector< double > &b) |
| The inner product of a row vector with a column held as a vector. | |
| Drift | line::smc::mg1_drift (const Blocks &A) |
| drift = theta * beta with beta = (Amax)e + (Amax+Amax-1)e + ..., the expected level increment per transition of the phase process. | |
| std::complex< double > | line::smc::max_eig (const Matrix< double > &M) |
| max(eig(M)) with MATLAB's semantics on a complex spectrum: the element of largest modulus, ties broken by the larger phase angle. | |
| Matrix< double > | line::smc::poly_at (const Blocks &A, double z) |
| A(z) = A0 + A1 z + ... + Amax z^max, by Horner as the reference writes it. | |
| std::vector< double > | line::smc::pf_vector (const Matrix< double > &M, bool left) |
| The Perron-Frobenius eigenvector of M, right (left == false) or left, scaled to unit sum: the null vector of M - PF(M) I (or its transpose), read off the SVD. | |
| double | line::smc::mg1_decay (const Blocks &A, std::vector< double > *uT=nullptr) |
| Decay rate of a recurrent M/G/1-type chain: the unique z > 1 with PF(A(z)) = z. | |
| double | line::smc::gim1_caudal (const Blocks &A, std::vector< double > *v=nullptr) |
| Caudal characteristic of a GI/M/1-type chain: the spectral radius of R, the unique z in (0,1) with PF(A(z)) = z. | |
| ShiftResult | line::smc::mg1_shifts (const Blocks &Ain, const std::string &shift_type) |
| Shift technique for the M/G/1-type sequence. | |
| void | line::smc::mg1_unshift (Matrix< double > &G, const ShiftResult &sh, const std::string &shift_type) |
| Put back on G the rank-one term a shift removed: ones/m for a 'one' shift at drift < 1, tau v e^T for a 'tau' shift at drift > 1, both for 'dbl'. | |
| Matrix< double > | line::smc::mg1_eg (const Blocks &Ain, bool &found) |
| G in closed form when A0 has rank one. | |
| Matrix< double > | line::smc::mg1_cr (const Blocks &Ain, const Mg1CrOptions &opts=Mg1CrOptions()) |
| Cyclic reduction for M/G/1-type Markov chains [Bini, Meini]. | |
| Matrix< double > | line::smc::mg1_fi (const Blocks &Ain, const Mg1FiOptions &opts=Mg1FiOptions()) |
| Functional iterations for M/G/1-type Markov chains [Neuts]. | |
| Matrix< double > | line::smc::mg1x_detail::sub (const Matrix< double > &A, std::size_t r0, std::size_t nr, std::size_t c0, std::size_t nc) |
| A(r0:r0+nr, c0:c0+nc), half-open, 0-based. | |
| void | line::smc::mg1x_detail::put (Matrix< double > &A, std::size_t r0, std::size_t c0, const Matrix< double > &S) |
| A(r0:, c0:) = S. | |
| Matrix< double > | line::smc::mg1x_detail::trans (const Matrix< double > &A) |
| Matrix< double > | line::smc::mg1x_detail::hjoin (const Matrix< double > &A, const Matrix< double > &B) |
| [A B], horizontally. | |
| Matrix< double > | line::smc::mg1x_detail::vjoin (const Matrix< double > &A, const Matrix< double > &B) |
| [A; B], vertically. | |
| std::vector< double > | line::smc::mg1x_detail::lsolve (const Matrix< double > &Z, const std::vector< double > &b) |
| Z \ b for a square Z. | |
| void | line::smc::mg1x_detail::qr_econ (const Matrix< double > &A, Matrix< double > &Q, Matrix< double > &R) |
| Economy QR, A = Q R with Q (r x k) orthonormal columns and R (k x c), k = min(r, c), by Householder reflections. | |
| Matrix< double > | line::smc::sylv_powers_direct (const Matrix< double > &A, const Blocks &B, const Matrix< double > &C) |
| Solve sum_{j=1}^N B_j Y A^{j-1} = C directly, through the Kronecker form of vec(Y). | |
| Matrix< double > | line::smc::sylv_powers_real_schur (const Matrix< double > &A, const Blocks &B, const Matrix< double > &C) |
| Solve sum_{j=1}^N B_j Y A^{j-1} = C by a real Schur form A = U T U': with Y = X U' the system becomes sum_j B_j X T^{j-1} = C U, and T quasi-triangular makes it a forward substitution over the columns of X, one n x n solve per 1 x 1 block and one 2n x 2n solve per 2 x 2 block. | |
| Matrix< double > | line::smc::mg1_ni (const Blocks &Ain, const Mg1NiOptions &opts=Mg1NiOptions()) |
| Newton iteration for M/G/1-type Markov chains. | |
| double | line::smc::rr_detail::sum_all (const Matrix< double > &A) |
| Matrix< double > | line::smc::mg1_rr (const Blocks &D, const Mg1RrOptions &opts=Mg1RrOptions()) |
| Ramaswami reduction for M/G/1-type Markov chains [Bini, Meini, Ramaswami]. | |
| Matrix< double > | line::smc::mg1_is (const Blocks &D, const Mg1IsOptions &opts=Mg1IsOptions()) |
| Invariant subspace method for M/G/1-type Markov chains [Akar, Sohraby]. | |
| Matrix< double > | line::smc::gim1_r (const Blocks &Ain, const std::string &dual, const std::string &algor) |
| R of a GI/M/1-type Markov chain, through the G of its dual. | |
The M/G/1-type and GI/M/1-type fundamental-matrix solvers of MAMSolver / SMCSolver, ported from matlab/lib/thirdparty/MG1files: stat.m, MG1_EG.m, MG1_Decay.m, GIM1_Caudal.m, MG1_Shifts.m, MG1_CR.m, MG1_FI.m, MG1_NI.m (with solveSylvPowersDirectSum.m and solveSylvPowersRealSchur_FW.m), MG1_RR.m (with MG1_RR_Btemp.m and MG1_RR_tempB.m), MG1_IS.m and GIM1_R.m.
These are the third-party numerics the ETAQA aggregation sits on: MG1_CR returns the minimal nonnegative G of an M/G/1-type chain and GIM1_R the minimal nonnegative R of a GI/M/1-type one, and without them solver_mam_bmap_map_1 and solver_mam_map_bmap_1 have nothing to aggregate. The port follows the MATLAB line by line, including the block index arithmetic, the stopping tests and the constants, so that a divergence against the reference is a bug here and not a design difference.
BLOCK LAYOUT. The reference passes a block sequence as one wide matrix A = [A0 A1 A2 ... Amax], m rows by m*(max+1) columns. This port carries the same sequence as a std::vector<Matrix<double>> of m x m blocks, which is the identical object with the index arithmetic done once, in blocks_of / hcat, instead of at every use. The GI/M/1 side stacks its blocks VERTICALLY in the reference; GIM1_R_ETAQA is what transposes that stack into the horizontal one, so everything below is horizontal.
DOUBLE ONLY, and not by preference. MG1_Decay and GIM1_Caudal bisect on the Perron-Frobenius eigenvalue of A(z), which needs LAPACK; MG1_CR evaluates its polynomials at complex roots of unity through an FFT, whose twiddle factors are cos/sin; MG1_pi_ETAQA drops a column chosen by a numerical rank test, which needs an SVD. None of the three has a multiprecision or exact counterpart in this tree, so the whole family is declared on Matrix<double> and the solvers that call it refuse at any other arithmetic rather than down-converting behind the caller's back.
MG1_Shifts IS PORTED IN FULL, all three ShiftTypes at either drift. Its reference writes the last block of a row shift as rowhatA(1,maxd*i:end) = uT*A(:,maxd*i:end) with i the value 1 left over from the beta loop, so the line addresses column maxd, not block maxd. It is nevertheless CORRECT: both sides take the same columns, so columns maxd*m+1:end receive exactly uT*A_maxd, and the stray columns before them are overwritten by the loop that follows. An earlier version of this header read that line as defective and refused the branches; the transient-chain Newton iteration under GIM1_R's Ramaswami dual is what reaches them. With maxd = 1 the beta loop never runs, i is the imaginary unit and the reference errors; this port computes the last block the line intends.
Definition in file mg1.h.
1.18.0