qsys_mmapgk1_lcfspr.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_QSYS_QSYS_MMAPGK1_LCFSPR_H
#define LINE_API_QSYS_QSYS_MMAPGK1_LCFSPR_H
 
/**
 * @file
 * @ingroup api_qsys
 * The MMAP[K]/G[K]/1 LCFS-PR queue: K customer types with class-dependent
 * GENERAL service under preemptive-resume last-come-first-served, fed by a
 * marked Markovian arrival process. Port of qsys_mmapgk1_lcfspr.m.
 *
 * THE METHOD, Takine (2002). With C = D0, D_k = D^(k) and F_k the service law
 * of type k, let Q be the minimal solution of
 *
 *     Q = C + sum_k D_k F_k(Q),      F_k(Q) = int_0^inf exp(Qx) dF_k(x),
 *
 * the generator of the arrival phase with the busy periods excised (He's U in
 * qsys_mmapgk1.h), and kappa its stationary vector. With beta_k = D_k e E[S_k],
 * Takine's Theorem 1 gives the rate matrices
 *
 *     R_k = (beta_k kappa + D_k (I - F_k(Q))) (e kappa - Q)^-1,  R = sum_k R_k,
 *
 * and the per-class numbers in system have the matrix product form
 * y(n) = sum_k y(n - e_k) R_k, y(0) = (1 - rho) kappa. From it
 *
 *     P(N = n)   = (1 - rho) kappa R^n e,
 *     P(Y_k = m) = (1 - rho) kappa (A_k R_k)^m A_k e,  A_k = (I - R + R_k)^-1,
 *     E[Y_k]     = pi R_k (I - R)^-1 e,
 *
 * with pi the stationary vector of C + D. A type-k sojourn under LCFS-PR is the
 * sub-busy period its own service opens, so
 *
 *     E[exp(-s W_k)] = pi D_k F_k(Q*(s)) e / lambda_k,
 *     Q*(s) = C - sI + sum_j D_j F_j(Q*(s)),
 *
 * evaluated at COMPLEX s for the Euler inversion of the CDF. That is why this
 * header carries a complex matrix transform: the real one in qsys_mmapgk1.h
 * returns real(V diag(g) V^-1), which is right only for a real argument.
 *
 * Reference:
 *   T. Takine, "Matrix product-form solution for an LCFS-PR single-server queue
 *   with multiple arrival streams governed by a Markov chain", Queueing Systems
 *   42(2):131-151, 2002.
 */
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstdio>
#include <functional>
#include <limits>
#include <string>
#include <vector>
 
#include "line/api/qsys/qsys_mmapgk1.h"
 
