docs: def rdl_decomposition_nrev(T, norm='standard'):
      r"""Decomposition into left and right eigenvectors.
  
      Parameters
      ----------
      T : (M, M) ndarray
          Transition matrix
      norm: {'standard', 'reversible'}
          standard: (L'R) = Id, L[:,0] is a probability distribution,
              the stationary distribution mu of T. Right eigenvectors
              R have a 2-norm of 1
          reversible: R and L are related via L=L[:,0]*R
  
      Returns
      -------
      R : (M, M) ndarray
          The normalized (with respect to L) right eigenvectors, such that the
          column R[:,i] is the right eigenvector corresponding to the eigenvalue
          w[i], dot(T,R[:,i])=w[i]*R[:,i]
      D : (M, M) ndarray
          A diagonal matrix containing the eigenvalues, each repeated
          according to its multiplicity
      L : (M, M) ndarray
          The normalized (with respect to `R`) left eigenvectors, such that the
          row ``L[i, :]`` is the left eigenvector corresponding to the eigenvalue
          ``w[i]``, ``dot(L[i, :], T)``=``w[i]*L[i, :]``
  
  
      """
      d = T.shape[0]
      w, R = eig(T)
  
      """Sort by decreasing magnitude of eigenvalue"""
      ind = np.argsort(np.abs(w))[::-1]
      w = w[ind]
      R = R[:, ind]
  
      """Diagonal matrix containing eigenvalues"""
      D = np.diag(w)
  
      # Standard norm: Euclidean norm is 1 for r and LR = I.
      if norm == 'standard':
          L = solve(np.transpose(R), np.eye(d))
  
          """l1- normalization of L[:, 0]"""
          R[:, 0] = R[:, 0] * np.sum(L[:, 0])
          L[:, 0] = L[:, 0] / np.sum(L[:, 0])
  
          return R, D, np.transpose(L)
  
      # Reversible norm:
      elif norm == 'reversible':
          b = np.zeros(d)
          b[0] = 1.0
  
          A = np.transpose(R)
          nu = solve(A, b)
          mu = nu / np.sum(nu)
  
          """Ensure that R[:,0] is positive"""
          R[:, 0] = R[:, 0] / np.sign(R[0, 0])
  
          """Use mu to connect L and R"""
          L = mu[:, np.newaxis] * R
  
          """Compute overlap"""
          s = np.diag(np.dot(np.transpose(L), R))
  
          """Renormalize left-and right eigenvectors to ensure L'R=Id"""
          R = R / np.sqrt(s[np.newaxis, :])
          L = L / np.sqrt(s[np.newaxis, :])
  
          return R, D, np.transpose(L)
  
      else:
          raise ValueError("Keyword 'norm' has to be either 'standard' or 'reversible'")