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   // chekc that linear acccessors works on blocks
 135   m1 = m1_copy;
 136   if((MatrixType::Flags&RowMajorBit)==0)
 137     VERIFY_IS_EQUAL(m1.leftCols(c1).coeff(r1+c1*rows), m1(r1,c1));
 138   else
 139     VERIFY_IS_EQUAL(m1.topRows(r1).coeff(c1+r1*cols), m1(r1,c1));
 140
 141
 142   // now test some block-inside-of-block.
 143
 144   // expressions with direct access
 145   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)) );
 146   VERIFY_IS_EQUAL( (m1.block(r1,c1,r2-r1+1,c2-c1+1).row(0)) , (m1.row(r1).segment(c1,c2-c1+1)) );
 147   VERIFY_IS_EQUAL( (m1.block(r1,c1,r2-r1+1,c2-c1+1).col(0)) , (m1.col(c1).segment(r1,r2-r1+1)) );
 148   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() );
 149   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() );
 150
 151   // expressions without direct access
 152   VERIFY_IS_APPROX( ((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)) );
 153   VERIFY_IS_APPROX( ((m1+m2).block(r1,c1,r2-r1+1,c2-c1+1).row(0)) , ((m1+m2).row(r1).segment(c1,c2-c1+1)) );
 154   VERIFY_IS_APPROX( ((m1+m2).block(r1,c1,r2-r1+1,c2-c1+1).col(0)) , ((m1+m2).col(c1).segment(r1,r2-r1+1)) );
 155   VERIFY_IS_APPROX( ((m1+m2).block(r1,c1,r2-r1+1,c2-c1+1).transpose().col(0)) , ((m1+m2).row(r1).segment(c1,c2-c1+1)).transpose() );
 156   VERIFY_IS_APPROX( ((m1+m2).transpose().block(c1,r1,c2-c1+1,r2-r1+1).col(0)) , ((m1+m2).row(r1).segment(c1,c2-c1+1)).transpose() );
 157
 158   // evaluation into plain matrices from expressions with direct access (stress MapBase)
 159   DynamicMatrixType dm;
 160   DynamicVectorType dv;
 161   dm.setZero();
 162   dm = m1.block(r1,c1,rows-r1,cols-c1).block(r2-r1,c2-c1,rows-r2,cols-c2);
 163   VERIFY_IS_EQUAL(dm, (m1.block(r2,c2,rows-r2,cols-c2)));
 164   dm.setZero();
 165   dv.setZero();
 166   dm = m1.block(r1,c1,r2-r1+1,c2-c1+1).row(0).transpose();
 167   dv = m1.row(r1).segment(c1,c2-c1+1);
 168   VERIFY_IS_EQUAL(dv, dm);
 169   dm.setZero();
 170   dv.setZero();
 171   dm = m1.col(c1).segment(r1,r2-r1+1);
 172   dv = m1.block(r1,c1,r2-r1+1,c2-c1+1).col(0);
 173   VERIFY_IS_EQUAL(dv, dm);
 174   dm.setZero();
 175   dv.setZero();
 176   dm = m1.block(r1,c1,r2-r1+1,c2-c1+1).transpose().col(0);
 177   dv = m1.row(r1).segment(c1,c2-c1+1);
 178   VERIFY_IS_EQUAL(dv, dm);
 179   dm.setZero();
 180   dv.setZero();
 181   dm = m1.row(r1).segment(c1,c2-c1+1).transpose();
 182   dv = m1.transpose().block(c1,r1,c2-c1+1,r2-r1+1).col(0);
 183   VERIFY_IS_EQUAL(dv, dm);
 184
 185   VERIFY_IS_EQUAL( (m1.template block<Dynamic,1>(1,0,0,1)), m1.block(1,0,0,1));
 186   VERIFY_IS_EQUAL( (m1.template block<1,Dynamic>(0,1,1,0)), m1.block(0,1,1,0));
 187   VERIFY_IS_EQUAL( ((m1*1).template block<Dynamic,1>(1,0,0,1)), m1.block(1,0,0,1));
 188   VERIFY_IS_EQUAL( ((m1*1).template block<1,Dynamic>(0,1,1,0)), m1.block(0,1,1,0));
 189
 190   if (rows>=2 && cols>=2)
 191   {
 192     VERIFY_RAISES_ASSERT( m1 += m1.col(0) );
 193     VERIFY_RAISES_ASSERT( m1 -= m1.col(0) );
 194     VERIFY_RAISES_ASSERT( m1.array() *= m1.col(0).array() );
 195     VERIFY_RAISES_ASSERT( m1.array() /= m1.col(0).array() );
 196   }
 197 }
 198
 199
 200 template<typename MatrixType>
 201 void compare_using_data_and_stride(const MatrixType& m)
 202 {
 203   typedef typename MatrixType::Index Index;
 204   Index rows = m.rows();
 205   Index cols = m.cols();
 206   Index size = m.size();
 207   Index innerStride = m.innerStride();
 208   Index outerStride = m.outerStride();
 209   Index rowStride = m.rowStride();
 210   Index colStride = m.colStride();
 211   const typename MatrixType::Scalar* data = m.data();
 212
 213   for(int j=0;j<cols;++j)
 214     for(int i=0;i<rows;++i)
 215       VERIFY(m.coeff(i,j) == data[i*rowStride + j*colStride]);
 216
 217   if(!MatrixType::IsVectorAtCompileTime)
 218   {
 219     for(int j=0;j<cols;++j)
 220       for(int i=0;i<rows;++i)
 221         VERIFY(m.coeff(i,j) == data[(MatrixType::Flags&RowMajorBit)
 222                                      ? i*outerStride + j*innerStride
 223                                      : j*outerStride + i*innerStride]);
 224   }
 225
 226   if(MatrixType::IsVectorAtCompileTime)
 227   {
 228     VERIFY(innerStride == int((&m.coeff(1))-(&m.coeff(0))));
 229     for (int i=0;i<size;++i)
 230       VERIFY(m.coeff(i) == data[i*innerStride]);
 231   }
 232 }
 233
 234 template<typename MatrixType>
 235 void data_and_stride(const MatrixType& m)
 236 {
 237   typedef typename MatrixType::Index Index;
 238   Index rows = m.rows();
 239   Index cols = m.cols();
 240
 241   Index r1 = internal::random<Index>(0,rows-1);
 242   Index r2 = internal::random<Index>(r1,rows-1);
 243   Index c1 = internal::random<Index>(0,cols-1);
 244   Index c2 = internal::random<Index>(c1,cols-1);
 245
 246   MatrixType m1 = MatrixType::Random(rows, cols);
 247   compare_using_data_and_stride(m1.block(r1, c1, r2-r1+1, c2-c1+1));
 248   compare_using_data_and_stride(m1.transpose().block(c1, r1, c2-c1+1, r2-r1+1));
 249   compare_using_data_and_stride(m1.row(r1));
 250   compare_using_data_and_stride(m1.col(c1));
 251   compare_using_data_and_stride(m1.row(r1).transpose());
 252   compare_using_data_and_stride(m1.col(c1).transpose());
 253 }
 254
 255 void test_block()
 256 {
 257   for(int i = 0; i < g_repeat; i++) {
 258     CALL_SUBTEST_1( block(Matrix<float, 1, 1>()) );
 259     CALL_SUBTEST_2( block(Matrix4d()) );
 260     CALL_SUBTEST_3( block(MatrixXcf(3, 3)) );
 261     CALL_SUBTEST_4( block(MatrixXi(8, 12)) );
 262     CALL_SUBTEST_5( block(MatrixXcd(20, 20)) );
 263     CALL_SUBTEST_6( block(MatrixXf(20, 20)) );
 264
 265     CALL_SUBTEST_8( block(Matrix<float,Dynamic,4>(3, 4)) );
 266
 267 #ifndef EIGEN_DEFAULT_TO_ROW_MAJOR
 268     CALL_SUBTEST_6( data_and_stride(MatrixXf(internal::random(5,50), internal::random(5,50))) );
 269     CALL_SUBTEST_7( data_and_stride(Matrix<int,Dynamic,Dynamic,RowMajor>(internal::random(5,50), internal::random(5,50))) );
 270 #endif
 271   }
 272 }