Refactor dense product evaluators

4 files changed, 185 insertions, 89 deletions

diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index 5d3789be7..52586e5c0 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -12,8 +12,7 @@
 
 namespace Eigen {
 
-template<typename Lhs, typename Rhs> class Product;
-template<typename Lhs, typename Rhs, typename StorageKind> class ProductImpl;
+template<typename Lhs, typename Rhs, int Option, int ProductTag, typename StorageKind> class ProductImpl;
 
 /** \class Product
   * \ingroup Core_Module
@@ -24,13 +23,17 @@ template<typename Lhs, typename Rhs, typename StorageKind> class ProductImpl;
   * \param Rhs the type of the right-hand side expression
   *
   * This class represents an expression of the product of two arbitrary matrices.
+  * 
+  * The other template parameters are:
+  * \tparam Option     can be DefaultProduct or LazyProduct
+  * \tparam ProductTag can be InnerProduct, OuterProduct, GemvProduct, GemmProduct. It is used to ease expression manipulations.
   *
   */
 
 // Use ProductReturnType to get correct traits, in particular vectorization flags
 namespace internal {
-template<typename Lhs, typename Rhs>
-struct traits<Product<Lhs, Rhs> >
+template<typename Lhs, typename Rhs, int Option, int ProductTag>
+struct traits<Product<Lhs, Rhs, Option, ProductTag> >
   : traits<typename ProductReturnType<Lhs, Rhs>::Type>
 { 
   // We want A+B*C to be of type Product<Matrix, Sum> and not Product<Matrix, Matrix>
@@ -42,14 +45,15 @@ struct traits<Product<Lhs, Rhs> >
 } // end namespace internal
 
 
-template<typename Lhs, typename Rhs>
-class Product : public ProductImpl<Lhs,Rhs,typename internal::promote_storage_type<typename internal::traits<Lhs>::StorageKind,
-                                                                            typename internal::traits<Rhs>::StorageKind>::ret>
+template<typename Lhs, typename Rhs, int Option, int ProductTag>
+class Product : public ProductImpl<Lhs,Rhs,Option,ProductTag,
+                                   typename internal::promote_storage_type<typename internal::traits<Lhs>::StorageKind,
+                                                                           typename internal::traits<Rhs>::StorageKind>::ret>
 {
   public:
     
     typedef typename ProductImpl<
-        Lhs, Rhs,
+        Lhs, Rhs, Option, ProductTag,
         typename internal::promote_storage_type<typename Lhs::StorageKind,
                                                 typename Rhs::StorageKind>::ret>::Base Base;
     EIGEN_GENERIC_PUBLIC_INTERFACE(Product)
@@ -78,13 +82,13 @@ class Product : public ProductImpl<Lhs,Rhs,typename internal::promote_storage_ty
     RhsNested m_rhs;
 };
 
-template<typename Lhs, typename Rhs>
-class ProductImpl<Lhs,Rhs,Dense> : public internal::dense_xpr_base<Product<Lhs,Rhs> >::type
+template<typename Lhs, typename Rhs, int Option, int ProductTag>
+class ProductImpl<Lhs,Rhs,Option,ProductTag,Dense> : public internal::dense_xpr_base<Product<Lhs,Rhs,Option,ProductTag> >::type
 {
     typedef Product<Lhs, Rhs> Derived;
   public:
 
-    typedef typename internal::dense_xpr_base<Product<Lhs, Rhs> >::type Base;
+    typedef typename internal::dense_xpr_base<Product<Lhs, Rhs, Option, ProductTag> >::type Base;
     EIGEN_DENSE_PUBLIC_INTERFACE(Derived)
 };
 
@@ -102,6 +106,15 @@ prod(const Lhs& lhs, const Rhs& rhs)
   return Product<Lhs,Rhs>(lhs,rhs);
 }
 
+/** \internal used to test the evaluator only
+  */
+template<typename Lhs,typename Rhs>
+const Product<Lhs,Rhs,LazyProduct>
+lazyprod(const Lhs& lhs, const Rhs& rhs)
+{
+  return Product<Lhs,Rhs,LazyProduct>(lhs,rhs);
+}
+
 } // end namespace Eigen
 
 #endif // EIGEN_PRODUCT_H
