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 +// Copyright (C) 2011 Gael Guennebaud +// Copyright (C) 2011 Jitse Niesen +// +// 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 . + +#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 +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::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::Traversal, + int Unrolling = copy_using_evaluator_traits::Unrolling> +struct copy_using_evaluator_impl; + +template +struct copy_using_evaluator_impl +{ + static void run(const LhsXprType& lhs, const RhsXprType& rhs) + { + typedef typename evaluator::type LhsEvaluatorType; + typedef typename evaluator::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 +void copy_using_evaluator(const LhsXprType& lhs, const RhsXprType& rhs) +{ + copy_using_evaluator_impl::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 +// Copyright (C) 2011 Gael Guennebaud +// Copyright (C) 2011 Jitse Niesen +// +// 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 . + + +#ifndef EIGEN_COREEVALUATORS_H +#define EIGEN_COREEVALUATORS_H + +namespace internal { + +template +struct evaluator_impl {}; + +template +struct evaluator +{ + typedef evaluator_impl type; +}; + +template +struct evaluator +{ + typedef evaluator_impl type; +}; + + +template +struct evaluator_impl > +{ + typedef Transpose 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::type m_argImpl; +}; + + +template +struct evaluator_impl > +{ + typedef Matrix 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 +struct evaluator_impl > +{ + typedef Array 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 +struct evaluator_impl > +{ + typedef CwiseNullaryOp 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 +struct evaluator_impl > +{ + typedef CwiseUnaryOp 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::type m_argImpl; +}; + + +template +struct evaluator_impl > +{ + typedef CwiseBinaryOp 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::type m_lhsImpl; + typename evaluator::type m_rhsImpl; +}; + +// products + +template +struct evaluator_impl > : public evaluator::PlainObject>::type +{ + typedef GeneralProduct XprType; + typedef typename XprType::PlainObject PlainObject; + typedef typename evaluator::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::type m_lhsImpl; + typename evaluator::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())); +}