namespace line {
namespace qsys {
 
/** Return value of qsys_mmapgk1_lcfspr, mirroring the MATLAB struct. */
template <class T>
struct MmapGk1LcfsprResult {
    std::vector<T> lambdas;                              ///< per-type arrival rates
    T arrivalRate;                                       ///< sum of the per-type rates
    T utilization;                                       ///< rho = sum_k lambda_k E[S_k]
    Matrix<T> Q;                                         ///< Takine's Q (He's U)
    std::vector<T> kappa;                                ///< stationary vector of Q
    std::vector<Matrix<T> > Rk;                          ///< the K rate matrices
    Matrix<T> R;                                         ///< sum of the R_k
    std::vector<T> classMeanQueueLength;                 ///< E[Y_k]
    T meanQueueLength;                                   ///< E[N]
    std::vector<T> meanSojournTime;                      ///< E[W_k]
    std::vector<T> queueLengthDistr;                     ///< P(N = n), n = 0..nMax
    std::vector<std::vector<T> > classQueueLengthDistr;  ///< [type][m], P(Y_k = m)
    /// s -> per-type E[exp(-s W_k)] at a complex argument
    std::function<std::vector<std::complex<double> >(std::complex<double>)> respLST;
    std::vector<std::vector<T> > respCDF;                ///< [type][point], P(W_k <= t)
    std::vector<T> respPoints;                           ///< the requested points
    std::string analyzer;                                ///< "LINE:MMAP[K]/G[K]/1-LCFSPR"
};
 
/** The MATLAB name-value options. nMax < 0 means "not requested". */
template <class T>
struct MmapGk1LcfsprOptions {
    long nMax = -1;
    std::vector<T> wPoints;
    std::size_t quadNodes = 120;
    double tol = 1e-12;
};
 
namespace mmapgk1lcfsprdetail {
 
typedef std::complex<double> cdbl;
typedef std::vector<std::vector<cdbl> > CMat;
 
inline CMat czeros(std::size_t n) { return CMat(n, std::vector<cdbl>(n, cdbl(0.0, 0.0))); }
 
inline CMat ceye(std::size_t n) {
    CMat I = czeros(n);
    for (std::size_t i = 0; i < n; ++i) I[i][i] = cdbl(1.0, 0.0);
    return I;
}
 
inline CMat cmatmul(const CMat& A, const CMat& B) {
    const std::size_t n = A.size(), m = B.empty() ? 0 : B[0].size(), p = B.size();
    CMat C(n, std::vector<cdbl>(m, cdbl(0.0, 0.0)));
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t q = 0; q < p; ++q) {
            const cdbl a = A[i][q];
            if (a == cdbl(0.0, 0.0)) continue;
            for (std::size_t j = 0; j < m; ++j) C[i][j] += a * B[q][j];
        }
    return C;
}
 
template <class T>
CMat to_cmat(const Matrix<T>& A) {
    CMat C(A.rows(), std::vector<cdbl>(A.cols(), cdbl(0.0, 0.0)));
    for (std::size_t i = 0; i < A.rows(); ++i)
        for (std::size_t j = 0; j < A.cols(); ++j)
            C[i][j] = cdbl(num_traits<T>::to_double(A(i, j)), 0.0);
    return C;
}
 
/** X = A^-1 B over the complex field, column by column. */
inline CMat csolve(const CMat& A, const CMat& B) {
    const std::size_t n = A.size(), m = B.empty() ? 0 : B[0].size();
    CMat X(n, std::vector<cdbl>(m, cdbl(0.0, 0.0)));
    for (std::size_t j = 0; j < m; ++j) {
        std::vector<cdbl> b(n);
        for (std::size_t i = 0; i < n; ++i) b[i] = B[i][j];
        const std::vector<cdbl> x = mmapgk1detail::complex_solve(A, b);
        for (std::size_t i = 0; i < n; ++i) X[i][j] = x[i];
    }
    return X;
}
 
/** Complex matrix exponential, scaling and squaring over a Taylor series. */
inline CMat cexpm(CMat A) {
    const std::size_t n = A.size();
    double nrm = 0.0;
    for (std::size_t j = 0; j < n; ++j) {
        double cs = 0.0;
        for (std::size_t i = 0; i < n; ++i) cs += std::abs(A[i][j]);
        nrm = std::max(nrm, cs);
    }
    int sq = 0;
    if (nrm > 0.5) sq = static_cast<int>(std::ceil(std::log2(nrm / 0.5)));
    const double scale = std::ldexp(1.0, -sq);
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) A[i][j] *= scale;
    CMat E = ceye(n), term = ceye(n);
    for (int k = 1; k <= 30; ++k) {
        term = cmatmul(term, A);
        double tn = 0.0;
        for (std::size_t i = 0; i < n; ++i)
            for (std::size_t j = 0; j < n; ++j) {
                term[i][j] /= static_cast<double>(k);
                E[i][j] += term[i][j];
                tn = std::max(tn, std::abs(term[i][j]));
            }
        if (tn < 1e-18) break;
    }
    for (int s = 0; s < sq; ++s) E = cmatmul(E, E);
    return E;
}
 
/**
 * U = V diag(ev) V^-1 for a COMPLEX matrix, the reference's `eig(U)` at a
 * complex argument. The same LAPACK-free route as mmapgk1detail::eig_decomp
 * (Faddeev-LeVerrier, Durand-Kerner, inverse iteration), with complex
 * coefficients throughout. `rcond` is the reciprocal 1-norm condition number
 * of V, 0 when V could not be inverted.
 */
inline void ceig_decomp(const CMat& U, std::vector<cdbl>& ev, CMat& V, CMat& Vinv, double& rcond) {
    const std::size_t n = U.size();
    ev.assign(n, cdbl(0.0, 0.0));
    V = czeros(n);
    Vinv = czeros(n);
    rcond = 0.0;
    if (n == 0) return;
    std::vector<cdbl> c(n + 1, cdbl(0.0, 0.0));
    c[0] = cdbl(1.0, 0.0);
    CMat Mk = ceye(n);
    for (std::size_t k = 1; k <= n; ++k) {
        const CMat TM = cmatmul(U, Mk);
        cdbl tr(0.0, 0.0);
        for (std::size_t i = 0; i < n; ++i) tr += TM[i][i];
        c[k] = -tr / static_cast<double>(k);
        Mk = TM;
        for (std::size_t i = 0; i < n; ++i) Mk[i][i] += c[k];
    }
    double scale = 0.0;
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) scale = std::max(scale, std::abs(U[i][j]));
    {
        const cdbl seed(0.4, 0.9);
        cdbl pw(1.0, 0.0);
        for (std::size_t i = 0; i < n; ++i) {
            ev[i] = pw * (scale > 0.0 ? scale : 1.0);
            pw *= seed;
        }
        for (int it = 0; it < 2000; ++it) {
            double move = 0.0;
            for (std::size_t i = 0; i < n; ++i) {
                cdbl den(1.0, 0.0);
                for (std::size_t j = 0; j < n; ++j)
                    if (j != i) den *= (ev[i] - ev[j]);
                if (std::abs(den) < 1e-300) continue;
                cdbl pv(0.0, 0.0);
                for (std::size_t q = 0; q <= n; ++q) pv = pv * ev[i] + c[q];
                const cdbl step = pv / den;
                ev[i] -= step;
                move = std::max(move, std::abs(step));
            }
            if (move < 1e-15 * std::max(1.0, scale)) break;
        }
    }
    double umax = 0.0;
    for (std::size_t i = 0; i < n; ++i) umax = std::max(umax, std::abs(ev[i]));
    for (std::size_t i = 0; i < n; ++i) {
        const double eps = 1e-10 * std::max(1.0, umax);
        const cdbl shift = ev[i] + cdbl(eps, eps);
        std::vector<cdbl> v(n, cdbl(1.0 / std::sqrt(static_cast<double>(n)), 0.0));
        for (int it = 0; it < 3; ++it) {
            CMat A = U;
            for (std::size_t p = 0; p < n; ++p) A[p][p] -= shift;
            const std::vector<cdbl> w = mmapgk1detail::complex_solve(A, v);
            double nrm = 0.0;
            for (std::size_t q = 0; q < n; ++q) nrm += std::norm(w[q]);
            nrm = std::sqrt(nrm);
            if (!(nrm > 0.0) || !std::isfinite(nrm)) break;
            for (std::size_t q = 0; q < n; ++q) v[q] = w[q] / nrm;
        }
        for (std::size_t q = 0; q < n; ++q) V[q][i] = v[q];
    }
    double normV = 0.0, normW = 0.0;
    for (std::size_t j = 0; j < n; ++j) {
        std::vector<cdbl> ej(n, cdbl(0.0, 0.0));
        ej[j] = cdbl(1.0, 0.0);
        const std::vector<cdbl> col = mmapgk1detail::complex_solve(V, ej);
        double cs = 0.0, vs = 0.0;
        for (std::size_t i = 0; i < n; ++i) {
            Vinv[i][j] = col[i];
            cs += std::abs(col[i]);
            vs += std::abs(V[i][j]);
        }
        normW = std::max(normW, cs);
        normV = std::max(normV, vs);
    }
    if (std::isfinite(normW) && normW > 0.0 && normV > 0.0) rcond = 1.0 / (normV * normW);
}
 
