cart-elc

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

evaluators.cpp (21038B)


      1 
      2 #include "main.h"
      3 
      4 namespace Eigen {
      5 
      6   template<typename Lhs,typename Rhs>
      7   const Product<Lhs,Rhs>
      8   prod(const Lhs& lhs, const Rhs& rhs)
      9   {
     10     return Product<Lhs,Rhs>(lhs,rhs);
     11   }
     12 
     13   template<typename Lhs,typename Rhs>
     14   const Product<Lhs,Rhs,LazyProduct>
     15   lazyprod(const Lhs& lhs, const Rhs& rhs)
     16   {
     17     return Product<Lhs,Rhs,LazyProduct>(lhs,rhs);
     18   }
     19   
     20   template<typename DstXprType, typename SrcXprType>
     21   EIGEN_STRONG_INLINE
     22   DstXprType& copy_using_evaluator(const EigenBase<DstXprType> &dst, const SrcXprType &src)
     23   {
     24     call_assignment(dst.const_cast_derived(), src.derived(), internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>());
     25     return dst.const_cast_derived();
     26   }
     27   
     28   template<typename DstXprType, template <typename> class StorageBase, typename SrcXprType>
     29   EIGEN_STRONG_INLINE
     30   const DstXprType& copy_using_evaluator(const NoAlias<DstXprType, StorageBase>& dst, const SrcXprType &src)
     31   {
     32     call_assignment(dst, src.derived(), internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>());
     33     return dst.expression();
     34   }
     35   
     36   template<typename DstXprType, typename SrcXprType>
     37   EIGEN_STRONG_INLINE
     38   DstXprType& copy_using_evaluator(const PlainObjectBase<DstXprType> &dst, const SrcXprType &src)
     39   {
     40     #ifdef EIGEN_NO_AUTOMATIC_RESIZING
     41     eigen_assert((dst.size()==0 || (IsVectorAtCompileTime ? (dst.size() == src.size())
     42                                                           : (dst.rows() == src.rows() && dst.cols() == src.cols())))
     43                 && "Size mismatch. Automatic resizing is disabled because EIGEN_NO_AUTOMATIC_RESIZING is defined");
     44   #else
     45     dst.const_cast_derived().resizeLike(src.derived());
     46   #endif
     47     
     48     call_assignment(dst.const_cast_derived(), src.derived(), internal::assign_op<typename DstXprType::Scalar,typename SrcXprType::Scalar>());
     49     return dst.const_cast_derived();
     50   }
     51 
     52   template<typename DstXprType, typename SrcXprType>
     53   void add_assign_using_evaluator(const DstXprType& dst, const SrcXprType& src)
     54   {
     55     typedef typename DstXprType::Scalar Scalar;
     56     call_assignment(const_cast<DstXprType&>(dst), src.derived(), internal::add_assign_op<Scalar,typename SrcXprType::Scalar>());
     57   }
     58 
     59   template<typename DstXprType, typename SrcXprType>
     60   void subtract_assign_using_evaluator(const DstXprType& dst, const SrcXprType& src)
     61   {
     62     typedef typename DstXprType::Scalar Scalar;
     63     call_assignment(const_cast<DstXprType&>(dst), src.derived(), internal::sub_assign_op<Scalar,typename SrcXprType::Scalar>());
     64   }
     65 
     66   template<typename DstXprType, typename SrcXprType>
     67   void multiply_assign_using_evaluator(const DstXprType& dst, const SrcXprType& src)
     68   {
     69     typedef typename DstXprType::Scalar Scalar;
     70     call_assignment(dst.const_cast_derived(), src.derived(), internal::mul_assign_op<Scalar,typename SrcXprType::Scalar>());
     71   }
     72 
     73   template<typename DstXprType, typename SrcXprType>
     74   void divide_assign_using_evaluator(const DstXprType& dst, const SrcXprType& src)
     75   {
     76     typedef typename DstXprType::Scalar Scalar;
     77     call_assignment(dst.const_cast_derived(), src.derived(), internal::div_assign_op<Scalar,typename SrcXprType::Scalar>());
     78   }
     79   
     80   template<typename DstXprType, typename SrcXprType>
     81   void swap_using_evaluator(const DstXprType& dst, const SrcXprType& src)
     82   {
     83     typedef typename DstXprType::Scalar Scalar;
     84     call_assignment(dst.const_cast_derived(), src.const_cast_derived(), internal::swap_assign_op<Scalar>());
     85   }
     86 
     87   namespace internal {
     88     template<typename Dst, template <typename> class StorageBase, typename Src, typename Func>
     89     EIGEN_DEVICE_FUNC void call_assignment(const NoAlias<Dst,StorageBase>& dst, const Src& src, const Func& func)
     90     {
     91       call_assignment_no_alias(dst.expression(), src, func);
     92     }
     93 
     94     template<typename Dst, template <typename> class StorageBase, typename Src, typename Func>
     95     EIGEN_DEVICE_FUNC void call_restricted_packet_assignment(const NoAlias<Dst,StorageBase>& dst, const Src& src, const Func& func)
     96     {
     97       call_restricted_packet_assignment_no_alias(dst.expression(), src, func);
     98     }
     99   }
    100   
    101 }
    102 
    103 template<typename XprType> long get_cost(const XprType& ) { return Eigen::internal::evaluator<XprType>::CoeffReadCost; }
    104 
    105 using namespace std;
    106 
    107 #define VERIFY_IS_APPROX_EVALUATOR(DEST,EXPR) VERIFY_IS_APPROX(copy_using_evaluator(DEST,(EXPR)), (EXPR).eval());
    108 #define VERIFY_IS_APPROX_EVALUATOR2(DEST,EXPR,REF) VERIFY_IS_APPROX(copy_using_evaluator(DEST,(EXPR)), (REF).eval());
    109 
    110 EIGEN_DECLARE_TEST(evaluators)
    111 {
    112   // Testing Matrix evaluator and Transpose
    113   Vector2d v = Vector2d::Random();
    114   const Vector2d v_const(v);
    115   Vector2d v2;
    116   RowVector2d w;
    117 
    118   VERIFY_IS_APPROX_EVALUATOR(v2, v);
    119   VERIFY_IS_APPROX_EVALUATOR(v2, v_const);
    120 
    121   // Testing Transpose
    122   VERIFY_IS_APPROX_EVALUATOR(w, v.transpose()); // Transpose as rvalue
    123   VERIFY_IS_APPROX_EVALUATOR(w, v_const.transpose());
    124 
    125   copy_using_evaluator(w.transpose(), v); // Transpose as lvalue
    126   VERIFY_IS_APPROX(w,v.transpose().eval());
    127 
    128   copy_using_evaluator(w.transpose(), v_const);
    129   VERIFY_IS_APPROX(w,v_const.transpose().eval());
    130 
    131   // Testing Array evaluator
    132   {
    133     ArrayXXf a(2,3);
    134     ArrayXXf b(3,2);
    135     a << 1,2,3, 4,5,6;
    136     const ArrayXXf a_const(a);
    137 
    138     VERIFY_IS_APPROX_EVALUATOR(b, a.transpose());
    139 
    140     VERIFY_IS_APPROX_EVALUATOR(b, a_const.transpose());
    141 
    142     // Testing CwiseNullaryOp evaluator
    143     copy_using_evaluator(w, RowVector2d::Random());
    144     VERIFY((w.array() >= -1).all() && (w.array() <= 1).all()); // not easy to test ...
    145 
    146     VERIFY_IS_APPROX_EVALUATOR(w, RowVector2d::Zero());
    147 
    148     VERIFY_IS_APPROX_EVALUATOR(w, RowVector2d::Constant(3));
    149     
    150     // mix CwiseNullaryOp and transpose
    151     VERIFY_IS_APPROX_EVALUATOR(w, Vector2d::Zero().transpose());
    152   }
    153 
    154   {
    155     // test product expressions
    156     int s = internal::random<int>(1,100);
    157     MatrixXf a(s,s), b(s,s), c(s,s), d(s,s);
    158     a.setRandom();
    159     b.setRandom();
    160     c.setRandom();
    161     d.setRandom();
    162     VERIFY_IS_APPROX_EVALUATOR(d, (a + b));
    163     VERIFY_IS_APPROX_EVALUATOR(d, (a + b).transpose());
    164     VERIFY_IS_APPROX_EVALUATOR2(d, prod(a,b), a*b);
    165     VERIFY_IS_APPROX_EVALUATOR2(d.noalias(), prod(a,b), a*b);
    166     VERIFY_IS_APPROX_EVALUATOR2(d, prod(a,b) + c, a*b + c);
    167     VERIFY_IS_APPROX_EVALUATOR2(d, s * prod(a,b), s * a*b);
    168     VERIFY_IS_APPROX_EVALUATOR2(d, prod(a,b).transpose(), (a*b).transpose());
    169     VERIFY_IS_APPROX_EVALUATOR2(d, prod(a,b) + prod(b,c), a*b + b*c);
    170 
    171     // check that prod works even with aliasing present
    172     c = a*a;
    173     copy_using_evaluator(a, prod(a,a));
    174     VERIFY_IS_APPROX(a,c);
    175 
    176     // check compound assignment of products
    177     d = c;
    178     add_assign_using_evaluator(c.noalias(), prod(a,b));
    179     d.noalias() += a*b;
    180     VERIFY_IS_APPROX(c, d);
    181 
    182     d = c;
    183     subtract_assign_using_evaluator(c.noalias(), prod(a,b));
    184     d.noalias() -= a*b;
    185     VERIFY_IS_APPROX(c, d);
    186   }
    187 
    188   {
    189     // test product with all possible sizes
    190     int s = internal::random<int>(1,100);
    191     Matrix<float,      1,      1> m11, res11;  m11.setRandom(1,1);
    192     Matrix<float,      1,      4> m14, res14;  m14.setRandom(1,4);
    193     Matrix<float,      1,Dynamic> m1X, res1X;  m1X.setRandom(1,s);
    194     Matrix<float,      4,      1> m41, res41;  m41.setRandom(4,1);
    195     Matrix<float,      4,      4> m44, res44;  m44.setRandom(4,4);
    196     Matrix<float,      4,Dynamic> m4X, res4X;  m4X.setRandom(4,s);
    197     Matrix<float,Dynamic,      1> mX1, resX1;  mX1.setRandom(s,1);
    198     Matrix<float,Dynamic,      4> mX4, resX4;  mX4.setRandom(s,4);
    199     Matrix<float,Dynamic,Dynamic> mXX, resXX;  mXX.setRandom(s,s);
    200 
    201     VERIFY_IS_APPROX_EVALUATOR2(res11, prod(m11,m11), m11*m11);
    202     VERIFY_IS_APPROX_EVALUATOR2(res11, prod(m14,m41), m14*m41);
    203     VERIFY_IS_APPROX_EVALUATOR2(res11, prod(m1X,mX1), m1X*mX1);
    204     VERIFY_IS_APPROX_EVALUATOR2(res14, prod(m11,m14), m11*m14);
    205     VERIFY_IS_APPROX_EVALUATOR2(res14, prod(m14,m44), m14*m44);
    206     VERIFY_IS_APPROX_EVALUATOR2(res14, prod(m1X,mX4), m1X*mX4);
    207     VERIFY_IS_APPROX_EVALUATOR2(res1X, prod(m11,m1X), m11*m1X);
    208     VERIFY_IS_APPROX_EVALUATOR2(res1X, prod(m14,m4X), m14*m4X);
    209     VERIFY_IS_APPROX_EVALUATOR2(res1X, prod(m1X,mXX), m1X*mXX);
    210     VERIFY_IS_APPROX_EVALUATOR2(res41, prod(m41,m11), m41*m11);
    211     VERIFY_IS_APPROX_EVALUATOR2(res41, prod(m44,m41), m44*m41);
    212     VERIFY_IS_APPROX_EVALUATOR2(res41, prod(m4X,mX1), m4X*mX1);
    213     VERIFY_IS_APPROX_EVALUATOR2(res44, prod(m41,m14), m41*m14);
    214     VERIFY_IS_APPROX_EVALUATOR2(res44, prod(m44,m44), m44*m44);
    215     VERIFY_IS_APPROX_EVALUATOR2(res44, prod(m4X,mX4), m4X*mX4);
    216     VERIFY_IS_APPROX_EVALUATOR2(res4X, prod(m41,m1X), m41*m1X);
    217     VERIFY_IS_APPROX_EVALUATOR2(res4X, prod(m44,m4X), m44*m4X);
    218     VERIFY_IS_APPROX_EVALUATOR2(res4X, prod(m4X,mXX), m4X*mXX);
    219     VERIFY_IS_APPROX_EVALUATOR2(resX1, prod(mX1,m11), mX1*m11);
    220     VERIFY_IS_APPROX_EVALUATOR2(resX1, prod(mX4,m41), mX4*m41);
    221     VERIFY_IS_APPROX_EVALUATOR2(resX1, prod(mXX,mX1), mXX*mX1);
    222     VERIFY_IS_APPROX_EVALUATOR2(resX4, prod(mX1,m14), mX1*m14);
    223     VERIFY_IS_APPROX_EVALUATOR2(resX4, prod(mX4,m44), mX4*m44);
    224     VERIFY_IS_APPROX_EVALUATOR2(resX4, prod(mXX,mX4), mXX*mX4);
    225     VERIFY_IS_APPROX_EVALUATOR2(resXX, prod(mX1,m1X), mX1*m1X);
    226     VERIFY_IS_APPROX_EVALUATOR2(resXX, prod(mX4,m4X), mX4*m4X);
    227     VERIFY_IS_APPROX_EVALUATOR2(resXX, prod(mXX,mXX), mXX*mXX);
    228   }
    229 
    230   {
    231     ArrayXXf a(2,3);
    232     ArrayXXf b(3,2);
    233     a << 1,2,3, 4,5,6;
    234     const ArrayXXf a_const(a);
    235     
    236     // this does not work because Random is eval-before-nested: 
    237     // copy_using_evaluator(w, Vector2d::Random().transpose());
    238 
    239     // test CwiseUnaryOp
    240     VERIFY_IS_APPROX_EVALUATOR(v2, 3 * v);
    241     VERIFY_IS_APPROX_EVALUATOR(w, (3 * v).transpose());
    242     VERIFY_IS_APPROX_EVALUATOR(b, (a + 3).transpose());
    243     VERIFY_IS_APPROX_EVALUATOR(b, (2 * a_const + 3).transpose());
    244 
    245     // test CwiseBinaryOp
    246     VERIFY_IS_APPROX_EVALUATOR(v2, v + Vector2d::Ones());
    247     VERIFY_IS_APPROX_EVALUATOR(w, (v + Vector2d::Ones()).transpose().cwiseProduct(RowVector2d::Constant(3)));
    248 
    249     // dynamic matrices and arrays
    250     MatrixXd mat1(6,6), mat2(6,6);
    251     VERIFY_IS_APPROX_EVALUATOR(mat1, MatrixXd::Identity(6,6));
    252     VERIFY_IS_APPROX_EVALUATOR(mat2, mat1);
    253     copy_using_evaluator(mat2.transpose(), mat1);
    254     VERIFY_IS_APPROX(mat2.transpose(), mat1);
    255 
    256     ArrayXXd arr1(6,6), arr2(6,6);
    257     VERIFY_IS_APPROX_EVALUATOR(arr1, ArrayXXd::Constant(6,6, 3.0));
    258     VERIFY_IS_APPROX_EVALUATOR(arr2, arr1);
    259     
    260     // test automatic resizing
    261     mat2.resize(3,3);
    262     VERIFY_IS_APPROX_EVALUATOR(mat2, mat1);
    263     arr2.resize(9,9);
    264     VERIFY_IS_APPROX_EVALUATOR(arr2, arr1);
    265 
    266     // test direct traversal
    267     Matrix3f m3;
    268     Array33f a3;
    269     VERIFY_IS_APPROX_EVALUATOR(m3, Matrix3f::Identity());  // matrix, nullary
    270     // TODO: find a way to test direct traversal with array
    271     VERIFY_IS_APPROX_EVALUATOR(m3.transpose(), Matrix3f::Identity().transpose());  // transpose
    272     VERIFY_IS_APPROX_EVALUATOR(m3, 2 * Matrix3f::Identity());  // unary
    273     VERIFY_IS_APPROX_EVALUATOR(m3, Matrix3f::Identity() + Matrix3f::Zero());  // binary
    274     VERIFY_IS_APPROX_EVALUATOR(m3.block(0,0,2,2), Matrix3f::Identity().block(1,1,2,2));  // block
    275 
    276     // test linear traversal
    277     VERIFY_IS_APPROX_EVALUATOR(m3, Matrix3f::Zero());  // matrix, nullary
    278     VERIFY_IS_APPROX_EVALUATOR(a3, Array33f::Zero());  // array
    279     VERIFY_IS_APPROX_EVALUATOR(m3.transpose(), Matrix3f::Zero().transpose());  // transpose
    280     VERIFY_IS_APPROX_EVALUATOR(m3, 2 * Matrix3f::Zero());  // unary
    281     VERIFY_IS_APPROX_EVALUATOR(m3, Matrix3f::Zero() + m3);  // binary  
    282 
    283     // test inner vectorization
    284     Matrix4f m4, m4src = Matrix4f::Random();
    285     Array44f a4, a4src = Matrix4f::Random();
    286     VERIFY_IS_APPROX_EVALUATOR(m4, m4src);  // matrix
    287     VERIFY_IS_APPROX_EVALUATOR(a4, a4src);  // array
    288     VERIFY_IS_APPROX_EVALUATOR(m4.transpose(), m4src.transpose());  // transpose
    289     // TODO: find out why Matrix4f::Zero() does not allow inner vectorization
    290     VERIFY_IS_APPROX_EVALUATOR(m4, 2 * m4src);  // unary
    291     VERIFY_IS_APPROX_EVALUATOR(m4, m4src + m4src);  // binary
    292 
    293     // test linear vectorization
    294     MatrixXf mX(6,6), mXsrc = MatrixXf::Random(6,6);
    295     ArrayXXf aX(6,6), aXsrc = ArrayXXf::Random(6,6);
    296     VERIFY_IS_APPROX_EVALUATOR(mX, mXsrc);  // matrix
    297     VERIFY_IS_APPROX_EVALUATOR(aX, aXsrc);  // array
    298     VERIFY_IS_APPROX_EVALUATOR(mX.transpose(), mXsrc.transpose());  // transpose
    299     VERIFY_IS_APPROX_EVALUATOR(mX, MatrixXf::Zero(6,6));  // nullary
    300     VERIFY_IS_APPROX_EVALUATOR(mX, 2 * mXsrc);  // unary
    301     VERIFY_IS_APPROX_EVALUATOR(mX, mXsrc + mXsrc);  // binary
    302 
    303     // test blocks and slice vectorization
    304     VERIFY_IS_APPROX_EVALUATOR(m4, (mXsrc.block<4,4>(1,0)));
    305     VERIFY_IS_APPROX_EVALUATOR(aX, ArrayXXf::Constant(10, 10, 3.0).block(2, 3, 6, 6));
    306 
    307     Matrix4f m4ref = m4;
    308     copy_using_evaluator(m4.block(1, 1, 2, 3), m3.bottomRows(2));
    309     m4ref.block(1, 1, 2, 3) = m3.bottomRows(2);
    310     VERIFY_IS_APPROX(m4, m4ref);
    311 
    312     mX.setIdentity(20,20);
    313     MatrixXf mXref = MatrixXf::Identity(20,20);
    314     mXsrc = MatrixXf::Random(9,12);
    315     copy_using_evaluator(mX.block(4, 4, 9, 12), mXsrc);
    316     mXref.block(4, 4, 9, 12) = mXsrc;
    317     VERIFY_IS_APPROX(mX, mXref);
    318 
    319     // test Map
    320     const float raw[3] = {1,2,3};
    321     float buffer[3] = {0,0,0};
    322     Vector3f v3;
    323     Array3f a3f;
    324     VERIFY_IS_APPROX_EVALUATOR(v3, Map<const Vector3f>(raw));
    325     VERIFY_IS_APPROX_EVALUATOR(a3f, Map<const Array3f>(raw));
    326     Vector3f::Map(buffer) = 2*v3;
    327     VERIFY(buffer[0] == 2);
    328     VERIFY(buffer[1] == 4);
    329     VERIFY(buffer[2] == 6);
    330 
    331     // test CwiseUnaryView
    332     mat1.setRandom();
    333     mat2.setIdentity();
    334     MatrixXcd matXcd(6,6), matXcd_ref(6,6);
    335     copy_using_evaluator(matXcd.real(), mat1);
    336     copy_using_evaluator(matXcd.imag(), mat2);
    337     matXcd_ref.real() = mat1;
    338     matXcd_ref.imag() = mat2;
    339     VERIFY_IS_APPROX(matXcd, matXcd_ref);
    340 
    341     // test Select
    342     VERIFY_IS_APPROX_EVALUATOR(aX, (aXsrc > 0).select(aXsrc, -aXsrc));
    343 
    344     // test Replicate
    345     mXsrc = MatrixXf::Random(6, 6);
    346     VectorXf vX = VectorXf::Random(6);
    347     mX.resize(6, 6);
    348     VERIFY_IS_APPROX_EVALUATOR(mX, mXsrc.colwise() + vX);
    349     matXcd.resize(12, 12);
    350     VERIFY_IS_APPROX_EVALUATOR(matXcd, matXcd_ref.replicate(2,2));
    351     VERIFY_IS_APPROX_EVALUATOR(matXcd, (matXcd_ref.replicate<2,2>()));
    352 
    353     // test partial reductions
    354     VectorXd vec1(6);
    355     VERIFY_IS_APPROX_EVALUATOR(vec1, mat1.rowwise().sum());
    356     VERIFY_IS_APPROX_EVALUATOR(vec1, mat1.colwise().sum().transpose());
    357 
    358     // test MatrixWrapper and ArrayWrapper
    359     mat1.setRandom(6,6);
    360     arr1.setRandom(6,6);
    361     VERIFY_IS_APPROX_EVALUATOR(mat2, arr1.matrix());
    362     VERIFY_IS_APPROX_EVALUATOR(arr2, mat1.array());
    363     VERIFY_IS_APPROX_EVALUATOR(mat2, (arr1 + 2).matrix());
    364     VERIFY_IS_APPROX_EVALUATOR(arr2, mat1.array() + 2);
    365     mat2.array() = arr1 * arr1;
    366     VERIFY_IS_APPROX(mat2, (arr1 * arr1).matrix());
    367     arr2.matrix() = MatrixXd::Identity(6,6);
    368     VERIFY_IS_APPROX(arr2, MatrixXd::Identity(6,6).array());
    369 
    370     // test Reverse
    371     VERIFY_IS_APPROX_EVALUATOR(arr2, arr1.reverse());
    372     VERIFY_IS_APPROX_EVALUATOR(arr2, arr1.colwise().reverse());
    373     VERIFY_IS_APPROX_EVALUATOR(arr2, arr1.rowwise().reverse());
    374     arr2.reverse() = arr1;
    375     VERIFY_IS_APPROX(arr2, arr1.reverse());
    376     mat2.array() = mat1.array().reverse();
    377     VERIFY_IS_APPROX(mat2.array(), mat1.array().reverse());
    378 
    379     // test Diagonal
    380     VERIFY_IS_APPROX_EVALUATOR(vec1, mat1.diagonal());
    381     vec1.resize(5);
    382     VERIFY_IS_APPROX_EVALUATOR(vec1, mat1.diagonal(1));
    383     VERIFY_IS_APPROX_EVALUATOR(vec1, mat1.diagonal<-1>());
    384     vec1.setRandom();
    385 
    386     mat2 = mat1;
    387     copy_using_evaluator(mat1.diagonal(1), vec1);
    388     mat2.diagonal(1) = vec1;
    389     VERIFY_IS_APPROX(mat1, mat2);
    390 
    391     copy_using_evaluator(mat1.diagonal<-1>(), mat1.diagonal(1));
    392     mat2.diagonal<-1>() = mat2.diagonal(1);
    393     VERIFY_IS_APPROX(mat1, mat2);
    394   }
    395   
    396   {
    397     // test swapping
    398     MatrixXd mat1, mat2, mat1ref, mat2ref;
    399     mat1ref = mat1 = MatrixXd::Random(6, 6);
    400     mat2ref = mat2 = 2 * mat1 + MatrixXd::Identity(6, 6);
    401     swap_using_evaluator(mat1, mat2);
    402     mat1ref.swap(mat2ref);
    403     VERIFY_IS_APPROX(mat1, mat1ref);
    404     VERIFY_IS_APPROX(mat2, mat2ref);
    405 
    406     swap_using_evaluator(mat1.block(0, 0, 3, 3), mat2.block(3, 3, 3, 3));
    407     mat1ref.block(0, 0, 3, 3).swap(mat2ref.block(3, 3, 3, 3));
    408     VERIFY_IS_APPROX(mat1, mat1ref);
    409     VERIFY_IS_APPROX(mat2, mat2ref);
    410 
    411     swap_using_evaluator(mat1.row(2), mat2.col(3).transpose());
    412     mat1.row(2).swap(mat2.col(3).transpose());
    413     VERIFY_IS_APPROX(mat1, mat1ref);
    414     VERIFY_IS_APPROX(mat2, mat2ref);
    415   }
    416 
    417   {
    418     // test compound assignment
    419     const Matrix4d mat_const = Matrix4d::Random(); 
    420     Matrix4d mat, mat_ref;
    421     mat = mat_ref = Matrix4d::Identity();
    422     add_assign_using_evaluator(mat, mat_const);
    423     mat_ref += mat_const;
    424     VERIFY_IS_APPROX(mat, mat_ref);
    425 
    426     subtract_assign_using_evaluator(mat.row(1), 2*mat.row(2));
    427     mat_ref.row(1) -= 2*mat_ref.row(2);
    428     VERIFY_IS_APPROX(mat, mat_ref);
    429 
    430     const ArrayXXf arr_const = ArrayXXf::Random(5,3); 
    431     ArrayXXf arr, arr_ref;
    432     arr = arr_ref = ArrayXXf::Constant(5, 3, 0.5);
    433     multiply_assign_using_evaluator(arr, arr_const);
    434     arr_ref *= arr_const;
    435     VERIFY_IS_APPROX(arr, arr_ref);
    436 
    437     divide_assign_using_evaluator(arr.row(1), arr.row(2) + 1);
    438     arr_ref.row(1) /= (arr_ref.row(2) + 1);
    439     VERIFY_IS_APPROX(arr, arr_ref);
    440   }
    441   
    442   {
    443     // test triangular shapes
    444     MatrixXd A = MatrixXd::Random(6,6), B(6,6), C(6,6), D(6,6);
    445     A.setRandom();B.setRandom();
    446     VERIFY_IS_APPROX_EVALUATOR2(B, A.triangularView<Upper>(), MatrixXd(A.triangularView<Upper>()));
    447     
    448     A.setRandom();B.setRandom();
    449     VERIFY_IS_APPROX_EVALUATOR2(B, A.triangularView<UnitLower>(), MatrixXd(A.triangularView<UnitLower>()));
    450     
    451     A.setRandom();B.setRandom();
    452     VERIFY_IS_APPROX_EVALUATOR2(B, A.triangularView<UnitUpper>(), MatrixXd(A.triangularView<UnitUpper>()));
    453     
    454     A.setRandom();B.setRandom();
    455     C = B; C.triangularView<Upper>() = A;
    456     copy_using_evaluator(B.triangularView<Upper>(), A);
    457     VERIFY(B.isApprox(C) && "copy_using_evaluator(B.triangularView<Upper>(), A)");
    458     
    459     A.setRandom();B.setRandom();
    460     C = B; C.triangularView<Lower>() = A.triangularView<Lower>();
    461     copy_using_evaluator(B.triangularView<Lower>(), A.triangularView<Lower>());
    462     VERIFY(B.isApprox(C) && "copy_using_evaluator(B.triangularView<Lower>(), A.triangularView<Lower>())");
    463     
    464     
    465     A.setRandom();B.setRandom();
    466     C = B; C.triangularView<Lower>() = A.triangularView<Upper>().transpose();
    467     copy_using_evaluator(B.triangularView<Lower>(), A.triangularView<Upper>().transpose());
    468     VERIFY(B.isApprox(C) && "copy_using_evaluator(B.triangularView<Lower>(), A.triangularView<Lower>().transpose())");
    469     
    470     
    471     A.setRandom();B.setRandom(); C = B; D = A;
    472     C.triangularView<Upper>().swap(D.triangularView<Upper>());
    473     swap_using_evaluator(B.triangularView<Upper>(), A.triangularView<Upper>());
    474     VERIFY(B.isApprox(C) && "swap_using_evaluator(B.triangularView<Upper>(), A.triangularView<Upper>())");
    475     
    476     
    477     VERIFY_IS_APPROX_EVALUATOR2(B, prod(A.triangularView<Upper>(),A), MatrixXd(A.triangularView<Upper>()*A));
    478     
    479     VERIFY_IS_APPROX_EVALUATOR2(B, prod(A.selfadjointView<Upper>(),A), MatrixXd(A.selfadjointView<Upper>()*A));
    480   }
    481 
    482   {
    483     // test diagonal shapes
    484     VectorXd d = VectorXd::Random(6);
    485     MatrixXd A = MatrixXd::Random(6,6), B(6,6);
    486     A.setRandom();B.setRandom();
    487     
    488     VERIFY_IS_APPROX_EVALUATOR2(B, lazyprod(d.asDiagonal(),A), MatrixXd(d.asDiagonal()*A));
    489     VERIFY_IS_APPROX_EVALUATOR2(B, lazyprod(A,d.asDiagonal()), MatrixXd(A*d.asDiagonal()));
    490   }
    491 
    492   {
    493     // test CoeffReadCost
    494     Matrix4d a, b;
    495     VERIFY_IS_EQUAL( get_cost(a), 1 );
    496     VERIFY_IS_EQUAL( get_cost(a+b), 3);
    497     VERIFY_IS_EQUAL( get_cost(2*a+b), 4);
    498     VERIFY_IS_EQUAL( get_cost(a*b), 1);
    499     VERIFY_IS_EQUAL( get_cost(a.lazyProduct(b)), 15);
    500     VERIFY_IS_EQUAL( get_cost(a*(a*b)), 1);
    501     VERIFY_IS_EQUAL( get_cost(a.lazyProduct(a*b)), 15);
    502     VERIFY_IS_EQUAL( get_cost(a*(a+b)), 1);
    503     VERIFY_IS_EQUAL( get_cost(a.lazyProduct(a+b)), 15);
    504   }
    505 
    506   // regression test for PR 544 and bug 1622 (introduced in #71609c4)
    507   {
    508     // test restricted_packet_assignment with an unaligned destination
    509     const size_t M = 2;
    510     const size_t K = 2;
    511     const size_t N = 5;
    512     float *destMem = new float[(M*N) + 1];
    513     float *dest = (internal::UIntPtr(destMem)%EIGEN_MAX_ALIGN_BYTES) == 0 ? destMem+1 : destMem;
    514 
    515     const Matrix<float, Dynamic, Dynamic, RowMajor> a = Matrix<float, Dynamic, Dynamic, RowMajor>::Random(M, K);
    516     const Matrix<float, Dynamic, Dynamic, RowMajor> b = Matrix<float, Dynamic, Dynamic, RowMajor>::Random(K, N);
    517     
    518     Map<Matrix<float, Dynamic, Dynamic, RowMajor> > z(dest, M, N);;
    519     Product<Matrix<float, Dynamic, Dynamic, RowMajor>, Matrix<float, Dynamic, Dynamic, RowMajor>, LazyProduct> tmp(a,b);
    520     internal::call_restricted_packet_assignment(z.noalias(), tmp.derived(), internal::assign_op<float, float>());
    521     
    522     VERIFY_IS_APPROX(z, a*b);
    523     delete[] destMem;
    524   }
    525 }