cart-elc

Source code for CART-ELC
git clone git://git.laack.co/cart-elc.git
Log | Files | Refs | README | LICENSE

ref.cpp (14363B)


      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 20013 Gael Guennebaud <gael.guennebaud@inria.fr>
      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 // This unit test cannot be easily written to work with EIGEN_DEFAULT_TO_ROW_MAJOR
     11 #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR
     12 #undef EIGEN_DEFAULT_TO_ROW_MAJOR
     13 #endif
     14 
     15 #define TEST_ENABLE_TEMPORARY_TRACKING
     16 #define TEST_CHECK_STATIC_ASSERTIONS
     17 #include "main.h"
     18 
     19 // test Ref.h
     20 
     21 // Deal with i387 extended precision
     22 #if EIGEN_ARCH_i386 && !(EIGEN_ARCH_x86_64)
     23 
     24 #if EIGEN_COMP_GNUC_STRICT && EIGEN_GNUC_AT_LEAST(4,4)
     25 #pragma GCC optimize ("-ffloat-store")
     26 #else
     27 #undef VERIFY_IS_EQUAL
     28 #define VERIFY_IS_EQUAL(X,Y) VERIFY_IS_APPROX(X,Y)
     29 #endif
     30 
     31 #endif
     32 
     33 template<typename MatrixType> void ref_matrix(const MatrixType& m)
     34 {
     35   typedef typename MatrixType::Scalar Scalar;
     36   typedef typename MatrixType::RealScalar RealScalar;
     37   typedef Matrix<Scalar,Dynamic,Dynamic,MatrixType::Options> DynMatrixType;
     38   typedef Matrix<RealScalar,Dynamic,Dynamic,MatrixType::Options> RealDynMatrixType;
     39   
     40   typedef Ref<MatrixType> RefMat;
     41   typedef Ref<DynMatrixType> RefDynMat;
     42   typedef Ref<const DynMatrixType> ConstRefDynMat;
     43   typedef Ref<RealDynMatrixType , 0, Stride<Dynamic,Dynamic> > RefRealMatWithStride;
     44 
     45   Index rows = m.rows(), cols = m.cols();
     46   
     47   MatrixType  m1 = MatrixType::Random(rows, cols),
     48               m2 = m1;
     49   
     50   Index i = internal::random<Index>(0,rows-1);
     51   Index j = internal::random<Index>(0,cols-1);
     52   Index brows = internal::random<Index>(1,rows-i);
     53   Index bcols = internal::random<Index>(1,cols-j);
     54   
     55   RefMat rm0 = m1;
     56   VERIFY_IS_EQUAL(rm0, m1);
     57   RefDynMat rm1 = m1;
     58   VERIFY_IS_EQUAL(rm1, m1);
     59   RefDynMat rm2 = m1.block(i,j,brows,bcols);
     60   VERIFY_IS_EQUAL(rm2, m1.block(i,j,brows,bcols));
     61   rm2.setOnes();
     62   m2.block(i,j,brows,bcols).setOnes();
     63   VERIFY_IS_EQUAL(m1, m2);
     64   
     65   m2.block(i,j,brows,bcols).setRandom();
     66   rm2 = m2.block(i,j,brows,bcols);
     67   VERIFY_IS_EQUAL(m1, m2);
     68   
     69   ConstRefDynMat rm3 = m1.block(i,j,brows,bcols);
     70   m1.block(i,j,brows,bcols) *= 2;
     71   m2.block(i,j,brows,bcols) *= 2;
     72   VERIFY_IS_EQUAL(rm3, m2.block(i,j,brows,bcols));
     73   RefRealMatWithStride rm4 = m1.real();
     74   VERIFY_IS_EQUAL(rm4, m2.real());
     75   rm4.array() += 1;
     76   m2.real().array() += 1;
     77   VERIFY_IS_EQUAL(m1, m2);
     78 }
     79 
     80 template<typename VectorType> void ref_vector(const VectorType& m)
     81 {
     82   typedef typename VectorType::Scalar Scalar;
     83   typedef typename VectorType::RealScalar RealScalar;
     84   typedef Matrix<Scalar,Dynamic,1,VectorType::Options> DynMatrixType;
     85   typedef Matrix<Scalar,Dynamic,Dynamic,ColMajor> MatrixType;
     86   typedef Matrix<RealScalar,Dynamic,1,VectorType::Options> RealDynMatrixType;
     87   
     88   typedef Ref<VectorType> RefMat;
     89   typedef Ref<DynMatrixType> RefDynMat;
     90   typedef Ref<const DynMatrixType> ConstRefDynMat;
     91   typedef Ref<RealDynMatrixType , 0, InnerStride<> > RefRealMatWithStride;
     92   typedef Ref<DynMatrixType , 0, InnerStride<> > RefMatWithStride;
     93 
     94   Index size = m.size();
     95   
     96   VectorType  v1 = VectorType::Random(size),
     97               v2 = v1;
     98   MatrixType mat1 = MatrixType::Random(size,size),
     99              mat2 = mat1,
    100              mat3 = MatrixType::Random(size,size);
    101   
    102   Index i = internal::random<Index>(0,size-1);
    103   Index bsize = internal::random<Index>(1,size-i);
    104   
    105   { RefMat    rm0 = v1;                   VERIFY_IS_EQUAL(rm0, v1); }
    106   { RefMat    rm0 = v1.block(0,0,size,1); VERIFY_IS_EQUAL(rm0, v1); }
    107   { RefDynMat rv1 = v1;                   VERIFY_IS_EQUAL(rv1, v1); }
    108   { RefDynMat rv1 = v1.block(0,0,size,1); VERIFY_IS_EQUAL(rv1, v1); }
    109   { VERIFY_RAISES_ASSERT( RefMat    rm0 = v1.block(0, 0, size, 0); EIGEN_UNUSED_VARIABLE(rm0); ); }
    110   if(VectorType::SizeAtCompileTime!=1)
    111   { VERIFY_RAISES_ASSERT( RefDynMat rv1 = v1.block(0, 0, size, 0); EIGEN_UNUSED_VARIABLE(rv1); ); }
    112 
    113   RefDynMat rv2 = v1.segment(i,bsize);
    114   VERIFY_IS_EQUAL(rv2, v1.segment(i,bsize));
    115   rv2.setOnes();
    116   v2.segment(i,bsize).setOnes();
    117   VERIFY_IS_EQUAL(v1, v2);
    118   
    119   v2.segment(i,bsize).setRandom();
    120   rv2 = v2.segment(i,bsize);
    121   VERIFY_IS_EQUAL(v1, v2);
    122   
    123   ConstRefDynMat rm3 = v1.segment(i,bsize);
    124   v1.segment(i,bsize) *= 2;
    125   v2.segment(i,bsize) *= 2;
    126   VERIFY_IS_EQUAL(rm3, v2.segment(i,bsize));
    127   
    128   RefRealMatWithStride rm4 = v1.real();
    129   VERIFY_IS_EQUAL(rm4, v2.real());
    130   rm4.array() += 1;
    131   v2.real().array() += 1;
    132   VERIFY_IS_EQUAL(v1, v2);
    133   
    134   RefMatWithStride rm5 = mat1.row(i).transpose();
    135   VERIFY_IS_EQUAL(rm5, mat1.row(i).transpose());
    136   rm5.array() += 1;
    137   mat2.row(i).array() += 1;
    138   VERIFY_IS_EQUAL(mat1, mat2);
    139   rm5.noalias() = rm4.transpose() * mat3;
    140   mat2.row(i) = v2.real().transpose() * mat3;
    141   VERIFY_IS_APPROX(mat1, mat2);
    142 }
    143 
    144 template<typename Scalar, int Rows, int Cols>
    145 void ref_vector_fixed_sizes()
    146 {
    147   typedef Matrix<Scalar,Rows,Cols,RowMajor> RowMajorMatrixType;
    148   typedef Matrix<Scalar,Rows,Cols,ColMajor> ColMajorMatrixType;
    149   typedef Matrix<Scalar,1,Cols> RowVectorType;
    150   typedef Matrix<Scalar,Rows,1> ColVectorType;
    151   typedef Matrix<Scalar,Cols,1> RowVectorTransposeType;
    152   typedef Matrix<Scalar,1,Rows> ColVectorTransposeType;
    153   typedef Stride<Dynamic, Dynamic> DynamicStride;
    154 
    155   RowMajorMatrixType mr = RowMajorMatrixType::Random();
    156   ColMajorMatrixType mc = ColMajorMatrixType::Random();
    157 
    158   Index i = internal::random<Index>(0,Rows-1);
    159   Index j = internal::random<Index>(0,Cols-1);
    160 
    161   // Reference ith row.
    162   Ref<RowVectorType, 0, DynamicStride> mr_ri = mr.row(i);
    163   VERIFY_IS_EQUAL(mr_ri, mr.row(i));
    164   Ref<RowVectorType, 0, DynamicStride> mc_ri = mc.row(i);
    165   VERIFY_IS_EQUAL(mc_ri, mc.row(i));
    166 
    167   // Reference jth col.
    168   Ref<ColVectorType, 0, DynamicStride> mr_cj = mr.col(j);
    169   VERIFY_IS_EQUAL(mr_cj, mr.col(j));
    170   Ref<ColVectorType, 0, DynamicStride> mc_cj = mc.col(j);
    171   VERIFY_IS_EQUAL(mc_cj, mc.col(j));
    172 
    173   // Reference the transpose of row i.
    174   Ref<RowVectorTransposeType, 0, DynamicStride> mr_rit = mr.row(i);
    175   VERIFY_IS_EQUAL(mr_rit, mr.row(i).transpose());
    176   Ref<RowVectorTransposeType, 0, DynamicStride> mc_rit = mc.row(i);
    177   VERIFY_IS_EQUAL(mc_rit, mc.row(i).transpose());
    178 
    179   // Reference the transpose of col j.
    180   Ref<ColVectorTransposeType, 0, DynamicStride> mr_cjt = mr.col(j);
    181   VERIFY_IS_EQUAL(mr_cjt, mr.col(j).transpose());
    182   Ref<ColVectorTransposeType, 0, DynamicStride> mc_cjt = mc.col(j);
    183   VERIFY_IS_EQUAL(mc_cjt, mc.col(j).transpose());
    184   
    185   // Const references without strides.
    186   Ref<const RowVectorType> cmr_ri = mr.row(i);
    187   VERIFY_IS_EQUAL(cmr_ri, mr.row(i));
    188   Ref<const RowVectorType> cmc_ri = mc.row(i);
    189   VERIFY_IS_EQUAL(cmc_ri, mc.row(i));
    190 
    191   Ref<const ColVectorType> cmr_cj = mr.col(j);
    192   VERIFY_IS_EQUAL(cmr_cj, mr.col(j));
    193   Ref<const ColVectorType> cmc_cj = mc.col(j);
    194   VERIFY_IS_EQUAL(cmc_cj, mc.col(j));
    195 
    196   Ref<const RowVectorTransposeType> cmr_rit = mr.row(i);
    197   VERIFY_IS_EQUAL(cmr_rit, mr.row(i).transpose());
    198   Ref<const RowVectorTransposeType> cmc_rit = mc.row(i);
    199   VERIFY_IS_EQUAL(cmc_rit, mc.row(i).transpose());
    200 
    201   Ref<const ColVectorTransposeType> cmr_cjt = mr.col(j);
    202   VERIFY_IS_EQUAL(cmr_cjt, mr.col(j).transpose());
    203   Ref<const ColVectorTransposeType> cmc_cjt = mc.col(j);
    204   VERIFY_IS_EQUAL(cmc_cjt, mc.col(j).transpose());
    205 }
    206 
    207 template<typename PlainObjectType> void check_const_correctness(const PlainObjectType&)
    208 {
    209   // verify that ref-to-const don't have LvalueBit
    210   typedef typename internal::add_const<PlainObjectType>::type ConstPlainObjectType;
    211   VERIFY( !(internal::traits<Ref<ConstPlainObjectType> >::Flags & LvalueBit) );
    212   VERIFY( !(internal::traits<Ref<ConstPlainObjectType, Aligned> >::Flags & LvalueBit) );
    213   VERIFY( !(Ref<ConstPlainObjectType>::Flags & LvalueBit) );
    214   VERIFY( !(Ref<ConstPlainObjectType, Aligned>::Flags & LvalueBit) );
    215 }
    216 
    217 template<typename B>
    218 EIGEN_DONT_INLINE void call_ref_1(Ref<VectorXf> a, const B &b) { VERIFY_IS_EQUAL(a,b); }
    219 template<typename B>
    220 EIGEN_DONT_INLINE void call_ref_2(const Ref<const VectorXf>& a, const B &b) { VERIFY_IS_EQUAL(a,b); }
    221 template<typename B>
    222 EIGEN_DONT_INLINE void call_ref_3(Ref<VectorXf,0,InnerStride<> > a, const B &b) { VERIFY_IS_EQUAL(a,b); }
    223 template<typename B>
    224 EIGEN_DONT_INLINE void call_ref_4(const Ref<const VectorXf,0,InnerStride<> >& a, const B &b) { VERIFY_IS_EQUAL(a,b); }
    225 template<typename B>
    226 EIGEN_DONT_INLINE void call_ref_5(Ref<MatrixXf,0,OuterStride<> > a, const B &b) { VERIFY_IS_EQUAL(a,b); }
    227 template<typename B>
    228 EIGEN_DONT_INLINE void call_ref_6(const Ref<const MatrixXf,0,OuterStride<> >& a, const B &b) { VERIFY_IS_EQUAL(a,b); }
    229 template<typename B>
    230 EIGEN_DONT_INLINE void call_ref_7(Ref<Matrix<float,Dynamic,3> > a, const B &b) { VERIFY_IS_EQUAL(a,b); }
    231 
    232 void call_ref()
    233 {
    234   VectorXcf ca  = VectorXcf::Random(10);
    235   VectorXf a    = VectorXf::Random(10);
    236   RowVectorXf b = RowVectorXf::Random(10);
    237   MatrixXf A    = MatrixXf::Random(10,10);
    238   RowVector3f c = RowVector3f::Random();
    239   const VectorXf& ac(a);
    240   VectorBlock<VectorXf> ab(a,0,3);
    241   const VectorBlock<VectorXf> abc(a,0,3);
    242   
    243 
    244   VERIFY_EVALUATION_COUNT( call_ref_1(a,a), 0);
    245   VERIFY_EVALUATION_COUNT( call_ref_1(b,b.transpose()), 0);
    246 //   call_ref_1(ac,a<c);           // does not compile because ac is const
    247   VERIFY_EVALUATION_COUNT( call_ref_1(ab,ab), 0);
    248   VERIFY_EVALUATION_COUNT( call_ref_1(a.head(4),a.head(4)), 0);
    249   VERIFY_EVALUATION_COUNT( call_ref_1(abc,abc), 0);
    250   VERIFY_EVALUATION_COUNT( call_ref_1(A.col(3),A.col(3)), 0);
    251 //   call_ref_1(A.row(3),A.row(3));    // does not compile because innerstride!=1
    252   VERIFY_EVALUATION_COUNT( call_ref_3(A.row(3),A.row(3).transpose()), 0);
    253   VERIFY_EVALUATION_COUNT( call_ref_4(A.row(3),A.row(3).transpose()), 0);
    254 //   call_ref_1(a+a, a+a);          // does not compile for obvious reason
    255 
    256   MatrixXf tmp = A*A.col(1);
    257   VERIFY_EVALUATION_COUNT( call_ref_2(A*A.col(1), tmp), 1);     // evaluated into a temp
    258   VERIFY_EVALUATION_COUNT( call_ref_2(ac.head(5),ac.head(5)), 0);
    259   VERIFY_EVALUATION_COUNT( call_ref_2(ac,ac), 0);
    260   VERIFY_EVALUATION_COUNT( call_ref_2(a,a), 0);
    261   VERIFY_EVALUATION_COUNT( call_ref_2(ab,ab), 0);
    262   VERIFY_EVALUATION_COUNT( call_ref_2(a.head(4),a.head(4)), 0);
    263   tmp = a+a;
    264   VERIFY_EVALUATION_COUNT( call_ref_2(a+a,tmp), 1);            // evaluated into a temp
    265   VERIFY_EVALUATION_COUNT( call_ref_2(ca.imag(),ca.imag()), 1);      // evaluated into a temp
    266 
    267   VERIFY_EVALUATION_COUNT( call_ref_4(ac.head(5),ac.head(5)), 0);
    268   tmp = a+a;
    269   VERIFY_EVALUATION_COUNT( call_ref_4(a+a,tmp), 1);           // evaluated into a temp
    270   VERIFY_EVALUATION_COUNT( call_ref_4(ca.imag(),ca.imag()), 0);
    271 
    272   VERIFY_EVALUATION_COUNT( call_ref_5(a,a), 0);
    273   VERIFY_EVALUATION_COUNT( call_ref_5(a.head(3),a.head(3)), 0);
    274   VERIFY_EVALUATION_COUNT( call_ref_5(A,A), 0);
    275 //   call_ref_5(A.transpose(),A.transpose());   // does not compile because storage order does not match
    276   VERIFY_EVALUATION_COUNT( call_ref_5(A.block(1,1,2,2),A.block(1,1,2,2)), 0);
    277   VERIFY_EVALUATION_COUNT( call_ref_5(b,b), 0);             // storage order do not match, but this is a degenerate case that should work
    278   VERIFY_EVALUATION_COUNT( call_ref_5(a.row(3),a.row(3)), 0);
    279 
    280   VERIFY_EVALUATION_COUNT( call_ref_6(a,a), 0);
    281   VERIFY_EVALUATION_COUNT( call_ref_6(a.head(3),a.head(3)), 0);
    282   VERIFY_EVALUATION_COUNT( call_ref_6(A.row(3),A.row(3)), 1);           // evaluated into a temp thouth it could be avoided by viewing it as a 1xn matrix
    283   tmp = A+A;
    284   VERIFY_EVALUATION_COUNT( call_ref_6(A+A,tmp), 1);                // evaluated into a temp
    285   VERIFY_EVALUATION_COUNT( call_ref_6(A,A), 0);
    286   VERIFY_EVALUATION_COUNT( call_ref_6(A.transpose(),A.transpose()), 1);      // evaluated into a temp because the storage orders do not match
    287   VERIFY_EVALUATION_COUNT( call_ref_6(A.block(1,1,2,2),A.block(1,1,2,2)), 0);
    288   
    289   VERIFY_EVALUATION_COUNT( call_ref_7(c,c), 0);
    290 }
    291 
    292 typedef Matrix<double,Dynamic,Dynamic,RowMajor> RowMatrixXd;
    293 int test_ref_overload_fun1(Ref<MatrixXd> )       { return 1; }
    294 int test_ref_overload_fun1(Ref<RowMatrixXd> )    { return 2; }
    295 int test_ref_overload_fun1(Ref<MatrixXf> )       { return 3; }
    296 
    297 int test_ref_overload_fun2(Ref<const MatrixXd> ) { return 4; }
    298 int test_ref_overload_fun2(Ref<const MatrixXf> ) { return 5; }
    299 
    300 void test_ref_ambiguous(const Ref<const ArrayXd> &A, Ref<ArrayXd> B)
    301 {
    302   B = A;
    303   B = A - A;
    304 }
    305 
    306 // See also bug 969
    307 void test_ref_overloads()
    308 {
    309   MatrixXd Ad, Bd;
    310   RowMatrixXd rAd, rBd;
    311   VERIFY( test_ref_overload_fun1(Ad)==1 );
    312   VERIFY( test_ref_overload_fun1(rAd)==2 );
    313   
    314   MatrixXf Af, Bf;
    315   VERIFY( test_ref_overload_fun2(Ad)==4 );
    316   VERIFY( test_ref_overload_fun2(Ad+Bd)==4 );
    317   VERIFY( test_ref_overload_fun2(Af+Bf)==5 );
    318   
    319   ArrayXd A, B;
    320   test_ref_ambiguous(A, B);
    321 }
    322 
    323 void test_ref_fixed_size_assert()
    324 {
    325   Vector4f v4 = Vector4f::Random();
    326   VectorXf vx = VectorXf::Random(10);
    327   VERIFY_RAISES_STATIC_ASSERT( Ref<Vector3f> y = v4; (void)y; );
    328   VERIFY_RAISES_STATIC_ASSERT( Ref<Vector3f> y = vx.head<4>(); (void)y; );
    329   VERIFY_RAISES_STATIC_ASSERT( Ref<const Vector3f> y = v4; (void)y; );
    330   VERIFY_RAISES_STATIC_ASSERT( Ref<const Vector3f> y = vx.head<4>(); (void)y; );
    331   VERIFY_RAISES_STATIC_ASSERT( Ref<const Vector3f> y = 2*v4; (void)y; );
    332 }
    333 
    334 EIGEN_DECLARE_TEST(ref)
    335 {
    336   for(int i = 0; i < g_repeat; i++) {
    337     CALL_SUBTEST_1( ref_vector(Matrix<float, 1, 1>()) );
    338     CALL_SUBTEST_1( check_const_correctness(Matrix<float, 1, 1>()) );
    339     CALL_SUBTEST_2( ref_vector(Vector4d()) );
    340     CALL_SUBTEST_2( check_const_correctness(Matrix4d()) );
    341     CALL_SUBTEST_3( ref_vector(Vector4cf()) );
    342     CALL_SUBTEST_4( ref_vector(VectorXcf(8)) );
    343     CALL_SUBTEST_5( ref_vector(VectorXi(12)) );
    344     CALL_SUBTEST_5( check_const_correctness(VectorXi(12)) );
    345 
    346     CALL_SUBTEST_1( ref_matrix(Matrix<float, 1, 1>()) );
    347     CALL_SUBTEST_2( ref_matrix(Matrix4d()) );
    348     CALL_SUBTEST_1( ref_matrix(Matrix<float,3,5>()) );
    349     CALL_SUBTEST_4( ref_matrix(MatrixXcf(internal::random<int>(1,10),internal::random<int>(1,10))) );
    350     CALL_SUBTEST_4( ref_matrix(Matrix<std::complex<double>,10,15>()) );
    351     CALL_SUBTEST_5( ref_matrix(MatrixXi(internal::random<int>(1,10),internal::random<int>(1,10))) );
    352     CALL_SUBTEST_6( call_ref() );
    353 
    354     CALL_SUBTEST_8( (ref_vector_fixed_sizes<float,3,5>()) );
    355     CALL_SUBTEST_8( (ref_vector_fixed_sizes<float,15,10>()) );
    356   }
    357   
    358   CALL_SUBTEST_7( test_ref_overloads() );
    359   CALL_SUBTEST_7( test_ref_fixed_size_assert() );
    360 }