/** V diag(g) V^-1, kept complex. */
inline CMat ceig_apply(const CMat& V, const CMat& Vinv, const std::vector<cdbl>& g) {
    const std::size_t n = V.size();
    CMat F = czeros(n);
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) {
            cdbl acc(0.0, 0.0);
            for (std::size_t q = 0; q < n; ++q) acc += V[i][q] * g[q] * Vinv[q][j];
            F[i][j] = acc;
        }
    return F;
}
 
/** int exp(Ut) dF(t) for a (D0, D1) phase-type law at a complex U, on the Kronecker sum. */
template <class T>
CMat cmatrix_lst_map(const mam::Map<T>& rep, const CMat& U) {
    const std::size_t n = U.size();
    const std::size_t ms = rep.D0.rows();
    const std::vector<T> beta = mam::map_pie(rep);
    const std::size_t N = n * ms;
    CMat KS = czeros(N), Rhs(N, std::vector<cdbl>(n, cdbl(0.0, 0.0)));
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j)
            for (std::size_t p = 0; p < ms; ++p) KS[i * ms + p][j * ms + p] += U[i][j];
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t p = 0; p < ms; ++p) {
            double s0 = 0.0;
            for (std::size_t q = 0; q < ms; ++q) {
                const double d = num_traits<T>::to_double(rep.D0(p, q));
                KS[i * ms + p][i * ms + q] += d;
                s0 -= d;
            }
            Rhs[i * ms + p][i] = cdbl(s0, 0.0);
        }
    const CMat X = csolve(KS, Rhs);
    CMat F = czeros(n);
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) {
            cdbl acc(0.0, 0.0);
            for (std::size_t p = 0; p < ms; ++p)
                acc += num_traits<T>::to_double(beta[p]) * X[i * ms + p][j];
            F[i][j] = -acc;
        }
    return F;
}
 
