From 2d53466fa936552f95b1c3609faa252caa464260 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 14 Jan 2009 21:27:54 +0000 Subject: [PATCH] Sparse module: * improved performance of mat*=scalar * bug fix in cwise* --- Eigen/src/Sparse/SparseCwise.h | 38 ++--- Eigen/src/Sparse/SparseCwiseBinaryOp.h | 197 +++++-------------------- Eigen/src/Sparse/SparseCwiseUnaryOp.h | 10 +- Eigen/src/Sparse/SparseMatrix.h | 1 + Eigen/src/Sparse/SparseVector.h | 1 + test/sparse_basic.cpp | 38 ++++- test/sparse_vector.cpp | 4 +- 7 files changed, 98 insertions(+), 191 deletions(-) diff --git a/Eigen/src/Sparse/SparseCwise.h b/Eigen/src/Sparse/SparseCwise.h index 0697439d3..2206883cc 100644 --- a/Eigen/src/Sparse/SparseCwise.h +++ b/Eigen/src/Sparse/SparseCwise.h @@ -73,8 +73,16 @@ template class SparseCwise operator*(const SparseMatrixBase &other) const; template - const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op) - operator/(const SparseMatrixBase &other) const; + const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE + operator*(const MatrixBase &other) const; + +// template +// const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op) +// operator/(const SparseMatrixBase &other) const; +// +// template +// const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op) +// operator/(const MatrixBase &other) const; template const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op) @@ -96,27 +104,13 @@ template class SparseCwise // const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_sin_op) sin() const; // const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_pow_op) pow(const Scalar& exponent) const; + template + inline ExpressionType& operator*=(const SparseMatrixBase &other); + +// template +// inline ExpressionType& operator/=(const SparseMatrixBase &other); + /* - const ScalarAddReturnType - operator+(const Scalar& scalar) const; - - friend const ScalarAddReturnType - operator+(const Scalar& scalar, const Cwise& mat) - { return mat + scalar; } - - ExpressionType& operator+=(const Scalar& scalar); - - const ScalarAddReturnType - operator-(const Scalar& scalar) const; - - ExpressionType& operator-=(const Scalar& scalar); - - template - inline ExpressionType& operator*=(const MatrixBase &other); - - template - inline ExpressionType& operator/=(const MatrixBase &other); - template const EIGEN_CWISE_BINOP_RETURN_TYPE(std::less) operator<(const MatrixBase& other) const; diff --git a/Eigen/src/Sparse/SparseCwiseBinaryOp.h b/Eigen/src/Sparse/SparseCwiseBinaryOp.h index 99bbae0cd..bdd19018b 100644 --- a/Eigen/src/Sparse/SparseCwiseBinaryOp.h +++ b/Eigen/src/Sparse/SparseCwiseBinaryOp.h @@ -42,40 +42,6 @@ // 4 - dense op dense product dense // generic dense -// template class ei_sparse_cwise_binary_op -// : ei_no_assignment_operator -// { -// typedef CwiseBinaryOp CwiseBinaryXpr; -// public: -// -// typedef typename ei_traits::LhsNested LhsNested; -// typedef typename ei_traits::RhsNested RhsNested; -// -// EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) -// : m_lhs(lhs), m_rhs(rhs), m_functor(func) -// { -// EIGEN_STATIC_ASSERT((ei_functor_allows_mixing_real_and_complex::ret -// ? int(ei_is_same_type::ret) -// : int(ei_is_same_type::ret)), -// YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) -// // require the sizes to match -// EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs) -// ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols()); -// } -// -// EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); } -// EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); } -// -// EIGEN_STRONG_INLINE const LhsNested& _lhs() const { return m_lhs; } -// EIGEN_STRONG_INLINE const RhsNested& _rhs() const { return m_rhs; } -// EIGEN_STRONG_INLINE const UnaryOp& _functor() const { return m_functor; } -// -// protected: -// const LhsNested m_lhs; -// const RhsNested m_rhs; -// const BinaryOp m_functor; -// } - template struct ei_traits > { @@ -132,103 +98,15 @@ class SparseCwiseBinaryOp : ei_no_assignment_operator, EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); } EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); } - EIGEN_STRONG_INLINE const _LhsNested& _lhs() const { return m_lhs; } - EIGEN_STRONG_INLINE const _RhsNested& _rhs() const { return m_rhs; } - EIGEN_STRONG_INLINE const BinaryOp& _functor() const { return m_functor; } + EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; } + EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; } + EIGEN_STRONG_INLINE const BinaryOp& functor() const { return m_functor; } protected: const LhsNested m_lhs; const RhsNested m_rhs; const BinaryOp m_functor; }; -/* -template -class CwiseBinaryOp, Rhs> - : public ei_sparse_cwise_binary_op, - public SparseMatrixBase, Rhs> > -{ - public: - EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp) - EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) - : ei_sparse_cwise_binary_op(lhs, rhs, func) - {} -}; -template -class CwiseBinaryOp, SparseMatrixBase > - : public ei_sparse_cwise_binary_op, - public SparseMatrixBase > -{ - public: - EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp) - EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) - : ei_sparse_cwise_binary_op(lhs, rhs, func) - {} -};*/ - -// template -// class CwiseBinaryOp > : ei_no_assignment_operator, -// public SparseMatrixBase, SparseMatrixBase > > -// { -// public: -// -// EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp) -// typedef typename ei_traits::LhsNested LhsNested; -// typedef typename ei_traits::RhsNested RhsNested; -// -// EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) -// : m_lhs(lhs), m_rhs(rhs), m_functor(func) -// { -// EIGEN_STATIC_ASSERT((ei_functor_allows_mixing_real_and_complex::ret -// ? int(ei_is_same_type::ret) -// : int(ei_is_same_type::ret)), -// YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) -// // require the sizes to match -// EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs) -// ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols()); -// } -// -// EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); } -// EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); } -// -// protected: -// const LhsNested m_lhs; -// const RhsNested m_rhs; -// const BinaryOp m_functor; -// }; -// -// template -// class CwiseBinaryOp, SparseMatrixBase > : ei_no_assignment_operator, -// public SparseMatrixBase, SparseMatrixBase > > -// { -// public: -// -// EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp) -// typedef typename ei_traits::LhsNested LhsNested; -// typedef typename ei_traits::RhsNested RhsNested; -// -// EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp()) -// : m_lhs(lhs), m_rhs(rhs), m_functor(func) -// { -// EIGEN_STATIC_ASSERT((ei_functor_allows_mixing_real_and_complex::ret -// ? int(ei_is_same_type::ret) -// : int(ei_is_same_type::ret)), -// YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) -// // require the sizes to match -// EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs) -// ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols()); -// } -// -// EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); } -// EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); } -// -// protected: -// const LhsNested m_lhs; -// const RhsNested m_rhs; -// const BinaryOp m_functor; -// }; - -// template struct ei_is_scalar_product { enum { ret = false }; }; -// template struct ei_is_scalar_product > { enum { ret = true }; }; templateoperator++(); } @@ -323,24 +201,6 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector -// class SparseCwiseBinaryOp::InnerIterator -// : public ei_sparse_cwise_binary_op_inner_iterator_selector::InnerIterator> -// { -// typedef CwiseBinaryOpInnerIterator< -// BinaryOp,Lhs,Rhs, typename CwiseBinaryOp::InnerIterator> Base; -// public: -// typedef typename CwiseBinaryOp::Scalar Scalar; -// typedef typename ei_traits::_LhsNested _LhsNested; -// typedef typename _LhsNested::InnerIterator LhsIterator; -// typedef typename ei_traits::_RhsNested _RhsNested; -// typedef typename _RhsNested::InnerIterator RhsIterator; -// // public: -// EIGEN_STRONG_INLINE InnerIterator(const CwiseBinaryOp& binOp, int outer) -// : Base(binOp.m_lhs,binOp.m_rhs,binOp.m_functor,outer) -// {} -// }; - // sparse - sparse (product) template class ei_sparse_cwise_binary_op_inner_iterator_selector, Lhs, Rhs, Derived, IsSparse, IsSparse> @@ -355,9 +215,9 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector, public: EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer) - : m_lhsIter(xpr._lhs(),outer), m_rhsIter(xpr._rhs()), m_functor(xpr._functor()) + : m_lhsIter(xpr.lhs(),outer), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor()) { - while (m_lhsIter && m_rhsIter && m_lhsIter.index() != m_rhsIter.index()) + while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index())) { if (m_lhsIter.index() < m_rhsIter.index()) ++m_lhsIter; @@ -368,15 +228,11 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector, EIGEN_STRONG_INLINE Derived& operator++() { - while (m_lhsIter && m_rhsIter) + ++m_lhsIter; + ++m_rhsIter; + while (m_lhsIter && m_rhsIter && (m_lhsIter.index() != m_rhsIter.index())) { - if (m_lhsIter.index() == m_rhsIter.index()) - { - ++m_lhsIter; - ++m_rhsIter; - break; - } - else if (m_lhsIter.index() < m_rhsIter.index()) + if (m_lhsIter.index() < m_rhsIter.index()) ++m_lhsIter; else ++m_rhsIter; @@ -388,7 +244,7 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector, EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); } - EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter && m_rhsIter; } + EIGEN_STRONG_INLINE operator bool() const { return (m_lhsIter && m_rhsIter); } protected: LhsIterator m_lhsIter; @@ -419,7 +275,8 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector, } EIGEN_STRONG_INLINE Scalar value() const - { return m_functor(m_lhsIter.value(), m_xpr.rhs().coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); } + { return m_functor(m_lhsIter.value(), + m_xpr.rhs().coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); } EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); } @@ -513,22 +370,38 @@ SparseCwise::operator*(const SparseMatrixBase &oth template template -EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op) -SparseCwise::operator/(const SparseMatrixBase &other) const +EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE +SparseCwise::operator*(const MatrixBase &other) const { - return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)(_expression(), other.derived()); + return EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE(_expression(), other.derived()); } // template // template -// inline ExpressionType& Cwise::operator*=(const SparseMatrixBase &other) +// EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op) +// SparseCwise::operator/(const SparseMatrixBase &other) const // { -// return m_matrix.const_cast_derived() = *this * other; +// return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)(_expression(), other.derived()); // } +// +// template +// template +// EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op) +// SparseCwise::operator/(const MatrixBase &other) const +// { +// return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)(_expression(), other.derived()); +// } + +template +template +inline ExpressionType& SparseCwise::operator*=(const SparseMatrixBase &other) +{ + return m_matrix.const_cast_derived() = *this * other; +} // template // template -// inline ExpressionType& Cwise::operator/=(const SparseMatrixBase &other) +// inline ExpressionType& SparseCwise::operator/=(const SparseMatrixBase &other) // { // return m_matrix.const_cast_derived() = *this / other; // } diff --git a/Eigen/src/Sparse/SparseCwiseUnaryOp.h b/Eigen/src/Sparse/SparseCwiseUnaryOp.h index 2b2050edd..e66baa0df 100644 --- a/Eigen/src/Sparse/SparseCwiseUnaryOp.h +++ b/Eigen/src/Sparse/SparseCwiseUnaryOp.h @@ -162,14 +162,20 @@ template EIGEN_STRONG_INLINE Derived& SparseMatrixBase::operator*=(const Scalar& other) { - return *this = *this * other; + for (int j=0; j EIGEN_STRONG_INLINE Derived& SparseMatrixBase::operator/=(const Scalar& other) { - return *this = *this / other; + for (int j=0; j::InnerIterator inline InnerIterator& operator++() { m_id++; return *this; } inline Scalar value() const { return m_matrix.m_data.value(m_id); } + inline Scalar& valueRef() { return const_cast(m_matrix.m_data.value(m_id)); } inline int index() const { return m_matrix.m_data.index(m_id); } diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h index b189cb5ee..49f5bde55 100644 --- a/Eigen/src/Sparse/SparseVector.h +++ b/Eigen/src/Sparse/SparseVector.h @@ -316,6 +316,7 @@ class SparseVector::InnerIterator inline InnerIterator& operator++() { m_id++; return *this; } inline Scalar value() const { return m_vector.m_data.value(m_id); } + inline Scalar& valueRef() { return const_cast(m_vector.m_data.value(m_id)); } inline int index() const { return m_vector.m_data.index(m_id); } diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index 07a38ddd8..84cf4cd16 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -52,6 +52,7 @@ template void sparse_basic(int rows, int cols) SparseMatrix m(rows, cols); DenseMatrix refMat = DenseMatrix::Zero(rows, cols); DenseVector vec1 = DenseVector::Random(rows); + Scalar s1 = ei_random(); std::vector zeroCoords; std::vector nonzeroCoords; @@ -191,6 +192,35 @@ template void sparse_basic(int rows, int cols) // std::cerr << m.transpose() << "\n\n" << refMat.transpose() << "\n\n"; // VERIFY_IS_APPROX(m, refMat); + // test basic computations + { + DenseMatrix refM1 = DenseMatrix::Zero(rows, rows); + DenseMatrix refM2 = DenseMatrix::Zero(rows, rows); + DenseMatrix refM3 = DenseMatrix::Zero(rows, rows); + DenseMatrix refM4 = DenseMatrix::Zero(rows, rows); + SparseMatrix m1(rows, rows); + SparseMatrix m2(rows, rows); + SparseMatrix m3(rows, rows); + SparseMatrix m4(rows, rows); + initSparse(density, refM1, m1); + initSparse(density, refM2, m2); + initSparse(density, refM3, m3); + initSparse(density, refM4, m4); + + VERIFY_IS_APPROX(m1+m2, refM1+refM2); + VERIFY_IS_APPROX(m1+m2+m3, refM1+refM2+refM3); + VERIFY_IS_APPROX(m3.cwise()*(m1+m2), refM3.cwise()*(refM1+refM2)); + VERIFY_IS_APPROX(m1*s1-m2, refM1*s1-refM2); + + VERIFY_IS_APPROX(m1*=s1, refM1*=s1); + VERIFY_IS_APPROX(m1/=s1, refM1/=s1); + + refM4.setRandom(); + // sparse cwise* dense + VERIFY_IS_APPROX(m3.cwise()*refM4, refM3.cwise()*refM4); +// VERIFY_IS_APPROX(m3.cwise()/refM4, refM3.cwise()/refM4); + } + // test innerVector() { DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows); @@ -198,8 +228,8 @@ template void sparse_basic(int rows, int cols) initSparse(density, refMat2, m2); int j0 = ei_random(0,rows-1); int j1 = ei_random(0,rows-1); -// VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0)); -// VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1)); + VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0)); + VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1)); } // test transpose @@ -227,13 +257,13 @@ template void sparse_basic(int rows, int cols) VERIFY_IS_APPROX(m4=m2.transpose()*m3, refMat4=refMat2.transpose()*refMat3); VERIFY_IS_APPROX(m4=m2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); VERIFY_IS_APPROX(m4=m2*m3.transpose(), refMat4=refMat2*refMat3.transpose()); - + // sparse * dense VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3); VERIFY_IS_APPROX(dm4=m2*refMat3.transpose(), refMat4=refMat2*refMat3.transpose()); VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3, refMat4=refMat2.transpose()*refMat3); VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose()); - + // dense * sparse VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3); VERIFY_IS_APPROX(dm4=refMat2*m3.transpose(), refMat4=refMat2*refMat3.transpose()); diff --git a/test/sparse_vector.cpp b/test/sparse_vector.cpp index 18ef74833..0a66af621 100644 --- a/test/sparse_vector.cpp +++ b/test/sparse_vector.cpp @@ -77,7 +77,9 @@ template void sparse_vector(int rows, int cols) VERIFY_IS_APPROX(v1*s1-v2, refV1*s1-refV2); -// std::cerr << v1.dot(v2) << " == " << refV1.dot(refV2) << "\n"; + VERIFY_IS_APPROX(v1*=s1, refV1*=s1); + VERIFY_IS_APPROX(v1/=s1, refV1/=s1); + VERIFY_IS_APPROX(v1.dot(v2), refV1.dot(refV2)); }