Evaluators: Implement LinearVectorizedTraversal, packet ops in evaluators.

Eigen/src/Core/AssignEvaluator.h

@@ -77,7 +77,7 @@ private:
 public:
   enum {
     Traversal = int(MayInnerVectorize)  ? int(DefaultTraversal)  // int(InnerVectorizedTraversal)
-              : int(MayLinearVectorize) ? int(DefaultTraversal)  // int(LinearVectorizedTraversal)
+              : int(MayLinearVectorize) ? int(LinearVectorizedTraversal)
               : int(MaySliceVectorize)  ? int(DefaultTraversal)  // int(SliceVectorizedTraversal)
               : int(MayLinearize)       ? int(DefaultTraversal)  // int(LinearTraversal)
                                         : int(DefaultTraversal),
@@ -140,40 +140,104 @@ public:
 
 // copy_using_evaluator_impl is based on assign_impl
 
-template<typename LhsXprType, typename RhsXprType,
+template<typename DstXprType, typename SrcXprType,
-         int Traversal = copy_using_evaluator_traits<LhsXprType, RhsXprType>::Traversal,
+         int Traversal = copy_using_evaluator_traits<DstXprType, SrcXprType>::Traversal,
-         int Unrolling = copy_using_evaluator_traits<LhsXprType, RhsXprType>::Unrolling>
+         int Unrolling = copy_using_evaluator_traits<DstXprType, SrcXprType>::Unrolling>
 struct copy_using_evaluator_impl;
 
-template<typename LhsXprType, typename RhsXprType>
+template<typename DstXprType, typename SrcXprType>
-struct copy_using_evaluator_impl<LhsXprType, RhsXprType, DefaultTraversal, NoUnrolling>
+struct copy_using_evaluator_impl<DstXprType, SrcXprType, DefaultTraversal, NoUnrolling>
 {
-  static void run(const LhsXprType& lhs, const RhsXprType& rhs)
+  static void run(const DstXprType& dst, const SrcXprType& src)
   {
-    typedef typename evaluator<LhsXprType>::type LhsEvaluatorType;
+    typedef typename evaluator<DstXprType>::type DstEvaluatorType;
-    typedef typename evaluator<RhsXprType>::type RhsEvaluatorType;
+    typedef typename evaluator<SrcXprType>::type SrcEvaluatorType;
-    typedef typename LhsXprType::Index Index;
+    typedef typename DstXprType::Index Index;
 
-    LhsEvaluatorType lhsEvaluator(lhs.const_cast_derived());
+    DstEvaluatorType dstEvaluator(dst.const_cast_derived());
-    RhsEvaluatorType rhsEvaluator(rhs);
+    SrcEvaluatorType srcEvaluator(src);
 
-    for(Index outer = 0; outer < lhs.outerSize(); ++outer) {
+    for(Index outer = 0; outer < dst.outerSize(); ++outer) {
-      for(Index inner = 0; inner < lhs.innerSize(); ++inner) {
+      for(Index inner = 0; inner < dst.innerSize(); ++inner) {
-	Index row = lhs.rowIndexByOuterInner(outer, inner);
+	Index row = dst.rowIndexByOuterInner(outer, inner);
-	Index col = lhs.colIndexByOuterInner(outer, inner);
+	Index col = dst.colIndexByOuterInner(outer, inner);
-	lhsEvaluator.coeffRef(row, col) = rhsEvaluator.coeff(row, col);
+	dstEvaluator.coeffRef(row, col) = srcEvaluator.coeff(row, col); // TODO: use copyCoeff ?
       }
     }
   }
 };
 
+template <bool IsAligned = false>
+struct unaligned_copy_using_evaluator_impl
+{
+  template <typename SrcEvaluatorType, typename DstEvaluatorType>
+  static EIGEN_STRONG_INLINE void run(const SrcEvaluatorType&, DstEvaluatorType&, 
+				      typename SrcEvaluatorType::Index, typename SrcEvaluatorType::Index) {}
+};
+
+// TODO: check why no ...<true> ????
+
+template <>
+struct unaligned_copy_using_evaluator_impl<false>
+{
+  // MSVC must not inline this functions. If it does, it fails to optimize the
+  // packet access path.
+#ifdef _MSC_VER
+  template <typename SrcEvaluatorType, typename DstEvaluatorType>
+  static EIGEN_DONT_INLINE void run(const SrcEvaluatorType& src, DstEvaluatorType& dst, 
+				    typename SrcEvaluatorType::Index start, typename SrcEvaluatorType::Index end)
+#else
+  template <typename SrcEvaluatorType, typename DstEvaluatorType>
+  static EIGEN_STRONG_INLINE void run(const SrcEvaluatorType& src, DstEvaluatorType& dst, 
+				      typename SrcEvaluatorType::Index start, typename SrcEvaluatorType::Index end)
+#endif
+  {
+    for (typename SrcEvaluatorType::Index index = start; index < end; ++index)
+      dst.copyCoeff(index, src);
+  }
+};
+
+template<typename DstXprType, typename SrcXprType>
+struct copy_using_evaluator_impl<DstXprType, SrcXprType, LinearVectorizedTraversal, NoUnrolling>
+{
+  EIGEN_STRONG_INLINE static void run(const DstXprType &dst, const SrcXprType &src)
+  {
+    typedef typename evaluator<DstXprType>::type DstEvaluatorType;
+    typedef typename evaluator<SrcXprType>::type SrcEvaluatorType;
+    typedef typename DstXprType::Index Index;
+
+    DstEvaluatorType dstEvaluator(dst.const_cast_derived());
+    SrcEvaluatorType srcEvaluator(src);
+
+    const Index size = dst.size();
+    typedef packet_traits<typename DstXprType::Scalar> PacketTraits;
+    enum {
+      packetSize = PacketTraits::size,
+      dstIsAligned = int(copy_using_evaluator_traits<DstXprType,SrcXprType>::DstIsAligned),
+      dstAlignment = PacketTraits::AlignedOnScalar ? Aligned : dstIsAligned,
+      srcAlignment = copy_using_evaluator_traits<DstXprType,SrcXprType>::JointAlignment
+    };
+    const Index alignedStart = dstIsAligned ? 0 : first_aligned(&dst.coeffRef(0), size);
+    const Index alignedEnd = alignedStart + ((size-alignedStart)/packetSize)*packetSize;
+
+    unaligned_copy_using_evaluator_impl<dstIsAligned!=0>::run(src,dst.const_cast_derived(),0,alignedStart);
+
+    for(Index index = alignedStart; index < alignedEnd; index += packetSize)
+    {
+      dstEvaluator.template writePacket<dstAlignment>(index, srcEvaluator.template packet<srcAlignment>(index));
+    }
+
+    unaligned_copy_using_evaluator_impl<>::run(src,dst.const_cast_derived(),alignedEnd,size);
+  }
+};
+
 // Based on DenseBase::LazyAssign()
 
-template<typename LhsXprType, typename RhsXprType>
+template<typename DstXprType, typename SrcXprType>
-const LhsXprType& copy_using_evaluator(const LhsXprType& lhs, const RhsXprType& rhs)
+const DstXprType& copy_using_evaluator(const DstXprType& dst, const SrcXprType& src)
 {
-  copy_using_evaluator_impl<LhsXprType, RhsXprType>::run(lhs, rhs);
+  copy_using_evaluator_impl<DstXprType, SrcXprType>::run(dst, src);
-  return lhs;
+  return dst;
 }
 
 } // namespace internal

Eigen/src/Core/CoreEvaluators.h

@@ -45,6 +45,7 @@ struct evaluator<const T>
   typedef evaluator_impl<T> type;
 };
 
+// -------------------- Transpose --------------------
 
 template<typename ExpressionType>
 struct evaluator_impl<Transpose<ExpressionType> >
@@ -64,10 +65,17 @@ struct evaluator_impl<Transpose<ExpressionType> >
     return m_argImpl.coeffRef(j, i);
   }
 
