import evaluator works

5 files changed, 502 insertions, 0 deletions

diff --git a/Eigen/Core b/Eigen/Core
index 7f384662e..aeaefbed6 100644
--- a/Eigen/Core
+++ b/Eigen/Core
@@ -350,6 +350,11 @@ using std::size_t;
 #include "src/Core/ArrayBase.h"
 #include "src/Core/ArrayWrapper.h"
 
+#ifdef EIGEN_ENABLE_EVALUATORS
+#include "src/Core/CoreEvaluators.h"
+#include "src/Core/AssignEvaluator.h"
+#endif 
+
 } // namespace Eigen
 
 #include "src/Core/GlobalFunctions.h"
diff --git a/Eigen/src/Core/AssignEvaluator.h b/Eigen/src/Core/AssignEvaluator.h
new file mode 100644
index 000000000..20c3b0911
--- /dev/null
+++ b/Eigen/src/Core/AssignEvaluator.h
@@ -0,0 +1,180 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2011 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_ASSIGN_EVALUATOR_H
+#define EIGEN_ASSIGN_EVALUATOR_H
+
+// This implementation is based on Assign.h
+
+// copy_using_evaluator_traits is based on assign_traits
+
+namespace internal {
+  
+template <typename Derived, typename OtherDerived>
+struct copy_using_evaluator_traits
+{
+public:
+  enum {
+    DstIsAligned = Derived::Flags & AlignedBit,
+    DstHasDirectAccess = Derived::Flags & DirectAccessBit,
+    SrcIsAligned = OtherDerived::Flags & AlignedBit,
+    JointAlignment = bool(DstIsAligned) && bool(SrcIsAligned) ? Aligned : Unaligned
+  };
+
+private:
+  enum {
+    InnerSize = int(Derived::IsVectorAtCompileTime) ? int(Derived::SizeAtCompileTime)
+              : int(Derived::Flags)&RowMajorBit ? int(Derived::ColsAtCompileTime)
+              : int(Derived::RowsAtCompileTime),
+    InnerMaxSize = int(Derived::IsVectorAtCompileTime) ? int(Derived::MaxSizeAtCompileTime)
+              : int(Derived::Flags)&RowMajorBit ? int(Derived::MaxColsAtCompileTime)
+              : int(Derived::MaxRowsAtCompileTime),
+    MaxSizeAtCompileTime = Derived::SizeAtCompileTime,
+    PacketSize = packet_traits<typename Derived::Scalar>::size
+  };
+
+  enum {
+    StorageOrdersAgree = (int(Derived::IsRowMajor) == int(OtherDerived::IsRowMajor)),
+    MightVectorize = StorageOrdersAgree
+                  && (int(Derived::Flags) & int(OtherDerived::Flags) & ActualPacketAccessBit),
+    MayInnerVectorize  = MightVectorize && int(InnerSize)!=Dynamic && int(InnerSize)%int(PacketSize)==0
+                       && int(DstIsAligned) && int(SrcIsAligned),
+    MayLinearize = StorageOrdersAgree && (int(Derived::Flags) & int(OtherDerived::Flags) & LinearAccessBit),
+    MayLinearVectorize = MightVectorize && MayLinearize && DstHasDirectAccess
+                       && (DstIsAligned || MaxSizeAtCompileTime == Dynamic),
+      /* If the destination isn't aligned, we have to do runtime checks and we don't unroll,
+         so it's only good for large enough sizes. */
+    MaySliceVectorize  = MightVectorize && DstHasDirectAccess
+                       && (int(InnerMaxSize)==Dynamic || int(InnerMaxSize)>=3*PacketSize)
+      /* slice vectorization can be slow, so we only want it if the slices are big, which is
+         indicated by InnerMaxSize rather than InnerSize, think of the case of a dynamic block
+         in a fixed-size matrix */
+  };
+
+public:
+  enum {
+    Traversal = int(MayInnerVectorize)  ? int(DefaultTraversal)  // int(InnerVectorizedTraversal)
+              : int(MayLinearVectorize) ? int(DefaultTraversal)  // int(LinearVectorizedTraversal)
+              : int(MaySliceVectorize)  ? int(DefaultTraversal)  // int(SliceVectorizedTraversal)
+              : int(MayLinearize)       ? int(DefaultTraversal)  // int(LinearTraversal)
+                                        : int(DefaultTraversal),
+    Vectorized = int(Traversal) == InnerVectorizedTraversal
+              || int(Traversal) == LinearVectorizedTraversal
+              || int(Traversal) == SliceVectorizedTraversal
+  };
+
+private:
+  enum {
+    UnrollingLimit      = EIGEN_UNROLLING_LIMIT * (Vectorized ? int(PacketSize) : 1),
+    MayUnrollCompletely = int(Derived::SizeAtCompileTime) != Dynamic
+                       && int(OtherDerived::CoeffReadCost) != Dynamic
+                       && int(Derived::SizeAtCompileTime) * int(OtherDerived::CoeffReadCost) <= int(UnrollingLimit),
+    MayUnrollInner      = int(InnerSize) != Dynamic
+                       && int(OtherDerived::CoeffReadCost) != Dynamic
+                       && int(InnerSize) * int(OtherDerived::CoeffReadCost) <= int(UnrollingLimit)
+  };
+
+public:
+  enum {
+    Unrolling = (int(Traversal) == int(InnerVectorizedTraversal) || int(Traversal) == int(DefaultTraversal))
+                ? (
+ 		    int(MayUnrollCompletely) ? int(NoUnrolling)    // int(CompleteUnrolling)
+                  : int(MayUnrollInner)      ? int(NoUnrolling)    // int(InnerUnrolling)
+                                             : int(NoUnrolling)
+                  )
+              : int(Traversal) == int(LinearVectorizedTraversal)
+                ? ( bool(MayUnrollCompletely) && bool(DstIsAligned) ? int(NoUnrolling)   // int(CompleteUnrolling) 
+                                                                    : int(NoUnrolling) )
+              : int(Traversal) == int(LinearTraversal)
+                ? ( bool(MayUnrollCompletely) ? int(NoUnrolling)   // int(CompleteUnrolling) 
+                                              : int(NoUnrolling) )
+              : int(NoUnrolling)
+  };
+
+#ifdef EIGEN_DEBUG_ASSIGN
+  static void debug()
+  {
+    EIGEN_DEBUG_VAR(DstIsAligned)
+    EIGEN_DEBUG_VAR(SrcIsAligned)
+    EIGEN_DEBUG_VAR(JointAlignment)
+    EIGEN_DEBUG_VAR(InnerSize)
+    EIGEN_DEBUG_VAR(InnerMaxSize)
+    EIGEN_DEBUG_VAR(PacketSize)
+    EIGEN_DEBUG_VAR(StorageOrdersAgree)
+    EIGEN_DEBUG_VAR(MightVectorize)
+    EIGEN_DEBUG_VAR(MayLinearize)
+    EIGEN_DEBUG_VAR(MayInnerVectorize)
+    EIGEN_DEBUG_VAR(MayLinearVectorize)
+    EIGEN_DEBUG_VAR(MaySliceVectorize)
+    EIGEN_DEBUG_VAR(Traversal)
+    EIGEN_DEBUG_VAR(UnrollingLimit)
+    EIGEN_DEBUG_VAR(MayUnrollCompletely)
+    EIGEN_DEBUG_VAR(MayUnrollInner)
+    EIGEN_DEBUG_VAR(Unrolling)
+  }
+#endif
+};
+
+// copy_using_evaluator_impl is based on assign_impl
+
+template<typename LhsXprType, typename RhsXprType,
+         int Traversal = copy_using_evaluator_traits<LhsXprType, RhsXprType>::Traversal,
+         int Unrolling = copy_using_evaluator_traits<LhsXprType, RhsXprType>::Unrolling>
+struct copy_using_evaluator_impl;
+
+template<typename LhsXprType, typename RhsXprType>
+struct copy_using_evaluator_impl<LhsXprType, RhsXprType, DefaultTraversal, NoUnrolling>
+{
+  static void run(const LhsXprType& lhs, const RhsXprType& rhs)
+  {
+    typedef typename evaluator<LhsXprType>::type LhsEvaluatorType;
+    typedef typename evaluator<RhsXprType>::type RhsEvaluatorType;
+    typedef typename LhsXprType::Index Index;
+
+    LhsEvaluatorType lhsEvaluator(lhs.const_cast_derived());
+    RhsEvaluatorType rhsEvaluator(rhs);
+
+    for(Index outer = 0; outer < lhs.outerSize(); ++outer) {
+      for(Index inner = 0; inner < lhs.innerSize(); ++inner) {
+	Index row = lhs.rowIndexByOuterInner(outer, inner);
+	Index col = lhs.colIndexByOuterInner(outer, inner);
+	lhsEvaluator.coeffRef(row, col) = rhsEvaluator.coeff(row, col);
+      }
+    }
+  }
+};
+
+// Based on DenseBase::LazyAssign()
+
+template<typename LhsXprType, typename RhsXprType>
+void copy_using_evaluator(const LhsXprType& lhs, const RhsXprType& rhs)
+{
+  copy_using_evaluator_impl<LhsXprType, RhsXprType>::run(lhs, rhs);
+}
+
+} // namespace internal
+
+#endif // EIGEN_ASSIGN_EVALUATOR_H
diff --git a/Eigen/src/Core/CoreEvaluators.h b/Eigen/src/Core/CoreEvaluators.h
new file mode 100644
index 000000000..64ec21a0f
--- /dev/null
+++ b/Eigen/src/Core/CoreEvaluators.h
@@ -0,0 +1,212 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2011 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+
+#ifndef EIGEN_COREEVALUATORS_H
+#define EIGEN_COREEVALUATORS_H
+
+namespace internal {
+  
+template<typename T>
+struct evaluator_impl {};
+
+template<typename T>
+struct evaluator
+{
+  typedef evaluator_impl<T> type;
+};
+
+template<typename T>
+struct evaluator<const T>
+{
+  typedef evaluator_impl<T> type;
+};
+
+
+template<typename ExpressionType>
+struct evaluator_impl<Transpose<ExpressionType> >
+{
+  typedef Transpose<ExpressionType> TransposeType;
+  evaluator_impl(const TransposeType& t) : m_argImpl(t.nestedExpression()) {}
+
+  typedef typename TransposeType::Index Index;
+
+  typename TransposeType::CoeffReturnType coeff(Index i, Index j) const
+  {
+    return m_argImpl.coeff(j, i);
+  }
+
+  typename TransposeType::Scalar& coeffRef(Index i, Index j)
+  {
+    return m_argImpl.coeffRef(j, i);
+  }
+
+protected:
+  typename evaluator<ExpressionType>::type m_argImpl;
+};
+
+
+template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
+struct evaluator_impl<Matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> >
+{
+  typedef Matrix<Scalar, Rows, Cols, Options, MaxRows, MaxCols> MatrixType;
+
+  evaluator_impl(const MatrixType& m) : m_matrix(m) {}
+
+  typedef typename MatrixType::Index Index;
+
+  typename MatrixType::CoeffReturnType coeff(Index i, Index j) const
+  {
+    return m_matrix.coeff(i, j);
+  }
+
+  typename MatrixType::Scalar& coeffRef(Index i, Index j)
+  {
+    return m_matrix.const_cast_derived().coeffRef(i, j);
+  }
+
+protected:
+  const MatrixType &m_matrix;
+};
+
+
+template<typename Scalar, int Rows, int Cols, int Options, int MaxRows, int MaxCols>
+struct evaluator_impl<Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols> >
+{
+  typedef Array<Scalar, Rows, Cols, Options, MaxRows, MaxCols> ArrayType;
+
+  evaluator_impl(const ArrayType& a) : m_array(a) {}
+
+  typedef typename ArrayType::Index Index;
+  
+  Index colIndexByOuterInner(Index outer, Index inner) const 
+  { 
+    return m_array.colIndexByOuterInner(outer, inner); 
+  }
+
+  typename ArrayType::CoeffReturnType coeff(Index i, Index j) const
+  {
+    return m_array.coeff(i, j);
+  }
+
+  typename ArrayType::Scalar& coeffRef(Index i, Index j)
+  {
+    return m_array.const_cast_derived().coeffRef(i, j);
+  }
+
+protected:
+  const ArrayType &m_array;
+};
+
+
+template<typename NullaryOp, typename PlainObjectType>
+struct evaluator_impl<CwiseNullaryOp<NullaryOp,PlainObjectType> >
+{
+  typedef CwiseNullaryOp<NullaryOp,PlainObjectType> NullaryOpType;
+
+  evaluator_impl(const NullaryOpType& n) : m_nullaryOp(n) {}
+
+  typedef typename NullaryOpType::Index Index;
+
+  typename NullaryOpType::CoeffReturnType coeff(Index i, Index j) const
+  {
+    return m_nullaryOp.coeff(i, j);
+  }
+
+protected:
+  const NullaryOpType& m_nullaryOp;
+};
+
+
+template<typename UnaryOp, typename ArgType>
+struct evaluator_impl<CwiseUnaryOp<UnaryOp, ArgType> >
+{
+  typedef CwiseUnaryOp<UnaryOp, ArgType> UnaryOpType;
+
+  evaluator_impl(const UnaryOpType& op) : m_unaryOp(op), m_argImpl(op.nestedExpression()) {}
+
+  typedef typename UnaryOpType::Index Index;
+
+  typename UnaryOpType::CoeffReturnType coeff(Index i, Index j) const
+  {
+    return m_unaryOp.functor()(m_argImpl.coeff(i, j));
+  }
+
+protected:
+  const UnaryOpType& m_unaryOp;
+  typename evaluator<ArgType>::type m_argImpl;
+};
+
+
+template<typename BinaryOp, typename Lhs, typename Rhs>
+struct evaluator_impl<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
+{
+  typedef CwiseBinaryOp<BinaryOp, Lhs, Rhs> BinaryOpType;
+
+  evaluator_impl(const BinaryOpType& xpr) : m_binaryOp(xpr), m_lhsImpl(xpr.lhs()), m_rhsImpl(xpr.rhs())  {}
+
+  typedef typename BinaryOpType::Index Index;
+
+  typename BinaryOpType::CoeffReturnType coeff(Index i, Index j) const
+  {
+    return m_binaryOp.functor()(m_lhsImpl.coeff(i, j),m_rhsImpl.coeff(i, j));
+  }
+
+protected:
+  const BinaryOpType& m_binaryOp;
+  typename evaluator<Lhs>::type m_lhsImpl;
+  typename evaluator<Rhs>::type m_rhsImpl;
+};
+
+// products
+
+template<typename Lhs, typename Rhs, int ProductType>
+struct evaluator_impl<GeneralProduct<Lhs,Rhs,ProductType> > : public evaluator<typename GeneralProduct<Lhs,Rhs,ProductType>::PlainObject>::type
+{
+  typedef GeneralProduct<Lhs,Rhs,ProductType> XprType;
+  typedef typename XprType::PlainObject PlainObject;
+  typedef typename evaluator<PlainObject>::type evaluator_base;
+  
+//   enum {
+//     EvaluateLhs = ;
+//     EvaluateRhs = ;
+//   };
+  
+  evaluator_impl(const XprType& product) : evaluator_base(m_result), m_lhsImpl(product.lhs()), m_rhsImpl(product.rhs())
+  {
+    m_result.resize(product.rows(), product.cols());
+    product.evalTo(m_result);
+  }
+  
+protected:  
+  PlainObject m_result;
+  typename evaluator<Lhs>::type m_lhsImpl;
+  typename evaluator<Rhs>::type m_rhsImpl;
+};
+
+} // namespace internal
+
+#endif // EIGEN_COREEVALUATORS_H
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index e9586f27b..14a765971 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -120,9 +120,11 @@ ei_add_test(nullary)
 ei_add_test(nesting_ops "${CMAKE_CXX_FLAGS_DEBUG}")
 ei_add_test(zerosized)
 ei_add_test(dontalign)
