From cc8c7cf0e6217778741e86ba278865c26b619984 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen <4643818-rmlarsen1@users.noreply.gitlab.com> Date: Sat, 21 Mar 2026 06:31:12 -0700 Subject: [PATCH] Fix bugs and clean up SparseCore module libeigen/eigen!2250 Co-authored-by: Rasmus Munk Larsen --- Eigen/src/SparseCore/AmbiVector.h | 4 +- .../ConservativeSparseSparseProduct.h | 10 +-- Eigen/src/SparseCore/SparseCwiseBinaryOp.h | 2 +- Eigen/src/SparseCore/SparsePermutation.h | 2 +- Eigen/src/SparseCore/SparseSelfAdjointView.h | 3 +- Eigen/src/SparseCore/TriangularSolver.h | 9 --- test/sparse_basic.cpp | 76 +++++++++++++++++++ 7 files changed, 82 insertions(+), 24 deletions(-) diff --git a/Eigen/src/SparseCore/AmbiVector.h b/Eigen/src/SparseCore/AmbiVector.h index 4af6f8825..99c896ab8 100644 --- a/Eigen/src/SparseCore/AmbiVector.h +++ b/Eigen/src/SparseCore/AmbiVector.h @@ -238,8 +238,8 @@ Scalar_& AmbiVector::coeff(Index i) { Index elid = m_llStart; while (elid >= 0 && llElements[elid].index < i) elid = llElements[elid].next; - if (llElements[elid].index == i) - return llElements[m_llCurrent].value; + if (elid >= 0 && llElements[elid].index == i) + return llElements[elid].value; else return m_zero; } diff --git a/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h b/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h index 2beadab9f..dc2866e81 100644 --- a/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h +++ b/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h @@ -33,7 +33,7 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r ei_declare_aligned_stack_constructed_variable(ResScalar, values, rows, 0); ei_declare_aligned_stack_constructed_variable(Index, indices, rows, 0); - std::memset(mask, 0, sizeof(bool) * rows); + std::fill_n(mask, rows, false); evaluator lhsEval(lhs); evaluator rhsEval(rhs); @@ -106,10 +106,6 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r res.finalize(); } -} // end namespace internal - -namespace internal { - // Helper template to generate new sparse matrix types template using WithStorageOrder = SparseMatrix; @@ -232,10 +228,6 @@ struct conservative_sparse_sparse_product_selector static void sparse_sparse_to_dense_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res) { typedef typename remove_all_t::Scalar LhsScalar; diff --git a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h index e3dcd57b2..b93194116 100644 --- a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h +++ b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h @@ -858,7 +858,7 @@ Derived& SparseMatrixBase::operator+=(const EigenBase& ot template template Derived& SparseMatrixBase::operator-=(const EigenBase& other) { - call_assignment(derived(), other.derived(), internal::assign_op()); + call_assignment(derived(), other.derived(), internal::sub_assign_op()); return derived(); } diff --git a/Eigen/src/SparseCore/SparsePermutation.h b/Eigen/src/SparseCore/SparsePermutation.h index 56f572d35..3fb370507 100644 --- a/Eigen/src/SparseCore/SparsePermutation.h +++ b/Eigen/src/SparseCore/SparsePermutation.h @@ -246,4 +246,4 @@ inline const Product, SparseDerived, AliasFreeProduct> } // end namespace Eigen -#endif // EIGEN_SPARSE_SELFADJOINTVIEW_H +#endif // EIGEN_SPARSE_PERMUTATION_H diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index 210462fc9..9b290ddd8 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -161,8 +161,7 @@ class SparseSelfAdjointView : public EigenBase void evalTo(Dest&) const; diff --git a/Eigen/src/SparseCore/TriangularSolver.h b/Eigen/src/SparseCore/TriangularSolver.h index b0054870a..f43f0c853 100644 --- a/Eigen/src/SparseCore/TriangularSolver.h +++ b/Eigen/src/SparseCore/TriangularSolver.h @@ -251,17 +251,8 @@ void TriangularViewImpl::solveInPlace(SparseMatrix eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows()); eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper | Lower))); - // enum { copy = internal::traits::Flags & RowMajorBit }; - - // typedef std::conditional_t::type, OtherDerived&> OtherCopy; - // OtherCopy otherCopy(other.derived()); - internal::sparse_solve_triangular_sparse_selector::run( derived().nestedExpression(), other.derived()); - - // if (copy) - // other = otherCopy; } #endif diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index f1fd6b771..4f1d4473e 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -1017,6 +1017,80 @@ void big_sparse_triplet(Index rows, Index cols, double density) { VERIFY_IS_APPROX(sum, m.sum()); } +// Regression test: SparseMatrixBase::operator-=(const EigenBase&) must subtract, not overwrite. +// SparseSelfAdjointView inherits EigenBase (not SparseMatrixBase), so it exercises this overload. +template +void sparse_sub_assign_eigenbase() { + Index n = 10; + SparseMatrix A(n, n), B(n, n); + A.setIdentity(); + B.setIdentity(); + B *= 3.0; + + // Make B symmetric so selfadjointView is valid + MatrixXd refA = MatrixXd(A); + MatrixXd refB = MatrixXd(B); + + // operator-= going through the EigenBase overload + A -= B.selfadjointView(); + refA -= refB.selfadjointView(); + VERIFY_IS_APPROX(MatrixXd(A), refA); + + // Also test operator+= for symmetry + A += B.selfadjointView(); + refA += refB.selfadjointView(); + VERIFY_IS_APPROX(MatrixXd(A), refA); +} + +// Regression test: AmbiVector::coeff() must return the correct linked-list element. +template +void ambivector_coeff() { + using namespace Eigen::internal; + const Index size = 20; + AmbiVector vec(size); + vec.setBounds(0, size); + + // Test in sparse (linked-list) mode + vec.init(IsSparse); + vec.restart(); + // Insert elements out of order to stress linked-list traversal: + // Insert at indices 2, 5, 10, 15 with distinct values + vec.coeffRef(2) = 2.0; + vec.coeffRef(5) = 5.0; + vec.coeffRef(10) = 10.0; + vec.coeffRef(15) = 15.0; + + // coeff() should return the correct value for each inserted element + VERIFY_IS_APPROX(vec.coeff(2), 2.0); + VERIFY_IS_APPROX(vec.coeff(5), 5.0); + VERIFY_IS_APPROX(vec.coeff(10), 10.0); + VERIFY_IS_APPROX(vec.coeff(15), 15.0); + + // coeff() should return zero for non-inserted elements + VERIFY_IS_EQUAL(vec.coeff(0), 0.0); + VERIFY_IS_EQUAL(vec.coeff(3), 0.0); + VERIFY_IS_EQUAL(vec.coeff(7), 0.0); + VERIFY_IS_EQUAL(vec.coeff(18), 0.0); + + // Verify coeff() still works after coeffRef() advances m_llCurrent. + // This is the specific scenario that triggered the bug: coeffRef(15) + // leaves m_llCurrent pointing at element 15, then coeff(2) should + // still return 2.0, not 15.0. + vec.restart(); + vec.coeffRef(15) += 0.0; // advances m_llCurrent to element at index 15 + VERIFY_IS_APPROX(vec.coeff(2), 2.0); + VERIFY_IS_APPROX(vec.coeff(5), 5.0); + + // Test in dense mode as well (simpler, but for completeness) + vec.init(IsDense); + vec.setZero(); + vec.coeffRef(3) = 3.0; + vec.coeffRef(7) = 7.0; + VERIFY_IS_APPROX(vec.coeff(3), 3.0); + VERIFY_IS_APPROX(vec.coeff(7), 7.0); + VERIFY_IS_EQUAL(vec.coeff(0), 0.0); +} + template void bug1105() { // Regression test for bug 1105 @@ -1066,5 +1140,7 @@ EIGEN_DECLARE_TEST(sparse_basic) { CALL_SUBTEST_5((big_sparse_triplet>(10000, 10000, 0.125))); CALL_SUBTEST_5(bug1105<0>()); + CALL_SUBTEST_1(sparse_sub_assign_eigenbase<0>()); + CALL_SUBTEST_1(ambivector_coeff<0>()); } #endif