Fix bugs and clean up SparseCore module

libeigen/eigen!2250

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>

Eigen/src/SparseCore/AmbiVector.h

@@ -238,8 +238,8 @@ Scalar_& AmbiVector<Scalar_, StorageIndex_>::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;
     }

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<Lhs> lhsEval(lhs);
   evaluator<Rhs> 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 <class Source, int Order>
 using WithStorageOrder = SparseMatrix<typename Source::Scalar, Order, typename Source::StorageIndex>;
@@ -232,10 +228,6 @@ struct conservative_sparse_sparse_product_selector<Lhs, Rhs, ResultType, RowMajo
   }
 };
 
-}  // end namespace internal
-
-namespace internal {
-
 template <typename Lhs, typename Rhs, typename ResultType>
 static void sparse_sparse_to_dense_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res) {
   typedef typename remove_all_t<Lhs>::Scalar LhsScalar;

Eigen/src/SparseCore/SparseCwiseBinaryOp.h

@@ -858,7 +858,7 @@ Derived& SparseMatrixBase<Derived>::operator+=(const EigenBase<OtherDerived>& ot
 template <typename Derived>
 template <typename OtherDerived>
 Derived& SparseMatrixBase<Derived>::operator-=(const EigenBase<OtherDerived>& other) {
-  call_assignment(derived(), other.derived(), internal::assign_op<Scalar, typename OtherDerived::Scalar>());
+  call_assignment(derived(), other.derived(), internal::sub_assign_op<Scalar, typename OtherDerived::Scalar>());
   return derived();
 }
 

Eigen/src/SparseCore/SparsePermutation.h

@@ -246,4 +246,4 @@ inline const Product<Inverse<PermutationType>, SparseDerived, AliasFreeProduct>
 
 }  // end namespace Eigen
 
-#endif  // EIGEN_SPARSE_SELFADJOINTVIEW_H
+#endif  // EIGEN_SPARSE_PERMUTATION_H

Eigen/src/SparseCore/SparseSelfAdjointView.h

@@ -161,8 +161,7 @@ class SparseSelfAdjointView : public EigenBase<SparseSelfAdjointView<MatrixType,
 
  protected:
   MatrixTypeNested m_matrix;
-  // mutable VectorI m_countPerRow;
-  // mutable VectorI m_countPerCol;
+
  private:
   template <typename Dest>
   void evalTo(Dest&) const;

Eigen/src/SparseCore/TriangularSolver.h

@@ -251,17 +251,8 @@ void TriangularViewImpl<ExpressionType, Mode, Sparse>::solveInPlace(SparseMatrix
   eigen_assert(derived().cols() == derived().rows() && derived().cols() == other.rows());
   eigen_assert((!(Mode & ZeroDiag)) && bool(Mode & (Upper | Lower)));
 
-  //   enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit };
-
-  //   typedef std::conditional_t<copy,
-  //     typename internal::plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&> OtherCopy;
-  //   OtherCopy otherCopy(other.derived());
-
   internal::sparse_solve_triangular_sparse_selector<ExpressionType, OtherDerived, Mode>::run(
       derived().nestedExpression(), other.derived());
-
-  //   if (copy)
-  //     other = otherCopy;
 }
 #endif
 

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 <int>
+void sparse_sub_assign_eigenbase() {
+  Index n = 10;
+  SparseMatrix<double> 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<OtherDerived> overload
+  A -= B.selfadjointView<Lower>();
+  refA -= refB.selfadjointView<Lower>();
+  VERIFY_IS_APPROX(MatrixXd(A), refA);
+
+  // Also test operator+= for symmetry
+  A += B.selfadjointView<Lower>();
+  refA += refB.selfadjointView<Lower>();
+  VERIFY_IS_APPROX(MatrixXd(A), refA);
+}
+
+// Regression test: AmbiVector::coeff() must return the correct linked-list element.
+template <int>
+void ambivector_coeff() {
+  using namespace Eigen::internal;
+  const Index size = 20;
+  AmbiVector<double, int> 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 <int>
 void bug1105() {
   // Regression test for bug 1105
@@ -1066,5 +1140,7 @@ EIGEN_DECLARE_TEST(sparse_basic) {
   CALL_SUBTEST_5((big_sparse_triplet<SparseMatrix<double, ColMajor, long int>>(10000, 10000, 0.125)));
 
   CALL_SUBTEST_5(bug1105<0>());
+  CALL_SUBTEST_1(sparse_sub_assign_eigenbase<0>());
+  CALL_SUBTEST_1(ambivector_coeff<0>());
 }
 #endif