JIT complied version of the main routine of gth_solve.

    Parameters
    ----------
    A : numpy.ndarray(float, ndim=2)
        Stochastic matrix or generator matrix. Must be of shape n x n.
        Data will be overwritten.

    out : numpy.ndarray(float, ndim=1)
        Output array in which to place the stationary distribution of A.

def _gth_solve_jit(A, out):
    """
    JIT complied version of the main routine of gth_solve.

    Parameters
    ----------
    A : numpy.ndarray(float, ndim=2)
        Stochastic matrix or generator matrix. Must be of shape n x n.
        Data will be overwritten.

    out : numpy.ndarray(float, ndim=1)
        Output array in which to place the stationary distribution of A.

    """
    n = A.shape[0]

    # === Reduction === #
    for k in range(n-1):
        scale = np.sum(A[k, k+1:n])
        if scale <= 0:
            # There is one (and only one) recurrent class contained in
            # {0, ..., k};
            # compute the solution associated with that recurrent class.
            n = k+1
            break
        for i in range(k+1, n):
            A[i, k] /= scale

            for j in range(k+1, n):
                A[i, j] += A[i, k] * A[k, j]

    # === Backward substitution === #
    out[n-1] = 1
    for k in range(n-2, -1, -1):
        for i in range(k+1, n):
            out[k] += out[i] * A[i, k]

    # === Normalization === #
    norm = np.sum(out)
    for k in range(n):
        out[k] /= norm