docs: 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