diff --git a/Eigen/src/Core/ProductEvaluators.h b/Eigen/src/Core/ProductEvaluators.h
index 855914f2e..42dd3c7ac 100644
--- a/Eigen/src/Core/ProductEvaluators.h
+++ b/Eigen/src/Core/ProductEvaluators.h
@@ -17,94 +17,172 @@ namespace Eigen {
   
 namespace internal {
   
-// We can evaluate the product either all at once, like GeneralProduct and its evalTo() function, or
-// traverse the matrix coefficient by coefficient, like CoeffBasedProduct.  Use the existing logic
-// in ProductReturnType to decide.
+  
+// Helper class to perform a dense product with the destination at hand.
+// Depending on the sizes of the factors, there are different evaluation strategies
+// as controlled by internal::product_type.
+template<typename Lhs, typename Rhs, int ProductType = internal::product_type<Lhs,Rhs>::value>
+struct dense_product_impl;
 
-template<typename XprType, typename ProductType>
-struct product_evaluator_dispatcher;
 
-template<typename Lhs, typename Rhs>
-struct evaluator_impl<Product<Lhs, Rhs> >
-  : product_evaluator_dispatcher<Product<Lhs, Rhs>, typename ProductReturnType<Lhs, Rhs>::Type> 
+// The evaluator for default dense products creates a temporary and call dense_product_impl
+template<typename Lhs, typename Rhs, int ProductTag>
+struct evaluator_impl<Product<Lhs, Rhs, DefaultProduct, ProductTag> > 
+  : public evaluator<typename Product<Lhs, Rhs, DefaultProduct, ProductTag>::PlainObject>::type
 {
-  typedef Product<Lhs, Rhs> XprType;
-  typedef product_evaluator_dispatcher<XprType, typename ProductReturnType<Lhs, Rhs>::Type> Base;
+  typedef Product<Lhs, Rhs, DefaultProduct, ProductTag> XprType;
+  typedef typename XprType::PlainObject PlainObject;
+  typedef typename evaluator<PlainObject>::type Base;
 
-  evaluator_impl(const XprType& xpr) : Base(xpr) 
-  { }
+  evaluator_impl(const XprType& xpr)
+    : m_result(xpr.rows(), xpr.cols())
+  {
+    ::new (static_cast<Base*>(this)) Base(m_result);
+    dense_product_impl<Lhs, Rhs>::evalTo(m_result, xpr.lhs(), xpr.rhs());
+  }
+  
+protected:  
+  PlainObject m_result;
 };
 
-template<typename XprType, typename ProductType>
-struct product_evaluator_traits_dispatcher;
 
 template<typename Lhs, typename Rhs>
-struct evaluator_traits<Product<Lhs, Rhs> >
-  : product_evaluator_traits_dispatcher<Product<Lhs, Rhs>, typename ProductReturnType<Lhs, Rhs>::Type> 
-{ 
-  static const int AssumeAliasing = 1;
+struct dense_product_impl<Lhs,Rhs,InnerProduct>
+{
+  template<typename Dst>
+  static inline void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
+  {
+    dst.coeffRef(0,0) = (lhs.transpose().cwiseProduct(rhs)).sum();
+  }
+  
+  template<typename Dst>
+  static inline void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
+  {
+    dst.coeffRef(0,0) += (lhs.transpose().cwiseProduct(rhs)).sum();
+  }
+  
+  template<typename Dst>
+  static void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
+  { dst.coeffRef(0,0) -= (lhs.transpose().cwiseProduct(rhs)).sum(); }
 };
 
-// Case 1: Evaluate all at once
-//
-// We can view the GeneralProduct class as a part of the product evaluator. 
-// Four sub-cases: InnerProduct, OuterProduct, GemmProduct and GemvProduct.
-// InnerProduct is special because GeneralProduct does not have an evalTo() method in this case.
+
 
 template<typename Lhs, typename Rhs>
-struct product_evaluator_traits_dispatcher<Product<Lhs, Rhs>, GeneralProduct<Lhs, Rhs, InnerProduct> > 
+struct dense_product_impl<Lhs,Rhs,OuterProduct>
+{
+  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
+  
+  template<typename Dst>
+  static inline void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
+  {
+    // TODO bypass GeneralProduct class
+    GeneralProduct<Lhs, Rhs, OuterProduct>(lhs,rhs).evalTo(dst);
+  }
+  
+  template<typename Dst>
+  static inline void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
+  {
+    // TODO bypass GeneralProduct class
+    GeneralProduct<Lhs, Rhs, OuterProduct>(lhs,rhs).addTo(dst);
+  }
+  
+  template<typename Dst>
+  static inline void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
+  {
+    // TODO bypass GeneralProduct class
+    GeneralProduct<Lhs, Rhs, OuterProduct>(lhs,rhs).subTo(dst);
+  }
+  
+  template<typename Dst>
+  static inline void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
+  {
+    // TODO bypass GeneralProduct class
+    GeneralProduct<Lhs, Rhs, OuterProduct>(lhs,rhs).scaleAndAddTo(dst, alpha);
+  }
+  
+};
+
+
+// This base class provides default implementations for evalTo, addTo, subTo, in terms of scaleAndAddTo
+template<typename Lhs, typename Rhs, typename Derived>
+struct dense_product_impl_base
 {
-  static const int HasEvalTo = 0;
+  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
+  
+  template<typename Dst>
+  static void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
+  { dst.setZero(); scaleAndAddTo(dst, lhs, rhs, Scalar(1)); }
+
+  template<typename Dst>
+  static void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
+  { scaleAndAddTo(dst,lhs, rhs, Scalar(1)); }
+
+  template<typename Dst>
+  static void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
+  { scaleAndAddTo(dst, lhs, rhs, Scalar(-1)); }
+  
+  template<typename Dst>
+  static void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
+  { Derived::scaleAndAddTo(dst,lhs,rhs,alpha); }
+
 };
 
 template<typename Lhs, typename Rhs>
-struct product_evaluator_dispatcher<Product<Lhs, Rhs>, GeneralProduct<Lhs, Rhs, InnerProduct> > 
-  : public evaluator<typename Product<Lhs, Rhs>::PlainObject>::type
+struct dense_product_impl<Lhs,Rhs,GemvProduct> : dense_product_impl_base<Lhs,Rhs,dense_product_impl<Lhs,Rhs,GemvProduct> >
 {
-  typedef Product<Lhs, Rhs> XprType;
-  typedef typename XprType::PlainObject PlainObject;
-  typedef typename evaluator<PlainObject>::type evaluator_base;
+  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
+  enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight };
+  typedef typename internal::conditional<int(Side)==OnTheRight,Lhs,Rhs>::type MatrixType;
 
-  // TODO: Computation is too early (?)
-  product_evaluator_dispatcher(const XprType& xpr) : evaluator_base(m_result)
+  template<typename Dest>
+  static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
   {
-    m_result.coeffRef(0,0) = (xpr.lhs().transpose().cwiseProduct(xpr.rhs())).sum();
+    internal::gemv_selector<Side,
+                            (int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor,
+                            bool(internal::blas_traits<MatrixType>::HasUsableDirectAccess)
+                           >::run(GeneralProduct<Lhs,Rhs,GemvProduct>(lhs,rhs), dst, alpha);
   }
-  
-protected:  
-  PlainObject m_result;
 };
 
-// For the other three subcases, simply call the evalTo() method of GeneralProduct
-// TODO: GeneralProduct should take evaluators, not expression objects.
-
-template<typename Lhs, typename Rhs, int ProductType>
-struct product_evaluator_traits_dispatcher<Product<Lhs, Rhs>, GeneralProduct<Lhs, Rhs, ProductType> > 
+template<typename Lhs, typename Rhs>
+struct dense_product_impl<Lhs,Rhs,GemmProduct> : dense_product_impl_base<Lhs,Rhs,dense_product_impl<Lhs,Rhs,GemmProduct> >
 {
-  static const int HasEvalTo = 1;
+  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
+  
+  template<typename Dest>
+  static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
+  {
+    // TODO bypass GeneralProduct class
+    GeneralProduct<Lhs, Rhs, GemmProduct>(lhs,rhs).scaleAndAddTo(dst, alpha);
+  }
 };
 
-template<typename Lhs, typename Rhs, int ProductType>
-struct product_evaluator_dispatcher<Product<Lhs, Rhs>, GeneralProduct<Lhs, Rhs, ProductType> > 
+template<typename Lhs, typename Rhs>
+struct dense_product_impl<Lhs,Rhs,CoeffBasedProductMode> 
 {
-  typedef Product<Lhs, Rhs> XprType;
-  typedef typename XprType::PlainObject PlainObject;
-  typedef typename evaluator<PlainObject>::type evaluator_base;
+  typedef typename Product<Lhs,Rhs>::Scalar Scalar;
   
-  product_evaluator_dispatcher(const XprType& xpr) : m_xpr(xpr)
-  { }
+  template<typename Dst>
+  static inline void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
+  { dst = lazyprod(lhs,rhs); }
   
-  template<typename DstEvaluatorType, typename DstXprType>
-  void evalTo(DstEvaluatorType /* not used */, DstXprType& dst) const
-  {
-    dst.resize(m_xpr.rows(), m_xpr.cols());
-    GeneralProduct<Lhs, Rhs, ProductType>(m_xpr.lhs(), m_xpr.rhs()).evalTo(dst);
-  }
+  template<typename Dst>
+  static inline void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
+  { dst += lazyprod(lhs,rhs); }
+  
+  template<typename Dst>
+  static inline void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs)
+  { dst -= lazyprod(lhs,rhs); }
   
-protected: 
-  const XprType& m_xpr;
+  template<typename Dst>
+  static inline void scaleAndAddTo(Dst& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha)
+  { dst += alpha * lazyprod(lhs,rhs); }
 };
 
+template<typename Lhs, typename Rhs>
+struct dense_product_impl<Lhs,Rhs,LazyCoeffBasedProductMode> : dense_product_impl<Lhs,Rhs,CoeffBasedProductMode> {};
+
 // Case 2: Evaluate coeff by coeff
 //
 // This is mostly taken from CoeffBasedProduct.h
@@ -117,20 +195,14 @@ struct etor_product_coeff_impl;
 template<int StorageOrder, int UnrollingIndex, typename Lhs, typename Rhs, typename Packet, int LoadMode>
 struct etor_product_packet_impl;
 
-template<typename Lhs, typename Rhs, typename LhsNested, typename RhsNested, int Flags>
-struct product_evaluator_traits_dispatcher<Product<Lhs, Rhs>, CoeffBasedProduct<LhsNested, RhsNested, Flags> >
+template<typename Lhs, typename Rhs, int ProductTag>
+struct evaluator_impl<Product<Lhs, Rhs, LazyProduct, ProductTag> > 
+    : evaluator_impl_base<Product<Lhs, Rhs, LazyProduct, ProductTag> >
 {
-  static const int HasEvalTo = 0;
-};
+  typedef Product<Lhs, Rhs, LazyProduct, ProductTag> XprType;
+  typedef CoeffBasedProduct<Lhs, Rhs, 0> CoeffBasedProductType;
 
-template<typename Lhs, typename Rhs, typename LhsNested, typename RhsNested, int Flags>
-struct product_evaluator_dispatcher<Product<Lhs, Rhs>, CoeffBasedProduct<LhsNested, RhsNested, Flags> >
-  : evaluator_impl_base<Product<Lhs, Rhs> >
-{
-  typedef Product<Lhs, Rhs> XprType;
-  typedef CoeffBasedProduct<LhsNested, RhsNested, Flags> CoeffBasedProductType;
-
-  product_evaluator_dispatcher(const XprType& xpr) 
+  evaluator_impl(const XprType& xpr) 
     : m_lhsImpl(xpr.lhs()), 
       m_rhsImpl(xpr.rhs()),  
       m_innerDim(xpr.lhs().cols())
@@ -150,11 +222,13 @@ struct product_evaluator_dispatcher<Product<Lhs, Rhs>, CoeffBasedProduct<LhsNest
     InnerSize  = traits<CoeffBasedProductType>::InnerSize,
     CoeffReadCost = traits<CoeffBasedProductType>::CoeffReadCost,
     Unroll = CoeffReadCost != Dynamic && CoeffReadCost <= EIGEN_UNROLLING_LIMIT,
-    CanVectorizeInner = traits<CoeffBasedProductType>::CanVectorizeInner
+    CanVectorizeInner = traits<CoeffBasedProductType>::CanVectorizeInner,
+    Flags = CoeffBasedProductType::Flags
   };
 
   typedef typename evaluator<Lhs>::type LhsEtorType;
   typedef typename evaluator<Rhs>::type RhsEtorType;
+  
   typedef etor_product_coeff_impl<CanVectorizeInner ? InnerVectorizedTraversal : DefaultTraversal,
                                   Unroll ? InnerSize-1 : Dynamic,
                                   LhsEtorType, RhsEtorType, Scalar> CoeffImpl;
@@ -183,8 +257,8 @@ struct product_evaluator_dispatcher<Product<Lhs, Rhs>, CoeffBasedProduct<LhsNest
   {
     PacketScalar res;
     typedef etor_product_packet_impl<Flags&RowMajorBit ? RowMajor : ColMajor,
-				     Unroll ? InnerSize-1 : Dynamic,
-				     LhsEtorType, RhsEtorType, PacketScalar, LoadMode> PacketImpl;
+                                     Unroll ? InnerSize-1 : Dynamic,
+                                     LhsEtorType, RhsEtorType, PacketScalar, LoadMode> PacketImpl;
     PacketImpl::run(row, col, m_lhsImpl, m_rhsImpl, m_innerDim, res);
     return res;
   }
@@ -197,6 +271,7 @@ protected:
   Index m_innerDim;
 };
 
+
 /***************************************************************************
 * Normal product .coeff() implementation (with meta-unrolling)
 ***************************************************************************/
@@ -275,7 +350,6 @@ struct etor_product_coeff_impl<InnerVectorizedTraversal, UnrollingIndex, Lhs, Rh
   {
     Packet pres;
     etor_product_coeff_vectorized_unroller<UnrollingIndex+1-PacketSize, Lhs, Rhs, Packet>::run(row, col, lhs, rhs, innerDim, pres);
-    etor_product_coeff_impl<DefaultTraversal,UnrollingIndex,Lhs,Rhs,RetScalar>::run(row, col, lhs, rhs, innerDim, res);
     res = predux(pres);
   }
 };
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h
index 05107fdfe..3178ff06e 100644
--- a/Eigen/src/Core/util/Constants.h
+++ b/Eigen/src/Core/util/Constants.h
@@ -425,7 +425,7 @@ namespace Architecture
 
 /** \internal \ingroup enums
   * Enum used as template parameter in GeneralProduct. */
