docs: def rdl_decomposition_rev(T, norm='reversible', mu=None):
      r"""Decomposition into left and right eigenvectors for reversible
      transition matrices.
  
      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.
      mu : (M,) ndarray, optional
          Stationary distribution of T
  
      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, :]``
  
      Notes
      -----
      The eigenvalues and eigenvectors of the similar symmetric matrix
      `\sqrt(\mu_i / \mu_j) p_{ij}` will be used to compute the
      eigenvalues and eigenvectors of T.
  
      The stationay distribution will be computed if no precomputed stationary
      distribution is given.
  
      """
      if mu is None:
          mu = stationary_distribution(T)
      """ symmetrize T """
      smu = np.sqrt(mu)
      S = smu[:,None] * T / smu
      val, eigvec = eigh(S)
      """Sort eigenvalues and eigenvectors"""
      perm = np.argsort(np.abs(val))[::-1]
      val = val[perm]
      eigvec = eigvec[:, perm]
  
      """Diagonal matrix of eigenvalues"""
      D = np.diag(val)
  
      """Right and left eigenvectors"""
      R = eigvec / smu[:, np.newaxis]
      L = eigvec * smu[:, np.newaxis]
  
      """Ensure that R[:,0] is positive and unity"""
      tmp = R[0, 0]
      R[:, 0] = R[:, 0] / tmp
  
      """Ensure that L[:, 0] is probability vector"""
      L[:, 0] = L[:, 0] *  tmp
  
      if norm == 'reversible':
          return R, D, L.T
      elif norm == 'standard':
          """Standard l2-norm of right eigenvectors"""
          w = np.diag(np.dot(R.T, R))
          sw = np.sqrt(w)
          """Don't change normalization of eigenvectors for dominant eigenvalue"""
          sw[0] = 1.0
  
          R = R / sw[np.newaxis, :]
          L = L * sw[np.newaxis, :]
          return R, D, L.T
      else:
          raise ValueError("Keyword 'norm' has to be either 'standard' or 'reversible'")