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.

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'")