From 2cd9bb7380b213cd370269ba072a332cba19f948 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Antonio=20S=C3=A1nchez?= <cantonios@google.com>
Date: Wed, 25 Feb 2026 00:27:06 +0000
Subject: [PATCH] Fix sparse product with entities that do not have direct
 access.

libeigen/eigen!2205
---
 Eigen/src/SparseCore/SparseDenseProduct.h | 42 ++++++++++++++++-------
 test/sparse_product.cpp                   |  5 +++
 2 files changed, 35 insertions(+), 12 deletions(-)

diff --git a/Eigen/src/SparseCore/SparseDenseProduct.h b/Eigen/src/SparseCore/SparseDenseProduct.h
index 4de497f42..98f81dd5e 100644
--- a/Eigen/src/SparseCore/SparseDenseProduct.h
+++ b/Eigen/src/SparseCore/SparseDenseProduct.h
@@ -61,6 +61,12 @@ struct sparse_time_dense_product_impl<SparseLhsType, DenseRhsType, DenseResType,
   // Direct pointer path: works for both compressed and non-compressed storage.
   static void runCol(const LhsEval& /*lhsEval*/, const SparseLhsType& lhs, const DenseRhsType& rhs, DenseResType& res,
                      const ResScalar& alpha, Index n, Index c, std::true_type /* has_compressed_storage */) {
+    runColImpl(lhs, rhs, res, alpha, n, c, std::integral_constant<bool, bool(DenseRhsType::Flags & DirectAccessBit)>());
+  }
+
+  template <typename RhsT>
+  static void runColImpl(const SparseLhsType& lhs, const RhsT& rhs, DenseResType& res, const ResScalar& alpha, Index n,
+                         Index c, std::true_type) {
     const Lhs& mat = lhs;
     const auto* vals = mat.valuePtr();
     const auto* inds = mat.innerIndexPtr();
@@ -105,20 +111,32 @@ struct sparse_time_dense_product_impl<SparseLhsType, DenseRhsType, DenseResType,
         }
       }
     } else {
-      // Non-unit rhs stride: use direct pointers for sparse side, coeff() for rhs
-      for (Index i = 0; i < n; ++i) {
-        Index k = outer[i];
-        const Index end = innerNnz ? outer[i] + innerNnz[i] : outer[i + 1];
-        ResScalar sum0(0), sum1(0);
-        for (; k < end; ++k) {
-          sum0 += vals[k] * rhs.coeff(inds[k], c);
-          ++k;
-          if (k < end) {
-            sum1 += vals[k] * rhs.coeff(inds[k], c);
-          }
+      runColImpl(lhs, rhs, res, alpha, n, c, std::false_type());
+    }
+  }
+
+  // Use fall-back path without direct access to rhs.
+  template <typename RhsT>
+  static void runColImpl(const SparseLhsType& lhs, const RhsT& rhs, DenseResType& res, const ResScalar& alpha, Index n,
+                         Index c, std::false_type) {
+    const Lhs& mat = lhs;
+    const auto* vals = mat.valuePtr();
+    const auto* inds = mat.innerIndexPtr();
+    const auto* outer = mat.outerIndexPtr();
+    const auto* innerNnz = mat.innerNonZeroPtr();
+    // Non-unit rhs stride (or no direct access): use direct pointers for sparse side, coeff() for rhs
+    for (Index i = 0; i < n; ++i) {
+      Index k = outer[i];
+      const Index end = innerNnz ? outer[i] + innerNnz[i] : outer[i + 1];
+      ResScalar sum0(0), sum1(0);
+      for (; k < end; ++k) {
+        sum0 += vals[k] * rhs.coeff(inds[k], c);
+        ++k;
+        if (k < end) {
+          sum1 += vals[k] * rhs.coeff(inds[k], c);
         }
-        res.coeffRef(i, c) += alpha * (sum0 + sum1);
       }
+      res.coeffRef(i, c) += alpha * (sum0 + sum1);
     }
   }
 
diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp
index 1083f3f1a..3f05f6613 100644
--- a/test/sparse_product.cpp
+++ b/test/sparse_product.cpp
@@ -158,6 +158,11 @@ void sparse_product() {
     VERIFY_IS_APPROX(dm4 = m2t.transpose() * (refMat3 + refMat5) * 0.5,
                      refMat4 = refMat2t.transpose() * (refMat3 + refMat5) * 0.5);
 
+    // sparse * dense expression without DirectAccessBit (e.g. CwiseNullaryOp)
+    VERIFY_IS_APPROX(dm4 = m2 * DenseMatrix::Constant(depth, cols, s1),
+                     refMat4 = refMat2 * DenseMatrix::Constant(depth, cols, s1));
+    VERIFY_IS_APPROX(dm4 = m2 * DenseMatrix::Zero(depth, cols), refMat4 = refMat2 * DenseMatrix::Zero(depth, cols));
+
     // sparse * dense vector
     VERIFY_IS_APPROX(dm4.col(0) = m2 * refMat3.col(0), refMat4.col(0) = refMat2 * refMat3.col(0));
     VERIFY_IS_APPROX(dm4.col(0) = m2 * refMat3t.transpose().col(0),