permutationmatrices.cpp revision c981c48f5bc9aefeffc0bcb0cc3934c2fae179dd 
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
12using namespace std;
13template<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
106void 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