+  template<int LoadMode>
+  const typename ExpressionType::PacketScalar packet(Index index) const
+  {
+    return m_argImpl.template packet<LoadMode>(index);
+  }
+
 protected:
   typename evaluator<ExpressionType>::type m_argImpl;
 };
 
+// -------------------- Matrix --------------------
 
 template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
 struct evaluator_impl<Matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> >
@@ -88,10 +96,27 @@ struct evaluator_impl<Matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> >
     return m_matrix.const_cast_derived().coeffRef(i, j);
   }
 
+  template<int LoadMode> 
+  typename MatrixType::PacketReturnType packet(Index index) const
+  {
+    // eigen_internal_assert(index >= 0 && index < size());
+    return m_matrix.template packet<LoadMode>(index);
+  }
+
+  template<int StoreMode> 
+  void writePacket(Index index, const typename MatrixType::PacketScalar& x)
+  {
+    // eigen_internal_assert(index >= 0 && index < size());
+    m_matrix.const_cast_derived().template writePacket<StoreMode>(index, x);
+  }
+
 protected:
   const MatrixType &m_matrix;
 };
 
+// -------------------- Array --------------------
+
+// TODO: should be sharing code with Matrix case
 
 template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
 struct evaluator_impl<Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols> >
@@ -117,10 +142,25 @@ struct evaluator_impl<Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols> >
     return m_array.const_cast_derived().coeffRef(i, j);
   }
 