-enum { CoeffBasedProductMode, LazyCoeffBasedProductMode, OuterProduct, InnerProduct, GemvProduct, GemmProduct };
+enum { DefaultProduct=0, CoeffBasedProductMode, LazyCoeffBasedProductMode, LazyProduct, OuterProduct, InnerProduct, GemvProduct, GemmProduct };
 
 /** \internal \ingroup enums
   * Enum used in experimental parallel implementation. */
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index 0a2144c69..459422524 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -87,11 +87,20 @@ template<typename NullaryOp, typename MatrixType>         class CwiseNullaryOp;
 template<typename UnaryOp,   typename MatrixType>         class CwiseUnaryOp;
 template<typename ViewOp,    typename MatrixType>         class CwiseUnaryView;
 template<typename BinaryOp,  typename Lhs, typename Rhs>  class CwiseBinaryOp;
-template<typename BinOp,     typename Lhs, typename Rhs>  class SelfCwiseBinaryOp;
+template<typename BinOp,     typename Lhs, typename Rhs>  class SelfCwiseBinaryOp;      // TODO deprecated
 template<typename Derived,   typename Lhs, typename Rhs>  class ProductBase;
-template<typename Lhs, typename Rhs>                      class Product;
-template<typename Lhs, typename Rhs, int Mode>            class GeneralProduct;
-template<typename Lhs, typename Rhs, int NestingFlags>    class CoeffBasedProduct;
+
+namespace internal {
+  template<typename Lhs, typename Rhs> struct product_tag;
+}
+
+template<typename Lhs, typename Rhs,
+         int Option     = DefaultProduct,
+         int ProductTag = internal::product_tag<Lhs,Rhs>::ret
+         > class Product;
+         
+template<typename Lhs, typename Rhs, int Mode>            class GeneralProduct;         // TODO deprecated
+template<typename Lhs, typename Rhs, int NestingFlags>    class CoeffBasedProduct;      // TODO deprecated
 
 template<typename Derived> class DiagonalBase;
 template<typename _DiagonalVectorType> class DiagonalWrapper;