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

    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, :]``

def rdl_decomposition_nrev(T, k, norm='standard', ncv=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

    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, :]``

    """
    # Standard norm: Euclidean norm is 1 for r and LR = I.
    if norm == 'standard':
        v, R = scipy.sparse.linalg.eigs(T, k=k, which='LM', ncv=ncv)
        r, L = scipy.sparse.linalg.eigs(T.transpose(), k=k, which='LM', ncv=ncv)

        """Sort right eigenvectors"""
        ind = np.argsort(np.abs(v))[::-1]
        v = v[ind]
        R = R[:, ind]

        """Sort left eigenvectors"""
        ind = np.argsort(np.abs(r))[::-1]
        r = r[ind]
        L = L[:, ind]

        """l1-normalization of L[:, 0]"""
        L[:, 0] = L[:, 0] / np.sum(L[:, 0])

        """Standard normalization L'R=Id"""
        ov = np.diag(np.dot(np.transpose(L), R))
        R = R / ov[np.newaxis, :]

        """Diagonal matrix with eigenvalues"""
        D = np.diag(v)

        return R, D, np.transpose(L)

    # Reversible norm:
    elif norm == 'reversible':
        v, R = scipy.sparse.linalg.eigs(T, k=k, which='LM', ncv=ncv)
        mu = stationary_distribution(T)

        """Sort right eigenvectors"""
        ind = np.argsort(np.abs(v))[::-1]
        v = v[ind]
        R = R[:, ind]

        """Ensure that R[:,0] is positive"""
        R[:, 0] = R[:, 0] / np.sign(R[0, 0])

        """Diagonal matrix with eigenvalues"""
        D = np.diag(v)

        """Compute left eigenvectors from right ones"""
        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'")