+  template<int LoadMode> 
+  typename ArrayType::PacketReturnType packet(Index index) const
+  {
+    // eigen_internal_assert(index >= 0 && index < size());
+    return m_array.template packet<LoadMode>(index);
+  }
+
+  template<int StoreMode> 
+  void writePacket(Index index, const typename ArrayType::PacketScalar& x)
+  {
+    // eigen_internal_assert(index >= 0 && index < size());
+    m_array.const_cast_derived().template writePacket<StoreMode>(index, x);
+  }
+
 protected:
   const ArrayType &m_array;
 };
 
+// -------------------- CwiseNullaryOp --------------------
 
 template<typename NullaryOp, typename PlainObjectType>
 struct evaluator_impl<CwiseNullaryOp<NullaryOp,PlainObjectType> >
@@ -136,10 +176,17 @@ struct evaluator_impl<CwiseNullaryOp<NullaryOp,PlainObjectType> >
     return m_nullaryOp.coeff(i, j);
   }
 
+  template<int LoadMode>
+  typename NullaryOpType::PacketScalar packet(Index index) const
+  {
+    return m_nullaryOp.template packet<LoadMode>(index);
+  }
+
 protected:
   const NullaryOpType& m_nullaryOp;
 };
 
+// -------------------- CwiseUnaryOp --------------------
 
 template<typename UnaryOp, typename ArgType>
 struct evaluator_impl<CwiseUnaryOp<UnaryOp, ArgType> >
@@ -155,11 +202,18 @@ struct evaluator_impl<CwiseUnaryOp<UnaryOp, ArgType> >
     return m_unaryOp.functor()(m_argImpl.coeff(i, j));
   }
 
+  template<int LoadMode>
+  typename UnaryOpType::PacketScalar packet(Index index) const
+  {
+    return m_unaryOp.functor().packetOp(m_argImpl.template packet<LoadMode>(index));
+  }
+
 protected:
   const UnaryOpType& m_unaryOp;
   typename evaluator<ArgType>::type m_argImpl;
 };
 
+// -------------------- CwiseBinaryOp --------------------
 
 template<typename BinaryOp, typename Lhs, typename Rhs>
 struct evaluator_impl<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
@@ -172,7 +226,14 @@ struct evaluator_impl<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
 
   typename BinaryOpType::CoeffReturnType coeff(Index i, Index j) const
   {
-    return m_binaryOp.functor()(m_lhsImpl.coeff(i, j),m_rhsImpl.coeff(i, j));
+    return m_binaryOp.functor()(m_lhsImpl.coeff(i, j), m_rhsImpl.coeff(i, j));
+  }
+
+  template<int LoadMode>
+  typename BinaryOpType::PacketScalar packet(Index index) const
+  {
+    return m_binaryOp.functor().packetOp(m_lhsImpl.template packet<LoadMode>(index),
+					 m_rhsImpl.template packet<LoadMode>(index));
   }
 
 protected:
@@ -181,7 +242,7 @@ protected:
   typename evaluator<Rhs>::type m_rhsImpl;
 };
 
-// product
+// -------------------- Product --------------------
 
 template<typename Lhs, typename Rhs>
 struct evaluator_impl<Product<Lhs,Rhs> > : public evaluator<typename Product<Lhs,Rhs>::PlainObject>::type

test/evaluators.cpp

@@ -78,4 +78,15 @@ void test_evaluators()
   // test CwiseBinaryOp
   VERIFY_IS_APPROX_EVALUATOR(v2, v + Vector2d::Ones());
   VERIFY_IS_APPROX_EVALUATOR(w, (v + Vector2d::Ones()).transpose().cwiseProduct(RowVector2d::Constant(3)));
+
+  // dynamic matrices and arrays
+  MatrixXd mat1(6,6), mat2(6,6);
+  VERIFY_IS_APPROX_EVALUATOR(mat1, MatrixXd::Identity(6,6));
+  VERIFY_IS_APPROX_EVALUATOR(mat2, mat1);
+  copy_using_evaluator(mat2.transpose(), mat1);
+  VERIFY_IS_APPROX(mat2.transpose(), mat1);
+
+  ArrayXXd arr1(6,6), arr2(6,6);
+  VERIFY_IS_APPROX_EVALUATOR(arr1, ArrayXXd::Constant(6,6, 3.0));
+  VERIFY_IS_APPROX_EVALUATOR(arr2, arr1);
 }