From ade01a1775f1911c0bf373190001923f452ad6cc Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Fri, 31 Jul 2026 10:12:45 -0700 Subject: [PATCH 1/2] Fix out-of-bounds read in expandSparsityPattern for empty block columns The diagonal-block detection in CSCMatrix::expandSparsityPattern reads Ai[Ap[block_j + 1] - 1], which for an empty block column reads the last entry of the previous column -- or Ai[-1] when the empty column precedes any nonzeros. When the stray value happens to equal block_j, the column size computation (numBlocks - 1) * N + 1 underflows with numBlocks == 0, corrupting the InOrderBuilder's column sizes and causing out-of-bounds writes (intermittent SIGSEGVs in BlockCSCHessian::toScalar/toEigen). FE Hessians always have diagonal blocks (every node belongs to an element), but patterns built from contact stencils are mostly empty columns -- only vertices currently in contact appear -- which is how this was found (AddressSanitizer repro: a SystemAssembler<3> blockSparsityPattern over a single 2-vertex stencil among 1000 block variables, followed by toScalar()). Guard the detection on numBlocks > 0; the filler loop below is already safe for empty columns. --- src/lib/MeshFEMSparse/SparseMatrices.hh | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/lib/MeshFEMSparse/SparseMatrices.hh b/src/lib/MeshFEMSparse/SparseMatrices.hh index e72b20d..18b98ae 100644 --- a/src/lib/MeshFEMSparse/SparseMatrices.hh +++ b/src/lib/MeshFEMSparse/SparseMatrices.hh @@ -1171,7 +1171,11 @@ struct CSCMatrix { for (_Index block_j = 0; block_j < blockHsp.n; ++block_j) { size_t gvar_j = block_j * N; size_t numBlocks = blockHsp.Ap[block_j + 1] - blockHsp.Ap[block_j]; - bool hasDiagonal = AssumeDiagonalExists || (blockHsp.Ai[blockHsp.Ap[block_j + 1] - 1] == block_j); + // Guard against empty block columns (possible, e.g., for contact + // Hessians, where most vertices are collision-free): reading + // Ai[Ap[block_j + 1] - 1] would be out of bounds, and a spurious + // `hasDiagonal` would underflow `colSize` below. + bool hasDiagonal = (numBlocks > 0) && (AssumeDiagonalExists || (blockHsp.Ai[blockHsp.Ap[block_j + 1] - 1] == block_j)); if (hasDiagonal) { size_t colSize = (numBlocks - 1) * N + 1; for (size_t c_j = 0; c_j < N; ++c_j) From 15a92834189ade69ed9e00373adc3036a72cd40a Mon Sep 17 00:00:00 2001 From: Zachary Ferguson Date: Fri, 31 Jul 2026 13:03:13 -0700 Subject: [PATCH 2/2] Skip empty block columns when visiting diagonal scalar entries visitDiagonalScalarEntries walked every block column and took diagBlockScalarLoc() for each. For an empty block column that offset is N^2 * Ap[bj] - N^2, which points into the preceding column's storage, or below the start of Ax entirely when the first block column is empty. trace() therefore summed unrelated values (and read out of bounds) on any matrix with empty block columns, and addDiag()/setDiag() wrote into the wrong entries. This is the same assumption behind the out-of-bounds read in expandSparsityPattern fixed in the previous commit: FE Hessians always have a diagonal block per column because every node belongs to an element, but Hessians assembled from contact stencils leave most columns empty, since only vertices currently in contact appear. Skip empty columns while still advancing the scalar column index, so trace() ignores their (structurally zero) diagonals. The mutating operations cannot be fixed by skipping, because there is no stored entry to write, so they now check for the missing blocks and throw instead of corrupting neighboring columns. Note that missingRequiredDiagonalBlocks(), and hence assertSupportsAssembly(), does not catch this: it excludes the StoreFullDiagonalBlocks case, where diagBlockScalarLoc() is equally invalid for an empty column. I left that alone rather than widen it, since assembly itself is unaffected -- the assembler only touches columns a stencil references, and those always contain their diagonal block. --- src/lib/MeshFEMSparse/BlockCSCHessian.hh | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/src/lib/MeshFEMSparse/BlockCSCHessian.hh b/src/lib/MeshFEMSparse/BlockCSCHessian.hh index 150cc15..0e0b68b 100644 --- a/src/lib/MeshFEMSparse/BlockCSCHessian.hh +++ b/src/lib/MeshFEMSparse/BlockCSCHessian.hh @@ -833,11 +833,18 @@ struct MESHFEM_EXPORT BlockCSCHessian final : public BlockToScalarPolicyDefault< // Call f(j, loc), passing the location in `Ax` of the diagonal entry // in scalar column j, for each j in 0..numScalarVars() - 1 + // Note: empty block columns are skipped. Their diagonal entries are not + // stored (so there is no location to hand to `f`), and + // `diagBlockScalarLoc()` would compute an offset into the preceding + // column's storage for them. `j` still advances so callers see the correct + // scalar column indices. Callers that must touch every diagonal entry + // should call `assertAllDiagonalBlocksPresent` first. template void visitDiagonalScalarEntries(F &&f) const { _Index j = 0; for (_Index bj = 0; bj < n; ++bj) { auto cs = columnScanner(bj); + if (col_nnz(bj) == 0) { j += cs.colBlockSize(); continue; } _Index loc = cs.diagBlockScalarLoc(); for (_Index c_j = 0; c_j < cs.colBlockSize(); ++c_j) { f(j++, loc); @@ -847,6 +854,14 @@ struct MESHFEM_EXPORT BlockCSCHessian final : public BlockToScalarPolicyDefault< } } + // Throw unless every block column has a diagonal block, which the + // diagonal-mutating operations below require: for a column without one + // there is no stored entry to write. + void assertAllDiagonalBlocksPresent(const char *op) const { + if (numDiagonalBlocks() < size_t(n)) + throw std::runtime_error(std::string("BlockCSCHessian::") + op + ": matrix is missing diagonal blocks; insert them first (see BorderedSparseHessian::insertSparsityPatternDiagonalBlocksIfNeeded)"); + } + virtual Real trace() const override { Real result = 0; visitDiagonalScalarEntries([&result, this](size_t /* j */, _Index loc) { result += Ax[loc]; }); @@ -1319,14 +1334,17 @@ private: VarStructure m_vars; virtual void m_addDiag(const _Real *d) override { + assertAllDiagonalBlocksPresent("addDiag"); visitDiagonalScalarEntries([d, this](size_t j, _Index loc) { Ax[loc] += d[j]; }); } virtual void m_addDiag(_Real d) override { + assertAllDiagonalBlocksPresent("addDiag"); visitDiagonalScalarEntries([d, this](size_t /* j */, _Index loc) { Ax[loc] += d; }); } virtual void m_setDiag(_Real d) override { + assertAllDiagonalBlocksPresent("setDiag"); visitDiagonalScalarEntries([d, this](size_t /* j */, _Index loc) { Ax[loc] = d; }); }