docs: def rdl_decomposition_rev(T, k, norm='reversible', ncv=None, mu=None):
      r"""Compute the decomposition into left and right eigenvectors.
  
      Parameters
      ----------
      T : sparse matrix
          Transition matrix
      k : int
          Number of eigenvector/eigenvalue pairs
      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.
      ncv : int (optional)
          The number of Lanczos vectors generated, `ncv` must be greater than k;
          it is recommended that ncv > 2*k
      mu : (M,) ndarray, optional
          Stationary distribution of T
  
      Returns
      -------
      R : (M, M) ndarray
          The normalized ("unit length") 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, :]``
  
      """
      if mu is None:
          mu = stationary_distribution(T)
      """ symmetrize T """
      smu = np.sqrt(mu)
      Dpi = diags(smu, 0)
      Dinv = diags(1.0/smu, 0)
      S = (Dpi.dot(T)).dot(Dinv)
      """Compute eigenvalues, eigenvecs using a solver for
      symmetric/hermititan eigenproblems"""
      val, eigvec = scipy.sparse.linalg.eigsh(S, k=k, ncv=ncv, which='LM',
                                              return_eigenvectors=True)
      """Sort eigenvalues and eigenvectors"""
      ind = np.argsort(np.abs(val))[::-1]
      val = val[ind]
      eigvec = eigvec[:, ind]
  
      """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'")