/** The matrix transform of a Distribution object at a complex U. */
template <class T>
CMat cmatrix_lst(const lang::Distrib<T>& law, const CMat& U, std::size_t quad_nodes) {
    const std::size_t n = U.size();
    if (lang::process_is_markovian(law.type)) return cmatrix_lst_map(lang::dist_to_map(law), U);
    if (law.type == lang::ProcessType::DET) {
        const double d = num_traits<T>::to_double(lang::dist_moment(law, 1u));
        CMat Ud = U;
        for (std::size_t i = 0; i < n; ++i)
            for (std::size_t j = 0; j < n; ++j) Ud[i][j] *= d;
        return cexpm(Ud);
    }
    {
        std::vector<cdbl> ev;
        CMat V, Vinv;
        double rcond = 0.0;
        ceig_decomp(U, ev, V, Vinv, rcond);
        bool ok = rcond > 1e-12;
        std::vector<cdbl> g(n, cdbl(0.0, 0.0));
        if (ok) {
            for (std::size_t i = 0; i < n; ++i) {
                try {
                    g[i] = lang::dist_lst(law, -ev[i]);
                } catch (const Error&) {
                    ok = false;
                    break;
                }
                if (!std::isfinite(g[i].real()) || !std::isfinite(g[i].imag())) {
                    ok = false;
                    break;
                }
            }
        }
        if (ok) return ceig_apply(V, Vinv, g);
    }
    std::vector<T> x, w;
    mmapgk1detail::stieltjes_nodes(law, x, w, quad_nodes);
    CMat F = czeros(n);
    for (std::size_t i = 0; i < x.size(); ++i) {
        const double xi = num_traits<T>::to_double(x[i]);
        const double wi = num_traits<T>::to_double(w[i]);
        if (wi == 0.0) continue;
        CMat Ut = U;
        for (std::size_t p = 0; p < n; ++p)
            for (std::size_t q = 0; q < n; ++q) Ut[p][q] *= xi;
        const CMat E = cexpm(Ut);
        for (std::size_t p = 0; p < n; ++p)
            for (std::size_t q = 0; q < n; ++q) F[p][q] += wi * E[p][q];
    }
    return F;
}
 
