From cd0ff253ec906f0b4604ec084df3d7e7cc6811d1 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 30 Jul 2014 15:22:50 +0200 Subject: Make permutation compatible with sparse matrices --- CMakeLists.txt | 2 +- Eigen/SparseCore | 3 - Eigen/src/Core/PermutationMatrix.h | 9 +- Eigen/src/Core/util/Constants.h | 3 + Eigen/src/Core/util/XprHelper.h | 24 +++-- Eigen/src/SparseCore/SparseCwiseUnaryOp.h | 92 +++++++++---------- Eigen/src/SparseCore/SparseMatrix.h | 8 +- Eigen/src/SparseCore/SparsePermutation.h | 131 ++++++++++++++++++++++++++- Eigen/src/SparseCore/SparseSelfAdjointView.h | 12 +-- test/CMakeLists.txt | 2 +- 10 files changed, 213 insertions(+), 73 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 470095680..cc3988196 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -143,7 +143,7 @@ if(NOT MSVC) ei_add_cxx_compiler_flag("-Wpointer-arith") ei_add_cxx_compiler_flag("-Wwrite-strings") ei_add_cxx_compiler_flag("-Wformat-security") - ei_add_cxx_compiler_flag("-Wshorten-64-to-32") +# ei_add_cxx_compiler_flag("-Wshorten-64-to-32") ei_add_cxx_compiler_flag("-Wenum-conversion") ei_add_cxx_compiler_flag("-Wc++11-extensions") diff --git a/Eigen/SparseCore b/Eigen/SparseCore index 578469c1c..41cae58b0 100644 --- a/Eigen/SparseCore +++ b/Eigen/SparseCore @@ -56,11 +56,8 @@ struct Sparse {}; #include "src/SparseCore/SparseSelfAdjointView.h" #include "src/SparseCore/SparseTriangularView.h" #include "src/SparseCore/TriangularSolver.h" - -#ifndef EIGEN_TEST_EVALUATORS #include "src/SparseCore/SparsePermutation.h" #include "src/SparseCore/SparseFuzzy.h" -#endif #include "src/Core/util/ReenableStupidWarnings.h" diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h index 5d648c68c..84d64edb3 100644 --- a/Eigen/src/Core/PermutationMatrix.h +++ b/Eigen/src/Core/PermutationMatrix.h @@ -13,7 +13,8 @@ namespace Eigen { -template class PermutedImpl; +// TODO: this does not seems to be needed at all: +// template class PermutedImpl; /** \class PermutationBase * \ingroup Core_Module @@ -276,6 +277,7 @@ template > : traits > { + typedef PermutationStorage StorageKind; typedef Matrix<_StorageIndexType, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType; typedef typename IndicesType::Index Index; typedef _StorageIndexType StorageIndexType; @@ -397,6 +399,7 @@ template,_PacketAccess> > : traits > { + typedef PermutationStorage StorageKind; typedef Map, _PacketAccess> IndicesType; typedef typename IndicesType::Index Index; typedef _StorageIndexType StorageIndexType; @@ -468,8 +471,6 @@ class Map class TranspositionsWrapper; namespace internal { template @@ -665,6 +666,8 @@ struct traits > > } // end namespace internal +// TODO: the specificties should be handled by the evaluator, +// at the very least we should only specialize TransposeImpl template class Transpose > : public EigenBase > > diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 3ab8d0ed3..c0463583d 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -443,6 +443,9 @@ enum Action {GetAction, SetAction}; /** The type used to identify a dense storage. */ struct Dense {}; +/** The type used to identify a permutation storage. */ +struct PermutationStorage {}; + /** The type used to identify a matrix expression */ struct MatrixXpr {}; diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 7779fb3fa..5ed3b39b7 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -557,18 +557,26 @@ template struct cwise_promote_s * K * dense -> dense * diag * K -> K * K * diag -> K + * Perm * K -> K + * K * Perm -> K * \endcode */ template struct product_promote_storage_type; -template struct product_promote_storage_type { typedef A ret;}; -template struct product_promote_storage_type { typedef Dense ret;}; -template struct product_promote_storage_type { typedef Dense ret; }; -template struct product_promote_storage_type { typedef Dense ret; }; -template struct product_promote_storage_type { typedef A ret; }; -template struct product_promote_storage_type { typedef B ret; }; -template struct product_promote_storage_type { typedef Dense ret; }; -template struct product_promote_storage_type { typedef Dense ret; }; +template struct product_promote_storage_type { typedef A ret;}; +template struct product_promote_storage_type { typedef Dense ret;}; +template struct product_promote_storage_type { typedef Dense ret; }; +template struct product_promote_storage_type { typedef Dense ret; }; + +template struct product_promote_storage_type { typedef A ret; }; +template struct product_promote_storage_type { typedef B ret; }; +template struct product_promote_storage_type { typedef Dense ret; }; +template struct product_promote_storage_type { typedef Dense ret; }; + +template struct product_promote_storage_type { typedef A ret; }; +template struct product_promote_storage_type { typedef B ret; }; +template struct product_promote_storage_type { typedef Dense ret; }; +template struct product_promote_storage_type { typedef Dense ret; }; /** \internal gives the plain matrix or array type to store a row/column/diagonal of a matrix type. * \param Scalar optional parameter allowing to pass a different scalar type than the one of the MatrixType. diff --git a/Eigen/src/SparseCore/SparseCwiseUnaryOp.h b/Eigen/src/SparseCore/SparseCwiseUnaryOp.h index 7a849c822..0e0f9f5f5 100644 --- a/Eigen/src/SparseCore/SparseCwiseUnaryOp.h +++ b/Eigen/src/SparseCore/SparseCwiseUnaryOp.h @@ -152,7 +152,7 @@ struct unary_evaluator, IteratorBased> typedef CwiseUnaryOp XprType; class InnerIterator; - class ReverseInnerIterator; +// class ReverseInnerIterator; enum { CoeffReadCost = evaluator::CoeffReadCost + functor_traits::Cost, @@ -163,7 +163,7 @@ struct unary_evaluator, IteratorBased> protected: typedef typename evaluator::InnerIterator EvalIterator; - typedef typename evaluator::ReverseInnerIterator EvalReverseIterator; +// typedef typename evaluator::ReverseInnerIterator EvalReverseIterator; const UnaryOp m_functor; typename evaluator::nestedType m_argImpl; @@ -192,28 +192,28 @@ class unary_evaluator, IteratorBased>::InnerIterat Scalar& valueRef(); }; -template -class unary_evaluator, IteratorBased>::ReverseInnerIterator - : public unary_evaluator, IteratorBased>::EvalReverseIterator -{ - typedef typename XprType::Scalar Scalar; - typedef typename unary_evaluator, IteratorBased>::EvalReverseIterator Base; - public: - - EIGEN_STRONG_INLINE ReverseInnerIterator(const XprType& unaryOp, typename XprType::Index outer) - : Base(unaryOp.derived().nestedExpression(),outer), m_functor(unaryOp.derived().functor()) - {} - - EIGEN_STRONG_INLINE ReverseInnerIterator& operator--() - { Base::operator--(); return *this; } - - EIGEN_STRONG_INLINE Scalar value() const { return m_functor(Base::value()); } - - protected: - const UnaryOp m_functor; - private: - Scalar& valueRef(); -}; +// template +// class unary_evaluator, IteratorBased>::ReverseInnerIterator +// : public unary_evaluator, IteratorBased>::EvalReverseIterator +// { +// typedef typename XprType::Scalar Scalar; +// typedef typename unary_evaluator, IteratorBased>::EvalReverseIterator Base; +// public: +// +// EIGEN_STRONG_INLINE ReverseInnerIterator(const XprType& unaryOp, typename XprType::Index outer) +// : Base(unaryOp.derived().nestedExpression(),outer), m_functor(unaryOp.derived().functor()) +// {} +// +// EIGEN_STRONG_INLINE ReverseInnerIterator& operator--() +// { Base::operator--(); return *this; } +// +// EIGEN_STRONG_INLINE Scalar value() const { return m_functor(Base::value()); } +// +// protected: +// const UnaryOp m_functor; +// private: +// Scalar& valueRef(); +// }; @@ -238,7 +238,7 @@ struct unary_evaluator, IteratorBased> protected: typedef typename evaluator::InnerIterator EvalIterator; - typedef typename evaluator::ReverseInnerIterator EvalReverseIterator; +// typedef typename evaluator::ReverseInnerIterator EvalReverseIterator; const ViewOp m_functor; typename evaluator::nestedType m_argImpl; @@ -266,27 +266,27 @@ class unary_evaluator, IteratorBased>::InnerItera const ViewOp m_functor; }; -template -class unary_evaluator, IteratorBased>::ReverseInnerIterator - : public unary_evaluator, IteratorBased>::EvalReverseIterator -{ - typedef typename XprType::Scalar Scalar; - typedef typename unary_evaluator, IteratorBased>::EvalReverseIterator Base; - public: - - EIGEN_STRONG_INLINE ReverseInnerIterator(const XprType& unaryOp, typename XprType::Index outer) - : Base(unaryOp.derived().nestedExpression(),outer), m_functor(unaryOp.derived().functor()) - {} - - EIGEN_STRONG_INLINE ReverseInnerIterator& operator--() - { Base::operator--(); return *this; } - - EIGEN_STRONG_INLINE Scalar value() const { return m_functor(Base::value()); } - EIGEN_STRONG_INLINE Scalar& valueRef() { return m_functor(Base::valueRef()); } - - protected: - const ViewOp m_functor; -}; +// template +// class unary_evaluator, IteratorBased>::ReverseInnerIterator +// : public unary_evaluator, IteratorBased>::EvalReverseIterator +// { +// typedef typename XprType::Scalar Scalar; +// typedef typename unary_evaluator, IteratorBased>::EvalReverseIterator Base; +// public: +// +// EIGEN_STRONG_INLINE ReverseInnerIterator(const XprType& unaryOp, typename XprType::Index outer) +// : Base(unaryOp.derived().nestedExpression(),outer), m_functor(unaryOp.derived().functor()) +// {} +// +// EIGEN_STRONG_INLINE ReverseInnerIterator& operator--() +// { Base::operator--(); return *this; } +// +// EIGEN_STRONG_INLINE Scalar value() const { return m_functor(Base::value()); } +// EIGEN_STRONG_INLINE Scalar& valueRef() { return m_functor(Base::valueRef()); } +// +// protected: +// const ViewOp m_functor; +// }; } // end namespace internal diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index 6080c272a..9e2dae1b3 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -736,8 +736,8 @@ class SparseMatrix return *this; } +#ifndef EIGEN_PARSED_BY_DOXYGEN #ifndef EIGEN_TEST_EVALUATORS - #ifndef EIGEN_PARSED_BY_DOXYGEN template inline SparseMatrix& operator=(const SparseSparseProduct& product) { return Base::operator=(product); } @@ -748,12 +748,12 @@ class SparseMatrix initAssignment(other); return Base::operator=(other.derived()); } - +#endif // EIGEN_TEST_EVALUATORS + template inline SparseMatrix& operator=(const EigenBase& other) { return Base::operator=(other.derived()); } - #endif -#endif // EIGEN_TEST_EVALUATORS +#endif // EIGEN_PARSED_BY_DOXYGEN template EIGEN_DONT_INLINE SparseMatrix& operator=(const SparseMatrixBase& other); diff --git a/Eigen/src/SparseCore/SparsePermutation.h b/Eigen/src/SparseCore/SparsePermutation.h index b85be93f6..ebfefab98 100644 --- a/Eigen/src/SparseCore/SparsePermutation.h +++ b/Eigen/src/SparseCore/SparsePermutation.h @@ -103,7 +103,7 @@ struct permut_sparsematrix_product_retval } - +#ifndef EIGEN_TEST_EVALUATORS /** \returns the matrix with the permutation applied to the columns */ @@ -143,6 +143,135 @@ operator*(const Transpose >& tperm, const SparseMat return internal::permut_sparsematrix_product_retval, SparseDerived, OnTheLeft, true>(tperm.nestedPermutation(), matrix.derived()); } +#else // EIGEN_TEST_EVALUATORS + +namespace internal { + +template struct product_promote_storage_type { typedef Sparse ret; }; +template struct product_promote_storage_type { typedef Sparse ret; }; + +// TODO, the following need cleaning, this is just a copy-past of the dense case + +template +struct generic_product_impl +{ + template + static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) + { + permut_sparsematrix_product_retval pmpr(lhs, rhs); + pmpr.evalTo(dst); + } +}; + +template +struct generic_product_impl +{ + template + static void evalTo(Dest& dst, const Lhs& lhs, const Rhs& rhs) + { + permut_sparsematrix_product_retval pmpr(rhs, lhs); + pmpr.evalTo(dst); + } +}; + +template +struct generic_product_impl, Rhs, PermutationShape, SparseShape, ProductTag> +{ + template + static void evalTo(Dest& dst, const Transpose& lhs, const Rhs& rhs) + { + permut_sparsematrix_product_retval pmpr(lhs.nestedPermutation(), rhs); + pmpr.evalTo(dst); + } +}; + +template +struct generic_product_impl, SparseShape, PermutationShape, ProductTag> +{ + template + static void evalTo(Dest& dst, const Lhs& lhs, const Transpose& rhs) + { + permut_sparsematrix_product_retval pmpr(rhs.nestedPermutation(), lhs); + pmpr.evalTo(dst); + } +}; + +template +struct product_evaluator, ProductTag, PermutationShape, SparseShape, typename traits::Scalar, typename traits::Scalar> + : public evaluator >::ReturnType>::type +{ + typedef Product XprType; + typedef typename traits >::ReturnType PlainObject; + typedef typename evaluator::type Base; + + product_evaluator(const XprType& xpr) + : m_result(xpr.rows(), xpr.cols()) + { + ::new (static_cast(this)) Base(m_result); + generic_product_impl::evalTo(m_result, xpr.lhs(), xpr.rhs()); + } + +protected: + PlainObject m_result; +}; + +template +struct product_evaluator, ProductTag, SparseShape, PermutationShape, typename traits::Scalar, typename traits::Scalar> + : public evaluator >::ReturnType>::type +{ + typedef Product XprType; + typedef typename traits >::ReturnType PlainObject; + typedef typename evaluator::type Base; + + product_evaluator(const XprType& xpr) + : m_result(xpr.rows(), xpr.cols()) + { + ::new (static_cast(this)) Base(m_result); + generic_product_impl::evalTo(m_result, xpr.lhs(), xpr.rhs()); + } + +protected: + PlainObject m_result; +}; + +} // end namespace internal + +/** \returns the matrix with the permutation applied to the columns + */ +template +inline const Product +operator*(const SparseMatrixBase& matrix, const PermutationBase& perm) +{ return Product(matrix.derived(), perm.derived()); } + +/** \returns the matrix with the permutation applied to the rows + */ +template +inline const Product +operator*( const PermutationBase& perm, const SparseMatrixBase& matrix) +{ return Product(perm.derived(), matrix.derived()); } + + +// TODO, the following specializations should not be needed as Transpose should be a PermutationBase. +/** \returns the matrix with the inverse permutation applied to the columns. + */ +template +inline const Product > > +operator*(const SparseMatrixBase& matrix, const Transpose >& tperm) +{ + return Product > >(matrix.derived(), tperm); +} + +/** \returns the matrix with the inverse permutation applied to the rows. + */ +template +inline const Product >, SparseDerived> +operator*(const Transpose >& tperm, const SparseMatrixBase& matrix) +{ + return Product >, SparseDerived>(tperm, matrix.derived()); +} + +#endif // EIGEN_TEST_EVALUATORS + } // end namespace Eigen #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index 10cd755a4..530ff27bf 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -20,7 +20,7 @@ class SparseSelfAdjointTimeDenseProduct; template class DenseTimeSparseSelfAdjointProduct; -#endif // #ifndef EIGEN_TEST_EVALUATORS +#endif // EIGEN_TEST_EVALUATORS /** \ingroup SparseCore_Module * \class SparseSelfAdjointView @@ -177,7 +177,7 @@ template class SparseSelfAdjointView } /** \returns an expression of P H P^-1 */ -#ifndef EIGEN_TEST_EVALUATORS +// #ifndef EIGEN_TEST_EVALUATORS SparseSymmetricPermutationProduct<_MatrixTypeNested,Mode> twistedBy(const PermutationMatrix& perm) const { return SparseSymmetricPermutationProduct<_MatrixTypeNested,Mode>(m_matrix, perm); @@ -189,7 +189,7 @@ template class SparseSelfAdjointView permutedMatrix.evalTo(*this); return *this; } -#endif // EIGEN_TEST_EVALUATORS +// #endif // EIGEN_TEST_EVALUATORS SparseSelfAdjointView& operator=(const SparseSelfAdjointView& src) { @@ -680,7 +680,7 @@ void permute_symm_to_symm(const MatrixType& mat, SparseMatrix