devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix sparse product with entities that do not have direct access.

Browse Source 
libeigen/eigen!2205
This commit is contained in:
Antonio Sánchez
2026-02-25 00:27:06 +00:00
committed by Rasmus Munk Larsen
parent 00cc497d32
commit 2cd9bb7380
2 changed files with 35 additions and 12 deletions
Download Patch File Download Diff File Expand all files Collapse all files

42
Eigen/src/SparseCore/SparseDenseProduct.h
View File

@@ -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);
}
}

5
test/sparse_product.cpp
View File

@@ -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),