/** The matrix transform of a service law in any of its three forms, at a complex U. */
template <class T>
CMat cmatrix_lst(const MmapGk1Service<T>& svc, const CMat& U, std::size_t quad_nodes) {
    switch (svc.kind) {
        case MmapGk1Service<T>::MAP:
            return cmatrix_lst_map(svc.map, U);
        case MmapGk1Service<T>::TRANSFORM: {
            const std::size_t n = U.size();
            std::vector<cdbl> ev;
            CMat V, Vinv;
            double rcond = 0.0;
            ceig_decomp(U, ev, V, Vinv, rcond);
            std::vector<cdbl> g(n);
            for (std::size_t i = 0; i < n; ++i) g[i] = svc.transform.lst(-ev[i]);
            return ceig_apply(V, Vinv, g);
        }
        default:
            return cmatrix_lst(svc.dist, U, quad_nodes);
    }
}
 
/**
 * Per-type sojourn LST pi D_k F_k(Q*(s)) e / lambda_k, Q*(s) the minimal
 * solution of Q* = C - sI + sum_j D_j F_j(Q*), iterated from C - sI.
 */
template <class T>
std::vector<cdbl> sojourn_lst(const Matrix<T>& D0, const std::vector<Matrix<T> >& Dk,
                              const std::vector<MmapGk1Service<T> >& svc,
                              const std::vector<T>& piArr, const std::vector<T>& lambdas, cdbl s,
                              std::size_t quad_nodes, double tol) {
    const std::size_t K = Dk.size(), ma = D0.rows();
    if (std::abs(s) < 1e-14) return std::vector<cdbl>(K, cdbl(1.0, 0.0));
    std::vector<CMat> Dc(K);
    for (std::size_t k = 0; k < K; ++k) Dc[k] = to_cmat(Dk[k]);
    CMat base = to_cmat(D0);
    for (std::size_t i = 0; i < ma; ++i) base[i][i] -= s;
    CMat Qs = base;
    for (std::size_t it = 0; it < 10000; ++it) {
        CMat Qnew = base;
        for (std::size_t j = 0; j < K; ++j) {
            const CMat P = cmatmul(Dc[j], cmatrix_lst(svc[j], Qs, quad_nodes));
            for (std::size_t a = 0; a < ma; ++a)
                for (std::size_t b = 0; b < ma; ++b) Qnew[a][b] += P[a][b];
        }
        double diff = 0.0;
        for (std::size_t a = 0; a < ma; ++a)
            for (std::size_t b = 0; b < ma; ++b) diff = std::max(diff, std::abs(Qnew[a][b] - Qs[a][b]));
        Qs = Qnew;
        if (diff <= tol) break;
    }
    std::vector<cdbl> w(K, cdbl(1.0, 0.0));
    for (std::size_t k = 0; k < K; ++k) {
        const double lam = num_traits<T>::to_double(lambdas[k]);
        if (!(lam > 0.0)) continue;
        const CMat F = cmatrix_lst(svc[k], Qs, quad_nodes);
        // pi D_k F e
        std::vector<cdbl> Fe(ma, cdbl(0.0, 0.0));
        for (std::size_t a = 0; a < ma; ++a)
            for (std::size_t b = 0; b < ma; ++b) Fe[a] += F[a][b];
        cdbl acc(0.0, 0.0);
        for (std::size_t a = 0; a < ma; ++a) {
            cdbl row(0.0, 0.0);
            for (std::size_t b = 0; b < ma; ++b) row += Dc[k][a][b] * Fe[b];
            acc += num_traits<T>::to_double(piArr[a]) * row;
        }
        w[k] = acc / lam;
    }
    return w;
}
 
