I  calculate two-site reduced density matrix   from  a MPS , just like 
//index linking i to i+1:
auto ir = commonIndex(psi.A(i),psi.A(i+1));
auto rho = psi.A(i)*dag(prime(psi.A(i),Site,ir));
for(int k = i+1; k < j; ++k)
    {
    rho *= psi.A(k);
    rho *= dag(prime(psi.A(k),Link));
    }
rho *= psi.A(j);
rho *= dag(prime(psi.A(j)));
for(int k = j+1; k <= psi.N(); ++k)
    {
    rho *= psi.A(k);
    rho *= dag(prime(psi.A(k),Link));
    }
How can I  obtain eigenvalue of  \rho