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

   1 // This file is part of Eigen, a lightweight C++ template library
   2 // for linear algebra.
   3 //
   4 // Copyright (C) 2006-2010 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 EIGEN_NO_STATIC_ASSERT // otherwise we fail at compile time on unused paths
  11 #include "main.h"
  12
  13 template<typename MatrixType, typename Index, typename Scalar>
  14 typename Eigen::internal::enable_if<!NumTraits<typename MatrixType::Scalar>::IsComplex,typename MatrixType::Scalar>::type
  15 block_real_only(const MatrixType &m1, Index r1, Index r2, Index c1, Index c2, const Scalar& s1) {
  16   // check cwise-Functions:
  17   VERIFY_IS_APPROX(m1.row(r1).cwiseMax(s1), m1.cwiseMax(s1).row(r1));
  18   VERIFY_IS_APPROX(m1.col(c1).cwiseMin(s1), m1.cwiseMin(s1).col(c1));
  19
  20   VERIFY_IS_APPROX(m1.block(r1,c1,r2-r1+1,c2-c1+1).cwiseMin(s1), m1.cwiseMin(s1).block(r1,c1,r2-r1+1,c2-c1+1));
  21   VERIFY_IS_APPROX(m1.block(r1,c1,r2-r1+1,c2-c1+1).cwiseMax(s1), m1.cwiseMax(s1).block(r1,c1,r2-r1+1,c2-c1+1));
  22
  23   return Scalar(0);
  24 }
  25
  26 template<typename MatrixType, typename Index, typename Scalar>
  27 typename Eigen::internal::enable_if<NumTraits<typename MatrixType::Scalar>::IsComplex,typename MatrixType::Scalar>::type
  28 block_real_only(const MatrixType &, Index, Index, Index, Index, const Scalar&) {
  29   return Scalar(0);
  30 }
  31
  32
  33 template<typename MatrixType> void block(const MatrixType& m)
  34 {
  35   typedef typename MatrixType::Index Index;
  36   typedef typename MatrixType::Scalar Scalar;
  37   typedef typename MatrixType::RealScalar RealScalar;
  38   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
  39   typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType;
  40   typedef Matrix<Scalar, Dynamic, Dynamic> DynamicMatrixType;
  41   typedef Matrix<Scalar, Dynamic, 1> DynamicVectorType;
  42
  43   Index rows = m.rows();
  44   Index cols = m.cols();
  45
  46   MatrixType m1 = MatrixType::Random(rows, cols),
  47              m1_copy = m1,
  48              m2 = MatrixType::Random(rows, cols),
  49              m3(rows, cols),
  50              ones = MatrixType::Ones(rows, cols);
  51   VectorType v1 = VectorType::Random(rows);
  52
  53   Scalar s1 = internal::random<Scalar>();
  54
  55   Index r1 = internal::random<Index>(0,rows-1);
  56   Index r2 = internal::random<Index>(r1,rows-1);
  57   Index c1 = internal::random<Index>(0,cols-1);
  58   Index c2 = internal::random<Index>(c1,cols-1);
  59
  60   block_real_only(m1, r1, r2, c1, c1, s1);
  61
  62   //check row() and col()
  63   VERIFY_IS_EQUAL(m1.col(c1).transpose(), m1.transpose().row(c1));
  64   //check operator(), both constant and non-constant, on row() and col()
  65   m1 = m1_copy;
  66   m1.row(r1) += s1 * m1_copy.row(r2);
  67   VERIFY_IS_APPROX(m1.row(r1), m1_copy.row(r1) + s1 * m1_copy.row(r2));
  68   // check nested block xpr on lhs
  69   m1.row(r1).row(0) += s1 * m1_copy.row(r2);
  70   VERIFY_IS_APPROX(m1.row(r1), m1_copy.row(r1) + Scalar(2) * s1 * m1_copy.row(r2));
  71   m1 = m1_copy;
  72   m1.col(c1) += s1 * m1_copy.col(c2);
  73   VERIFY_IS_APPROX(m1.col(c1), m1_copy.col(c1) + s1 * m1_copy.col(c2));
  74   m1.col(c1).col(0) += s1 * m1_copy.col(c2);
  75   VERIFY_IS_APPROX(m1.col(c1), m1_copy.col(c1) + Scalar(2) * s1 * m1_copy.col(c2));
  76
  77
  78   //check block()
  79   Matrix<Scalar,Dynamic,Dynamic> b1(1,1); b1(0,0) = m1(r1,c1);
  80
  81   RowVectorType br1(m1.block(r1,0,1,cols));
  82   VectorType bc1(m1.block(0,c1,rows,1));
  83   VERIFY_IS_EQUAL(b1, m1.block(r1,c1,1,1));
  84   VERIFY_IS_EQUAL(m1.row(r1), br1);
  85   VERIFY_IS_EQUAL(m1.col(c1), bc1);
  86   //check operator(), both constant and non-constant, on block()
  87   m1.block(r1,c1,r2-r1+1,c2-c1+1) = s1 * m2.block(0, 0, r2-r1+1,c2-c1+1);
  88   m1.block(r1,c1,r2-r1+1,c2-c1+1)(r2-r1,c2-c1) = m2.block(0, 0, r2-r1+1,c2-c1+1)(0,0);
  89
  90   enum {
  91     BlockRows = 2,
  92     BlockCols = 5
  93   };
  94   if (rows>=5 && cols>=8)
  95   {
  96     // test fixed block() as lvalue
  97     m1.template block<BlockRows,BlockCols>(1,1) *= s1;
  98     // test operator() on fixed block() both as constant and non-constant
  99     m1.template block<BlockRows,BlockCols>(1,1)(0, 3) = m1.template block<2,5>(1,1)(1,2);
 100     // check that fixed block() and block() agree
 101     Matrix<Scalar,Dynamic,Dynamic> b = m1.template block<BlockRows,BlockCols>(3,3);
 102     VERIFY_IS_EQUAL(b, m1.block(3,3,BlockRows,BlockCols));
 103
 104     // same tests with mixed fixed/dynamic size
 105     m1.template block<BlockRows,Dynamic>(1,1,BlockRows,BlockCols) *= s1;
 106     m1.template block<BlockRows,Dynamic>(1,1,BlockRows,BlockCols)(0,3) = m1.template block<2,5>(1,1)(1,2);
 107     Matrix<Scalar,Dynamic,Dynamic> b2 = m1.template block<Dynamic,BlockCols>(3,3,2,5);
 108     VERIFY_IS_EQUAL(b2, m1.block(3,3,BlockRows,BlockCols));
 109   }
 110
 111   if (rows>2)
 112   {
 113     // test sub vectors
 114     VERIFY_IS_EQUAL(v1.template head<2>(), v1.block(0,0,2,1));
 115     VERIFY_IS_EQUAL(v1.template head<2>(), v1.head(2));
 116     VERIFY_IS_EQUAL(v1.template head<2>(), v1.segment(0,2));
 117     VERIFY_IS_EQUAL(v1.template head<2>(), v1.template segment<2>(0));
 118     Index i = rows-2;
 119     VERIFY_IS_EQUAL(v1.template tail<2>(), v1.block(i,0,2,1));
 120     VERIFY_IS_EQUAL(v1.template tail<2>(), v1.tail(2));
 121     VERIFY_IS_EQUAL(v1.template tail<2>(), v1.segment(i,2));
 122     VERIFY_IS_EQUAL(v1.template tail<2>(), v1.template segment<2>(i));
 123     i = internal::random<Index>(0,rows-2);
 124     VERIFY_IS_EQUAL(v1.segment(i,2), v1.template segment<2>(i));
 125   }
 126
 127   // stress some basic stuffs with block matrices
 128   VERIFY(numext::real(ones.col(c1).sum()) == RealScalar(rows));
 129   VERIFY(numext::real(ones.row(r1).sum()) == RealScalar(cols));
 130
 131   VERIFY(numext::real(ones.col(c1).dot(ones.col(c2))) == RealScalar(rows));
 132   VERIFY(numext::real(ones.row(r1).dot(ones.row(r2))) == RealScalar(cols));
 133
 134   // now test some block-inside-of-block.
 135
 136   // expressions with direct access
 137   VERIFY_IS_EQUAL( (m1.block(r1,c1,rows-r1,cols-c1).block(r2-r1,c2-c1,rows-r2,cols-c2)) , (m1.block(r2,c2,rows-r2,cols-c2)) );
 138   VERIFY_IS_EQUAL( (m1.block(r1,c1,r2-r1+1,c2-c1+1).row(0)) , (m1.row(r1).segment(c1,c2-c1+1)) );
 139   VERIFY_IS_EQUAL( (m1.block(r1,c1,r2-r1+1,c2-c1+1).col(0)) , (m1.col(c1).segment(r1,r2-r1+1)) );
 140   VERIFY_IS_EQUAL( (m1.block(r1,c1,r2-r1+1,c2-c1+1).transpose().col(0)) , (m1.row(r1).segment(c1,c2-c1+1)).transpose() );
 141   VERIFY_IS_EQUAL( (m1.transpose().block(c1,r1,c2-c1+1,r2-r1+1).col(0)) , (m1.row(r1).segment(c1,c2-c1+1)).transpose() );
 142
 143   // expressions without direct access
 144   VERIFY_IS_EQUAL( ((m1+m2).block(r1,c1,rows-r1,cols-c1).block(r2-r1,c2-c1,rows-r2,cols-c2)) , ((m1+m2).block(r2,c2,rows-r2,cols-c2)) );
 145   VERIFY_IS_EQUAL( ((m1+m2).block(r1,c1,r2-r1+1,c2-c1+1).row(0)) , ((m1+m2).row(r1).segment(c1,c2-c1+1)) );
 146   VERIFY_IS_EQUAL( ((m1+m2).block(r1,c1,r2-r1+1,c2-c1+1).col(0)) , ((m1+m2).col(c1).segment(r1,r2-r1+1)) );
 147   VERIFY_IS_EQUAL( ((m1+m2).block(r1,c1,r2-r1+1,c2-c1+1).transpose().col(0)) , ((m1+m2).row(r1).segment(c1,c2-c1+1)).transpose() );
 148   VERIFY_IS_EQUAL( ((m1+m2).transpose().block(c1,r1,c2-c1+1,r2-r1+1).col(0)) , ((m1+m2).row(r1).segment(c1,c2-c1+1)).transpose() );
 149
 150   // evaluation into plain matrices from expressions with direct access (stress MapBase)
 151   DynamicMatrixType dm;
 152   DynamicVectorType dv;
 153   dm.setZero();
 154   dm = m1.block(r1,c1,rows-r1,cols-c1).block(r2-r1,c2-c1,rows-r2,cols-c2);
 155   VERIFY_IS_EQUAL(dm, (m1.block(r2,c2,rows-r2,cols-c2)));
 156   dm.setZero();
 157   dv.setZero();
 158   dm = m1.block(r1,c1,r2-r1+1,c2-c1+1).row(0).transpose();
 159   dv = m1.row(r1).segment(c1,c2-c1+1);
 160   VERIFY_IS_EQUAL(dv, dm);
 161   dm.setZero();
 162   dv.setZero();
 163   dm = m1.col(c1).segment(r1,r2-r1+1);
 164   dv = m1.block(r1,c1,r2-r1+1,c2-c1+1).col(0);
 165   VERIFY_IS_EQUAL(dv, dm);
 166   dm.setZero();
 167   dv.setZero();
 168   dm = m1.block(r1,c1,r2-r1+1,c2-c1+1).transpose().col(0);
 169   dv = m1.row(r1).segment(c1,c2-c1+1);
 170   VERIFY_IS_EQUAL(dv, dm);
 171   dm.setZero();
 172   dv.setZero();
 173   dm = m1.row(r1).segment(c1,c2-c1+1).transpose();
 174   dv = m1.transpose().block(c1,r1,c2-c1+1,r2-r1+1).col(0);
 175   VERIFY_IS_EQUAL(dv, dm);
 176 }
 177
 178
 179 template<typename MatrixType>
 180 void compare_using_data_and_stride(const MatrixType& m)
 181 {
 182   typedef typename MatrixType::Index Index;
 183   Index rows = m.rows();
 184   Index cols = m.cols();
 185   Index size = m.size();
 186   Index innerStride = m.innerStride();
 187   Index outerStride = m.outerStride();
 188   Index rowStride = m.rowStride();
 189   Index colStride = m.colStride();
 190   const typename MatrixType::Scalar* data = m.data();
 191
 192   for(int j=0;j<cols;++j)
 193     for(int i=0;i<rows;++i)
 194       VERIFY(m.coeff(i,j) == data[i*rowStride + j*colStride]);
 195
 196   if(!MatrixType::IsVectorAtCompileTime)
 197   {
 198     for(int j=0;j<cols;++j)
 199       for(int i=0;i<rows;++i)
 200         VERIFY(m.coeff(i,j) == data[(MatrixType::Flags&RowMajorBit)
 201                                      ? i*outerStride + j*innerStride
 202                                      : j*outerStride + i*innerStride]);
 203   }
 204
 205   if(MatrixType::IsVectorAtCompileTime)
 206   {
 207     VERIFY(innerStride == int((&m.coeff(1))-(&m.coeff(0))));
 208     for (int i=0;i<size;++i)
 209       VERIFY(m.coeff(i) == data[i*innerStride]);
 210   }
 211 }
 212
 213 template<typename MatrixType>
 214 void data_and_stride(const MatrixType& m)
 215 {
 216   typedef typename MatrixType::Index Index;
 217   Index rows = m.rows();
 218   Index cols = m.cols();
 219
 220   Index r1 = internal::random<Index>(0,rows-1);
 221   Index r2 = internal::random<Index>(r1,rows-1);
 222   Index c1 = internal::random<Index>(0,cols-1);
 223   Index c2 = internal::random<Index>(c1,cols-1);
 224
 225   MatrixType m1 = MatrixType::Random(rows, cols);
 226   compare_using_data_and_stride(m1.block(r1, c1, r2-r1+1, c2-c1+1));
 227   compare_using_data_and_stride(m1.transpose().block(c1, r1, c2-c1+1, r2-r1+1));
 228   compare_using_data_and_stride(m1.row(r1));
 229   compare_using_data_and_stride(m1.col(c1));
 230   compare_using_data_and_stride(m1.row(r1).transpose());
 231   compare_using_data_and_stride(m1.col(c1).transpose());
 232 }
 233
 234 void test_block()
 235 {
 236   for(int i = 0; i < g_repeat; i++) {
 237     CALL_SUBTEST_1( block(Matrix<float, 1, 1>()) );
 238     CALL_SUBTEST_2( block(Matrix4d()) );
 239     CALL_SUBTEST_3( block(MatrixXcf(3, 3)) );
 240     CALL_SUBTEST_4( block(MatrixXi(8, 12)) );
 241     CALL_SUBTEST_5( block(MatrixXcd(20, 20)) );
 242     CALL_SUBTEST_6( block(MatrixXf(20, 20)) );
 243
 244     CALL_SUBTEST_8( block(Matrix<float,Dynamic,4>(3, 4)) );
 245
 246 #ifndef EIGEN_DEFAULT_TO_ROW_MAJOR
 247     CALL_SUBTEST_6( data_and_stride(MatrixXf(internal::random(5,50), internal::random(5,50))) );
 248     CALL_SUBTEST_7( data_and_stride(Matrix<int,Dynamic,Dynamic,RowMajor>(internal::random(5,50), internal::random(5,50))) );
 249 #endif
 250   }
 251 }