| 13 | |
| 14 | using namespace std; |
| 15 | template<typename MatrixType> void permutationmatrices(const MatrixType& m) |
| 16 | { |
| 17 | typedef typename MatrixType::Index Index; |
| 18 | typedef typename MatrixType::Scalar Scalar; |
| 19 | enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime, |
| 20 | Options = MatrixType::Options }; |
| 21 | typedef PermutationMatrix<Rows> LeftPermutationType; |
| 22 | typedef Matrix<int, Rows, 1> LeftPermutationVectorType; |
| 23 | typedef Map<LeftPermutationType> MapLeftPerm; |
| 24 | typedef PermutationMatrix<Cols> RightPermutationType; |
| 25 | typedef Matrix<int, Cols, 1> RightPermutationVectorType; |
| 26 | typedef Map<RightPermutationType> MapRightPerm; |
| 27 | |
| 28 | Index rows = m.rows(); |
| 29 | Index cols = m.cols(); |
| 30 | |
| 31 | MatrixType m_original = MatrixType::Random(rows,cols); |
| 32 | LeftPermutationVectorType lv; |
| 33 | randomPermutationVector(lv, rows); |
| 34 | LeftPermutationType lp(lv); |
| 35 | RightPermutationVectorType rv; |
| 36 | randomPermutationVector(rv, cols); |
| 37 | RightPermutationType rp(rv); |
| 38 | MatrixType m_permuted = MatrixType::Random(rows,cols); |
| 39 | |
| 40 | VERIFY_EVALUATION_COUNT(m_permuted = lp * m_original * rp, 1); // 1 temp for sub expression "lp * m_original" |
| 41 | |
| 42 | for (int i=0; i<rows; i++) |
| 43 | for (int j=0; j<cols; j++) |
| 44 | VERIFY_IS_APPROX(m_permuted(lv(i),j), m_original(i,rv(j))); |
| 45 | |
| 46 | Matrix<Scalar,Rows,Rows> lm(lp); |
| 47 | Matrix<Scalar,Cols,Cols> rm(rp); |
| 48 | |
| 49 | VERIFY_IS_APPROX(m_permuted, lm*m_original*rm); |
| 50 | |
| 51 | m_permuted = m_original; |
| 52 | VERIFY_EVALUATION_COUNT(m_permuted = lp * m_permuted * rp, 1); |
| 53 | VERIFY_IS_APPROX(m_permuted, lm*m_original*rm); |
| 54 | |
| 55 | VERIFY_IS_APPROX(lp.inverse()*m_permuted*rp.inverse(), m_original); |
| 56 | VERIFY_IS_APPROX(lv.asPermutation().inverse()*m_permuted*rv.asPermutation().inverse(), m_original); |
| 57 | VERIFY_IS_APPROX(MapLeftPerm(lv.data(),lv.size()).inverse()*m_permuted*MapRightPerm(rv.data(),rv.size()).inverse(), m_original); |
| 58 | |
| 59 | VERIFY((lp*lp.inverse()).toDenseMatrix().isIdentity()); |
| 60 | VERIFY((lv.asPermutation()*lv.asPermutation().inverse()).toDenseMatrix().isIdentity()); |
| 61 | VERIFY((MapLeftPerm(lv.data(),lv.size())*MapLeftPerm(lv.data(),lv.size()).inverse()).toDenseMatrix().isIdentity()); |
| 62 | |
| 63 | LeftPermutationVectorType lv2; |
| 64 | randomPermutationVector(lv2, rows); |
| 65 | LeftPermutationType lp2(lv2); |
| 66 | Matrix<Scalar,Rows,Rows> lm2(lp2); |
| 67 | VERIFY_IS_APPROX((lp*lp2).toDenseMatrix().template cast<Scalar>(), lm*lm2); |
| 68 | VERIFY_IS_APPROX((lv.asPermutation()*lv2.asPermutation()).toDenseMatrix().template cast<Scalar>(), lm*lm2); |
| 69 | VERIFY_IS_APPROX((MapLeftPerm(lv.data(),lv.size())*MapLeftPerm(lv2.data(),lv2.size())).toDenseMatrix().template cast<Scalar>(), lm*lm2); |
| 70 | |
| 71 | LeftPermutationType identityp; |
| 72 | identityp.setIdentity(rows); |