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