diff --git a/Eigen/src/Sparse/DynamicSparseMatrix.h b/Eigen/src/Sparse/DynamicSparseMatrix.h index 88ad85f84..49a78fd26 100644 --- a/Eigen/src/Sparse/DynamicSparseMatrix.h +++ b/Eigen/src/Sparse/DynamicSparseMatrix.h @@ -108,7 +108,7 @@ class DynamicSparseMatrix /** \returns a reference to the coefficient value at given position \a row, \a col * This operation involes a log(rho*outer_size) binary search. If the coefficient does not - * exist yet, then a sorted insertion Indexo a sequential buffer is performed. + * exist yet, then a sorted insertion into a sequential buffer is performed. */ inline Scalar& coeffRef(Index row, Index col) { diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index 38e810725..4fa5d5e33 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -69,6 +69,7 @@ class SparseMatrix { public: EIGEN_SPARSE_PUBLIC_INTERFACE(SparseMatrix) + using Base::operator=; EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseMatrix, +=) EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseMatrix, -=) // FIXME: why are these operator already alvailable ??? diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index ba6f64ab7..bc16be955 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -46,6 +46,16 @@ template class SparseMatrixBase : public EigenBase typedef typename internal::traits::Index Index; typedef SparseMatrixBase StorageBaseType; + typedef EigenBase Base; + + template + Derived& operator=(const EigenBase &other) + { + other.derived().evalTo(derived()); + return derived(); + } + +// using Base::operator=; enum { @@ -173,7 +183,7 @@ template class SparseMatrixBase : public EigenBase Derived& markAsRValue() { m_isRValue = true; return derived(); } SparseMatrixBase() : m_isRValue(false) { /* TODO check flags */ } - + inline Derived& operator=(const Derived& other) { // std::cout << "Derived& operator=(const Derived& other)\n"; diff --git a/Eigen/src/Sparse/SparseSelfAdjointView.h b/Eigen/src/Sparse/SparseSelfAdjointView.h index 2a7101082..b4ece8dd4 100644 --- a/Eigen/src/Sparse/SparseSelfAdjointView.h +++ b/Eigen/src/Sparse/SparseSelfAdjointView.h @@ -45,12 +45,24 @@ class SparseSelfAdjointTimeDenseProduct; template class DenseTimeSparseSelfAdjointProduct; +namespace internal { + +template +struct traits > : traits { +}; + +} + template class SparseSelfAdjointView + : public EigenBase > { public: 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) { @@ -77,7 +89,7 @@ template class SparseSelfAdjointView DenseTimeSparseSelfAdjointProduct operator*(const MatrixBase& lhs, const SparseSelfAdjointView& rhs) { - return DenseTimeSparseSelfAdjointProduct(lhs.derived(), rhs.m_matrix); + return DenseTimeSparseSelfAdjointProduct(lhs.derived(), rhs.m_matrix); } /** Perform a symmetric rank K update of the selfadjoint matrix \c *this: @@ -93,6 +105,70 @@ template class SparseSelfAdjointView */ template SparseSelfAdjointView& rankUpdate(const SparseMatrixBase& u, Scalar alpha = Scalar(1)); + + /** \internal triggered by sparse_matrix = SparseSelfadjointView; */ + template void evalTo(SparseMatrix& _dest) const + { + typedef SparseMatrix Dest; + Dest& dest(_dest.derived()); + enum { + StorageOrderMatch = Dest::IsRowMajor == _MatrixTypeNested::IsRowMajor + }; + VectorI count = Dest::IsRowMajor ? m_countPerRow : m_countPerCol; + Index size = m_matrix.rows(); + count.resize(size); + count.setZero(); + dest.resize(size,size); + for(Index j = 0; jj) || (UpLo==Upper && ij) || (UpLo==Upper && i llt() const; // const SparseLDLT ldlt() const; @@ -100,6 +176,8 @@ template class SparseSelfAdjointView protected: const typename MatrixType::Nested m_matrix; + VectorI m_countPerRow; + VectorI m_countPerCol; }; /*************************************************************************** @@ -133,7 +211,7 @@ SparseSelfAdjointView::rankUpdate(const SparseMatrixBase(); else - m_matrix.const_cast_derived() /*+*/= alpha * tmp.template triangularView(); + m_matrix.const_cast_derived() += alpha * tmp.template triangularView(); return *this; } diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index 92fff4057..9d79ca740 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -251,10 +251,21 @@ template void sparse_basic(const SparseMatrixType& re VERIFY_IS_APPROX(m2, refM2); } + // test selfadjointView + { + DenseMatrix refMat2(rows, rows), refMat3(rows, rows); + SparseMatrixType m2(rows, rows), m3(rows, rows); + initSparse(density, refMat2, m2); + refMat3 = refMat2.template selfadjointView(); + m3 = m2.template selfadjointView(); + VERIFY_IS_APPROX(m3, refMat3); + } + // test sparseView { DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); SparseMatrixType m2(rows, rows); + initSparse(density, refMat2, m2); VERIFY_IS_APPROX(m2.eval(), refMat2.sparseView().eval()); } }