ctmc_stochcomp.m

function [S,Q11,Q12,Q21,Q22,T]=ctmc_stochcomp(Q,I)
% [S,Q11,Q12,Q21,Q22,T] = CTMC_STOCHCOMP(Q,I)
% Copyright (c) 2012-2026, Imperial College London
% All rights reserved.
if nargin==1
    I=1:ceil(length(Q)/2);
end
isSelected = false(1, length(Q));
for idx = 1:length(I)
    isSelected(I(idx)) = true;
end
Ic = zeros(1, length(Q) - length(I));
icPos = 0;
for idx = 1:length(Q)
    if ~isSelected(idx)
        icPos = icPos + 1;
        Ic(icPos) = idx;
    end
end
Q11 = Q(I,I);
Q12 = Q(I,Ic);
Q21 = Q(Ic,I);
Q22 = Q(Ic,Ic);
%I = eye(size(Q22));
% Iterative path. The backslash below factorizes the whole complement block, so
% its fill-in is what limits the model above the dispatch threshold rather than
% the arithmetic. One ILUT factorization serves every column of Q21. The direct
% solve stays the default and remains the fallback when a column fails.
GMRES_MIN_STATES = 6000;
T = [];
if size(Q22,1) > GMRES_MIN_STATES
    [T,gflag] = ctmc_gmres_multi(-Q22, Q21);
    if gflag ~= 0
        T = [];
    end
end
if isempty(T)
    T = (-Q22) \ Q21;
end
T = Q12*T;
S = Q11+T;
end