 |
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 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"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... | |
Namespaces | |
| namespace | line |
| namespace | line::smc |
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. | |
| double | line::smc::mg1_decay (const Blocks &A) |
| 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) |
| 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. | |
| 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::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 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.
TWO REFERENCE BRANCHES ARE REFUSED BY NAME rather than transcribed, because they cannot run in MATLAB either. MG1_Shifts writes rowhatA(1,maxd*i:end) = uT*A(:,maxd*i:end) in three places (the drift < 1 'tau' branch and the drift > 1 'one' branch), where i is not a loop variable at that point: it is either MATLAB's imaginary unit, which makes the index complex and errors, or the leftover value 1 from the beta loop above, which addresses column maxd of a matrix whose blocks start every m columns. Either way the line does not compute "the last block", which is what the surrounding code needs. ETAQA reaches neither branch: it shifts a positive recurrent chain (drift < 1) with the default ShiftType 'one', which is the branch that is correct and is ported.
Definition in file mg1.h.
1.18.0