From 2a251ffab0f3abd1fcfe340a4bee61e352d83424 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 22 Jul 2014 09:32:40 +0200 Subject: [PATCH] Implement evaluator for sparse-selfadjoint products --- Eigen/SparseCore | 3 +- Eigen/src/Core/GeneralProduct.h | 16 +- Eigen/src/Core/Product.h | 24 +- Eigen/src/Core/SelfAdjointView.h | 1 + Eigen/src/SparseCore/SparseAssign.h | 3 + Eigen/src/SparseCore/SparseDenseProduct.h | 39 -- Eigen/src/SparseCore/SparseMatrix.h | 4 + Eigen/src/SparseCore/SparseSelfAdjointView.h | 397 +++++++++++++++---- test/sparse_product.cpp | 10 +- 9 files changed, 352 insertions(+), 145 deletions(-) diff --git a/Eigen/SparseCore b/Eigen/SparseCore index 7cbfb47f2..69b413ec2 100644 --- a/Eigen/SparseCore +++ b/Eigen/SparseCore @@ -53,11 +53,12 @@ struct Sparse {}; #include "src/SparseCore/SparseSparseProductWithPruning.h" #include "src/SparseCore/SparseProduct.h" #include "src/SparseCore/SparseDenseProduct.h" +#include "src/SparseCore/SparseSelfAdjointView.h" #ifndef EIGEN_TEST_EVALUATORS #include "src/SparseCore/SparsePermutation.h" #include "src/SparseCore/SparseFuzzy.h" #include "src/SparseCore/SparseTriangularView.h" -#include "src/SparseCore/SparseSelfAdjointView.h" + #include "src/SparseCore/TriangularSolver.h" #endif diff --git a/Eigen/src/Core/GeneralProduct.h b/Eigen/src/Core/GeneralProduct.h index 4c1922741..76271a87d 100644 --- a/Eigen/src/Core/GeneralProduct.h +++ b/Eigen/src/Core/GeneralProduct.h @@ -62,14 +62,14 @@ template struct product_type typedef typename remove_all::type _Lhs; typedef typename remove_all::type _Rhs; enum { - MaxRows = _Lhs::MaxRowsAtCompileTime, - Rows = _Lhs::RowsAtCompileTime, - MaxCols = _Rhs::MaxColsAtCompileTime, - Cols = _Rhs::ColsAtCompileTime, - MaxDepth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::MaxColsAtCompileTime, - _Rhs::MaxRowsAtCompileTime), - Depth = EIGEN_SIZE_MIN_PREFER_FIXED(_Lhs::ColsAtCompileTime, - _Rhs::RowsAtCompileTime) + MaxRows = traits<_Lhs>::MaxRowsAtCompileTime, + Rows = traits<_Lhs>::RowsAtCompileTime, + MaxCols = traits<_Rhs>::MaxColsAtCompileTime, + Cols = traits<_Rhs>::ColsAtCompileTime, + MaxDepth = EIGEN_SIZE_MIN_PREFER_FIXED(traits<_Lhs>::MaxColsAtCompileTime, + traits<_Rhs>::MaxRowsAtCompileTime), + Depth = EIGEN_SIZE_MIN_PREFER_FIXED(traits<_Lhs>::ColsAtCompileTime, + traits<_Rhs>::RowsAtCompileTime) }; // the splitting into different lines of code here, introducing the _select enums and the typedef below, diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index cb4d4c924..9e5e47d13 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -58,24 +58,26 @@ struct traits > { typedef typename remove_all::type LhsCleaned; typedef typename remove_all::type RhsCleaned; + typedef traits LhsTraits; + typedef traits RhsTraits; typedef MatrixXpr XprKind; typedef typename product_result_scalar::Scalar Scalar; - typedef typename product_promote_storage_type::StorageKind, - typename traits::StorageKind, + typedef typename product_promote_storage_type::ret>::ret StorageKind; - typedef typename promote_index_type::Index, - typename traits::Index>::type Index; + typedef typename promote_index_type::type Index; enum { - RowsAtCompileTime = LhsCleaned::RowsAtCompileTime, - ColsAtCompileTime = RhsCleaned::ColsAtCompileTime, - MaxRowsAtCompileTime = LhsCleaned::MaxRowsAtCompileTime, - MaxColsAtCompileTime = RhsCleaned::MaxColsAtCompileTime, + RowsAtCompileTime = LhsTraits::RowsAtCompileTime, + ColsAtCompileTime = RhsTraits::ColsAtCompileTime, + MaxRowsAtCompileTime = LhsTraits::MaxRowsAtCompileTime, + MaxColsAtCompileTime = RhsTraits::MaxColsAtCompileTime, // FIXME: only needed by GeneralMatrixMatrixTriangular - InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(LhsCleaned::ColsAtCompileTime, RhsCleaned::RowsAtCompileTime), + InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(LhsTraits::ColsAtCompileTime, RhsTraits::RowsAtCompileTime), #ifndef EIGEN_TEST_EVALUATORS // dummy, for evaluators unit test only @@ -84,8 +86,8 @@ struct traits > // The storage order is somewhat arbitrary here. The correct one will be determined through the evaluator. Flags = ( MaxRowsAtCompileTime==1 - || ((LhsCleaned::Flags&NoPreferredStorageOrderBit) && (RhsCleaned::Flags&RowMajorBit)) - || ((RhsCleaned::Flags&NoPreferredStorageOrderBit) && (LhsCleaned::Flags&RowMajorBit)) ) + || ((LhsTraits::Flags&NoPreferredStorageOrderBit) && (RhsTraits::Flags&RowMajorBit)) + || ((RhsTraits::Flags&NoPreferredStorageOrderBit) && (LhsTraits::Flags&RowMajorBit)) ) ? RowMajorBit : (MaxColsAtCompileTime==1 ? 0 : NoPreferredStorageOrderBit) }; }; diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index 9319ad87c..c8bdec5d4 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -328,6 +328,7 @@ struct triangular_assignment_selector // in the future selfadjoint-ness should be defined by the expression traits // such that Transpose > is valid. (currently TriangularBase::transpose() is overloaded to make it work) diff --git a/Eigen/src/SparseCore/SparseAssign.h b/Eigen/src/SparseCore/SparseAssign.h index 528d63e35..f3da8c6c4 100644 --- a/Eigen/src/SparseCore/SparseAssign.h +++ b/Eigen/src/SparseCore/SparseAssign.h @@ -127,6 +127,7 @@ template template Derived& SparseMatrixBase::operator=(const EigenBase &other) { + // TODO use the evaluator mechanism other.derived().evalTo(derived()); return derived(); } @@ -135,6 +136,7 @@ template template Derived& SparseMatrixBase::operator=(const ReturnByValue& other) { + // TODO use the evaluator mechanism other.evalTo(derived()); return derived(); } @@ -143,6 +145,7 @@ template template inline Derived& SparseMatrixBase::operator=(const SparseMatrixBase& other) { + // FIXME, by default sparse evaluation do not alias, so we should be able to bypass the generic call_assignment internal::call_assignment/*_no_alias*/(derived(), other.derived()); return derived(); } diff --git a/Eigen/src/SparseCore/SparseDenseProduct.h b/Eigen/src/SparseCore/SparseDenseProduct.h index 883e24acb..8864b7308 100644 --- a/Eigen/src/SparseCore/SparseDenseProduct.h +++ b/Eigen/src/SparseCore/SparseDenseProduct.h @@ -448,45 +448,6 @@ protected: PlainObject m_result; }; - -// template -// class sparse_dense_outer_product_iterator : public LhsIterator -// { -// typedef typename SparseDenseOuterProduct::Index Index; -// public: -// template -// EIGEN_STRONG_INLINE InnerIterator(const XprEval& prod, Index outer) -// : LhsIterator(prod.lhs(), 0), -// m_outer(outer), m_empty(false), m_factor(get(prod.rhs(), outer, typename internal::traits::StorageKind() )) -// {} -// -// inline Index outer() const { return m_outer; } -// inline Index row() const { return Transpose ? m_outer : Base::index(); } -// inline Index col() const { return Transpose ? Base::index() : m_outer; } -// -// inline Scalar value() const { return Base::value() * m_factor; } -// inline operator bool() const { return Base::operator bool() && !m_empty; } -// -// protected: -// Scalar get(const _RhsNested &rhs, Index outer, Dense = Dense()) const -// { -// return rhs.coeff(outer); -// } -// -// Scalar get(const _RhsNested &rhs, Index outer, Sparse = Sparse()) -// { -// typename Traits::_RhsNested::InnerIterator it(rhs, outer); -// if (it && it.index()==0 && it.value()!=Scalar(0)) -// return it.value(); -// m_empty = true; -// return Scalar(0); -// } -// -// Index m_outer; -// bool m_empty; -// Scalar m_factor; -// }; - template struct sparse_dense_outer_product_evaluator { diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index 36a024cc7..5f0c3d0a7 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -664,7 +664,11 @@ class SparseMatrix : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0) { check_template_parameters(); +#ifndef EIGEN_TEST_EVALUATORS *this = other; +#else + Base::operator=(other); +#endif } /** Copy constructor (it performs a deep copy) */ diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index 56c922929..8bd836883 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2009 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 @@ -12,13 +12,23 @@ namespace Eigen { +#ifndef EIGEN_TEST_EVALUATORS + +template +class SparseSelfAdjointTimeDenseProduct; + +template +class DenseTimeSparseSelfAdjointProduct; + +#endif // #ifndef EIGEN_TEST_EVALUATORS + /** \ingroup SparseCore_Module * \class SparseSelfAdjointView * * \brief Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix. * * \param MatrixType the type of the dense matrix storing the coefficients - * \param UpLo can be either \c #Lower or \c #Upper + * \param Mode can be either \c #Lower or \c #Upper * * This class is an expression of a sefladjoint matrix from a triangular part of a matrix * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView() @@ -26,37 +36,33 @@ namespace Eigen { * * \sa SparseMatrixBase::selfadjointView() */ -template -class SparseSelfAdjointTimeDenseProduct; - -template -class DenseTimeSparseSelfAdjointProduct; - namespace internal { -template -struct traits > : traits { +template +struct traits > : traits { }; -template +template void permute_symm_to_symm(const MatrixType& mat, SparseMatrix& _dest, const typename MatrixType::Index* perm = 0); -template +template void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix& _dest, const typename MatrixType::Index* perm = 0); } -template class SparseSelfAdjointView - : public EigenBase > +template class SparseSelfAdjointView + : public EigenBase > { public: + + enum { Mode = _Mode }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Index Index; typedef Matrix VectorI; typedef typename MatrixType::Nested MatrixTypeNested; typedef typename internal::remove_all::type _MatrixTypeNested; - + inline SparseSelfAdjointView(const MatrixType& matrix) : m_matrix(matrix) { eigen_assert(rows()==cols() && "SelfAdjointView is only for squared matrices"); @@ -74,40 +80,76 @@ template class SparseSelfAdjointView * Note that there is no algorithmic advantage of performing such a product compared to a general sparse-sparse matrix product. * Indeed, the SparseSelfadjointView operand is first copied into a temporary SparseMatrix before computing the product. */ +#ifndef EIGEN_TEST_EVALUATORS template SparseSparseProduct operator*(const SparseMatrixBase& rhs) const { return SparseSparseProduct(*this, rhs.derived()); } +#else + template + Product + operator*(const SparseMatrixBase& rhs) const + { + return Product(*this, rhs.derived()); + } +#endif /** \returns an expression of the matrix product between a sparse matrix \a lhs and a sparse self-adjoint matrix \a rhs. * * Note that there is no algorithmic advantage of performing such a product compared to a general sparse-sparse matrix product. * Indeed, the SparseSelfadjointView operand is first copied into a temporary SparseMatrix before computing the product. */ +#ifndef EIGEN_TEST_EVALUATORS template friend SparseSparseProduct operator*(const SparseMatrixBase& lhs, const SparseSelfAdjointView& rhs) { return SparseSparseProduct(lhs.derived(), rhs); } +#else // EIGEN_TEST_EVALUATORS + template friend + Product + operator*(const SparseMatrixBase& lhs, const SparseSelfAdjointView& rhs) + { + return Product(lhs.derived(), rhs); + } +#endif // EIGEN_TEST_EVALUATORS /** Efficient sparse self-adjoint matrix times dense vector/matrix product */ +#ifndef EIGEN_TEST_EVALUATORS template - SparseSelfAdjointTimeDenseProduct + SparseSelfAdjointTimeDenseProduct operator*(const MatrixBase& rhs) const { - return SparseSelfAdjointTimeDenseProduct(m_matrix, rhs.derived()); + return SparseSelfAdjointTimeDenseProduct(m_matrix, rhs.derived()); } +#else + template + Product + operator*(const MatrixBase& rhs) const + { + return Product(*this, rhs.derived()); + } +#endif /** Efficient dense vector/matrix times sparse self-adjoint matrix product */ +#ifndef EIGEN_TEST_EVALUATORS template friend - DenseTimeSparseSelfAdjointProduct + DenseTimeSparseSelfAdjointProduct operator*(const MatrixBase& lhs, const SparseSelfAdjointView& rhs) { - return DenseTimeSparseSelfAdjointProduct(lhs.derived(), rhs.m_matrix); + return DenseTimeSparseSelfAdjointProduct(lhs.derived(), rhs.m_matrix); } +#else + template friend + Product + operator*(const MatrixBase& lhs, const SparseSelfAdjointView& rhs) + { + return Product(lhs.derived(), rhs); + } +#endif /** Perform a symmetric rank K update of the selfadjoint matrix \c *this: * \f$ this = this + \alpha ( u u^* ) \f$ where \a u is a vector or matrix. @@ -123,30 +165,31 @@ template class SparseSelfAdjointView /** \internal triggered by sparse_matrix = SparseSelfadjointView; */ template void evalTo(SparseMatrix& _dest) const { - internal::permute_symm_to_fullsymm(m_matrix, _dest); + internal::permute_symm_to_fullsymm(m_matrix, _dest); } template void evalTo(DynamicSparseMatrix& _dest) const { // TODO directly evaluate into _dest; SparseMatrix tmp(_dest.rows(),_dest.cols()); - internal::permute_symm_to_fullsymm(m_matrix, tmp); + internal::permute_symm_to_fullsymm(m_matrix, tmp); _dest = tmp; } /** \returns an expression of P H P^-1 */ - SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo> twistedBy(const PermutationMatrix& perm) const +#ifndef EIGEN_TEST_EVALUATORS + SparseSymmetricPermutationProduct<_MatrixTypeNested,Mode> twistedBy(const PermutationMatrix& perm) const { - return SparseSymmetricPermutationProduct<_MatrixTypeNested,UpLo>(m_matrix, perm); + return SparseSymmetricPermutationProduct<_MatrixTypeNested,Mode>(m_matrix, perm); } - - template - SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct& permutedMatrix) + + template + SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct& permutedMatrix) { permutedMatrix.evalTo(*this); return *this; } - +#endif // EIGEN_TEST_EVALUATORS SparseSelfAdjointView& operator=(const SparseSelfAdjointView& src) { @@ -154,22 +197,18 @@ template class SparseSelfAdjointView return *this = src.twistedBy(pnull); } - template - SparseSelfAdjointView& operator=(const SparseSelfAdjointView& src) + template + SparseSelfAdjointView& operator=(const SparseSelfAdjointView& src) { PermutationMatrix pnull; return *this = src.twistedBy(pnull); } - - // const SparseLLT llt() const; - // const SparseLDLT ldlt() const; - protected: typename MatrixType::Nested m_matrix; - mutable VectorI m_countPerRow; - mutable VectorI m_countPerCol; + //mutable VectorI m_countPerRow; + //mutable VectorI m_countPerCol; }; /*************************************************************************** @@ -177,15 +216,15 @@ template class SparseSelfAdjointView ***************************************************************************/ template -template -const SparseSelfAdjointView SparseMatrixBase::selfadjointView() const +template +const SparseSelfAdjointView SparseMatrixBase::selfadjointView() const { return derived(); } template -template -SparseSelfAdjointView SparseMatrixBase::selfadjointView() +template +SparseSelfAdjointView SparseMatrixBase::selfadjointView() { return derived(); } @@ -194,16 +233,16 @@ SparseSelfAdjointView SparseMatrixBase::selfadjointView( * Implementation of SparseSelfAdjointView methods ***************************************************************************/ -template +template template -SparseSelfAdjointView& -SparseSelfAdjointView::rankUpdate(const SparseMatrixBase& u, const Scalar& alpha) +SparseSelfAdjointView& +SparseSelfAdjointView::rankUpdate(const SparseMatrixBase& u, const Scalar& alpha) { SparseMatrix tmp = u * u.adjoint(); if(alpha==Scalar(0)) - m_matrix.const_cast_derived() = tmp.template triangularView(); + m_matrix.const_cast_derived() = tmp.template triangularView(); else - m_matrix.const_cast_derived() += alpha * tmp.template triangularView(); + m_matrix.const_cast_derived() += alpha * tmp.template triangularView(); return *this; } @@ -212,18 +251,19 @@ SparseSelfAdjointView::rankUpdate(const SparseMatrixBase -struct traits > - : traits, Lhs, Rhs> > +template +struct traits > + : traits, Lhs, Rhs> > { typedef Dense StorageKind; }; } -template +template class SparseSelfAdjointTimeDenseProduct - : public ProductBase, Lhs, Rhs> + : public ProductBase, Lhs, Rhs> { public: EIGEN_PRODUCT_PUBLIC_INTERFACE(SparseSelfAdjointTimeDenseProduct) @@ -241,9 +281,9 @@ class SparseSelfAdjointTimeDenseProduct enum { LhsIsRowMajor = (_Lhs::Flags&RowMajorBit)==RowMajorBit, ProcessFirstHalf = - ((UpLo&(Upper|Lower))==(Upper|Lower)) - || ( (UpLo&Upper) && !LhsIsRowMajor) - || ( (UpLo&Lower) && LhsIsRowMajor), + ((Mode&(Upper|Lower))==(Upper|Lower)) + || ( (Mode&Upper) && !LhsIsRowMajor) + || ( (Mode&Lower) && LhsIsRowMajor), ProcessSecondHalf = !ProcessFirstHalf }; for (typename _Lhs::Index j=0; j -struct traits > - : traits, Lhs, Rhs> > +template +struct traits > + : traits, Lhs, Rhs> > {}; } -template +template class DenseTimeSparseSelfAdjointProduct - : public ProductBase, Lhs, Rhs> + : public ProductBase, Lhs, Rhs> { public: EIGEN_PRODUCT_PUBLIC_INTERFACE(DenseTimeSparseSelfAdjointProduct) @@ -301,16 +341,197 @@ class DenseTimeSparseSelfAdjointProduct DenseTimeSparseSelfAdjointProduct& operator=(const DenseTimeSparseSelfAdjointProduct&); }; +#else // EIGEN_TEST_EVALUATORS + +namespace internal { + +template +inline void sparse_selfadjoint_time_dense_product(const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res, const AlphaType& alpha) +{ + EIGEN_ONLY_USED_FOR_DEBUG(alpha); + // TODO use alpha + eigen_assert(alpha==AlphaType(1) && "alpha != 1 is not implemented yet, sorry"); + + typedef typename evaluator::type LhsEval; + typedef typename evaluator::InnerIterator LhsIterator; + typedef typename SparseLhsType::Index Index; + typedef typename SparseLhsType::Scalar LhsScalar; + + enum { + LhsIsRowMajor = (LhsEval::Flags&RowMajorBit)==RowMajorBit, + ProcessFirstHalf = + ((Mode&(Upper|Lower))==(Upper|Lower)) + || ( (Mode&Upper) && !LhsIsRowMajor) + || ( (Mode&Lower) && LhsIsRowMajor), + ProcessSecondHalf = !ProcessFirstHalf + }; + + LhsEval lhsEval(lhs); + + for (Index j=0; j +// in the future selfadjoint-ness should be defined by the expression traits +// such that Transpose > is valid. (currently TriangularBase::transpose() is overloaded to make it work) +template +struct evaluator_traits > +{ + typedef typename storage_kind_to_evaluator_kind::Kind Kind; + typedef SparseSelfAdjointShape Shape; + + static const int AssumeAliasing = 0; +}; + +template +struct generic_product_impl +{ + template + static void evalTo(Dest& dst, const LhsView& lhsView, const Rhs& rhs) + { + typedef typename LhsView::_MatrixTypeNested Lhs; + typedef typename nested_eval::type LhsNested; + typedef typename nested_eval::type RhsNested; + LhsNested lhsNested(lhsView.matrix()); + RhsNested rhsNested(rhs); + + dst.setZero(); + internal::sparse_selfadjoint_time_dense_product(lhsNested, rhsNested, dst, typename Dest::Scalar(1)); + } +}; + +template +struct generic_product_impl +{ + template + static void evalTo(Dest& dst, const Lhs& lhs, const RhsView& rhsView) + { + typedef typename RhsView::_MatrixTypeNested Rhs; + typedef typename nested_eval::type LhsNested; + typedef typename nested_eval::type RhsNested; + LhsNested lhsNested(lhs); + RhsNested rhsNested(rhsView.matrix()); + + dst.setZero(); + // transpoe everything + Transpose dstT(dst); + internal::sparse_selfadjoint_time_dense_product(rhsNested.transpose(), lhsNested.transpose(), dstT, typename Dest::Scalar(1)); + } +}; + +template +struct product_evaluator, ProductTag, SparseSelfAdjointShape, DenseShape, typename Lhs::Scalar, typename Rhs::Scalar> + : public evaluator::PlainObject>::type +{ + typedef Product XprType; + typedef typename XprType::PlainObject 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, DenseShape, SparseSelfAdjointShape, typename Lhs::Scalar, typename Rhs::Scalar> + : public evaluator::PlainObject>::type +{ + typedef Product XprType; + typedef typename XprType::PlainObject 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, SparseSelfAdjointShape, SparseShape, typename LhsView::Scalar, typename Rhs::Scalar> + : public evaluator::PlainObject>::type +{ + typedef Product XprType; + typedef typename XprType::PlainObject PlainObject; + typedef typename evaluator::type Base; + + product_evaluator(const XprType& xpr) + : /*m_lhs(xpr.lhs()),*/ m_result(xpr.rows(), xpr.cols()) + { + m_lhs = xpr.lhs(); + ::new (static_cast(this)) Base(m_result); + generic_product_impl::evalTo(m_result, m_lhs, xpr.rhs()); + } + +protected: + typename Rhs::PlainObject m_lhs; + PlainObject m_result; +}; + +template +struct product_evaluator, ProductTag, SparseShape, SparseSelfAdjointShape, typename Lhs::Scalar, typename RhsView::Scalar> + : public evaluator::PlainObject>::type +{ + typedef Product XprType; + typedef typename XprType::PlainObject PlainObject; + typedef typename evaluator::type Base; + + product_evaluator(const XprType& xpr) + : m_rhs(xpr.rhs()), m_result(xpr.rows(), xpr.cols()) + { + ::new (static_cast(this)) Base(m_result); + generic_product_impl::evalTo(m_result, xpr.lhs(), m_rhs); + } + +protected: + typename Lhs::PlainObject m_rhs; + PlainObject m_result; +}; + +} // namespace internal + +#endif // EIGEN_TEST_EVALUATORS + /*************************************************************************** * Implementation of symmetric copies and permutations ***************************************************************************/ namespace internal { - -template -struct traits > : traits { -}; -template +template void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrix& _dest, const typename MatrixType::Index* perm) { typedef typename MatrixType::Index Index; @@ -337,11 +558,11 @@ void permute_symm_to_fullsymm(const MatrixType& mat, SparseMatrixc) || ( UpLo==Upper && rc) || ( Mode==Upper && rc) || ( (UpLo&Upper)==Upper && rc) || ( (Mode&Upper)==Upper && r +template void permute_symm_to_symm(const MatrixType& mat, SparseMatrix& _dest, const typename MatrixType::Index* perm) { typedef typename MatrixType::Index Index; @@ -407,8 +628,8 @@ void permute_symm_to_symm(const MatrixType& mat, SparseMatrixj)) + if((int(SrcMode)==int(Lower) && ij)) continue; Index ip = perm ? perm[i] : i; - count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; + count[int(DstMode)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; } } dest.outerIndexPtr()[0] = 0; @@ -441,17 +662,17 @@ void permute_symm_to_symm(const MatrixType& mat, SparseMatrixj)) + if((int(SrcMode)==int(Lower) && ij)) continue; Index jp = perm ? perm[j] : j; Index ip = perm? perm[i] : i; - Index k = count[int(DstUpLo)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; - dest.innerIndexPtr()[k] = int(DstUpLo)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp); + Index k = count[int(DstMode)==int(Lower) ? (std::min)(ip,jp) : (std::max)(ip,jp)]++; + dest.innerIndexPtr()[k] = int(DstMode)==int(Lower) ? (std::max)(ip,jp) : (std::min)(ip,jp); if(!StorageOrderMatch) std::swap(ip,jp); - if( ((int(DstUpLo)==int(Lower) && ipjp))) + if( ((int(DstMode)==int(Lower) && ipjp))) dest.valuePtr()[k] = numext::conj(it.value()); else dest.valuePtr()[k] = it.value(); @@ -461,9 +682,19 @@ void permute_symm_to_symm(const MatrixType& mat, SparseMatrix +#ifndef EIGEN_TEST_EVALUATORS + +namespace internal { + +template +struct traits > : traits { +}; + +} + +template class SparseSymmetricPermutationProduct - : public EigenBase > + : public EigenBase > { public: typedef typename MatrixType::Scalar Scalar; @@ -485,15 +716,15 @@ class SparseSymmetricPermutationProduct template void evalTo(SparseMatrix& _dest) const { -// internal::permute_symm_to_fullsymm(m_matrix,_dest,m_perm.indices().data()); +// internal::permute_symm_to_fullsymm(m_matrix,_dest,m_perm.indices().data()); SparseMatrix tmp; - internal::permute_symm_to_fullsymm(m_matrix,tmp,m_perm.indices().data()); + internal::permute_symm_to_fullsymm(m_matrix,tmp,m_perm.indices().data()); _dest = tmp; } - template void evalTo(SparseSelfAdjointView& dest) const + template void evalTo(SparseSelfAdjointView& dest) const { - internal::permute_symm_to_symm(m_matrix,dest.matrix(),m_perm.indices().data()); + internal::permute_symm_to_symm(m_matrix,dest.matrix(),m_perm.indices().data()); } protected: @@ -502,6 +733,10 @@ class SparseSymmetricPermutationProduct }; +#else // EIGEN_TEST_EVALUATORS + +#endif // EIGEN_TEST_EVALUATORS + } // end namespace Eigen #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp index 27bc548f8..fa9be5440 100644 --- a/test/sparse_product.cpp +++ b/test/sparse_product.cpp @@ -19,7 +19,7 @@ template void sparse_product() typedef typename SparseMatrixType::Scalar Scalar; enum { Flags = SparseMatrixType::Flags }; - double density = (std::max)(8./(rows*cols), 0.1); + double density = (std::max)(8./(rows*cols), 0.2); typedef Matrix DenseMatrix; typedef Matrix DenseVector; typedef Matrix RowDenseVector; @@ -109,7 +109,7 @@ template void sparse_product() Index c1 = internal::random(0,cols-1); Index r1 = internal::random(0,depth-1); DenseMatrix dm5 = DenseMatrix::Random(depth, cols); - + VERIFY_IS_APPROX( m4=m2.col(c)*dm5.col(c1).transpose(), refMat4=refMat2.col(c)*dm5.col(c1).transpose()); VERIFY_IS_EQUAL(m4.nonZeros(), (refMat4.array()!=0).count()); VERIFY_IS_APPROX( m4=m2.middleCols(c,1)*dm5.col(c1).transpose(), refMat4=refMat2.col(c)*dm5.col(c1).transpose()); @@ -153,11 +153,11 @@ template void sparse_product() RowSpVector rv0(depth), rv1; RowDenseVector drv0(depth), drv1(rv1); initSparse(2*density,drv0, rv0); - - VERIFY_IS_APPROX(cv1=rv0*m3, dcv1=drv0*refMat3); + + VERIFY_IS_APPROX(cv1=m3*cv0, dcv1=refMat3*dcv0); VERIFY_IS_APPROX(rv1=rv0*m3, drv1=drv0*refMat3); - VERIFY_IS_APPROX(cv1=m3*cv0, dcv1=refMat3*dcv0); VERIFY_IS_APPROX(cv1=m3t.adjoint()*cv0, dcv1=refMat3t.adjoint()*dcv0); + VERIFY_IS_APPROX(cv1=rv0*m3, dcv1=drv0*refMat3); VERIFY_IS_APPROX(rv1=m3*cv0, drv1=refMat3*dcv0); }