src/Core/ProductEvaluators.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
0005 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
0006 // Copyright (C) 2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
0007 //
0008 // This Source Code Form is subject to the terms of the Mozilla
0009 // Public License v. 2.0. If a copy of the MPL was not distributed
0010 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0011
0012
0013 #ifndef EIGEN_PRODUCTEVALUATORS_H
0014 #define EIGEN_PRODUCTEVALUATORS_H
0015
0016 namespace Eigen {
0017
0018 namespace internal {
0019
0020 /** \internal
0021   * Evaluator of a product expression.
0022   * Since products require special treatments to handle all possible cases,
0023   * we simply defer the evaluation logic to a product_evaluator class
0024   * which offers more partial specialization possibilities.
0025   *
0026   * \sa class product_evaluator
0027   */
0028 template<typename Lhs, typename Rhs, int Options>
0029 struct evaluator<Product<Lhs, Rhs, Options> >
0030  : public product_evaluator<Product<Lhs, Rhs, Options> >
0031 {
0032   typedef Product<Lhs, Rhs, Options> XprType;
0033   typedef product_evaluator<XprType> Base;
0034
0035   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit evaluator(const XprType& xpr) : Base(xpr) {}
0036 };
0037
0038 // Catch "scalar * ( A * B )" and transform it to "(A*scalar) * B"
0039 // TODO we should apply that rule only if that's really helpful
0040 template<typename Lhs, typename Rhs, typename Scalar1, typename Scalar2, typename Plain1>
0041 struct evaluator_assume_aliasing<CwiseBinaryOp<internal::scalar_product_op<Scalar1,Scalar2>,
0042                                                const CwiseNullaryOp<internal::scalar_constant_op<Scalar1>, Plain1>,
0043                                                const Product<Lhs, Rhs, DefaultProduct> > >
0044 {
0045   static const bool value = true;
0046 };
0047 template<typename Lhs, typename Rhs, typename Scalar1, typename Scalar2, typename Plain1>
0048 struct evaluator<CwiseBinaryOp<internal::scalar_product_op<Scalar1,Scalar2>,
0049                                const CwiseNullaryOp<internal::scalar_constant_op<Scalar1>, Plain1>,
0050                                const Product<Lhs, Rhs, DefaultProduct> > >
0051  : public evaluator<Product<EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar1,Lhs,product), Rhs, DefaultProduct> >
0052 {
0053   typedef CwiseBinaryOp<internal::scalar_product_op<Scalar1,Scalar2>,
0054                                const CwiseNullaryOp<internal::scalar_constant_op<Scalar1>, Plain1>,
0055                                const Product<Lhs, Rhs, DefaultProduct> > XprType;
0056   typedef evaluator<Product<EIGEN_SCALAR_BINARYOP_EXPR_RETURN_TYPE(Scalar1,Lhs,product), Rhs, DefaultProduct> > Base;
0057
0058   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit evaluator(const XprType& xpr)
0059     : Base(xpr.lhs().functor().m_other * xpr.rhs().lhs() * xpr.rhs().rhs())
0060   {}
0061 };
0062
0063
0064 template<typename Lhs, typename Rhs, int DiagIndex>
0065 struct evaluator<Diagonal<const Product<Lhs, Rhs, DefaultProduct>, DiagIndex> >
0066  : public evaluator<Diagonal<const Product<Lhs, Rhs, LazyProduct>, DiagIndex> >
0067 {
0068   typedef Diagonal<const Product<Lhs, Rhs, DefaultProduct>, DiagIndex> XprType;
0069   typedef evaluator<Diagonal<const Product<Lhs, Rhs, LazyProduct>, DiagIndex> > Base;
0070
0071   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit evaluator(const XprType& xpr)
0072     : Base(Diagonal<const Product<Lhs, Rhs, LazyProduct>, DiagIndex>(
0073         Product<Lhs, Rhs, LazyProduct>(xpr.nestedExpression().lhs(), xpr.nestedExpression().rhs()),
0074         xpr.index() ))
0075   {}
0076 };
0077
0078
0079 // Helper class to perform a matrix product with the destination at hand.
0080 // Depending on the sizes of the factors, there are different evaluation strategies
0081 // as controlled by internal::product_type.
0082 template< typename Lhs, typename Rhs,
0083           typename LhsShape = typename evaluator_traits<Lhs>::Shape,
0084           typename RhsShape = typename evaluator_traits<Rhs>::Shape,
0085           int ProductType = internal::product_type<Lhs,Rhs>::value>
0086 struct generic_product_impl;
0087
0088 template<typename Lhs, typename Rhs>
0089 struct evaluator_assume_aliasing<Product<Lhs, Rhs, DefaultProduct> > {
0090   static const bool value = true;
0091 };
0092
0093 // This is the default evaluator implementation for products:
0094 // It creates a temporary and call generic_product_impl
0095 template<typename Lhs, typename Rhs, int Options, int ProductTag, typename LhsShape, typename RhsShape>
0096 struct product_evaluator<Product<Lhs, Rhs, Options>, ProductTag, LhsShape, RhsShape>
0097   : public evaluator<typename Product<Lhs, Rhs, Options>::PlainObject>
0098 {
0099   typedef Product<Lhs, Rhs, Options> XprType;
0100   typedef typename XprType::PlainObject PlainObject;
0101   typedef evaluator<PlainObject> Base;
0102   enum {
0103     Flags = Base::Flags | EvalBeforeNestingBit
0104   };
0105
0106   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0107   explicit product_evaluator(const XprType& xpr)
0108     : m_result(xpr.rows(), xpr.cols())
0109   {
0110     ::new (static_cast<Base*>(this)) Base(m_result);
0111
0112 // FIXME shall we handle nested_eval here?,
0113 // if so, then we must take care at removing the call to nested_eval in the specializations (e.g., in permutation_matrix_product, transposition_matrix_product, etc.)
0114 //     typedef typename internal::nested_eval<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
0115 //     typedef typename internal::nested_eval<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
0116 //     typedef typename internal::remove_all<LhsNested>::type LhsNestedCleaned;
0117 //     typedef typename internal::remove_all<RhsNested>::type RhsNestedCleaned;
0118 //
0119 //     const LhsNested lhs(xpr.lhs());
0120 //     const RhsNested rhs(xpr.rhs());
0121 //
0122 //     generic_product_impl<LhsNestedCleaned, RhsNestedCleaned>::evalTo(m_result, lhs, rhs);
0123
0124     generic_product_impl<Lhs, Rhs, LhsShape, RhsShape, ProductTag>::evalTo(m_result, xpr.lhs(), xpr.rhs());
0125   }
0126
0127 protected:
0128   PlainObject m_result;
0129 };
0130
0131 // The following three shortcuts are enabled only if the scalar types match exactly.
0132 // TODO: we could enable them for different scalar types when the product is not vectorized.
0133
0134 // Dense = Product
0135 template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar>
0136 struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::assign_op<Scalar,Scalar>, Dense2Dense,
0137   typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct)>::type>
0138 {
0139   typedef Product<Lhs,Rhs,Options> SrcXprType;
0140   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0141   void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &)
0142   {
0143     Index dstRows = src.rows();
0144     Index dstCols = src.cols();
0145     if((dst.rows()!=dstRows) || (dst.cols()!=dstCols))
0146       dst.resize(dstRows, dstCols);
0147     // FIXME shall we handle nested_eval here?
0148     generic_product_impl<Lhs, Rhs>::evalTo(dst, src.lhs(), src.rhs());
0149   }
0150 };
0151
0152 // Dense += Product
0153 template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar>
0154 struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::add_assign_op<Scalar,Scalar>, Dense2Dense,
0155   typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct)>::type>
0156 {
0157   typedef Product<Lhs,Rhs,Options> SrcXprType;
0158   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0159   void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op<Scalar,Scalar> &)
0160   {
0161     eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
0162     // FIXME shall we handle nested_eval here?
0163     generic_product_impl<Lhs, Rhs>::addTo(dst, src.lhs(), src.rhs());
0164   }
0165 };
0166
0167 // Dense -= Product
0168 template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar>
0169 struct Assignment<DstXprType, Product<Lhs,Rhs,Options>, internal::sub_assign_op<Scalar,Scalar>, Dense2Dense,
0170   typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct)>::type>
0171 {
0172   typedef Product<Lhs,Rhs,Options> SrcXprType;
0173   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0174   void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op<Scalar,Scalar> &)
0175   {
0176     eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols());
0177     // FIXME shall we handle nested_eval here?
0178     generic_product_impl<Lhs, Rhs>::subTo(dst, src.lhs(), src.rhs());
0179   }
0180 };
0181
0182
0183 // Dense ?= scalar * Product
0184 // TODO we should apply that rule if that's really helpful
0185 // for instance, this is not good for inner products
0186 template< typename DstXprType, typename Lhs, typename Rhs, typename AssignFunc, typename Scalar, typename ScalarBis, typename Plain>
0187 struct Assignment<DstXprType, CwiseBinaryOp<internal::scalar_product_op<ScalarBis,Scalar>, const CwiseNullaryOp<internal::scalar_constant_op<ScalarBis>,Plain>,
0188                                            const Product<Lhs,Rhs,DefaultProduct> >, AssignFunc, Dense2Dense>
0189 {
0190   typedef CwiseBinaryOp<internal::scalar_product_op<ScalarBis,Scalar>,
0191                         const CwiseNullaryOp<internal::scalar_constant_op<ScalarBis>,Plain>,
0192                         const Product<Lhs,Rhs,DefaultProduct> > SrcXprType;
0193   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0194   void run(DstXprType &dst, const SrcXprType &src, const AssignFunc& func)
0195   {
0196     call_assignment_no_alias(dst, (src.lhs().functor().m_other * src.rhs().lhs())*src.rhs().rhs(), func);
0197   }
0198 };
0199
0200 //----------------------------------------
0201 // Catch "Dense ?= xpr + Product<>" expression to save one temporary
0202 // FIXME we could probably enable these rules for any product, i.e., not only Dense and DefaultProduct
0203
0204 template<typename OtherXpr, typename Lhs, typename Rhs>
0205 struct evaluator_assume_aliasing<CwiseBinaryOp<internal::scalar_sum_op<typename OtherXpr::Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, const OtherXpr,
0206                                                const Product<Lhs,Rhs,DefaultProduct> >, DenseShape > {
0207   static const bool value = true;
0208 };
0209
0210 template<typename OtherXpr, typename Lhs, typename Rhs>
0211 struct evaluator_assume_aliasing<CwiseBinaryOp<internal::scalar_difference_op<typename OtherXpr::Scalar,typename Product<Lhs,Rhs,DefaultProduct>::Scalar>, const OtherXpr,
0212                                                const Product<Lhs,Rhs,DefaultProduct> >, DenseShape > {
0213   static const bool value = true;
0214 };
0215
0216 template<typename DstXprType, typename OtherXpr, typename ProductType, typename Func1, typename Func2>
0217 struct assignment_from_xpr_op_product
0218 {
0219   template<typename SrcXprType, typename InitialFunc>
0220   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0221   void run(DstXprType &dst, const SrcXprType &src, const InitialFunc& /*func*/)
0222   {
0223     call_assignment_no_alias(dst, src.lhs(), Func1());
0224     call_assignment_no_alias(dst, src.rhs(), Func2());
0225   }
0226 };
0227
0228 #define EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(ASSIGN_OP,BINOP,ASSIGN_OP2) \
0229   template< typename DstXprType, typename OtherXpr, typename Lhs, typename Rhs, typename DstScalar, typename SrcScalar, typename OtherScalar,typename ProdScalar> \
0230   struct Assignment<DstXprType, CwiseBinaryOp<internal::BINOP<OtherScalar,ProdScalar>, const OtherXpr, \
0231                                             const Product<Lhs,Rhs,DefaultProduct> >, internal::ASSIGN_OP<DstScalar,SrcScalar>, Dense2Dense> \
0232     : assignment_from_xpr_op_product<DstXprType, OtherXpr, Product<Lhs,Rhs,DefaultProduct>, internal::ASSIGN_OP<DstScalar,OtherScalar>, internal::ASSIGN_OP2<DstScalar,ProdScalar> > \
0233   {}
0234
0235 EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(assign_op,    scalar_sum_op,add_assign_op);
0236 EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(add_assign_op,scalar_sum_op,add_assign_op);
0237 EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(sub_assign_op,scalar_sum_op,sub_assign_op);
0238
0239 EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(assign_op,    scalar_difference_op,sub_assign_op);
0240 EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(add_assign_op,scalar_difference_op,sub_assign_op);
0241 EIGEN_CATCH_ASSIGN_XPR_OP_PRODUCT(sub_assign_op,scalar_difference_op,add_assign_op);
0242
0243 //----------------------------------------
0244
0245 template<typename Lhs, typename Rhs>
0246 struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,InnerProduct>
0247 {
0248   template<typename Dst>
0249   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
0250   {
0251     dst.coeffRef(0,0) = (lhs.transpose().cwiseProduct(rhs)).sum();
0252   }
0253
0254   template<typename Dst>
0255   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
0256   {
0257     dst.coeffRef(0,0) += (lhs.transpose().cwiseProduct(rhs)).sum();
0258   }
0259
0260   template<typename Dst>
0261   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
0262   { dst.coeffRef(0,0) -= (lhs.transpose().cwiseProduct(rhs)).sum(); }
0263 };
0264
0265
0266 /***********************************************************************
0267 *  Implementation of outer dense * dense vector product
0268 ***********************************************************************/
0269
0270 // Column major result
0271 template<typename Dst, typename Lhs, typename Rhs, typename Func>
0272 void EIGEN_DEVICE_FUNC outer_product_selector_run(Dst& dst, const Lhs &lhs, const Rhs &rhs, const Func& func, const false_type&)
0273 {
0274   evaluator<Rhs> rhsEval(rhs);
0275   ei_declare_local_nested_eval(Lhs,lhs,Rhs::SizeAtCompileTime,actual_lhs);
0276   // FIXME if cols is large enough, then it might be useful to make sure that lhs is sequentially stored
0277   // FIXME not very good if rhs is real and lhs complex while alpha is real too
0278   const Index cols = dst.cols();
0279   for (Index j=0; j<cols; ++j)
0280     func(dst.col(j), rhsEval.coeff(Index(0),j) * actual_lhs);
0281 }
0282
0283 // Row major result
0284 template<typename Dst, typename Lhs, typename Rhs, typename Func>
0285 void EIGEN_DEVICE_FUNC outer_product_selector_run(Dst& dst, const Lhs &lhs, const Rhs &rhs, const Func& func, const true_type&)
0286 {
0287   evaluator<Lhs> lhsEval(lhs);
0288   ei_declare_local_nested_eval(Rhs,rhs,Lhs::SizeAtCompileTime,actual_rhs);
0289   // FIXME if rows is large enough, then it might be useful to make sure that rhs is sequentially stored
0290   // FIXME not very good if lhs is real and rhs complex while alpha is real too
0291   const Index rows = dst.rows();
0292   for (Index i=0; i<rows; ++i)
0293     func(dst.row(i), lhsEval.coeff(i,Index(0)) * actual_rhs);
0294 }
0295
0296 template<typename Lhs, typename Rhs>
0297 struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,OuterProduct>
0298 {
0299   template<typename T> struct is_row_major : internal::conditional<(int(T::Flags)&RowMajorBit), internal::true_type, internal::false_type>::type {};
0300   typedef typename Product<Lhs,Rhs>::Scalar Scalar;
0301
0302   // TODO it would be nice to be able to exploit our *_assign_op functors for that purpose
0303   struct set  { template<typename Dst, typename Src> EIGEN_DEVICE_FUNC void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived()  = src; } };
0304   struct add  { template<typename Dst, typename Src> EIGEN_DEVICE_FUNC void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() += src; } };
0305   struct sub  { template<typename Dst, typename Src> EIGEN_DEVICE_FUNC void operator()(const Dst& dst, const Src& src) const { dst.const_cast_derived() -= src; } };
0306   struct adds {
0307     Scalar m_scale;
0308     explicit adds(const Scalar& s) : m_scale(s) {}
0309     template<typename Dst, typename Src> void EIGEN_DEVICE_FUNC operator()(const Dst& dst, const Src& src) const {
0310       dst.const_cast_derived() += m_scale * src;
0311     }
0312   };
0313
0314   template<typename Dst>
0315   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
0316   {
0317     internal::outer_product_selector_run(dst, lhs, rhs, set(), is_row_major<Dst>());
0318   }
0319
0320   template<typename Dst>
0321   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
0322   {
0323     internal::outer_product_selector_run(dst, lhs, rhs, add(), is_row_major<Dst>());
0324   }
0325
0326   template<typename Dst>
0327   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
0328   {
0329     internal::outer_product_selector_run(dst, lhs, rhs, sub(), is_row_major<Dst>());
0330   }
0331
0332   template<typename Dst>
0333   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
0334   {
0335     internal::outer_product_selector_run(dst, lhs, rhs, adds(alpha), is_row_major<Dst>());
0336   }
0337
0338 };
0339
0340
0341 // This base class provides default implementations for evalTo, addTo, subTo, in terms of scaleAndAddTo
0342 template<typename Lhs, typename Rhs, typename Derived>
0343 struct generic_product_impl_base
0344 {
0345   typedef typename Product<Lhs,Rhs>::Scalar Scalar;
0346
0347   template<typename Dst>
0348   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
0349   { dst.setZero(); scaleAndAddTo(dst, lhs, rhs, Scalar(1)); }
0350
0351   template<typename Dst>
0352   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
0353   { scaleAndAddTo(dst,lhs, rhs, Scalar(1)); }
0354
0355   template<typename Dst>
0356   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
0357   { scaleAndAddTo(dst, lhs, rhs, Scalar(-1)); }
0358
0359   template<typename Dst>
0360   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
0361   { Derived::scaleAndAddTo(dst,lhs,rhs,alpha); }
0362
0363 };
0364
0365 template<typename Lhs, typename Rhs>
0366 struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,GemvProduct>
0367   : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,GemvProduct> >
0368 {
0369   typedef typename nested_eval<Lhs,1>::type LhsNested;
0370   typedef typename nested_eval<Rhs,1>::type RhsNested;
0371   typedef typename Product<Lhs,Rhs>::Scalar Scalar;
0372   enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight };
0373   typedef typename internal::remove_all<typename internal::conditional<int(Side)==OnTheRight,LhsNested,RhsNested>::type>::type MatrixType;
0374
0375   template<typename Dest>
0376   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
0377   {
0378     // Fallback to inner product if both the lhs and rhs is a runtime vector.
0379     if (lhs.rows() == 1 && rhs.cols() == 1) {
0380       dst.coeffRef(0,0) += alpha * lhs.row(0).conjugate().dot(rhs.col(0));
0381       return;
0382     }
0383     LhsNested actual_lhs(lhs);
0384     RhsNested actual_rhs(rhs);
0385     internal::gemv_dense_selector<Side,
0386                             (int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor,
0387                             bool(internal::blas_traits<MatrixType>::HasUsableDirectAccess)
0388                            >::run(actual_lhs, actual_rhs, dst, alpha);
0389   }
0390 };
0391
0392 template<typename Lhs, typename Rhs>
0393 struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,CoeffBasedProductMode>
0394 {
0395   typedef typename Product<Lhs,Rhs>::Scalar Scalar;
0396
0397   template<typename Dst>
0398   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
0399   {
0400     // Same as: dst.noalias() = lhs.lazyProduct(rhs);
0401     // but easier on the compiler side
0402     call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::assign_op<typename Dst::Scalar,Scalar>());
0403   }
0404
0405   template<typename Dst>
0406   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
0407   {
0408     // dst.noalias() += lhs.lazyProduct(rhs);
0409     call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::add_assign_op<typename Dst::Scalar,Scalar>());
0410   }
0411
0412   template<typename Dst>
0413   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
0414   {
0415     // dst.noalias() -= lhs.lazyProduct(rhs);
0416     call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::sub_assign_op<typename Dst::Scalar,Scalar>());
0417   }
0418
0419   // This is a special evaluation path called from generic_product_impl<...,GemmProduct> in file GeneralMatrixMatrix.h
0420   // This variant tries to extract scalar multiples from both the LHS and RHS and factor them out. For instance:
0421   //   dst {,+,-}= (s1*A)*(B*s2)
0422   // will be rewritten as:
0423   //   dst {,+,-}= (s1*s2) * (A.lazyProduct(B))
0424   // There are at least four benefits of doing so:
0425   //  1 - huge performance gain for heap-allocated matrix types as it save costly allocations.
0426   //  2 - it is faster than simply by-passing the heap allocation through stack allocation.
0427   //  3 - it makes this fallback consistent with the heavy GEMM routine.
0428   //  4 - it fully by-passes huge stack allocation attempts when multiplying huge fixed-size matrices.
0429   //      (see https://stackoverflow.com/questions/54738495)
0430   // For small fixed sizes matrices, howver, the gains are less obvious, it is sometimes x2 faster, but sometimes x3 slower,
0431   // and the behavior depends also a lot on the compiler... This is why this re-writting strategy is currently
0432   // enabled only when falling back from the main GEMM.
0433   template<typename Dst, typename Func>
0434   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0435   void eval_dynamic(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Func &func)
0436   {
0437     enum {
0438       HasScalarFactor = blas_traits<Lhs>::HasScalarFactor || blas_traits<Rhs>::HasScalarFactor,
0439       ConjLhs = blas_traits<Lhs>::NeedToConjugate,
0440       ConjRhs = blas_traits<Rhs>::NeedToConjugate
0441     };
0442     // FIXME: in c++11 this should be auto, and extractScalarFactor should also return auto
0443     //        this is important for real*complex_mat
0444     Scalar actualAlpha = combine_scalar_factors<Scalar>(lhs, rhs);
0445
0446     eval_dynamic_impl(dst,
0447                       blas_traits<Lhs>::extract(lhs).template conjugateIf<ConjLhs>(),
0448                       blas_traits<Rhs>::extract(rhs).template conjugateIf<ConjRhs>(),
0449                       func,
0450                       actualAlpha,
0451                       typename conditional<HasScalarFactor,true_type,false_type>::type());
0452   }
0453
0454 protected:
0455
0456   template<typename Dst, typename LhsT, typename RhsT, typename Func, typename Scalar>
0457   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0458   void eval_dynamic_impl(Dst& dst, const LhsT& lhs, const RhsT& rhs, const Func &func, const Scalar&  s /* == 1 */, false_type)
0459   {
0460     EIGEN_UNUSED_VARIABLE(s);
0461     eigen_internal_assert(s==Scalar(1));
0462     call_restricted_packet_assignment_no_alias(dst, lhs.lazyProduct(rhs), func);
0463   }
0464
0465   template<typename Dst, typename LhsT, typename RhsT, typename Func, typename Scalar>
0466   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0467   void eval_dynamic_impl(Dst& dst, const LhsT& lhs, const RhsT& rhs, const Func &func, const Scalar& s, true_type)
0468   {
0469     call_restricted_packet_assignment_no_alias(dst, s * lhs.lazyProduct(rhs), func);
0470   }
0471 };
0472
0473 // This specialization enforces the use of a coefficient-based evaluation strategy
0474 template<typename Lhs, typename Rhs>
0475 struct generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,LazyCoeffBasedProductMode>
0476   : generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,CoeffBasedProductMode> {};
0477
0478 // Case 2: Evaluate coeff by coeff
0479 //
0480 // This is mostly taken from CoeffBasedProduct.h
0481 // The main difference is that we add an extra argument to the etor_product_*_impl::run() function
0482 // for the inner dimension of the product, because evaluator object do not know their size.
0483
0484 template<int Traversal, int UnrollingIndex, typename Lhs, typename Rhs, typename RetScalar>
0485 struct etor_product_coeff_impl;
0486
0487 template<int StorageOrder, int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode>
0488 struct etor_product_packet_impl;
0489
0490 template<typename Lhs, typename Rhs, int ProductTag>
0491 struct product_evaluator<Product<Lhs, Rhs, LazyProduct>, ProductTag, DenseShape, DenseShape>
0492     : evaluator_base<Product<Lhs, Rhs, LazyProduct> >
0493 {
0494   typedef Product<Lhs, Rhs, LazyProduct> XprType;
0495   typedef typename XprType::Scalar Scalar;
0496   typedef typename XprType::CoeffReturnType CoeffReturnType;
0497
0498   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0499   explicit product_evaluator(const XprType& xpr)
0500     : m_lhs(xpr.lhs()),
0501       m_rhs(xpr.rhs()),
0502       m_lhsImpl(m_lhs),     // FIXME the creation of the evaluator objects should result in a no-op, but check that!
0503       m_rhsImpl(m_rhs),     //       Moreover, they are only useful for the packet path, so we could completely disable them when not needed,
0504                             //       or perhaps declare them on the fly on the packet method... We have experiment to check what's best.
0505       m_innerDim(xpr.lhs().cols())
0506   {
0507     EIGEN_INTERNAL_CHECK_COST_VALUE(NumTraits<Scalar>::MulCost);
0508     EIGEN_INTERNAL_CHECK_COST_VALUE(NumTraits<Scalar>::AddCost);
0509     EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
0510 #if 0
0511     std::cerr << "LhsOuterStrideBytes=  " << LhsOuterStrideBytes << "\n";
0512     std::cerr << "RhsOuterStrideBytes=  " << RhsOuterStrideBytes << "\n";
0513     std::cerr << "LhsAlignment=         " << LhsAlignment << "\n";
0514     std::cerr << "RhsAlignment=         " << RhsAlignment << "\n";
0515     std::cerr << "CanVectorizeLhs=      " << CanVectorizeLhs << "\n";
0516     std::cerr << "CanVectorizeRhs=      " << CanVectorizeRhs << "\n";
0517     std::cerr << "CanVectorizeInner=    " << CanVectorizeInner << "\n";
0518     std::cerr << "EvalToRowMajor=       " << EvalToRowMajor << "\n";
0519     std::cerr << "Alignment=            " << Alignment << "\n";
0520     std::cerr << "Flags=                " << Flags << "\n";
0521 #endif
0522   }
0523
0524   // Everything below here is taken from CoeffBasedProduct.h
0525
0526   typedef typename internal::nested_eval<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
0527   typedef typename internal::nested_eval<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
0528
0529   typedef typename internal::remove_all<LhsNested>::type LhsNestedCleaned;
0530   typedef typename internal::remove_all<RhsNested>::type RhsNestedCleaned;
0531
0532   typedef evaluator<LhsNestedCleaned> LhsEtorType;
0533   typedef evaluator<RhsNestedCleaned> RhsEtorType;
0534
0535   enum {
0536     RowsAtCompileTime = LhsNestedCleaned::RowsAtCompileTime,
0537     ColsAtCompileTime = RhsNestedCleaned::ColsAtCompileTime,
0538     InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(LhsNestedCleaned::ColsAtCompileTime, RhsNestedCleaned::RowsAtCompileTime),
0539     MaxRowsAtCompileTime = LhsNestedCleaned::MaxRowsAtCompileTime,
0540     MaxColsAtCompileTime = RhsNestedCleaned::MaxColsAtCompileTime
0541   };
0542
0543   typedef typename find_best_packet<Scalar,RowsAtCompileTime>::type LhsVecPacketType;
0544   typedef typename find_best_packet<Scalar,ColsAtCompileTime>::type RhsVecPacketType;
0545
0546   enum {
0547
0548     LhsCoeffReadCost = LhsEtorType::CoeffReadCost,
0549     RhsCoeffReadCost = RhsEtorType::CoeffReadCost,
0550     CoeffReadCost = InnerSize==0 ? NumTraits<Scalar>::ReadCost
0551                   : InnerSize == Dynamic ? HugeCost
0552                     : InnerSize * (NumTraits<Scalar>::MulCost + int(LhsCoeffReadCost) + int(RhsCoeffReadCost))
0553                     + (InnerSize - 1) * NumTraits<Scalar>::AddCost,
0554
0555     Unroll = CoeffReadCost <= EIGEN_UNROLLING_LIMIT,
0556
0557     LhsFlags = LhsEtorType::Flags,
0558     RhsFlags = RhsEtorType::Flags,
0559
0560     LhsRowMajor = LhsFlags & RowMajorBit,
0561     RhsRowMajor = RhsFlags & RowMajorBit,
0562
0563     LhsVecPacketSize = unpacket_traits<LhsVecPacketType>::size,
0564     RhsVecPacketSize = unpacket_traits<RhsVecPacketType>::size,
0565
0566     // Here, we don't care about alignment larger than the usable packet size.
0567     LhsAlignment = EIGEN_PLAIN_ENUM_MIN(LhsEtorType::Alignment,LhsVecPacketSize*int(sizeof(typename LhsNestedCleaned::Scalar))),
0568     RhsAlignment = EIGEN_PLAIN_ENUM_MIN(RhsEtorType::Alignment,RhsVecPacketSize*int(sizeof(typename RhsNestedCleaned::Scalar))),
0569
0570     SameType = is_same<typename LhsNestedCleaned::Scalar,typename RhsNestedCleaned::Scalar>::value,
0571
0572     CanVectorizeRhs = bool(RhsRowMajor) && (RhsFlags & PacketAccessBit) && (ColsAtCompileTime!=1),
0573     CanVectorizeLhs = (!LhsRowMajor) && (LhsFlags & PacketAccessBit) && (RowsAtCompileTime!=1),
0574
0575     EvalToRowMajor = (MaxRowsAtCompileTime==1&&MaxColsAtCompileTime!=1) ? 1
0576                     : (MaxColsAtCompileTime==1&&MaxRowsAtCompileTime!=1) ? 0
0577                     : (bool(RhsRowMajor) && !CanVectorizeLhs),
0578
0579     Flags = ((int(LhsFlags) | int(RhsFlags)) & HereditaryBits & ~RowMajorBit)
0580           | (EvalToRowMajor ? RowMajorBit : 0)
0581           // TODO enable vectorization for mixed types
0582           | (SameType && (CanVectorizeLhs || CanVectorizeRhs) ? PacketAccessBit : 0)
0583           | (XprType::IsVectorAtCompileTime ? LinearAccessBit : 0),
0584
0585     LhsOuterStrideBytes = int(LhsNestedCleaned::OuterStrideAtCompileTime) * int(sizeof(typename LhsNestedCleaned::Scalar)),
0586     RhsOuterStrideBytes = int(RhsNestedCleaned::OuterStrideAtCompileTime) * int(sizeof(typename RhsNestedCleaned::Scalar)),
0587
0588     Alignment = bool(CanVectorizeLhs) ? (LhsOuterStrideBytes<=0 || (int(LhsOuterStrideBytes) % EIGEN_PLAIN_ENUM_MAX(1,LhsAlignment))!=0 ? 0 : LhsAlignment)
0589               : bool(CanVectorizeRhs) ? (RhsOuterStrideBytes<=0 || (int(RhsOuterStrideBytes) % EIGEN_PLAIN_ENUM_MAX(1,RhsAlignment))!=0 ? 0 : RhsAlignment)
0590               : 0,
0591
0592     /* CanVectorizeInner deserves special explanation. It does not affect the product flags. It is not used outside
0593      * of Product. If the Product itself is not a packet-access expression, there is still a chance that the inner
0594      * loop of the product might be vectorized. This is the meaning of CanVectorizeInner. Since it doesn't affect
0595      * the Flags, it is safe to make this value depend on ActualPacketAccessBit, that doesn't affect the ABI.
0596      */
0597     CanVectorizeInner =    SameType
0598                         && LhsRowMajor
0599                         && (!RhsRowMajor)
0600                         && (int(LhsFlags) & int(RhsFlags) & ActualPacketAccessBit)
0601                         && (int(InnerSize) % packet_traits<Scalar>::size == 0)
0602   };
0603
0604   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const CoeffReturnType coeff(Index row, Index col) const
0605   {
0606     return (m_lhs.row(row).transpose().cwiseProduct( m_rhs.col(col) )).sum();
0607   }
0608
0609   /* Allow index-based non-packet access. It is impossible though to allow index-based packed access,
0610    * which is why we don't set the LinearAccessBit.
0611    * TODO: this seems possible when the result is a vector
0612    */
0613   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0614   const CoeffReturnType coeff(Index index) const
0615   {
0616     const Index row = (RowsAtCompileTime == 1 || MaxRowsAtCompileTime==1) ? 0 : index;
0617     const Index col = (RowsAtCompileTime == 1 || MaxRowsAtCompileTime==1) ? index : 0;
0618     return (m_lhs.row(row).transpose().cwiseProduct( m_rhs.col(col) )).sum();
0619   }
0620
0621   template<int LoadMode, typename PacketType>
0622   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0623   const PacketType packet(Index row, Index col) const
0624   {
0625     PacketType res;
0626     typedef etor_product_packet_impl<bool(int(Flags)&RowMajorBit) ? RowMajor : ColMajor,
0627                                      Unroll ? int(InnerSize) : Dynamic,
0628                                      LhsEtorType, RhsEtorType, PacketType, LoadMode> PacketImpl;
0629     PacketImpl::run(row, col, m_lhsImpl, m_rhsImpl, m_innerDim, res);
0630     return res;
0631   }
0632
0633   template<int LoadMode, typename PacketType>
0634   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0635   const PacketType packet(Index index) const
0636   {
0637     const Index row = (RowsAtCompileTime == 1 || MaxRowsAtCompileTime==1) ? 0 : index;
0638     const Index col = (RowsAtCompileTime == 1 || MaxRowsAtCompileTime==1) ? index : 0;
0639     return packet<LoadMode,PacketType>(row,col);
0640   }
0641
0642 protected:
0643   typename internal::add_const_on_value_type<LhsNested>::type m_lhs;
0644   typename internal::add_const_on_value_type<RhsNested>::type m_rhs;
0645
0646   LhsEtorType m_lhsImpl;
0647   RhsEtorType m_rhsImpl;
0648
0649   // TODO: Get rid of m_innerDim if known at compile time
0650   Index m_innerDim;
0651 };
0652
0653 template<typename Lhs, typename Rhs>
0654 struct product_evaluator<Product<Lhs, Rhs, DefaultProduct>, LazyCoeffBasedProductMode, DenseShape, DenseShape>
0655   : product_evaluator<Product<Lhs, Rhs, LazyProduct>, CoeffBasedProductMode, DenseShape, DenseShape>
0656 {
0657   typedef Product<Lhs, Rhs, DefaultProduct> XprType;
0658   typedef Product<Lhs, Rhs, LazyProduct> BaseProduct;
0659   typedef product_evaluator<BaseProduct, CoeffBasedProductMode, DenseShape, DenseShape> Base;
0660   enum {
0661     Flags = Base::Flags | EvalBeforeNestingBit
0662   };
0663   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
0664   explicit product_evaluator(const XprType& xpr)
0665     : Base(BaseProduct(xpr.lhs(),xpr.rhs()))
0666   {}
0667 };
0668
0669 /****************************************
0670 *** Coeff based product, Packet path  ***
0671 ****************************************/
0672
0673 template<int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode>
0674 struct etor_product_packet_impl<RowMajor, UnrollingIndex, Lhs, Rhs, Packet, LoadMode>
0675 {
0676   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet &res)
0677   {
0678     etor_product_packet_impl<RowMajor, UnrollingIndex-1, Lhs, Rhs, Packet, LoadMode>::run(row, col, lhs, rhs, innerDim, res);
0679     res =  pmadd(pset1<Packet>(lhs.coeff(row, Index(UnrollingIndex-1))), rhs.template packet<LoadMode,Packet>(Index(UnrollingIndex-1), col), res);
0680   }
0681 };
0682
0683 template<int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode>
0684 struct etor_product_packet_impl<ColMajor, UnrollingIndex, Lhs, Rhs, Packet, LoadMode>
0685 {
0686   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet &res)
0687   {
0688     etor_product_packet_impl<ColMajor, UnrollingIndex-1, Lhs, Rhs, Packet, LoadMode>::run(row, col, lhs, rhs, innerDim, res);
0689     res =  pmadd(lhs.template packet<LoadMode,Packet>(row, Index(UnrollingIndex-1)), pset1<Packet>(rhs.coeff(Index(UnrollingIndex-1), col)), res);
0690   }
0691 };
0692
0693 template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
0694 struct etor_product_packet_impl<RowMajor, 1, Lhs, Rhs, Packet, LoadMode>
0695 {
0696   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, Packet &res)
0697   {
0698     res = pmul(pset1<Packet>(lhs.coeff(row, Index(0))),rhs.template packet<LoadMode,Packet>(Index(0), col));
0699   }
0700 };
0701
0702 template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
0703 struct etor_product_packet_impl<ColMajor, 1, Lhs, Rhs, Packet, LoadMode>
0704 {
0705   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index /*innerDim*/, Packet &res)
0706   {
0707     res = pmul(lhs.template packet<LoadMode,Packet>(row, Index(0)), pset1<Packet>(rhs.coeff(Index(0), col)));
0708   }
0709 };
0710
0711 template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
0712 struct etor_product_packet_impl<RowMajor, 0, Lhs, Rhs, Packet, LoadMode>
0713 {
0714   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index /*row*/, Index /*col*/, const Lhs& /*lhs*/, const Rhs& /*rhs*/, Index /*innerDim*/, Packet &res)
0715   {
0716     res = pset1<Packet>(typename unpacket_traits<Packet>::type(0));
0717   }
0718 };
0719
0720 template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
0721 struct etor_product_packet_impl<ColMajor, 0, Lhs, Rhs, Packet, LoadMode>
0722 {
0723   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index /*row*/, Index /*col*/, const Lhs& /*lhs*/, const Rhs& /*rhs*/, Index /*innerDim*/, Packet &res)
0724   {
0725     res = pset1<Packet>(typename unpacket_traits<Packet>::type(0));
0726   }
0727 };
0728
0729 template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
0730 struct etor_product_packet_impl<RowMajor, Dynamic, Lhs, Rhs, Packet, LoadMode>
0731 {
0732   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet& res)
0733   {
0734     res = pset1<Packet>(typename unpacket_traits<Packet>::type(0));
0735     for(Index i = 0; i < innerDim; ++i)
0736       res =  pmadd(pset1<Packet>(lhs.coeff(row, i)), rhs.template packet<LoadMode,Packet>(i, col), res);
0737   }
0738 };
0739
0740 template<typename Lhs, typename Rhs, typename Packet, int LoadMode>
0741 struct etor_product_packet_impl<ColMajor, Dynamic, Lhs, Rhs, Packet, LoadMode>
0742 {
0743   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, Index innerDim, Packet& res)
0744   {
0745     res = pset1<Packet>(typename unpacket_traits<Packet>::type(0));
0746     for(Index i = 0; i < innerDim; ++i)
0747       res =  pmadd(lhs.template packet<LoadMode,Packet>(row, i), pset1<Packet>(rhs.coeff(i, col)), res);
0748   }
0749 };
0750
0751
0752 /***************************************************************************
0753 * Triangular products
0754 ***************************************************************************/
0755 template<int Mode, bool LhsIsTriangular,
0756          typename Lhs, bool LhsIsVector,
0757          typename Rhs, bool RhsIsVector>
0758 struct triangular_product_impl;
0759
0760 template<typename Lhs, typename Rhs, int ProductTag>
0761 struct generic_product_impl<Lhs,Rhs,TriangularShape,DenseShape,ProductTag>
0762   : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,TriangularShape,DenseShape,ProductTag> >
0763 {
0764   typedef typename Product<Lhs,Rhs>::Scalar Scalar;
0765
0766   template<typename Dest>
0767   static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
0768   {
0769     triangular_product_impl<Lhs::Mode,true,typename Lhs::MatrixType,false,Rhs, Rhs::ColsAtCompileTime==1>
0770         ::run(dst, lhs.nestedExpression(), rhs, alpha);
0771   }
0772 };
0773
0774 template<typename Lhs, typename Rhs, int ProductTag>
0775 struct generic_product_impl<Lhs,Rhs,DenseShape,TriangularShape,ProductTag>
0776 : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,TriangularShape,ProductTag> >
0777 {
0778   typedef typename Product<Lhs,Rhs>::Scalar Scalar;
0779
0780   template<typename Dest>
0781   static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
0782   {
0783     triangular_product_impl<Rhs::Mode,false,Lhs,Lhs::RowsAtCompileTime==1, typename Rhs::MatrixType, false>::run(dst, lhs, rhs.nestedExpression(), alpha);
0784   }
0785 };
0786
0787
0788 /***************************************************************************
0789 * SelfAdjoint products
0790 ***************************************************************************/
0791 template <typename Lhs, int LhsMode, bool LhsIsVector,
0792           typename Rhs, int RhsMode, bool RhsIsVector>
0793 struct selfadjoint_product_impl;
0794
0795 template<typename Lhs, typename Rhs, int ProductTag>
0796 struct generic_product_impl<Lhs,Rhs,SelfAdjointShape,DenseShape,ProductTag>
0797   : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,SelfAdjointShape,DenseShape,ProductTag> >
0798 {
0799   typedef typename Product<Lhs,Rhs>::Scalar Scalar;
0800
0801   template<typename Dest>
0802   static EIGEN_DEVICE_FUNC
0803   void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
0804   {
0805     selfadjoint_product_impl<typename Lhs::MatrixType,Lhs::Mode,false,Rhs,0,Rhs::IsVectorAtCompileTime>::run(dst, lhs.nestedExpression(), rhs, alpha);
0806   }
0807 };
0808
0809 template<typename Lhs, typename Rhs, int ProductTag>
0810 struct generic_product_impl<Lhs,Rhs,DenseShape,SelfAdjointShape,ProductTag>
0811 : generic_product_impl_base<Lhs,Rhs,generic_product_impl<Lhs,Rhs,DenseShape,SelfAdjointShape,ProductTag> >
0812 {
0813   typedef typename Product<Lhs,Rhs>::Scalar Scalar;
0814
0815   template<typename Dest>
0816   static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
0817   {
0818     selfadjoint_product_impl<Lhs,0,Lhs::IsVectorAtCompileTime,typename Rhs::MatrixType,Rhs::Mode,false>::run(dst, lhs, rhs.nestedExpression(), alpha);
0819   }
0820 };
0821
0822
0823 /***************************************************************************
0824 * Diagonal products
0825 ***************************************************************************/
0826
0827 template<typename MatrixType, typename DiagonalType, typename Derived, int ProductOrder>
0828 struct diagonal_product_evaluator_base
0829   : evaluator_base<Derived>
0830 {
0831    typedef typename ScalarBinaryOpTraits<typename MatrixType::Scalar, typename DiagonalType::Scalar>::ReturnType Scalar;
0832 public:
0833   enum {
0834     CoeffReadCost = int(NumTraits<Scalar>::MulCost) + int(evaluator<MatrixType>::CoeffReadCost) + int(evaluator<DiagonalType>::CoeffReadCost),
0835
0836     MatrixFlags = evaluator<MatrixType>::Flags,
0837     DiagFlags = evaluator<DiagonalType>::Flags,
0838
0839     _StorageOrder = (Derived::MaxRowsAtCompileTime==1 && Derived::MaxColsAtCompileTime!=1) ? RowMajor
0840                   : (Derived::MaxColsAtCompileTime==1 && Derived::MaxRowsAtCompileTime!=1) ? ColMajor
0841                   : MatrixFlags & RowMajorBit ? RowMajor : ColMajor,
0842     _SameStorageOrder = _StorageOrder == (MatrixFlags & RowMajorBit ? RowMajor : ColMajor),
0843
0844     _ScalarAccessOnDiag =  !((int(_StorageOrder) == ColMajor && int(ProductOrder) == OnTheLeft)
0845                            ||(int(_StorageOrder) == RowMajor && int(ProductOrder) == OnTheRight)),
0846     _SameTypes = is_same<typename MatrixType::Scalar, typename DiagonalType::Scalar>::value,
0847     // FIXME currently we need same types, but in the future the next rule should be the one
0848     //_Vectorizable = bool(int(MatrixFlags)&PacketAccessBit) && ((!_PacketOnDiag) || (_SameTypes && bool(int(DiagFlags)&PacketAccessBit))),
0849     _Vectorizable =   bool(int(MatrixFlags)&PacketAccessBit)
0850                   &&  _SameTypes
0851                   && (_SameStorageOrder || (MatrixFlags&LinearAccessBit)==LinearAccessBit)
0852                   && (_ScalarAccessOnDiag || (bool(int(DiagFlags)&PacketAccessBit))),
0853     _LinearAccessMask = (MatrixType::RowsAtCompileTime==1 || MatrixType::ColsAtCompileTime==1) ? LinearAccessBit : 0,
0854     Flags = ((HereditaryBits|_LinearAccessMask) & (unsigned int)(MatrixFlags)) | (_Vectorizable ? PacketAccessBit : 0),
0855     Alignment = evaluator<MatrixType>::Alignment,
0856
0857     AsScalarProduct =     (DiagonalType::SizeAtCompileTime==1)
0858                       ||  (DiagonalType::SizeAtCompileTime==Dynamic && MatrixType::RowsAtCompileTime==1 && ProductOrder==OnTheLeft)
0859                       ||  (DiagonalType::SizeAtCompileTime==Dynamic && MatrixType::ColsAtCompileTime==1 && ProductOrder==OnTheRight)
0860   };
0861
0862   EIGEN_DEVICE_FUNC diagonal_product_evaluator_base(const MatrixType &mat, const DiagonalType &diag)
0863     : m_diagImpl(diag), m_matImpl(mat)
0864   {
0865     EIGEN_INTERNAL_CHECK_COST_VALUE(NumTraits<Scalar>::MulCost);
0866     EIGEN_INTERNAL_CHECK_COST_VALUE(CoeffReadCost);
0867   }
0868
0869   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar coeff(Index idx) const
0870   {
0871     if(AsScalarProduct)
0872       return m_diagImpl.coeff(0) * m_matImpl.coeff(idx);
0873     else
0874       return m_diagImpl.coeff(idx) * m_matImpl.coeff(idx);
0875   }
0876
0877 protected:
0878   template<int LoadMode,typename PacketType>
0879   EIGEN_STRONG_INLINE PacketType packet_impl(Index row, Index col, Index id, internal::true_type) const
0880   {
0881     return internal::pmul(m_matImpl.template packet<LoadMode,PacketType>(row, col),
0882                           internal::pset1<PacketType>(m_diagImpl.coeff(id)));
0883   }
0884
0885   template<int LoadMode,typename PacketType>
0886   EIGEN_STRONG_INLINE PacketType packet_impl(Index row, Index col, Index id, internal::false_type) const
0887   {
0888     enum {
0889       InnerSize = (MatrixType::Flags & RowMajorBit) ? MatrixType::ColsAtCompileTime : MatrixType::RowsAtCompileTime,
0890       DiagonalPacketLoadMode = EIGEN_PLAIN_ENUM_MIN(LoadMode,((InnerSize%16) == 0) ? int(Aligned16) : int(evaluator<DiagonalType>::Alignment)) // FIXME hardcoded 16!!
0891     };
0892     return internal::pmul(m_matImpl.template packet<LoadMode,PacketType>(row, col),
0893                           m_diagImpl.template packet<DiagonalPacketLoadMode,PacketType>(id));
0894   }
0895
0896   evaluator<DiagonalType> m_diagImpl;
0897   evaluator<MatrixType>   m_matImpl;
0898 };
0899
0900 // diagonal * dense
0901 template<typename Lhs, typename Rhs, int ProductKind, int ProductTag>
0902 struct product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, DiagonalShape, DenseShape>
0903   : diagonal_product_evaluator_base<Rhs, typename Lhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheLeft>
0904 {
0905   typedef diagonal_product_evaluator_base<Rhs, typename Lhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheLeft> Base;
0906   using Base::m_diagImpl;
0907   using Base::m_matImpl;
0908   using Base::coeff;
0909   typedef typename Base::Scalar Scalar;
0910
0911   typedef Product<Lhs, Rhs, ProductKind> XprType;
0912   typedef typename XprType::PlainObject PlainObject;
0913   typedef typename Lhs::DiagonalVectorType DiagonalType;
0914
0915
0916   enum { StorageOrder = Base::_StorageOrder };
0917
0918   EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr)
0919     : Base(xpr.rhs(), xpr.lhs().diagonal())
0920   {
0921   }
0922
0923   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar coeff(Index row, Index col) const
0924   {
0925     return m_diagImpl.coeff(row) * m_matImpl.coeff(row, col);
0926   }
0927
0928 #ifndef EIGEN_GPUCC
0929   template<int LoadMode,typename PacketType>
0930   EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const
0931   {
0932     // FIXME: NVCC used to complain about the template keyword, but we have to check whether this is still the case.
0933     // See also similar calls below.
0934     return this->template packet_impl<LoadMode,PacketType>(row,col, row,
0935                                  typename internal::conditional<int(StorageOrder)==RowMajor, internal::true_type, internal::false_type>::type());
0936   }
0937
0938   template<int LoadMode,typename PacketType>
0939   EIGEN_STRONG_INLINE PacketType packet(Index idx) const
0940   {
0941     return packet<LoadMode,PacketType>(int(StorageOrder)==ColMajor?idx:0,int(StorageOrder)==ColMajor?0:idx);
0942   }
0943 #endif
0944 };
0945
0946 // dense * diagonal
0947 template<typename Lhs, typename Rhs, int ProductKind, int ProductTag>
0948 struct product_evaluator<Product<Lhs, Rhs, ProductKind>, ProductTag, DenseShape, DiagonalShape>
0949   : diagonal_product_evaluator_base<Lhs, typename Rhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheRight>
0950 {
0951   typedef diagonal_product_evaluator_base<Lhs, typename Rhs::DiagonalVectorType, Product<Lhs, Rhs, LazyProduct>, OnTheRight> Base;
0952   using Base::m_diagImpl;
0953   using Base::m_matImpl;
0954   using Base::coeff;
0955   typedef typename Base::Scalar Scalar;
0956
0957   typedef Product<Lhs, Rhs, ProductKind> XprType;
0958   typedef typename XprType::PlainObject PlainObject;
0959
0960   enum { StorageOrder = Base::_StorageOrder };
0961
0962   EIGEN_DEVICE_FUNC explicit product_evaluator(const XprType& xpr)
0963     : Base(xpr.lhs(), xpr.rhs().diagonal())
0964   {
0965   }
0966
0967   EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar coeff(Index row, Index col) const
0968   {
0969     return m_matImpl.coeff(row, col) * m_diagImpl.coeff(col);
0970   }
0971
0972 #ifndef EIGEN_GPUCC
0973   template<int LoadMode,typename PacketType>
0974   EIGEN_STRONG_INLINE PacketType packet(Index row, Index col) const
0975   {
0976     return this->template packet_impl<LoadMode,PacketType>(row,col, col,
0977                                  typename internal::conditional<int(StorageOrder)==ColMajor, internal::true_type, internal::false_type>::type());
0978   }
0979
0980   template<int LoadMode,typename PacketType>
0981   EIGEN_STRONG_INLINE PacketType packet(Index idx) const
0982   {
0983     return packet<LoadMode,PacketType>(int(StorageOrder)==ColMajor?idx:0,int(StorageOrder)==ColMajor?0:idx);
0984   }
0985 #endif
0986 };
0987
0988 /***************************************************************************
0989 * Products with permutation matrices
0990 ***************************************************************************/
0991
0992 /** \internal
0993   * \class permutation_matrix_product
0994   * Internal helper class implementing the product between a permutation matrix and a matrix.
0995   * This class is specialized for DenseShape below and for SparseShape in SparseCore/SparsePermutation.h
0996   */
0997 template<typename ExpressionType, int Side, bool Transposed, typename ExpressionShape>
0998 struct permutation_matrix_product;
0999
1000 template<typename ExpressionType, int Side, bool Transposed>
1001 struct permutation_matrix_product<ExpressionType, Side, Transposed, DenseShape>
1002 {
1003     typedef typename nested_eval<ExpressionType, 1>::type MatrixType;
1004     typedef typename remove_all<MatrixType>::type MatrixTypeCleaned;
1005
1006     template<typename Dest, typename PermutationType>
1007     static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Dest& dst, const PermutationType& perm, const ExpressionType& xpr)
1008     {
1009       MatrixType mat(xpr);
1010       const Index n = Side==OnTheLeft ? mat.rows() : mat.cols();
1011       // FIXME we need an is_same for expression that is not sensitive to constness. For instance
1012       // is_same_xpr<Block<const Matrix>, Block<Matrix> >::value should be true.
1013       //if(is_same<MatrixTypeCleaned,Dest>::value && extract_data(dst) == extract_data(mat))
1014       if(is_same_dense(dst, mat))
1015       {
1016         // apply the permutation inplace
1017         Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(perm.size());
1018         mask.fill(false);
1019         Index r = 0;
1020         while(r < perm.size())
1021         {
1022           // search for the next seed
1023           while(r<perm.size() && mask[r]) r++;
1024           if(r>=perm.size())
1025             break;
1026           // we got one, let's follow it until we are back to the seed
1027           Index k0 = r++;
1028           Index kPrev = k0;
1029           mask.coeffRef(k0) = true;
1030           for(Index k=perm.indices().coeff(k0); k!=k0; k=perm.indices().coeff(k))
1031           {
1032                   Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>(dst, k)
1033             .swap(Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
1034                        (dst,((Side==OnTheLeft) ^ Transposed) ? k0 : kPrev));
1035
1036             mask.coeffRef(k) = true;
1037             kPrev = k;
1038           }
1039         }
1040       }
1041       else
1042       {
1043         for(Index i = 0; i < n; ++i)
1044         {
1045           Block<Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime>
1046                (dst, ((Side==OnTheLeft) ^ Transposed) ? perm.indices().coeff(i) : i)
1047
1048           =
1049
1050           Block<const MatrixTypeCleaned,Side==OnTheLeft ? 1 : MatrixTypeCleaned::RowsAtCompileTime,Side==OnTheRight ? 1 : MatrixTypeCleaned::ColsAtCompileTime>
1051                (mat, ((Side==OnTheRight) ^ Transposed) ? perm.indices().coeff(i) : i);
1052         }
1053       }
1054     }
1055 };
1056
1057 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1058 struct generic_product_impl<Lhs, Rhs, PermutationShape, MatrixShape, ProductTag>
1059 {
1060   template<typename Dest>
1061   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
1062   {
1063     permutation_matrix_product<Rhs, OnTheLeft, false, MatrixShape>::run(dst, lhs, rhs);
1064   }
1065 };
1066
1067 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1068 struct generic_product_impl<Lhs, Rhs, MatrixShape, PermutationShape, ProductTag>
1069 {
1070   template<typename Dest>
1071   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
1072   {
1073     permutation_matrix_product<Lhs, OnTheRight, false, MatrixShape>::run(dst, rhs, lhs);
1074   }
1075 };
1076
1077 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1078 struct generic_product_impl<Inverse<Lhs>, Rhs, PermutationShape, MatrixShape, ProductTag>
1079 {
1080   template<typename Dest>
1081   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Inverse<Lhs>& lhs, const Rhs& rhs)
1082   {
1083     permutation_matrix_product<Rhs, OnTheLeft, true, MatrixShape>::run(dst, lhs.nestedExpression(), rhs);
1084   }
1085 };
1086
1087 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1088 struct generic_product_impl<Lhs, Inverse<Rhs>, MatrixShape, PermutationShape, ProductTag>
1089 {
1090   template<typename Dest>
1091   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Inverse<Rhs>& rhs)
1092   {
1093     permutation_matrix_product<Lhs, OnTheRight, true, MatrixShape>::run(dst, rhs.nestedExpression(), lhs);
1094   }
1095 };
1096
1097
1098 /***************************************************************************
1099 * Products with transpositions matrices
1100 ***************************************************************************/
1101
1102 // FIXME could we unify Transpositions and Permutation into a single "shape"??
1103
1104 /** \internal
1105   * \class transposition_matrix_product
1106   * Internal helper class implementing the product between a permutation matrix and a matrix.
1107   */
1108 template<typename ExpressionType, int Side, bool Transposed, typename ExpressionShape>
1109 struct transposition_matrix_product
1110 {
1111   typedef typename nested_eval<ExpressionType, 1>::type MatrixType;
1112   typedef typename remove_all<MatrixType>::type MatrixTypeCleaned;
1113
1114   template<typename Dest, typename TranspositionType>
1115   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void run(Dest& dst, const TranspositionType& tr, const ExpressionType& xpr)
1116   {
1117     MatrixType mat(xpr);
1118     typedef typename TranspositionType::StorageIndex StorageIndex;
1119     const Index size = tr.size();
1120     StorageIndex j = 0;
1121
1122     if(!is_same_dense(dst,mat))
1123       dst = mat;
1124
1125     for(Index k=(Transposed?size-1:0) ; Transposed?k>=0:k<size ; Transposed?--k:++k)
1126       if(Index(j=tr.coeff(k))!=k)
1127       {
1128         if(Side==OnTheLeft)        dst.row(k).swap(dst.row(j));
1129         else if(Side==OnTheRight)  dst.col(k).swap(dst.col(j));
1130       }
1131   }
1132 };
1133
1134 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1135 struct generic_product_impl<Lhs, Rhs, TranspositionsShape, MatrixShape, ProductTag>
1136 {
1137   template<typename Dest>
1138   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
1139   {
1140     transposition_matrix_product<Rhs, OnTheLeft, false, MatrixShape>::run(dst, lhs, rhs);
1141   }
1142 };
1143
1144 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1145 struct generic_product_impl<Lhs, Rhs, MatrixShape, TranspositionsShape, ProductTag>
1146 {
1147   template<typename Dest>
1148   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs)
1149   {
1150     transposition_matrix_product<Lhs, OnTheRight, false, MatrixShape>::run(dst, rhs, lhs);
1151   }
1152 };
1153
1154
1155 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1156 struct generic_product_impl<Transpose<Lhs>, Rhs, TranspositionsShape, MatrixShape, ProductTag>
1157 {
1158   template<typename Dest>
1159   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Transpose<Lhs>& lhs, const Rhs& rhs)
1160   {
1161     transposition_matrix_product<Rhs, OnTheLeft, true, MatrixShape>::run(dst, lhs.nestedExpression(), rhs);
1162   }
1163 };
1164
1165 template<typename Lhs, typename Rhs, int ProductTag, typename MatrixShape>
1166 struct generic_product_impl<Lhs, Transpose<Rhs>, MatrixShape, TranspositionsShape, ProductTag>
1167 {
1168   template<typename Dest>
1169   static EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void evalTo(Dest& dst, const Lhs& lhs, const Transpose<Rhs>& rhs)
1170   {
1171     transposition_matrix_product<Lhs, OnTheRight, true, MatrixShape>::run(dst, rhs.nestedExpression(), lhs);
1172   }
1173 };
1174
1175 } // end namespace internal
1176
1177 } // end namespace Eigen
1178
1179 #endif // EIGEN_PRODUCT_EVALUATORS_H