/** Abate-Whitt Euler inversion of a CDF from its LST, qsys_mmapgk1_invert.m. */
inline double euler_invert(const std::function<cdbl(cdbl)>& lst, double t) {
    if (t <= 0.0) return lst(cdbl(1e12, 0.0)).real();
    const double A = 18.4;
    const std::size_t nE = 15, mE = 11;
    const double u = std::exp(A / 2) / t;
    const double x = A / (2 * t);
    std::vector<double> terms(nE + mE + 1, 0.0);
    terms[0] = lst(cdbl(x, 0.0)).real() / x / 2.0;
    for (std::size_t j = 1; j <= nE + mE; ++j) {
        const cdbl s(x, M_PI * static_cast<double>(j) / t);
        terms[j] = ((j % 2 == 0) ? 1.0 : -1.0) * (lst(s) / s).real();
    }
    double run = 0.0, F = 0.0;
    std::vector<double> partial(terms.size(), 0.0);
    for (std::size_t j = 0; j < terms.size(); ++j) {
        run += terms[j];
        partial[j] = run;
    }
    for (std::size_t j = 0; j <= mE; ++j)
        F += mmapgk1detail::binom(mE, j) / std::pow(2.0, static_cast<double>(mE)) * partial[nE + j];
    F = u * F;
    return std::min(std::max(F, 0.0), 1.0);
}
 
}  // namespace mmapgk1lcfsprdetail
 
/**
 * MMAP[K]/G[K]/1 LCFS-PR, per type, by Takine's matrix product form.
 *
 * @param MMAP LINE convention {D0, D1, D^(1), ..., D^(K)}
 * @param svc  K service laws, one per marked type, each a Distribution object,
 *             a (D0, D1) phase-type pair or a transform struct
 * @param opt  nMax (distributions on 0..nMax when >= 0), wPoints (sojourn CDF
 *             points), quadNodes (default 120), tol (default 1e-12)
 */
