ground/gcs/src/libs/eigen/test/permutationmatrices.cpp

   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;
  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
 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 }