+ei_add_test(evaluators)
 
 ei_add_test(prec_inverse_4x4)
 
+
 string(TOLOWER "${CMAKE_CXX_COMPILER}" cmake_cxx_compiler_tolower)
 if(cmake_cxx_compiler_tolower MATCHES "qcc")
   set(CXX_IS_QCC "ON")
diff --git a/test/evaluators.cpp b/test/evaluators.cpp
new file mode 100644
index 000000000..ea65439a3
--- /dev/null
+++ b/test/evaluators.cpp
@@ -0,0 +1,103 @@
+
+#define EIGEN_ENABLE_EVALUATORS
+#include "main.h"
+
+using internal::copy_using_evaluator;
+using namespace std;
+
+void test_evaluators()
+{
+  // Testing Matrix evaluator and Transpose
+  Vector2d v(1,2);
+  const Vector2d v_const(v);
+  Vector2d v2;
+  RowVector2d w;
+
+  copy_using_evaluator(v2, v);
+  assert(v2.isApprox((Vector2d() << 1,2).finished()));
+
+  copy_using_evaluator(v2, v_const);
+  assert(v2.isApprox((Vector2d() << 1,2).finished()));
+
+  // Testing Transpose
+  copy_using_evaluator(w, v.transpose()); // Transpose as rvalue
+  assert(w.isApprox((RowVector2d() << 1,2).finished()));
+
+  copy_using_evaluator(w, v_const.transpose());
+  assert(w.isApprox((RowVector2d() << 1,2).finished()));
+
+  copy_using_evaluator(w.transpose(), v); // Transpose as lvalue
+  assert(w.isApprox((RowVector2d() << 1,2).finished()));
+
+  copy_using_evaluator(w.transpose(), v_const);
+  assert(w.isApprox((RowVector2d() << 1,2).finished()));
+
+  // Testing Array evaluator
+  ArrayXXf a(2,3);
+  ArrayXXf b(3,2);
+  a << 1,2,3, 4,5,6;
+  const ArrayXXf a_const(a);
+
+  ArrayXXf b_expected(3,2);
+  b_expected << 1,4, 2,5, 3,6;
+  copy_using_evaluator(b, a.transpose());
+  assert(b.isApprox(b_expected));
+
+  copy_using_evaluator(b, a_const.transpose());
+  assert(b.isApprox(b_expected));
+
+  // Testing CwiseNullaryOp evaluator
+  copy_using_evaluator(w, RowVector2d::Random());
+  assert((w.array() >= -1).all() && (w.array() <= 1).all()); // not easy to test ...
+
+  copy_using_evaluator(w, RowVector2d::Zero());
+  assert(w.isApprox((RowVector2d() << 0,0).finished()));
+
+  copy_using_evaluator(w, RowVector2d::Constant(3));
+  assert(w.isApprox((RowVector2d() << 3,3).finished()));
+  
+  // mix CwiseNullaryOp and transpose
+  copy_using_evaluator(w, Vector2d::Zero().transpose());
+  assert(w.isApprox((RowVector2d() << 0,0).finished()));
+
+  {
+    MatrixXf a(2,2), b(2,2), c(2,2), d(2,2);
+    a << 1, 2, 3, 4; b << 5, 6, 7, 8; c << 9, 10, 11, 12;
+    copy_using_evaluator(d, (a + b));
+    cout << d << endl;
+    
+    copy_using_evaluator(d, (a + b).transpose());
+    cout << d << endl;
+    
+//     copy_using_evaluator(d, (a * b).transpose());
+//     cout << d << endl;
+    
+//     copy_using_evaluator(d, a.transpose() + (a.transpose() * (b+b)));
+//     cout << d << endl;
+  }
+  
+  // this does not work because Random is eval-before-nested: 
+  // copy_using_evaluator(w, Vector2d::Random().transpose());
+  
+  // test CwiseUnaryOp
+  copy_using_evaluator(v2, 3 * v);
+  assert(v2.isApprox((Vector2d() << 3,6).finished()));
+
+  copy_using_evaluator(w, (3 * v).transpose());
+  assert(w.isApprox((RowVector2d() << 3,6).finished()));
+
+  copy_using_evaluator(b, (a + 3).transpose());
+  b_expected << 4,7, 5,8, 6,9;
+  assert(b.isApprox(b_expected));
+
+  copy_using_evaluator(b, (2 * a_const + 3).transpose());
+  b_expected << 5,11, 7,13, 9,15;
+  assert(b.isApprox(b_expected));
+
+  // test CwiseBinaryOp
+  copy_using_evaluator(v2, v + Vector2d::Ones());
+  assert(v2.isApprox((Vector2d() << 2,3).finished()));
+
+  copy_using_evaluator(w, (v + Vector2d::Ones()).transpose().cwiseProduct(RowVector2d::Constant(3)));
+  assert(w.isApprox((RowVector2d() << 6,9).finished()));
+}