permutationmatrices.cpp
Go to the documentation of this file.
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 #define TEST_ENABLE_TEMPORARY_TRACKING
11 
12 #include "main.h"
13 
14 using namespace std;
15 template<typename MatrixType> void permutationmatrices(const MatrixType& m)
16 {
17  typedef typename MatrixType::Scalar Scalar;
18  enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime,
19  Options = MatrixType::Options };
20  typedef PermutationMatrix<Rows> LeftPermutationType;
21  typedef Transpositions<Rows> LeftTranspositionsType;
22  typedef Matrix<int, Rows, 1> LeftPermutationVectorType;
23  typedef Map<LeftPermutationType> MapLeftPerm;
24  typedef PermutationMatrix<Cols> RightPermutationType;
25  typedef Transpositions<Cols> RightTranspositionsType;
26  typedef Matrix<int, Cols, 1> RightPermutationVectorType;
27  typedef Map<RightPermutationType> MapRightPerm;
28 
29  Index rows = m.rows();
30  Index cols = m.cols();
31 
32  MatrixType m_original = MatrixType::Random(rows,cols);
33  LeftPermutationVectorType lv;
34  randomPermutationVector(lv, rows);
35  LeftPermutationType lp(lv);
36  RightPermutationVectorType rv;
37  randomPermutationVector(rv, cols);
38  RightPermutationType rp(rv);
39  LeftTranspositionsType lt(lv);
40  RightTranspositionsType rt(rv);
41  MatrixType m_permuted = MatrixType::Random(rows,cols);
42 
43  VERIFY_EVALUATION_COUNT(m_permuted = lp * m_original * rp, 1); // 1 temp for sub expression "lp * m_original"
44 
45  for (int i=0; i<rows; i++)
46  for (int j=0; j<cols; j++)
47  VERIFY_IS_APPROX(m_permuted(lv(i),j), m_original(i,rv(j)));
48 
49  Matrix<Scalar,Rows,Rows> lm(lp);
50  Matrix<Scalar,Cols,Cols> rm(rp);
51 
52  VERIFY_IS_APPROX(m_permuted, lm*m_original*rm);
53 
54  m_permuted = m_original;
55  VERIFY_EVALUATION_COUNT(m_permuted = lp * m_permuted * rp, 1);
56  VERIFY_IS_APPROX(m_permuted, lm*m_original*rm);
57 
58  LeftPermutationType lpi;
59  lpi = lp.inverse();
60  VERIFY_IS_APPROX(lpi*m_permuted,lp.inverse()*m_permuted);
61 
62  VERIFY_IS_APPROX(lp.inverse()*m_permuted*rp.inverse(), m_original);
63  VERIFY_IS_APPROX(lv.asPermutation().inverse()*m_permuted*rv.asPermutation().inverse(), m_original);
64  VERIFY_IS_APPROX(MapLeftPerm(lv.data(),lv.size()).inverse()*m_permuted*MapRightPerm(rv.data(),rv.size()).inverse(), m_original);
65 
66  VERIFY((lp*lp.inverse()).toDenseMatrix().isIdentity());
67  VERIFY((lv.asPermutation()*lv.asPermutation().inverse()).toDenseMatrix().isIdentity());
68  VERIFY((MapLeftPerm(lv.data(),lv.size())*MapLeftPerm(lv.data(),lv.size()).inverse()).toDenseMatrix().isIdentity());
69 
70  LeftPermutationVectorType lv2;
71  randomPermutationVector(lv2, rows);
72  LeftPermutationType lp2(lv2);
73  Matrix<Scalar,Rows,Rows> lm2(lp2);
74  VERIFY_IS_APPROX((lp*lp2).toDenseMatrix().template cast<Scalar>(), lm*lm2);
75  VERIFY_IS_APPROX((lv.asPermutation()*lv2.asPermutation()).toDenseMatrix().template cast<Scalar>(), lm*lm2);
76  VERIFY_IS_APPROX((MapLeftPerm(lv.data(),lv.size())*MapLeftPerm(lv2.data(),lv2.size())).toDenseMatrix().template cast<Scalar>(), lm*lm2);
77 
78  LeftPermutationType identityp;
79  identityp.setIdentity(rows);
80  VERIFY_IS_APPROX(m_original, identityp*m_original);
81 
82  // check inplace permutations
83  m_permuted = m_original;
84  VERIFY_EVALUATION_COUNT(m_permuted.noalias()= lp.inverse() * m_permuted, 1); // 1 temp to allocate the mask
85  VERIFY_IS_APPROX(m_permuted, lp.inverse()*m_original);
86 
87  m_permuted = m_original;
88  VERIFY_EVALUATION_COUNT(m_permuted.noalias() = m_permuted * rp.inverse(), 1); // 1 temp to allocate the mask
89  VERIFY_IS_APPROX(m_permuted, m_original*rp.inverse());
90 
91  m_permuted = m_original;
92  VERIFY_EVALUATION_COUNT(m_permuted.noalias() = lp * m_permuted, 1); // 1 temp to allocate the mask
93  VERIFY_IS_APPROX(m_permuted, lp*m_original);
94 
95  m_permuted = m_original;
96  VERIFY_EVALUATION_COUNT(m_permuted.noalias() = m_permuted * rp, 1); // 1 temp to allocate the mask
97  VERIFY_IS_APPROX(m_permuted, m_original*rp);
98 
99  if(rows>1 && cols>1)
100  {
101  lp2 = lp;
102  Index i = internal::random<Index>(0, rows-1);
103  Index j;
104  do j = internal::random<Index>(0, rows-1); while(j==i);
105  lp2.applyTranspositionOnTheLeft(i, j);
106  lm = lp;
107  lm.row(i).swap(lm.row(j));
108  VERIFY_IS_APPROX(lm, lp2.toDenseMatrix().template cast<Scalar>());
109 
110  RightPermutationType rp2 = rp;
111  i = internal::random<Index>(0, cols-1);
112  do j = internal::random<Index>(0, cols-1); while(j==i);
113  rp2.applyTranspositionOnTheRight(i, j);
114  rm = rp;
115  rm.col(i).swap(rm.col(j));
116  VERIFY_IS_APPROX(rm, rp2.toDenseMatrix().template cast<Scalar>());
117  }
118 
119  {
120  // simple compilation check
121  Matrix<Scalar, Cols, Cols> A = rp;
122  Matrix<Scalar, Cols, Cols> B = rp.transpose();
123  VERIFY_IS_APPROX(A, B.transpose());
124  }
125 
126  m_permuted = m_original;
127  lp = lt;
128  rp = rt;
129  VERIFY_EVALUATION_COUNT(m_permuted = lt * m_permuted * rt, 1);
130  VERIFY_IS_APPROX(m_permuted, lp*m_original*rp.transpose());
131 
132  VERIFY_IS_APPROX(lt.inverse()*m_permuted*rt.inverse(), m_original);
133 
134  // Check inplace transpositions
135  m_permuted = m_original;
136  VERIFY_IS_APPROX(m_permuted = lt * m_permuted, lp * m_original);
137  m_permuted = m_original;
138  VERIFY_IS_APPROX(m_permuted = lt.inverse() * m_permuted, lp.inverse() * m_original);
139  m_permuted = m_original;
140  VERIFY_IS_APPROX(m_permuted = m_permuted * rt, m_original * rt);
141  m_permuted = m_original;
142  VERIFY_IS_APPROX(m_permuted = m_permuted * rt.inverse(), m_original * rt.inverse());
143 }
144 
145 template<typename T>
146 void bug890()
147 {
148  typedef Matrix<T, Dynamic, Dynamic> MatrixType;
149  typedef Matrix<T, Dynamic, 1> VectorType;
150  typedef Stride<Dynamic,Dynamic> S;
151  typedef Map<MatrixType, Aligned, S> MapType;
152  typedef PermutationMatrix<Dynamic> Perm;
153 
154  VectorType v1(2), v2(2), op(4), rhs(2);
155  v1 << 666,667;
156  op << 1,0,0,1;
157  rhs << 42,42;
158 
159  Perm P(2);
160  P.indices() << 1, 0;
161 
162  MapType(v1.data(),2,1,S(1,1)) = P * MapType(rhs.data(),2,1,S(1,1));
163  VERIFY_IS_APPROX(v1, (P * rhs).eval());
164 
165  MapType(v1.data(),2,1,S(1,1)) = P.inverse() * MapType(rhs.data(),2,1,S(1,1));
166  VERIFY_IS_APPROX(v1, (P.inverse() * rhs).eval());
167 }
168 
169 EIGEN_DECLARE_TEST(permutationmatrices)
170 {
171  for(int i = 0; i < g_repeat; i++) {
172  CALL_SUBTEST_1( permutationmatrices(Matrix<float, 1, 1>()) );
173  CALL_SUBTEST_2( permutationmatrices(Matrix3f()) );
174  CALL_SUBTEST_3( permutationmatrices(Matrix<double,3,3,RowMajor>()) );
175  CALL_SUBTEST_4( permutationmatrices(Matrix4d()) );
176  CALL_SUBTEST_5( permutationmatrices(Matrix<double,40,60>()) );
177  CALL_SUBTEST_6( permutationmatrices(Matrix<double,Dynamic,Dynamic,RowMajor>(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
178  CALL_SUBTEST_7( permutationmatrices(MatrixXcf(internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
179  }
180  CALL_SUBTEST_5( bug890<double>() );
181 }
B
Matrix< SCALARB, Dynamic, Dynamic, opt_B > B
Definition: bench_gemm.cpp:49
Eigen::Stride
Holds strides information for Map.
Definition: Stride.h:48
MatrixType
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::swap
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void swap(DenseBase< OtherDerived > &other)
Definition: PlainObjectBase.h:953
Eigen::Transpositions
Represents a sequence of transpositions (row/column interchange)
Definition: Transpositions.h:155
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
A
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
CALL_SUBTEST_4
#define CALL_SUBTEST_4(FUNC)
Definition: split_test_helper.h:22
CALL_SUBTEST_3
#define CALL_SUBTEST_3(FUNC)
Definition: split_test_helper.h:16
CALL_SUBTEST_1
#define CALL_SUBTEST_1(FUNC)
Definition: split_test_helper.h:4
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
EIGEN_DECLARE_TEST
EIGEN_DECLARE_TEST(permutationmatrices)
Definition: permutationmatrices.cpp:169
CALL_SUBTEST_5
#define CALL_SUBTEST_5(FUNC)
Definition: split_test_helper.h:28
Eigen::g_repeat
static int g_repeat
Definition: main.h:169
m
Matrix3f m
Definition: AngleAxis_mimic_euler.cpp:1
CALL_SUBTEST_6
#define CALL_SUBTEST_6(FUNC)
Definition: split_test_helper.h:34
CALL_SUBTEST_2
#define CALL_SUBTEST_2(FUNC)
Definition: split_test_helper.h:10
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
VERIFY_IS_APPROX
#define VERIFY_IS_APPROX(a, b)
Definition: integer_types.cpp:15
main.h
std
Definition: BFloat16.h:88
Eigen::PermutationMatrix
Permutation matrix.
Definition: PermutationMatrix.h:297
EIGEN_TEST_MAX_SIZE
#define EIGEN_TEST_MAX_SIZE
Definition: boostmultiprec.cpp:16
VERIFY_EVALUATION_COUNT
#define VERIFY_EVALUATION_COUNT(XPR, N)
Definition: test/sparse_product.cpp:27
v2
Vector v2
Definition: testSerializationBase.cpp:39
P
static double P[]
Definition: ellpe.c:68
Eigen::randomPermutationVector
void randomPermutationVector(PermutationVectorType &v, Index size)
Definition: main.h:693
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
VectorType
Definition: FFTW.cpp:65
bug890
void bug890()
Definition: permutationmatrices.cpp:146
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
CALL_SUBTEST_7
#define CALL_SUBTEST_7(FUNC)
Definition: split_test_helper.h:40
permutationmatrices
void permutationmatrices(const MatrixType &m)
Definition: permutationmatrices.cpp:15
MapType
Map< MatrixType > MapType
Definition: Tutorial_Map_using.cpp:2
eval
internal::nested_eval< T, 1 >::type eval(const T &xpr)
Definition: sparse_permutations.cpp:38
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
S
DiscreteKey S(1, 2)
v1
Vector v1
Definition: testSerializationBase.cpp:38
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
VERIFY
#define VERIFY(a)
Definition: main.h:380
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Cols
static const int Cols
Definition: testCallRecord.cpp:32

gtsam
Author(s):
autogenerated on Tue Jun 25 2024 03:02:05