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 typedef typename MatrixType::RealScalar RealScalar;
18 enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime,
19 Options = MatrixType::Options };
20 typedef PermutationMatrix<Rows> LeftPermutationType;
21 typedef Matrix<int, Rows, 1> LeftPermutationVectorType;
22 typedef Map<LeftPermutationType> MapLeftPerm;
23 typedef PermutationMatrix<Cols> RightPermutationType;
24 typedef Matrix<int, Cols, 1> RightPermutationVectorType;
25 typedef Map<RightPermutationType> MapRightPerm;
26
27 Index rows = m.rows();
28 Index cols = m.cols();
29
30 MatrixType m_original = MatrixType::Random(rows,cols);
31 LeftPermutationVectorType lv;
32 randomPermutationVector(lv, rows);
33 LeftPermutationType lp(lv);
34 RightPermutationVectorType rv;
35 randomPermutationVector(rv, cols);
36 RightPermutationType rp(rv);
37 MatrixType m_permuted = lp * m_original * rp;
38
39 for (int i=0; i<rows; i++)
40 for (int j=0; j<cols; j++)
41 VERIFY_IS_APPROX(m_permuted(lv(i),j), m_original(i,rv(j)));
42
43 Matrix<Scalar,Rows,Rows> lm(lp);
44 Matrix<Scalar,Cols,Cols> rm(rp);
45
46 VERIFY_IS_APPROX(m_permuted, lm*m_original*rm);
47
48 VERIFY_IS_APPROX(lp.inverse()*m_permuted*rp.inverse(), m_original);
49 VERIFY_IS_APPROX(lv.asPermutation().inverse()*m_permuted*rv.asPermutation().inverse(), m_original);
50 VERIFY_IS_APPROX(MapLeftPerm(lv.data(),lv.size()).inverse()*m_permuted*MapRightPerm(rv.data(),rv.size()).inverse(), m_original);
51
52 VERIFY((lp*lp.inverse()).toDenseMatrix().isIdentity());
53 VERIFY((lv.asPermutation()*lv.asPermutation().inverse()).toDenseMatrix().isIdentity());
54 VERIFY((MapLeftPerm(lv.data(),lv.size())*MapLeftPerm(lv.data(),lv.size()).inverse()).toDenseMatrix().isIdentity());
55
56 LeftPermutationVectorType lv2;
57 randomPermutationVector(lv2, rows);
58 LeftPermutationType lp2(lv2);
59 Matrix<Scalar,Rows,Rows> lm2(lp2);
60 VERIFY_IS_APPROX((lp*lp2).toDenseMatrix().template cast<Scalar>(), lm*lm2);
61 VERIFY_IS_APPROX((lv.asPermutation()*lv2.asPermutation()).toDenseMatrix().template cast<Scalar>(), lm*lm2);
62 VERIFY_IS_APPROX((MapLeftPerm(lv.data(),lv.size())*MapLeftPerm(lv2.data(),lv2.size())).toDenseMatrix().template cast<Scalar>(), lm*lm2);
63
64 LeftPermutationType identityp;
65 identityp.setIdentity(rows);
66 VERIFY_IS_APPROX(m_original, identityp*m_original);
67
68 // check inplace permutations
69 m_permuted = m_original;
70 m_permuted = lp.inverse() * m_permuted;
71 VERIFY_IS_APPROX(m_permuted, lp.inverse()*m_original);
72
73 m_permuted = m_original;
74 m_permuted = m_permuted * rp.inverse();
75 VERIFY_IS_APPROX(m_permuted, m_original*rp.inverse());
76
77 m_permuted = m_original;
78 m_permuted = lp * m_permuted;
79 VERIFY_IS_APPROX(m_permuted, lp*m_original);
80
81 m_permuted = m_original;
82 m_permuted = m_permuted * rp;
83 VERIFY_IS_APPROX(m_permuted, m_original*rp);
84
85 if(rows>1 && cols>1)
86 {
87 lp2 = lp;
88 Index i = internal::random<Index>(0, rows-1);
89 Index j;
90 do j = internal::random<Index>(0, rows-1); while(j==i);
91 lp2.applyTranspositionOnTheLeft(i, j);
92 lm = lp;
93 lm.row(i).swap(lm.row(j));
94 VERIFY_IS_APPROX(lm, lp2.toDenseMatrix().template cast<Scalar>());
95
96 RightPermutationType rp2 = rp;
97 i = internal::random<Index>(0, cols-1);
98 do j = internal::random<Index>(0, cols-1); while(j==i);
99 rp2.applyTranspositionOnTheRight(i, j);
100 rm = rp;
101 rm.col(i).swap(rm.col(j));
102 VERIFY_IS_APPROX(rm, rp2.toDenseMatrix().template cast<Scalar>());
103 }
104 }
105
test_permutationmatrices()106 void test_permutationmatrices()
107 {
108 for(int i = 0; i < g_repeat; i++) {
109 CALL_SUBTEST_1( permutationmatrices(Matrix<float, 1, 1>()) );
110 CALL_SUBTEST_2( permutationmatrices(Matrix3f()) );
111 CALL_SUBTEST_3( permutationmatrices(Matrix<double,3,3,RowMajor>()) );
112 CALL_SUBTEST_4( permutationmatrices(Matrix4d()) );
113 CALL_SUBTEST_5( permutationmatrices(Matrix<double,40,60>()) );
114 CALL_SUBTEST_6( permutationmatrices(Matrix<double,Dynamic,Dynamic,RowMajor>(20, 30)) );
115 CALL_SUBTEST_7( permutationmatrices(MatrixXcf(15, 10)) );
116 }
117 }
118