template <class T>
MmapGk1LcfsprResult<T> qsys_mmapgk1_lcfspr(const std::vector<Matrix<T> >& MMAP,
                                           const std::vector<MmapGk1Service<T> >& svc,
                                           const MmapGk1LcfsprOptions<T>& opt) {
    static_assert(num_traits<T>::has_transcendental,
                  "qsys_mmapgk1_lcfspr requires transcendental arithmetic");
    using namespace mmapgk1detail;
    using mmapgk1lcfsprdetail::cdbl;
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    if (MMAP.size() < 3)
        throw InputError("qsys_mmapgk1_lcfspr: The MMAP must carry at least one marked arrival block.");
    const std::size_t K = MMAP.size() - 2;
    if (svc.size() != K)
        throw InputError("qsys_mmapgk1_lcfspr: One service law per marked type is required (" +
                         std::to_string(svc.size()) + " given, " + std::to_string(K) + " types).");
    const std::size_t quad = opt.quadNodes;
    const double tol = opt.tol;
    const Matrix<T>& D0 = MMAP[0];
    const std::size_t ma = D0.rows();
    const Matrix<T> I = eye<T>(ma);
    const std::vector<T> e(ma, one);
    std::vector<Matrix<T> > Dk;
    Matrix<T> Dsum = D0;
    for (std::size_t k = 0; k < K; ++k) {
        Dk.push_back(MMAP[k + 2]);
        for (std::size_t i = 0; i < ma; ++i)
            for (std::size_t j = 0; j < ma; ++j) Dsum(i, j) += Dk[k](i, j);
    }
 
    const std::vector<T> piArr = stat_left_null(Dsum);
    std::vector<T> lambdas(K, zero), meanS(K, zero);
    T rho = zero;
    for (std::size_t k = 0; k < K; ++k) {
        const std::vector<T> De = mulvec(Dk[k], e);
        T lam = zero;
        for (std::size_t i = 0; i < ma; ++i) lam += piArr[i] * De[i];
        lambdas[k] = lam;
        meanS[k] = raw_moment(svc[k], 1);
        rho += lam * meanS[k];
    }
    if (num_traits<T>::to_double(rho) >= 1.0 - std::numeric_limits<double>::epsilon()) {
        char buf[64];
        std::snprintf(buf, sizeof(buf), "%g", num_traits<T>::to_double(rho));
        throw InputError(std::string("qsys_mmapgk1_lcfspr: The load ") + buf +
                         " of the system is not below one.");
    }
 
    // Q = C + sum_k D_k F_k(Q), iterated from C
    Matrix<T> Q = D0;
    std::vector<Matrix<T> > FQ(K);
    for (std::size_t it = 0; it < 10000; ++it) {
        Matrix<T> Qnew = D0;
        for (std::size_t k = 0; k < K; ++k) {
            FQ[k] = matrix_lst(svc[k], Q, quad);
            const Matrix<T> P = matmul(Dk[k], FQ[k]);
            for (std::size_t i = 0; i < ma; ++i)
                for (std::size_t j = 0; j < ma; ++j) Qnew(i, j) += P(i, j);
        }
        double diff = 0.0;
        for (std::size_t i = 0; i < ma; ++i)
            for (std::size_t j = 0; j < ma; ++j)
                diff = std::max(diff, std::fabs(num_traits<T>::to_double(Qnew(i, j) - Q(i, j))));
        Q = Qnew;
        if (diff <= tol) break;
    }
    // Q e = 0 at the fixed point but only to tolerance in the iterate: restore it
    // so kappa and (e kappa - Q)^-1 are those of a proper generator.
    for (std::size_t i = 0; i < ma; ++i) {
        T rowsum = zero;
        for (std::size_t j = 0; j < ma; ++j) rowsum += Q(i, j);
        Q(i, i) -= rowsum;
    }
    for (std::size_t k = 0; k < K; ++k) FQ[k] = matrix_lst(svc[k], Q, quad);
    const std::vector<T> kappa = stat_left_null(Q);
 
    // Rate matrices (Takine, Theorem 1)
    Matrix<T> ekQ(ma, ma, zero);
    for (std::size_t i = 0; i < ma; ++i)
        for (std::size_t j = 0; j < ma; ++j) ekQ(i, j) = kappa[j] - Q(i, j);
    const Matrix<T> Z = inverse(ekQ);
    std::vector<Matrix<T> > Rk(K);
    Matrix<T> R(ma, ma, zero);
    for (std::size_t k = 0; k < K; ++k) {
        const std::vector<T> De = mulvec(Dk[k], e);
        Matrix<T> IF = I;
        for (std::size_t i = 0; i < ma; ++i)
            for (std::size_t j = 0; j < ma; ++j) IF(i, j) -= FQ[k](i, j);
        Matrix<T> Bk = matmul(Dk[k], IF);
        for (std::size_t i = 0; i < ma; ++i)
            for (std::size_t j = 0; j < ma; ++j) Bk(i, j) += De[i] * meanS[k] * kappa[j];
        Rk[k] = matmul(Bk, Z);
        for (std::size_t i = 0; i < ma; ++i)
            for (std::size_t j = 0; j < ma; ++j) R(i, j) += Rk[k](i, j);
    }
    Matrix<T> IR = I;
    for (std::size_t i = 0; i < ma; ++i)
        for (std::size_t j = 0; j < ma; ++j) IR(i, j) -= R(i, j);
    const std::vector<T> IRinve = mulvec(inverse(IR), e);
 
    MmapGk1LcfsprResult<T> r;
    r.classMeanQueueLength.assign(K, zero);
    r.meanSojournTime.assign(K, zero);
    r.meanQueueLength = zero;
    for (std::size_t k = 0; k < K; ++k) {
        const std::vector<T> v = vecmul(piArr, Rk[k]);
        T q = zero;
        for (std::size_t i = 0; i < ma; ++i) q += v[i] * IRinve[i];
        if (q < zero) q = zero;
        r.classMeanQueueLength[k] = q;
        r.meanSojournTime[k] = (lambdas[k] > zero) ? T(q / lambdas[k]) : zero;
        r.meanQueueLength += q;
    }
 
    if (opt.nMax >= 0) {
        const std::size_t nMax = static_cast<std::size_t>(opt.nMax);
        std::vector<T> y0(ma, zero);
        for (std::size_t i = 0; i < ma; ++i) y0[i] = (one - rho) * kappa[i];
        r.queueLengthDistr.assign(nMax + 1, zero);
        std::vector<T> v = y0;
        for (std::size_t n = 0; n <= nMax; ++n) {
            T s = zero;
            for (std::size_t i = 0; i < ma; ++i) s += v[i];
            r.queueLengthDistr[n] = (s > zero) ? s : zero;
            v = vecmul(v, R);
        }
        r.classQueueLengthDistr.assign(K, std::vector<T>(nMax + 1, zero));
        for (std::size_t k = 0; k < K; ++k) {
            Matrix<T> M = IR;
            for (std::size_t i = 0; i < ma; ++i)
                for (std::size_t j = 0; j < ma; ++j) M(i, j) += Rk[k](i, j);
            const Matrix<T> Ak = inverse(M);
            const std::vector<T> Ake = mulvec(Ak, e);
            const Matrix<T> ARk = matmul(Ak, Rk[k]);
            v = y0;
            for (std::size_t m = 0; m <= nMax; ++m) {
                T s = zero;
                for (std::size_t i = 0; i < ma; ++i) s += v[i] * Ake[i];
                r.classQueueLengthDistr[k][m] = (s > zero) ? s : zero;
                v = vecmul(v, ARk);
            }
        }
    }
 
    const std::vector<MmapGk1Service<T> > svcCopy = svc;
    r.respLST = [D0, Dk, svcCopy, piArr, lambdas, quad, tol](cdbl s) {
        return mmapgk1lcfsprdetail::sojourn_lst(D0, Dk, svcCopy, piArr, lambdas, s, quad, tol);
    };
 
    r.respPoints = opt.wPoints;
    if (!opt.wPoints.empty()) {
        r.respCDF.assign(K, std::vector<T>(opt.wPoints.size(), zero));
        for (std::size_t it = 0; it < opt.wPoints.size(); ++it) {
            const double t = num_traits<T>::to_double(opt.wPoints[it]);
            for (std::size_t k = 0; k < K; ++k) {
                const auto& lstAll = r.respLST;
                const std::function<cdbl(cdbl)> f = [&lstAll, k](cdbl s) { return lstAll(s)[k]; };
                r.respCDF[k][it] =
                    num_traits<T>::from_double(mmapgk1lcfsprdetail::euler_invert(f, t));
            }
        }
    }
 
    r.lambdas = lambdas;
    r.arrivalRate = zero;
    for (std::size_t k = 0; k < K; ++k) r.arrivalRate += lambdas[k];
    r.utilization = rho;
    r.Q = Q;
    r.kappa = kappa;
    r.Rk = Rk;
    r.R = R;
    r.analyzer = "LINE:MMAP[" + std::to_string(K) + "]/G[" + std::to_string(K) + "]/1-LCFSPR";
    return r;
}
 
