From 763c833637d3918c32dc9c7ce5c9fcf647c7479b Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 24 Jun 2015 17:54:09 +0200 Subject: [PATCH] Make SparseSelfAdjointView, twists, and SparseQR more evaluator friendly --- Eigen/src/SparseCore/SparseAssign.h | 13 +- Eigen/src/SparseCore/SparseCwiseBinaryOp.h | 16 +++ Eigen/src/SparseCore/SparseSelfAdjointView.h | 135 +++++++++++++------ Eigen/src/SparseQR/SparseQR.h | 68 +++++++--- 4 files changed, 166 insertions(+), 66 deletions(-) diff --git a/Eigen/src/SparseCore/SparseAssign.h b/Eigen/src/SparseCore/SparseAssign.h index 469c2b188..93e0adbff 100644 --- a/Eigen/src/SparseCore/SparseAssign.h +++ b/Eigen/src/SparseCore/SparseAssign.h @@ -16,8 +16,7 @@ template template Derived& SparseMatrixBase::operator=(const EigenBase &other) { - // TODO use the evaluator mechanism - other.derived().evalTo(derived()); + internal::call_assignment_no_alias(derived(), other.derived()); return derived(); } @@ -137,8 +136,8 @@ struct Assignment }; // Sparse to Dense assignment -template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar> -struct Assignment +template< typename DstXprType, typename SrcXprType, typename Functor> +struct Assignment { static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) { @@ -153,8 +152,8 @@ struct Assignment } }; -template< typename DstXprType, typename SrcXprType, typename Scalar> -struct Assignment, Sparse2Dense, Scalar> +template< typename DstXprType, typename SrcXprType> +struct Assignment, Sparse2Dense> { static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { @@ -173,7 +172,7 @@ struct Assignment -struct Assignment, internal::assign_op, Sparse2Sparse, Scalar> +struct Assignment, internal::assign_op, Sparse2Sparse> { typedef Solve SrcXprType; static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) diff --git a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h index 096af7fb0..78d16892d 100644 --- a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h +++ b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h @@ -390,6 +390,22 @@ SparseMatrixBase::operator+=(const SparseMatrixBase& othe return derived() = derived() + other.derived(); } +template +template +Derived& SparseMatrixBase::operator+=(const DiagonalBase& other) +{ + call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op()); + return derived(); +} + +template +template +Derived& SparseMatrixBase::operator-=(const DiagonalBase& other) +{ + call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op()); + return derived(); +} + template template EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index 3da856799..dec96efb5 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -45,8 +45,13 @@ template class SparseSelfAdjointView { public: - enum { Mode = _Mode }; + enum { + Mode = _Mode, + RowsAtCompileTime = internal::traits::RowsAtCompileTime, + ColsAtCompileTime = internal::traits::ColsAtCompileTime + }; + typedef EigenBase Base; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::StorageIndex StorageIndex; typedef Matrix VectorI; @@ -116,20 +121,6 @@ template class SparseSelfAdjointView template SparseSelfAdjointView& rankUpdate(const SparseMatrixBase& u, const Scalar& alpha = Scalar(1)); - /** \internal triggered by sparse_matrix = SparseSelfadjointView; */ - template void evalTo(SparseMatrix& _dest) const - { - 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); - _dest = tmp; - } - /** \returns an expression of P H P^-1 */ // TODO implement twists in a more evaluator friendly fashion SparseSymmetricPermutationProduct<_MatrixTypeNested,Mode> twistedBy(const PermutationMatrix& perm) const @@ -140,7 +131,7 @@ template class SparseSelfAdjointView template SparseSelfAdjointView& operator=(const SparseSymmetricPermutationProduct& permutedMatrix) { - permutedMatrix.evalTo(*this); + internal::call_assignment_no_alias_no_transpose(*this, permutedMatrix); return *this; } @@ -157,11 +148,21 @@ template class SparseSelfAdjointView return *this = src.twistedBy(pnull); } + void resize(Index rows, Index cols) + { + EIGEN_ONLY_USED_FOR_DEBUG(rows); + EIGEN_ONLY_USED_FOR_DEBUG(cols); + eigen_assert(rows == this->rows() && cols == this->cols() + && "SparseSelfadjointView::resize() does not actually allow to resize."); + } + protected: typename MatrixType::Nested m_matrix; //mutable VectorI m_countPerRow; //mutable VectorI m_countPerCol; + private: + template void evalTo(Dest &) const; }; /*************************************************************************** @@ -200,6 +201,47 @@ SparseSelfAdjointView::rankUpdate(const SparseMatrixBase +// 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; +}; + +struct SparseSelfAdjoint2Sparse {}; + +template<> struct AssignmentKind { typedef SparseSelfAdjoint2Sparse Kind; }; +template<> struct AssignmentKind { typedef Sparse2Sparse Kind; }; + +template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar> +struct Assignment +{ + typedef typename DstXprType::StorageIndex StorageIndex; + template + static void run(SparseMatrix &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + { + internal::permute_symm_to_fullsymm(src.matrix(), dst); + } + + template + static void run(DynamicSparseMatrix& dst, const SrcXprType &src, const internal::assign_op &/*func*/) + { + // TODO directly evaluate into dst; + SparseMatrix tmp(dst.rows(),dst.cols()); + internal::permute_symm_to_fullsymm(src.matrix(), tmp); + dst = tmp; + } +}; + +} // end namespace internal + /*************************************************************************** * Implementation of sparse self-adjoint time dense matrix ***************************************************************************/ @@ -252,18 +294,7 @@ inline void sparse_selfadjoint_time_dense_product(const SparseLhsType& lhs, cons res.row(j) += i.value() * rhs.row(j); } } - -// TODO currently a selfadjoint expression has the form SelfAdjointView<.,.> -// 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 @@ -519,12 +550,16 @@ class SparseSymmetricPermutationProduct public: typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::StorageIndex StorageIndex; + enum { + RowsAtCompileTime = internal::traits::RowsAtCompileTime, + ColsAtCompileTime = internal::traits::ColsAtCompileTime + }; protected: typedef PermutationMatrix Perm; public: typedef Matrix VectorI; typedef typename MatrixType::Nested MatrixTypeNested; - typedef typename internal::remove_all::type _MatrixTypeNested; + typedef typename internal::remove_all::type NestedExpression; SparseSymmetricPermutationProduct(const MatrixType& mat, const Perm& perm) : m_matrix(mat), m_perm(perm) @@ -532,20 +567,9 @@ class SparseSymmetricPermutationProduct inline Index rows() const { return m_matrix.rows(); } inline Index cols() const { return m_matrix.cols(); } - - template - void evalTo(SparseMatrix& _dest) const - { -// 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()); - _dest = tmp; - } - - template void evalTo(SparseSelfAdjointView& dest) const - { - internal::permute_symm_to_symm(m_matrix,dest.matrix(),m_perm.indices().data()); - } + + const NestedExpression& matrix() const { return m_matrix; } + const Perm& perm() const { return m_perm; } protected: MatrixTypeNested m_matrix; @@ -553,6 +577,31 @@ class SparseSymmetricPermutationProduct }; +namespace internal { + +template +struct Assignment, internal::assign_op, Sparse2Sparse> +{ + typedef SparseSymmetricPermutationProduct SrcXprType; + typedef typename DstXprType::StorageIndex DstIndex; + template + static void run(SparseMatrix &dst, const SrcXprType &src, const internal::assign_op &) + { + // internal::permute_symm_to_fullsymm(m_matrix,_dest,m_perm.indices().data()); + SparseMatrix tmp; + internal::permute_symm_to_fullsymm(src.matrix(),tmp,src.perm().indices().data()); + dst = tmp; + } + + template + static void run(SparseSelfAdjointView& dst, const SrcXprType &src, const internal::assign_op &) + { + internal::permute_symm_to_symm(src.matrix(),dst.matrix(),src.perm().indices().data()); + } +}; + +} // end namespace internal + } // end namespace Eigen #endif // EIGEN_SPARSE_SELFADJOINTVIEW_H diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h index ce4a70454..548b3f9b0 100644 --- a/Eigen/src/SparseQR/SparseQR.h +++ b/Eigen/src/SparseQR/SparseQR.h @@ -23,6 +23,10 @@ namespace internal { typedef typename SparseQRType::MatrixType ReturnType; typedef typename ReturnType::StorageIndex StorageIndex; typedef typename ReturnType::StorageKind StorageKind; + enum { + RowsAtCompileTime = Dynamic, + ColsAtCompileTime = Dynamic + }; }; template struct traits > { @@ -235,8 +239,9 @@ class SparseQR : public SparseSolverBase > return m_info; } - protected: - inline void sort_matrix_Q() + + /** \internal */ + inline void _sort_matrix_Q() { if(this->m_isQSorted) return; // The matrix Q is sorted during the transposition @@ -267,7 +272,6 @@ class SparseQR : public SparseSolverBase > bool m_isEtreeOk; // whether the elimination tree match the initial input matrix template friend struct SparseQR_QProduct; - template friend struct SparseQRMatrixQReturnType; }; @@ -635,6 +639,10 @@ struct SparseQRMatrixQReturnType : public EigenBase DenseMatrix; + enum { + RowsAtCompileTime = Dynamic, + ColsAtCompileTime = Dynamic + }; explicit SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {} template SparseQR_QProduct operator*(const MatrixBase& other) @@ -652,19 +660,6 @@ struct SparseQRMatrixQReturnType : public EigenBase(m_qr); } - template void evalTo(MatrixBase& dest) const - { - dest.derived() = m_qr.matrixQ() * Dest::Identity(m_qr.rows(), m_qr.rows()); - } - template void evalTo(SparseMatrixBase& dest) const - { - Dest idMat(m_qr.rows(), m_qr.rows()); - idMat.setIdentity(); - // Sort the sparse householder reflectors if needed - const_cast(&m_qr)->sort_matrix_Q(); - dest.derived() = SparseQR_QProduct(m_qr, idMat, false); - } - const SparseQRType& m_qr; }; @@ -680,6 +675,47 @@ struct SparseQRMatrixQTransposeReturnType const SparseQRType& m_qr; }; +namespace internal { + +template +struct evaluator_traits > +{ + typedef typename SparseQRType::MatrixType MatrixType; + typedef typename storage_kind_to_evaluator_kind::Kind Kind; + typedef SparseShape Shape; + static const int AssumeAliasing = 0; +}; + +template< typename DstXprType, typename SparseQRType> +struct Assignment, internal::assign_op, Sparse2Sparse> +{ + typedef SparseQRMatrixQReturnType SrcXprType; + typedef typename DstXprType::Scalar Scalar; + typedef typename DstXprType::StorageIndex StorageIndex; + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + { + typename DstXprType::PlainObject idMat(src.m_qr.rows(), src.m_qr.rows()); + idMat.setIdentity(); + // Sort the sparse householder reflectors if needed + const_cast(&src.m_qr)->_sort_matrix_Q(); + dst = SparseQR_QProduct(src.m_qr, idMat, false); + } +}; + +template< typename DstXprType, typename SparseQRType> +struct Assignment, internal::assign_op, Sparse2Dense> +{ + typedef SparseQRMatrixQReturnType SrcXprType; + typedef typename DstXprType::Scalar Scalar; + typedef typename DstXprType::StorageIndex StorageIndex; + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + { + dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows()); + } +}; + +} // end namespace internal + } // end namespace Eigen #endif