def rdl_decomposition(T, k=None, reversible=False, norm='standard', mu=None):
    r"""Compute the decomposition into left and right eigenvectors.

    Parameters
    ----------
    T : (M, M) ndarray
        Transition matrix
    k : int (optional)
        Number of eigenvector/eigenvalue pairs
    norm: {'standard', 'reversible', 'auto'}
        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.
        auto: will be reversible if T is reversible, otherwise standard
    reversible : bool, optional
        Indicate that transition matrix is reversible
    mu : (d,) 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
    -----
    If reversible=True the 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 precomputed stationary distribution will only be used if
    reversible=True.

    """
    # auto-set norm
    if norm == 'auto':
        if is_reversible(T):
            norm = 'reversible'
        else:
            norm = 'standard'

    if reversible:
        R, D, L = rdl_decomposition_rev(T, norm=norm, mu=mu)
    else:
        R, D, L = rdl_decomposition_nrev(T, norm=norm)

    if k is None:
        return R, D, L
    else:
        return R[:, 0:k], D[0:k, 0:k], L[0:k, :]