From 3849cc65ee2192ddd9d63cdc2e65cb74128bcbb3 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 23 Jun 2014 10:40:03 +0200 Subject: [PATCH] Implement binaryop and transpose evaluators for sparse matrices --- Eigen/SparseCore | 4 +- Eigen/src/Core/AssignEvaluator.h | 10 +- Eigen/src/Core/CoreEvaluators.h | 6 +- Eigen/src/Core/CwiseBinaryOp.h | 2 +- Eigen/src/Core/CwiseUnaryOp.h | 6 +- Eigen/src/Core/Transpose.h | 13 +- Eigen/src/SparseCore/SparseCwiseBinaryOp.h | 312 ++++++++++++++++++++- Eigen/src/SparseCore/SparseCwiseUnaryOp.h | 2 +- Eigen/src/SparseCore/SparseMatrix.h | 2 +- Eigen/src/SparseCore/SparseMatrixBase.h | 4 +- Eigen/src/SparseCore/SparseTranspose.h | 54 +++- 11 files changed, 396 insertions(+), 19 deletions(-) diff --git a/Eigen/SparseCore b/Eigen/SparseCore index c0f0a4121..799cb1ba1 100644 --- a/Eigen/SparseCore +++ b/Eigen/SparseCore @@ -42,10 +42,10 @@ struct Sparse {}; #include "src/SparseCore/MappedSparseMatrix.h" #include "src/SparseCore/SparseVector.h" #include "src/SparseCore/SparseCwiseUnaryOp.h" +#include "src/SparseCore/SparseCwiseBinaryOp.h" +#include "src/SparseCore/SparseTranspose.h" #ifndef EIGEN_TEST_EVALUATORS #include "src/SparseCore/SparseBlock.h" -#include "src/SparseCore/SparseTranspose.h" -#include "src/SparseCore/SparseCwiseBinaryOp.h" #include "src/SparseCore/SparseDot.h" #include "src/SparseCore/SparsePermutation.h" #include "src/SparseCore/SparseRedux.h" diff --git a/Eigen/src/Core/AssignEvaluator.h b/Eigen/src/Core/AssignEvaluator.h index a083e340e..77b4dd9bc 100644 --- a/Eigen/src/Core/AssignEvaluator.h +++ b/Eigen/src/Core/AssignEvaluator.h @@ -2,7 +2,7 @@ // for linear algebra. // // Copyright (C) 2011 Benoit Jacob -// Copyright (C) 2011-2013 Gael Guennebaud +// Copyright (C) 2011-2014 Gael Guennebaud // Copyright (C) 2011-2012 Jitse Niesen // // This Source Code Form is subject to the terms of the Mozilla @@ -738,8 +738,10 @@ void call_assignment_no_alias(Dst& dst, const Src& src, const Func& func) && int(Dst::SizeAtCompileTime) != 1 }; - dst.resize(NeedToTranspose ? src.cols() : src.rows(), - NeedToTranspose ? src.rows() : src.cols()); + typename Dst::Index dstRows = NeedToTranspose ? src.cols() : src.rows(); + typename Dst::Index dstCols = NeedToTranspose ? src.rows() : src.cols(); + if((dst.rows()!=dstRows) || (dst.cols()!=dstCols)) + dst.resize(dstRows, dstCols); typedef typename internal::conditional, Dst>::type ActualDstTypeCleaned; typedef typename internal::conditional, Dst&>::type ActualDstType; @@ -749,7 +751,7 @@ void call_assignment_no_alias(Dst& dst, const Src& src, const Func& func) EIGEN_STATIC_ASSERT_LVALUE(Dst) EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(ActualDstTypeCleaned,Src) EIGEN_CHECK_BINARY_COMPATIBILIY(Func,typename ActualDstTypeCleaned::Scalar,typename Src::Scalar); - + Assignment::run(actualDst, src, func); } template diff --git a/Eigen/src/Core/CoreEvaluators.h b/Eigen/src/Core/CoreEvaluators.h index 5c4da9978..b56c3c635 100644 --- a/Eigen/src/Core/CoreEvaluators.h +++ b/Eigen/src/Core/CoreEvaluators.h @@ -2,7 +2,7 @@ // for linear algebra. // // Copyright (C) 2011 Benoit Jacob -// Copyright (C) 2011 Gael Guennebaud +// Copyright (C) 2011-2014 Gael Guennebaud // Copyright (C) 2011-2012 Jitse Niesen // // This Source Code Form is subject to the terms of the Mozilla @@ -253,7 +253,7 @@ struct evaluator > // -------------------- Transpose -------------------- template -struct unary_evaluator > +struct unary_evaluator, IndexBased> : evaluator_base > { typedef Transpose XprType; @@ -440,7 +440,7 @@ struct evaluator > }; template -struct binary_evaluator > +struct binary_evaluator, IndexBased, IndexBased> : evaluator_base > { typedef CwiseBinaryOp XprType; diff --git a/Eigen/src/Core/CwiseBinaryOp.h b/Eigen/src/Core/CwiseBinaryOp.h index 5e4fb147b..355f3bdc8 100644 --- a/Eigen/src/Core/CwiseBinaryOp.h +++ b/Eigen/src/Core/CwiseBinaryOp.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2009 Gael Guennebaud +// Copyright (C) 2008-2014 Gael Guennebaud // Copyright (C) 2006-2008 Benoit Jacob // // This Source Code Form is subject to the terms of the Mozilla diff --git a/Eigen/src/Core/CwiseUnaryOp.h b/Eigen/src/Core/CwiseUnaryOp.h index 098034a1e..c2bc47c93 100644 --- a/Eigen/src/Core/CwiseUnaryOp.h +++ b/Eigen/src/Core/CwiseUnaryOp.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2010 Gael Guennebaud +// Copyright (C) 2008-2014 Gael Guennebaud // Copyright (C) 2006-2008 Benoit Jacob // // This Source Code Form is subject to the terms of the Mozilla @@ -135,7 +135,7 @@ class CwiseUnaryOpImpl return derived().functor().packetOp(derived().nestedExpression().template packet(index)); } }; -#else +#else // EIGEN_TEST_EVALUATORS // Generic API dispatcher template class CwiseUnaryOpImpl @@ -144,7 +144,7 @@ class CwiseUnaryOpImpl public: typedef typename internal::generic_xpr_base >::type Base; }; -#endif +#endif // EIGEN_TEST_EVALUATORS } // end namespace Eigen diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 4b605dfd6..f5148221d 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -2,7 +2,7 @@ // for linear algebra. // // Copyright (C) 2006-2008 Benoit Jacob -// Copyright (C) 2009-2010 Gael Guennebaud +// Copyright (C) 2009-2014 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -108,6 +108,17 @@ struct TransposeImpl_base } // end namespace internal +#ifdef EIGEN_TEST_EVALUATORS +// Generic API dispatcher +template +class TransposeImpl + : public internal::generic_xpr_base >::type +{ +public: + typedef typename internal::generic_xpr_base >::type Base; +}; +#endif + template class TransposeImpl : public internal::TransposeImpl_base::type { diff --git a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h index ec86ca933..f43976641 100644 --- a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h +++ b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Gael Guennebaud +// Copyright (C) 2014 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -44,6 +44,8 @@ class sparse_cwise_binary_op_inner_iterator_selector; } // end namespace internal +#ifndef EIGEN_TEST_EVALUATORS + template class CwiseBinaryOpImpl : public SparseMatrixBase > @@ -291,6 +293,314 @@ class sparse_cwise_binary_op_inner_iterator_selector, Lhs, } // end namespace internal +#else // EIGEN_TEST_EVALUATORS + +namespace internal { + + +// Generic "sparse OP sparse" +template +struct binary_evaluator, IteratorBased, IteratorBased> + : evaluator_base > +{ +protected: + typedef typename evaluator::InnerIterator LhsIterator; + typedef typename evaluator::InnerIterator RhsIterator; +public: + typedef CwiseBinaryOp XprType; + + class ReverseInnerIterator; + class InnerIterator + { + typedef typename traits::Scalar Scalar; + typedef typename XprType::Index Index; + + public: + + EIGEN_STRONG_INLINE InnerIterator(const binary_evaluator& aEval, Index outer) + : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor) + { + this->operator++(); + } + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { + if (m_lhsIter && m_rhsIter && (m_lhsIter.index() == m_rhsIter.index())) + { + m_id = m_lhsIter.index(); + m_value = m_functor(m_lhsIter.value(), m_rhsIter.value()); + ++m_lhsIter; + ++m_rhsIter; + } + else if (m_lhsIter && (!m_rhsIter || (m_lhsIter.index() < m_rhsIter.index()))) + { + m_id = m_lhsIter.index(); + m_value = m_functor(m_lhsIter.value(), Scalar(0)); + ++m_lhsIter; + } + else if (m_rhsIter && (!m_lhsIter || (m_lhsIter.index() > m_rhsIter.index()))) + { + m_id = m_rhsIter.index(); + m_value = m_functor(Scalar(0), m_rhsIter.value()); + ++m_rhsIter; + } + else + { + m_value = 0; // this is to avoid a compilation warning + m_id = -1; + } + return *this; + } + + EIGEN_STRONG_INLINE Scalar value() const { return m_value; } + + EIGEN_STRONG_INLINE Index index() const { return m_id; } + EIGEN_STRONG_INLINE Index row() const { return Lhs::IsRowMajor ? m_lhsIter.row() : index(); } + EIGEN_STRONG_INLINE Index col() const { return Lhs::IsRowMajor ? index() : m_lhsIter.col(); } + + EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; } + + protected: + LhsIterator m_lhsIter; + RhsIterator m_rhsIter; + const BinaryOp& m_functor; + Scalar m_value; + Index m_id; + }; + + + enum { + CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, + Flags = XprType::Flags + }; + + binary_evaluator(const XprType& xpr) + : m_functor(xpr.functor()), + m_lhsImpl(xpr.lhs()), + m_rhsImpl(xpr.rhs()) + { } + +protected: + const BinaryOp m_functor; + typename evaluator::nestedType m_lhsImpl; + typename evaluator::nestedType m_rhsImpl; +}; + +// "sparse .* sparse" +template +struct binary_evaluator, Lhs, Rhs>, IteratorBased, IteratorBased> + : evaluator_base, Lhs, Rhs> > +{ +protected: + typedef scalar_product_op BinaryOp; + typedef typename evaluator::InnerIterator LhsIterator; + typedef typename evaluator::InnerIterator RhsIterator; +public: + typedef CwiseBinaryOp XprType; + + class ReverseInnerIterator; + class InnerIterator + { + typedef typename traits::Scalar Scalar; + typedef typename XprType::Index Index; + + public: + + EIGEN_STRONG_INLINE InnerIterator(const binary_evaluator& aEval, Index outer) + : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor) + { + while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index())) + { + if (m_lhsIter.index() < m_rhsIter.index()) + ++m_lhsIter; + else + ++m_rhsIter; + } + } + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { + ++m_lhsIter; + ++m_rhsIter; + while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index())) + { + if (m_lhsIter.index() < m_rhsIter.index()) + ++m_lhsIter; + else + ++m_rhsIter; + } + return *this; + } + + EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_lhsIter.value(), m_rhsIter.value()); } + + EIGEN_STRONG_INLINE Index index() const { return m_lhsIter.index(); } + EIGEN_STRONG_INLINE Index row() const { return m_lhsIter.row(); } + EIGEN_STRONG_INLINE Index col() const { return m_lhsIter.col(); } + + EIGEN_STRONG_INLINE operator bool() const { return (m_lhsIter && m_rhsIter); } + + protected: + LhsIterator m_lhsIter; + RhsIterator m_rhsIter; + const BinaryOp& m_functor; + }; + + + enum { + CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, + Flags = XprType::Flags + }; + + binary_evaluator(const XprType& xpr) + : m_functor(xpr.functor()), + m_lhsImpl(xpr.lhs()), + m_rhsImpl(xpr.rhs()) + { } + +protected: + const BinaryOp m_functor; + typename evaluator::nestedType m_lhsImpl; + typename evaluator::nestedType m_rhsImpl; +}; + +// "dense .* sparse" +template +struct binary_evaluator, Lhs, Rhs>, IndexBased, IteratorBased> + : evaluator_base, Lhs, Rhs> > +{ +protected: + typedef scalar_product_op BinaryOp; + typedef typename evaluator::type LhsEvaluator; + typedef typename evaluator::InnerIterator RhsIterator; +public: + typedef CwiseBinaryOp XprType; + + class ReverseInnerIterator; + class InnerIterator + { + typedef typename traits::Scalar Scalar; + typedef typename XprType::Index Index; + enum { IsRowMajor = (int(Rhs::Flags)&RowMajorBit)==RowMajorBit }; + + public: + + EIGEN_STRONG_INLINE InnerIterator(const binary_evaluator& aEval, Index outer) + : m_lhsEval(aEval.m_lhsImpl), m_rhsIter(aEval.m_rhsImpl,outer), m_functor(aEval.m_functor), m_outer(outer) + {} + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { + ++m_rhsIter; + return *this; + } + + EIGEN_STRONG_INLINE Scalar value() const + { return m_functor(m_lhsEval.coeff(IsRowMajor?m_outer:m_rhsIter.index(),IsRowMajor?m_rhsIter.index():m_outer), m_rhsIter.value()); } + + EIGEN_STRONG_INLINE Index index() const { return m_rhsIter.index(); } + EIGEN_STRONG_INLINE Index row() const { return m_rhsIter.row(); } + EIGEN_STRONG_INLINE Index col() const { return m_rhsIter.col(); } + + EIGEN_STRONG_INLINE operator bool() const { return m_rhsIter; } + + protected: + const LhsEvaluator &m_lhsEval; + RhsIterator m_rhsIter; + const BinaryOp& m_functor; + const Index m_outer; + }; + + + enum { + CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, + Flags = XprType::Flags + }; + + binary_evaluator(const XprType& xpr) + : m_functor(xpr.functor()), + m_lhsImpl(xpr.lhs()), + m_rhsImpl(xpr.rhs()) + { } + +protected: + const BinaryOp m_functor; + typename evaluator::nestedType m_lhsImpl; + typename evaluator::nestedType m_rhsImpl; +}; + +// "sparse .* dense" +template +struct binary_evaluator, Lhs, Rhs>, IteratorBased, IndexBased> + : evaluator_base, Lhs, Rhs> > +{ +protected: + typedef scalar_product_op BinaryOp; + typedef typename evaluator::InnerIterator LhsIterator; + typedef typename evaluator::type RhsEvaluator; +public: + typedef CwiseBinaryOp XprType; + + class ReverseInnerIterator; + class InnerIterator + { + typedef typename traits::Scalar Scalar; + typedef typename XprType::Index Index; + enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit }; + + public: + + EIGEN_STRONG_INLINE InnerIterator(const binary_evaluator& aEval, Index outer) + : m_lhsIter(aEval.m_lhsImpl,outer), m_rhsEval(aEval.m_rhsImpl), m_functor(aEval.m_functor), m_outer(outer) + {} + + EIGEN_STRONG_INLINE InnerIterator& operator++() + { + ++m_lhsIter; + return *this; + } + + EIGEN_STRONG_INLINE Scalar value() const + { return m_functor(m_lhsIter.value(), + m_rhsEval.coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); } + + EIGEN_STRONG_INLINE Index index() const { return m_lhsIter.index(); } + EIGEN_STRONG_INLINE Index row() const { return m_lhsIter.row(); } + EIGEN_STRONG_INLINE Index col() const { return m_lhsIter.col(); } + + EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter; } + + protected: + LhsIterator m_lhsIter; + const RhsEvaluator &m_rhsEval; + const BinaryOp& m_functor; + const Index m_outer; + }; + + + enum { + CoeffReadCost = evaluator::CoeffReadCost + evaluator::CoeffReadCost + functor_traits::Cost, + Flags = XprType::Flags + }; + + binary_evaluator(const XprType& xpr) + : m_functor(xpr.functor()), + m_lhsImpl(xpr.lhs()), + m_rhsImpl(xpr.rhs()) + { } + +protected: + const BinaryOp m_functor; + typename evaluator::nestedType m_lhsImpl; + typename evaluator::nestedType m_rhsImpl; +}; + +} + +#endif + + + /*************************************************************************** * Implementation of SparseMatrixBase and SparseCwise functions/operators ***************************************************************************/ diff --git a/Eigen/src/SparseCore/SparseCwiseUnaryOp.h b/Eigen/src/SparseCore/SparseCwiseUnaryOp.h index c5042504b..ead3af59d 100644 --- a/Eigen/src/SparseCore/SparseCwiseUnaryOp.h +++ b/Eigen/src/SparseCore/SparseCwiseUnaryOp.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2010 Gael Guennebaud +// Copyright (C) 2008-2014 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index aa8bfc0f5..8aa0ed7e6 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2010 Gael Guennebaud +// Copyright (C) 2008-2014 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed diff --git a/Eigen/src/SparseCore/SparseMatrixBase.h b/Eigen/src/SparseCore/SparseMatrixBase.h index a87407fe9..e4960a445 100644 --- a/Eigen/src/SparseCore/SparseMatrixBase.h +++ b/Eigen/src/SparseCore/SparseMatrixBase.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2011 Gael Guennebaud +// Copyright (C) 2008-2014 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -79,10 +79,12 @@ template class SparseMatrixBase : public EigenBase * constructed from this one. See the \ref flags "list of flags". */ +#ifndef EIGEN_TEST_EVALUATORS CoeffReadCost = internal::traits::CoeffReadCost, /**< This is a rough measure of how expensive it is to read one coefficient from * this expression. */ +#endif IsRowMajor = Flags&RowMajorBit ? 1 : 0, diff --git a/Eigen/src/SparseCore/SparseTranspose.h b/Eigen/src/SparseCore/SparseTranspose.h index 7c300ee8d..b20843dd3 100644 --- a/Eigen/src/SparseCore/SparseTranspose.h +++ b/Eigen/src/SparseCore/SparseTranspose.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2009 Gael Guennebaud +// Copyright (C) 2008-2014 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -12,6 +12,7 @@ namespace Eigen { +#ifndef EIGEN_TEST_EVALUATORS template class TransposeImpl : public SparseMatrixBase > { @@ -58,6 +59,57 @@ template class TransposeImpl::ReverseInn Index col() const { return Base::row(); } }; +#else // EIGEN_TEST_EVALUATORS + +namespace internal { + +template +struct unary_evaluator, IteratorBased> + : public evaluator_base > +{ + typedef typename evaluator::InnerIterator EvalIterator; + typedef typename evaluator::ReverseInnerIterator EvalReverseIterator; + public: + typedef Transpose XprType; + typedef typename XprType::Index Index; + + class InnerIterator : public EvalIterator + { + public: + EIGEN_STRONG_INLINE InnerIterator(const unary_evaluator& unaryOp, typename XprType::Index outer) + : EvalIterator(unaryOp.m_argImpl,outer) + {} + + Index row() const { return EvalIterator::col(); } + Index col() const { return EvalIterator::row(); } + }; + + class ReverseInnerIterator : public EvalReverseIterator + { + public: + EIGEN_STRONG_INLINE ReverseInnerIterator(const unary_evaluator& unaryOp, typename XprType::Index outer) + : EvalReverseIterator(unaryOp.m_argImpl,outer) + {} + + Index row() const { return EvalReverseIterator::col(); } + Index col() const { return EvalReverseIterator::row(); } + }; + + enum { + CoeffReadCost = evaluator::CoeffReadCost, + Flags = XprType::Flags + }; + + unary_evaluator(const XprType& op) :m_argImpl(op.nestedExpression()) {} + + protected: + typename evaluator::nestedType m_argImpl; +}; + +} // end namespace internal + +#endif // EIGEN_TEST_EVALUATORS + } // end namespace Eigen #endif // EIGEN_SPARSETRANSPOSE_H