1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10 #include "main.h"
11
12 using namespace std;
permutationmatrices(const MatrixType & m)13 template<typename MatrixType> void permutationmatrices(const MatrixType& m)
14 {
15 typedef typename MatrixType::Index Index;
16 typedef typename MatrixType::Scalar Scalar;
17 enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime,
18 Options = MatrixType::Options };
19 typedef PermutationMatrix<Rows> LeftPermutationType;
20 typedef Matrix<int, Rows, 1> LeftPermutationVectorType;
21 typedef Map<LeftPermutationType> MapLeftPerm;
22 typedef PermutationMatrix<Cols> RightPermutationType;
23 typedef Matrix<int, Cols, 1> RightPermutationVectorType;
24 typedef Map<RightPermutationType> MapRightPerm;
25
26 Index rows = m.rows();
27 Index cols = m.cols();
28
29 MatrixType m_original = MatrixType::Random(rows,cols);
30 LeftPermutationVectorType lv;
31 randomPermutationVector(lv, rows);
32 LeftPermutationType lp(lv);
33 RightPermutationVectorType rv;
34 randomPermutationVector(rv, cols);
35 RightPermutationType rp(rv);
36 MatrixType m_permuted = lp * m_original * rp;
37
38 for (int i=0; i<rows; i++)
39 for (int j=0; j<cols; j++)
40 VERIFY_IS_APPROX(m_permuted(lv(i),j), m_original(i,rv(j)));
41
42 Matrix<Scalar,Rows,Rows> lm(lp);
43 Matrix<Scalar,Cols,Cols> rm(rp);
44
45 VERIFY_IS_APPROX(m_permuted, lm*m_original*rm);
46
47 VERIFY_IS_APPROX(lp.inverse()*m_permuted*rp.inverse(), m_original);
48 VERIFY_IS_APPROX(lv.asPermutation().inverse()*m_permuted*rv.asPermutation().inverse(), m_original);
49 VERIFY_IS_APPROX(MapLeftPerm(lv.data(),lv.size()).inverse()*m_permuted*MapRightPerm(rv.data(),rv.size()).inverse(), m_original);
50
51 VERIFY((lp*lp.inverse()).toDenseMatrix().isIdentity());
52 VERIFY((lv.asPermutation()*lv.asPermutation().inverse()).toDenseMatrix().isIdentity());
53 VERIFY((MapLeftPerm(lv.data(),lv.size())*MapLeftPerm(lv.data(),lv.size()).inverse()).toDenseMatrix().isIdentity());
54
55 LeftPermutationVectorType lv2;
56 randomPermutationVector(lv2, rows);
57 LeftPermutationType lp2(lv2);
58 Matrix<Scalar,Rows,Rows> lm2(lp2);
59 VERIFY_IS_APPROX((lp*lp2).toDenseMatrix().template cast<Scalar>(), lm*lm2);
60 VERIFY_IS_APPROX((lv.asPermutation()*lv2.asPermutation()).toDenseMatrix().template cast<Scalar>(), lm*lm2);
61 VERIFY_IS_APPROX((MapLeftPerm(lv.data(),lv.size())*MapLeftPerm(lv2.data(),lv2.size())).toDenseMatrix().template cast<Scalar>(), lm*lm2);
62
63 LeftPermutationType identityp;
64 identityp.setIdentity(rows);
65 VERIFY_IS_APPROX(m_original, identityp*m_original);
66
67 // check inplace permutations
68 m_permuted = m_original;
69 m_permuted = lp.inverse() * m_permuted;
70 VERIFY_IS_APPROX(m_permuted, lp.inverse()*m_original);
71
72 m_permuted = m_original;
73 m_permuted = m_permuted * rp.inverse();
74 VERIFY_IS_APPROX(m_permuted, m_original*rp.inverse());
75
76 m_permuted = m_original;
77 m_permuted = lp * m_permuted;
78 VERIFY_IS_APPROX(m_permuted, lp*m_original);
79
80 m_permuted = m_original;
81 m_permuted = m_permuted * rp;
82 VERIFY_IS_APPROX(m_permuted, m_original*rp);
83
84 if(rows>1 && cols>1)
85 {
86 lp2 = lp;
87 Index i = internal::random<Index>(0, rows-1);
88 Index j;
89 do j = internal::random<Index>(0, rows-1); while(j==i);
90 lp2.applyTranspositionOnTheLeft(i, j);
91 lm = lp;
92 lm.row(i).swap(lm.row(j));
93 VERIFY_IS_APPROX(lm, lp2.toDenseMatrix().template cast<Scalar>());
94
95 RightPermutationType rp2 = rp;
96 i = internal::random<Index>(0, cols-1);
97 do j = internal::random<Index>(0, cols-1); while(j==i);
98 rp2.applyTranspositionOnTheRight(i, j);
99 rm = rp;
100 rm.col(i).swap(rm.col(j));
101 VERIFY_IS_APPROX(rm, rp2.toDenseMatrix().template cast<Scalar>());
102 }
103 }
104
test_permutationmatrices()105 void test_permutationmatrices()
106 {
107 for(int i = 0; i < g_repeat; i++) {
108 CALL_SUBTEST_1( permutationmatrices(Matrix<float, 1, 1>()) );
109 CALL_SUBTEST_2( permutationmatrices(Matrix3f()) );
110 CALL_SUBTEST_3( permutationmatrices(Matrix<double,3,3,RowMajor>()) );
111 CALL_SUBTEST_4( permutationmatrices(Matrix4d()) );
112 CALL_SUBTEST_5( permutationmatrices(Matrix<double,40,60>()) );
113 CALL_SUBTEST_6( permutationmatrices(Matrix<double,Dynamic,Dynamic,RowMajor>(20, 30)) );
114 CALL_SUBTEST_7( permutationmatrices(MatrixXcf(15, 10)) );
115 }
116 }
117