 |
LINE Solver (C++)
Templated C++ port of the LINE queueing solver
|
|
LINE Solver (C++)
Templated C++ port of the LINE queueing solver
|
Namespaces | |
| namespace | mg1x_detail |
| namespace | rr_detail |
Classes | |
| struct | Drift |
| The drift of an M/G/1-type sequence, and the invariant vector it uses. More... | |
| struct | ShiftResult |
| What MG1_Shifts returns: the shifted sequence and the drift it measured. More... | |
| struct | Mg1CrOptions |
| Options of MG1_CR, with the reference's defaults. More... | |
| struct | Mg1FiOptions |
| Options of MG1_FI, with the reference's defaults. More... | |
| struct | Mg1NiOptions |
| Options of MG1_NI, with the reference's defaults. More... | |
| struct | Mg1RrOptions |
| Options of MG1_RR, with the reference's defaults. More... | |
| struct | Mg1IsOptions |
| Options of MG1_IS, with the reference's defaults. More... | |
Typedefs | |
| using | Blocks = std::vector<Matrix<double>> |
Functions | |
| Matrix< double > | mg1_g_etaqa (const Matrix< double > &A) |
| G of an M/G/1-type chain, uniformized first. | |
| std::vector< double > | mg1_pi_etaqa (const Matrix< double > &Bin, const Matrix< double > &Ain, const Matrix< double > &G, const Matrix< double > &C0in=Matrix< double >()) |
| Aggregated stationary vector [pi0, pi1, pi2+pi3+...] of an M/G/1-type chain. | |
| double | mg1_qlen_etaqa (const Matrix< double > &Bin, const Matrix< double > &Ain, const std::vector< double > &pi, std::size_t n, const Matrix< double > &C0in=Matrix< double >()) |
| n-th moment of the level (the queue length) of an M/G/1-type chain from the ETAQA aggregates. | |
| Matrix< double > | gim1_r_etaqa (const Matrix< double > &A) |
| R of a GI/M/1-type chain, uniformized first. | |
| std::vector< double > | gim1_pi_etaqa (const Matrix< double > &Bin, const Matrix< double > &Ain, const Matrix< double > &R, const Matrix< double > &B0in=Matrix< double >()) |
| Aggregated stationary vector [pi0, pi1, pi2+pi3+...] of a GI/M/1-type chain. | |
| double | gim1_qlen_etaqa (const Matrix< double > &Bin, const Matrix< double > &Ain, const Matrix< double > &R, const std::vector< double > &pi, std::size_t n, const Matrix< double > &B0in=Matrix< double >()) |
| n-th moment of the level of a GI/M/1-type chain from the ETAQA aggregates. | |
| template<class T> | |
| Matrix< T > | madd (const Matrix< T > &A, const Matrix< T > &B) |
| A + B. | |
| template<class T> | |
| Matrix< T > | msub (const Matrix< T > &A, const Matrix< T > &B) |
| A - B. | |
| Matrix< double > | mscale (const Matrix< double > &A, double c) |
| c * A. | |
| std::vector< double > | rowsums (const Matrix< double > &A) |
| sum(A,2), the row sums, as a column held in a vector. | |
| double | inf_norm (const Matrix< double > &A) |
| norm(A,inf), the largest absolute row sum. | |
| double | inf_norm (const Blocks &blk, std::size_t from) |
| norm(A,inf) over a whole block sequence stacked vertically. | |
| double | max_abs_diff (const Matrix< double > &A, const Matrix< double > &B) |
| max(max(abs(A-B))). | |
| double | max_col_sum (const Matrix< double > &A) |
| max(sum(A)), the largest column sum WITHOUT absolute values, as in MATLAB. | |
| Blocks | blocks_of (const Matrix< double > &A, std::size_t m) |
| Splits the wide [A0 A1 ... Amax] into its m x m blocks. | |
| Matrix< double > | hcat (const Blocks &blk) |
| Re-assembles a block sequence into the wide [A0 A1 ... Amax]. | |
| Matrix< double > | vcat (const Blocks &blk) |
| Stacks a block sequence vertically, [A0; A1; ...; Amax]. | |
| Blocks | vblocks_of (const Matrix< double > &A, std::size_t r) |
| Splits a vertical stack into its blocks of r rows. | |
| std::vector< double > | 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 > | rowvec_times (const std::vector< double > &v, const Matrix< double > &A) |
| theta A, the row vector times matrix product used throughout. | |
| double | 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 | 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 > | 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 > | 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 > | 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 | 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 | 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 | mg1_shifts (const Blocks &Ain, const std::string &shift_type) |
| Shift technique for the M/G/1-type sequence. | |
| void | 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 > | mg1_eg (const Blocks &Ain, bool &found) |
| G in closed form when A0 has rank one. | |
| Matrix< double > | mg1_cr (const Blocks &Ain, const Mg1CrOptions &opts=Mg1CrOptions()) |
| Cyclic reduction for M/G/1-type Markov chains [Bini, Meini]. | |
| Matrix< double > | mg1_fi (const Blocks &Ain, const Mg1FiOptions &opts=Mg1FiOptions()) |
| Functional iterations for M/G/1-type Markov chains [Neuts]. | |
| Matrix< double > | 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 > | 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 > | mg1_ni (const Blocks &Ain, const Mg1NiOptions &opts=Mg1NiOptions()) |
| Newton iteration for M/G/1-type Markov chains. | |
| Matrix< double > | mg1_rr (const Blocks &D, const Mg1RrOptions &opts=Mg1RrOptions()) |
| Ramaswami reduction for M/G/1-type Markov chains [Bini, Meini, Ramaswami]. | |
| Matrix< double > | mg1_is (const Blocks &D, const Mg1IsOptions &opts=Mg1IsOptions()) |
| Invariant subspace method for M/G/1-type Markov chains [Akar, Sohraby]. | |
| Matrix< double > | 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. | |
| using line::smc::Blocks = std::vector<Matrix<double>> |
Splits the wide [A0 A1 ... Amax] into its m x m blocks.
Definition at line 162 of file mg1.h.
References blocks_of(), line::Matrix< T >::cols(), line::InputError::InputError(), and line::Matrix< T >::rows().
Referenced by blocks_of(), mg1_g_etaqa(), mg1_pi_etaqa(), and mg1_qlen_etaqa().
|
inline |
The inner product of a row vector with a column held as a vector.
Definition at line 240 of file mg1.h.
References dot(), and line::InputError::InputError().
Referenced by dot(), and mg1_drift().
|
inline |
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.
Port of GIM1_Caudal.m. When v is given it receives the right PF eigenvector of A(z) at the last bisection point, scaled to unit sum.
Definition at line 363 of file mg1.h.
References gim1_caudal(), max_eig(), pf_vector(), and poly_at().
Referenced by gim1_caudal(), gim1_r(), mg1_eg(), and mg1_shifts().
|
inline |
Aggregated stationary vector [pi0, pi1, pi2+pi3+...] of a GI/M/1-type chain.
Port of GIM1_pi_ETAQA.m. B and A are vertical stacks; B0 is the reference's 'Boundary' option (pass empty for the default A0).
Definition at line 559 of file etaqa.h.
References line::Matrix< T >::cols(), line::eye(), gim1_pi_etaqa(), line::InputError::InputError(), line::inverse(), madd(), line::matmul(), msub(), line::NumericError::NumericError(), line::Matrix< T >::rows(), rowsums(), and vblocks_of().
Referenced by gim1_pi_etaqa(), and line::mam::solver_mam_map_bmap_1().
|
inline |
n-th moment of the level of a GI/M/1-type chain from the ETAQA aggregates.
Port of GIM1_qlen_ETAQA.m, including the scalar A(3) of defect 3.
Definition at line 677 of file etaqa.h.
References line::Matrix< T >::cols(), line::eye(), gim1_qlen_etaqa(), line::InputError::InputError(), line::inverse(), madd(), line::matmul(), mscale(), msub(), line::NumericError::NumericError(), line::Matrix< T >::rows(), rowsums(), and vblocks_of().
Referenced by gim1_qlen_etaqa(), and line::mam::solver_mam_map_bmap_1().
|
inline |
R of a GI/M/1-type Markov chain, through the G of its dual.
Port of GIM1_R.m for Dual 'A', 'R' and 'B' and every Algor: 'FI', 'CR', 'NI', 'RR', 'IS'.
THE DUAL IS THE WHOLE IDEA. There is no cyclic reduction for R directly, so the chain is transposed into an M/G/1-type one whose G carries the same information: the Ramaswami dual diag(theta)^-1 A_i' diag(theta) for a transient chain, and the Bright dual, which additionally rescales block i by eta^(i-1) with eta the caudal characteristic, for a positive recurrent one. 'A' picks between them by the drift, which is what makes it the fastest default. R is then read back off G by the inverse similarity, times eta in the Bright case.
Definition at line 1621 of file mg1.h.
References gim1_caudal(), gim1_r(), line::InputError::InputError(), madd(), mg1_cr(), mg1_decay(), mg1_drift(), mg1_fi(), mg1_is(), mg1_ni(), mg1_rr(), mscale(), stat(), line::smc::Drift::theta, and line::smc::Drift::value.
Referenced by gim1_r(), and gim1_r_etaqa().
R of a GI/M/1-type chain, uniformized first.
Port of GIM1_R_ETAQA.m.
A is the VERTICAL stack [A0; A1; ...; Amax], which is how a GI/M/1-type sequence is written; the reference transposes that stack into the horizontal one GIM1_R wants, and asks for the automatic dual with functional iterations.
Definition at line 531 of file etaqa.h.
References line::Matrix< T >::cols(), gim1_r(), gim1_r_etaqa(), line::InputError::InputError(), line::Matrix< T >::rows(), and vblocks_of().
Referenced by gim1_r_etaqa(), and line::mam::solver_mam_map_bmap_1().
Re-assembles a block sequence into the wide [A0 A1 ... Amax].
Definition at line 174 of file mg1.h.
References line::Matrix< T >::cols(), line::Matrix< T >::empty(), hcat(), line::Matrix< T >::rows(), and line::Matrix< T >::size().
Referenced by hcat().
|
inline |
norm(A,inf) over a whole block sequence stacked vertically.
Definition at line 135 of file mg1.h.
References inf_norm(), and line::Matrix< T >::size().
|
inline |
norm(A,inf), the largest absolute row sum.
Definition at line 124 of file mg1.h.
References line::Matrix< T >::cols(), inf_norm(), and line::Matrix< T >::rows().
Referenced by inf_norm(), inf_norm(), mg1_cr(), mg1_fi(), mg1_ni(), and mg1_rr().
A + B.
Templated so the point-wise CR step can add complex blocks.
Definition at line 87 of file mg1.h.
References line::Matrix< T >::cols(), line::InputError::InputError(), madd(), and line::Matrix< T >::rows().
Referenced by gim1_pi_etaqa(), gim1_qlen_etaqa(), gim1_r(), line::smc::rr_detail::Disp::left(), madd(), mg1_cr(), mg1_drift(), mg1_eg(), mg1_fi(), mg1_is(), mg1_ni(), mg1_pi_etaqa(), mg1_qlen_etaqa(), mg1_rr(), poly_at(), and line::smc::rr_detail::Disp::right().
max(max(abs(A-B))).
Definition at line 142 of file mg1.h.
References line::Matrix< T >::cols(), max_abs_diff(), and line::Matrix< T >::rows().
Referenced by max_abs_diff(), and mg1_cr().
|
inline |
max(sum(A)), the largest column sum WITHOUT absolute values, as in MATLAB.
Definition at line 151 of file mg1.h.
References line::Matrix< T >::cols(), max_col_sum(), and line::Matrix< T >::rows().
Referenced by max_col_sum(), mg1_cr(), and mg1_is().
|
inline |
max(eig(M)) with MATLAB's semantics on a complex spectrum: the element of largest modulus, ties broken by the larger phase angle.
For the nonnegative A(z) both callers evaluate, this is the Perron-Frobenius eigenvalue and is real; the comparisons below then take its real part, which is what MATLAB's relational operators do on a complex value.
Definition at line 285 of file mg1.h.
References line::eig_values(), max_eig(), and line::NumericError::NumericError().
Referenced by gim1_caudal(), max_eig(), mg1_decay(), and pf_vector().
|
inline |
Cyclic reduction for M/G/1-type Markov chains [Bini, Meini].
Port of MG1_CR.m, default mode 'ShiftPWCR' with ShiftType 'one'.
WHAT THE ALGORITHM DOES, since the transcription is otherwise opaque. One step of cyclic reduction eliminates every odd level of the chain and leaves a chain of the same M/G/1-type shape on the even ones, so the level distance halves per iteration and the iteration converges quadratically. Doing that on the block sequences directly is a polynomial composition; the reference instead evaluates the four sequences at the (nj+1)-th roots of unity, does the composition POINT-WISE (a small dense inverse per root), and interpolates back with an inverse transform – the "point-wise" in PWCR. The number of roots doubles until the interpolated tail is below (nj+1) eps, which is the reference's own accuracy control and the reason MaxNumRoot exists.
Everything runs on the TRANSPOSED blocks, as the reference does after D=D', and the final G is transposed back.
Definition at line 672 of file mg1.h.
References line::eye(), line::smc::ShiftResult::hatA, inf_norm(), line::InputError::InputError(), line::inverse(), madd(), line::matmul(), max_abs_diff(), max_col_sum(), mg1_cr(), mg1_eg(), mg1_shifts(), mg1_unshift(), msub(), line::svd_values(), line::Matrix< T >::transpose(), and line::UnsupportedError::UnsupportedError().
Referenced by gim1_r(), mg1_cr(), and mg1_g_etaqa().
|
inline |
Decay rate of a recurrent M/G/1-type chain: the unique z > 1 with PF(A(z)) = z.
Port of MG1_Decay.m. When uT is given it receives the left PF eigenvector of A(z) at the last bisection point, as the reference's second output, scaled to unit sum.
Definition at line 333 of file mg1.h.
References max_eig(), mg1_decay(), pf_vector(), and poly_at().
Referenced by gim1_r(), mg1_decay(), and mg1_shifts().
drift = theta * beta with beta = (Amax)e + (Amax+Amax-1)e + ..., the expected level increment per transition of the phase process.
Repeated verbatim in MG1_EG, MG1_Shifts, MG1_pi_ETAQA and GIM1_R, so it lives here.
Definition at line 258 of file mg1.h.
References dot(), madd(), mg1_drift(), rowsums(), stat(), line::smc::Drift::theta, and line::smc::Drift::value.
Referenced by gim1_r(), mg1_drift(), mg1_eg(), mg1_is(), mg1_pi_etaqa(), and mg1_shifts().
G in closed form when A0 has rank one.
Port of MG1_EG.m.
found is false when the shortcut does not apply, which is the reference's empty return. A rank-one A0 means every down-transition forgets the phase it came from, so G is the same rank-one matrix e beta in the recurrent case; this is not an approximation and it is why an M/M/1-shaped input never enters cyclic reduction at all.
Definition at line 517 of file mg1.h.
References line::eye(), gim1_caudal(), line::inverse(), madd(), line::matmul(), line::matrix_rank(), mg1_drift(), mg1_eg(), mscale(), msub(), rowsums(), line::smc::Drift::theta, and line::smc::Drift::value.
Referenced by mg1_cr(), mg1_eg(), mg1_fi(), mg1_is(), mg1_ni(), and mg1_rr().
|
inline |
Functional iterations for M/G/1-type Markov chains [Neuts].
Port of MG1_FI.m for the three modes and the shift variants; the reference's NonZeroBlocks option is not exposed, because it changes only which products are skipped when some A_i vanish and converges to the same G.
'U-Based' is the default and the one GIM1_R(...,'FI') uses: it solves G = (I - sum_{j>=1} A_j G^{j-1})^{-1} A0, which is the U-based iteration and converges monotonically from below to the minimal nonnegative solution.
Definition at line 863 of file mg1.h.
References line::eye(), line::smc::ShiftResult::hatA, inf_norm(), line::inverse(), madd(), line::matmul(), mg1_eg(), mg1_fi(), mg1_shifts(), mg1_unshift(), msub(), and line::UnsupportedError::UnsupportedError().
G of an M/G/1-type chain, uniformized first.
Port of MG1_G_ETAQA.m.
A is the wide [A0 A1 ... Amax]. The generator is turned into the transition matrix of the uniformized chain by dividing through by -min(diag(A1)) and adding the identity back onto A1, which is what cyclic reduction expects; see defect 1 in the header for why that branch is unconditional.
Definition at line 199 of file etaqa.h.
References blocks_of(), line::Matrix< T >::cols(), line::InputError::InputError(), mg1_cr(), mg1_g_etaqa(), and line::Matrix< T >::rows().
Referenced by mg1_g_etaqa(), and line::mam::solver_mam_bmap_map_1().
|
inline |
Invariant subspace method for M/G/1-type Markov chains [Akar, Sohraby].
Port of MG1_IS.m.
F(z) = zI - A(z) is mapped by the Moebius transform z = (1+s)/(1-s) onto a matrix polynomial H(s) whose companion matrix, after a rank-one correction that moves the unit root off the imaginary axis, has exactly m eigenvalues in the open left half plane. A basis T of that invariant subspace gives G = (T1 + T2)(T1 - T2)^{-1}, invariant to the choice of basis. The matrix-sign modes find the subspace as the range of sign(Z) - I; 'Schur' orders a real Schur form with the left half plane first, as the reference's ordschur 'lhp' does.
Definition at line 1487 of file mg1.h.
References line::eye(), line::InputError::InputError(), line::inverse(), line::lstsq(), line::lu_det(), madd(), line::matmul(), max_col_sum(), mg1_drift(), mg1_eg(), mg1_is(), mscale(), msub(), line::NumericError::NumericError(), line::SvdFactors::s, line::schur_decomposition(), line::schur_reorder(), line::svd_full(), line::RealSchur::T, line::SvdFactors::U, line::smc::Drift::value, line::vecmul(), and line::RealSchur::Z.
|
inline |
Newton iteration for M/G/1-type Markov chains.
Port of MG1_NI.m.
Each step linearizes G = sum_i A_i G^i at the current G and solves the resulting Sylvester-power equation sum_j B_j Y G^{j-1} = G - B_0 for the correction Y, so the iteration converges quadratically; the shift variants run it on the sequence with the unit root moved to zero (mg1_shifts) and put the removed rank-one term back at the end.
'ComplexSchur' asks the reference for a complex Schur form of G. It solves the SAME linear system as 'RealSchur', whose real quasi-triangular form this port uses for both, so the two modes differ only in rounding: there is no complex Schur factorization in this tree.
Definition at line 1146 of file mg1.h.
References line::smc::ShiftResult::hatA, inf_norm(), line::InputError::InputError(), madd(), line::matmul(), mg1_eg(), mg1_ni(), mg1_shifts(), mg1_unshift(), msub(), sylv_powers_direct(), and sylv_powers_real_schur().
|
inline |
Aggregated stationary vector [pi0, pi1, pi2+pi3+...] of an M/G/1-type chain.
Port of MG1_pi_ETAQA.m.
B may be empty, in which case the boundary repeats the repetitive blocks. C0 is the reference's 'Boundary' option, the block that takes level 1 back to a boundary of a different size; pass an empty matrix for the default A0.
Definition at line 230 of file etaqa.h.
References blocks_of(), line::Matrix< T >::cols(), line::Matrix< T >::empty(), line::InputError::InputError(), madd(), line::matmul(), line::matrix_rank(), mg1_drift(), mg1_pi_etaqa(), line::NumericError::NumericError(), line::Matrix< T >::rows(), rowsums(), and line::smc::Drift::value.
Referenced by mg1_pi_etaqa(), and line::mam::solver_mam_bmap_map_1().
|
inline |
n-th moment of the level (the queue length) of an M/G/1-type chain from the ETAQA aggregates.
Port of MG1_qlen_ETAQA.m.
Definition at line 382 of file etaqa.h.
References blocks_of(), line::Matrix< T >::cols(), line::Matrix< T >::empty(), line::InputError::InputError(), madd(), mg1_qlen_etaqa(), mscale(), and line::Matrix< T >::rows().
Referenced by mg1_qlen_etaqa(), and line::mam::solver_mam_bmap_map_1().
|
inline |
Ramaswami reduction for M/G/1-type Markov chains [Bini, Meini, Ramaswami].
Port of MG1_RR.m.
The chain is reduced to a QBD whose level is N blocks wide, and cyclic reduction on that QBD is run in the Sherman-Morrison-Woodbury form of the reference, carrying only the rank-m corrections uhat and vT; G is then read off by the formula at the end of Section 3 of the paper, G = (I - vT uhat)^{-1} A0.
'Direct' (the default) stores the N*m x N*m matrix B outright. 'DispStruct' stores it through its displacement generators (rr_detail::Disp), compressed to rank 2m by a QR of each side and an SVD of the small core at every step. 'DispStructFFT' is the reference's FFT evaluation of the same block-Toeplitz products; this port evaluates those products directly, so it returns the 'DispStruct' answer, which the FFT reproduces up to rounding.
Definition at line 1308 of file mg1.h.
References line::smc::rr_detail::Disp::b, line::smc::rr_detail::Disp::c1, line::smc::rr_detail::Disp::c2, line::Matrix< T >::cols(), line::eye(), inf_norm(), line::InputError::InputError(), line::inverse(), line::smc::rr_detail::Disp::left(), madd(), line::matmul(), mg1_eg(), mg1_rr(), mscale(), msub(), line::smc::rr_detail::Disp::r1, line::smc::rr_detail::Disp::r2, line::smc::rr_detail::Disp::right(), line::Matrix< T >::rows(), line::SvdFactors::s, line::smc::rr_detail::sum_all(), line::svd_full(), line::SvdFactors::U, and line::SvdFactors::Vt.
|
inline |
Shift technique for the M/G/1-type sequence.
Port of MG1_Shifts.m, ShiftType 'one' (the default), 'tau' and 'dbl'.
For a positive recurrent chain (drift < 1) 'one' shifts the eigenvalue 1 of A(z) to zero by subtracting (A0+...+Ai)e u^T from block i with u^T = e^T/m, and 'tau' shifts the decay rate tau to infinity by subtracting e rowhatA_i; for a transient one (drift >= 1) 'one' shifts 1 to infinity through the row theta(Amax+...+Ai) and 'tau' shifts the caudal value to zero through the column built on its right eigenvector v. 'dbl' applies both, in the reference's order. The solver converges in the next root instead of stalling on the removed one and puts the removed rank-one term back on G (mg1_unshift).
Definition at line 404 of file mg1.h.
References line::smc::ShiftResult::drift, gim1_caudal(), line::smc::ShiftResult::hatA, line::InputError::InputError(), mg1_decay(), mg1_drift(), mg1_shifts(), line::mulvec(), rowsums(), line::smc::ShiftResult::tau, line::smc::Drift::theta, line::smc::ShiftResult::v, line::smc::Drift::value, and line::vecmul().
Referenced by mg1_cr(), mg1_fi(), mg1_ni(), and mg1_shifts().
|
inline |
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'.
The reference repeats this switch verbatim in MG1_CR, MG1_FI and MG1_NI.
Definition at line 492 of file mg1.h.
References line::smc::ShiftResult::drift, mg1_unshift(), line::Matrix< T >::rows(), line::smc::ShiftResult::tau, and line::smc::ShiftResult::v.
Referenced by mg1_cr(), mg1_fi(), mg1_ni(), and mg1_unshift().
c * A.
Definition at line 108 of file mg1.h.
References line::Matrix< T >::cols(), mscale(), and line::Matrix< T >::rows().
Referenced by gim1_qlen_etaqa(), gim1_r(), mg1_eg(), mg1_is(), mg1_qlen_etaqa(), mg1_rr(), mscale(), and poly_at().
A - B.
Definition at line 98 of file mg1.h.
References line::Matrix< T >::cols(), line::InputError::InputError(), msub(), and line::Matrix< T >::rows().
Referenced by gim1_pi_etaqa(), gim1_qlen_etaqa(), mg1_cr(), mg1_eg(), mg1_fi(), mg1_is(), mg1_ni(), mg1_rr(), and msub().
|
inline |
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.
The reference picks the column of eig at the largest eigenvalue; both callers normalize it by its sum, so the scale and sign eig happens to return do not matter.
Definition at line 309 of file mg1.h.
References max_eig(), pf_vector(), line::Matrix< T >::rows(), line::svd_full(), and line::SvdFactors::Vt.
Referenced by gim1_caudal(), mg1_decay(), and pf_vector().
A(z) = A0 + A1 z + ... + Amax z^max, by Horner as the reference writes it.
Definition at line 297 of file mg1.h.
References madd(), mscale(), and poly_at().
Referenced by gim1_caudal(), mg1_decay(), and poly_at().
|
inline |
sum(A,2), the row sums, as a column held in a vector.
Definition at line 116 of file mg1.h.
References line::Matrix< T >::cols(), line::Matrix< T >::rows(), and rowsums().
Referenced by gim1_pi_etaqa(), gim1_qlen_etaqa(), mg1_drift(), mg1_eg(), mg1_pi_etaqa(), mg1_shifts(), and rowsums().
|
inline |
theta A, the row vector times matrix product used throughout.
Definition at line 235 of file mg1.h.
References rowvec_times(), and line::vecmul().
Referenced by rowvec_times().
|
inline |
Stationary distribution of a stochastic matrix: the left eigenvector for eigenvalue 1, nonnegative and summing to one.
Port of stat.m, including its shape: [A - I, e] is S x (S+1), so the reference's y / B is a least-squares solve of a consistent overdetermined system and not a square solve. The normalization is IN the system (the appended column of ones against the appended 1 on the right), which is why the result needs no rescaling afterwards.
Definition at line 221 of file mg1.h.
References line::Matrix< T >::cols(), line::InputError::InputError(), line::lstsq(), line::Matrix< T >::rows(), and stat().
Referenced by gim1_r(), mg1_drift(), and stat().
|
inline |
Solve sum_{j=1}^N B_j Y A^{j-1} = C directly, through the Kronecker form of vec(Y).
Port of solveSylvPowersDirectSum.m; B holds the N blocks B_1..B_N (n x n), A is m x m.
Definition at line 1024 of file mg1.h.
References line::eye(), line::smc::mg1x_detail::lsolve(), line::matmul(), line::Matrix< T >::rows(), sylv_powers_direct(), and line::smc::mg1x_detail::trans().
Referenced by mg1_ni(), and sylv_powers_direct().
|
inline |
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.
Port of solveSylvPowersRealSchur_FW.m, including its 1e-13 test on the subdiagonal that tells the two block shapes apart.
Definition at line 1057 of file mg1.h.
References line::eye(), line::smc::mg1x_detail::lsolve(), line::matmul(), line::smc::mg1x_detail::put(), line::Matrix< T >::rows(), line::schur_decomposition(), sylv_powers_real_schur(), line::RealSchur::T, line::smc::mg1x_detail::trans(), and line::RealSchur::Z.
Referenced by mg1_ni(), and sylv_powers_real_schur().
Splits a vertical stack into its blocks of r rows.
Definition at line 196 of file mg1.h.
References line::Matrix< T >::cols(), line::InputError::InputError(), line::Matrix< T >::rows(), and vblocks_of().
Referenced by gim1_pi_etaqa(), gim1_qlen_etaqa(), gim1_r_etaqa(), and vblocks_of().
Stacks a block sequence vertically, [A0; A1; ...; Amax].
Definition at line 185 of file mg1.h.
References line::Matrix< T >::cols(), line::Matrix< T >::empty(), line::Matrix< T >::rows(), line::Matrix< T >::size(), and vcat().
Referenced by vcat().
1.18.0