/** qsys_mmapgk1_lcfspr with the reference's defaults (no distributions, no CDF). */
template <class T>
MmapGk1LcfsprResult<T> qsys_mmapgk1_lcfspr(const std::vector<Matrix<T> >& MMAP,
                                           const std::vector<MmapGk1Service<T> >& svc) {
    return qsys_mmapgk1_lcfspr(MMAP, svc, MmapGk1LcfsprOptions<T>());
}
 
/** qsys_mmapgk1_lcfspr on Distribution objects. */
template <class T>
MmapGk1LcfsprResult<T> qsys_mmapgk1_lcfspr(const std::vector<Matrix<T> >& MMAP,
                                           const std::vector<lang::Distrib<T> >& svc,
                                           const MmapGk1LcfsprOptions<T>& opt) {
    std::vector<MmapGk1Service<T> > laws;
    laws.reserve(svc.size());
    for (std::size_t k = 0; k < svc.size(); ++k) laws.push_back(MmapGk1Service<T>(svc[k]));
    return qsys_mmapgk1_lcfspr(MMAP, laws, opt);
}
 
/** qsys_mmapgk1_lcfspr on Distribution objects with the reference's defaults. */
template <class T>
MmapGk1LcfsprResult<T> qsys_mmapgk1_lcfspr(const std::vector<Matrix<T> >& MMAP,
                                           const std::vector<lang::Distrib<T> >& svc) {
    return qsys_mmapgk1_lcfspr(MMAP, svc, MmapGk1LcfsprOptions<T>());
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MMAPGK1_LCFSPR_H