Fix out-of-bounds read in expandSparsityPattern for empty block columns

[zfergus]

CSCMatrix::expandSparsityPattern<UniformBlockSize>() detects the presence of a diagonal block with

bool hasDiagonal = AssumeDiagonalExists || (blockHsp.Ai[blockHsp.Ap[block_j + 1] - 1] == block_j);

For a block column with zero entries (Ap[block_j] == Ap[block_j + 1]), this reads the last row index of the previous non-empty column — or Ai[-1] when the empty column precedes any nonzeros. When the stray value happens to equal block_j, hasDiagonal is spuriously true and

size_t colSize = (numBlocks - 1) * N + 1; // numBlocks == 0 → wraps around

underflows to a huge size_t, corrupting the InOrderBuilder's column-size array and causing out-of-bounds writes (we observed intermittent SIGSEGVs inside BlockCSCHessian::toScalar/toEigen).

Why FEM never hits this: every node belongs to an element, so every block column has at least its diagonal block. We hit it using SystemAssembler::blockSparsityPattern for IPC contact Hessians, where the overwhelming majority of block columns are empty (only vertices currently in contact appear).

AddressSanitizer repro (fails before this change, clean after):

SystemAssembler<3> assembler(1000);
auto stencil = [](size_t) {
    ElementBlockVarsWithSizeRange<1, 4> s(2);
    s[0] = 500; s[1] = 501;
    return s;
};
auto H = assembler.blockSparsityPattern(1, stencil); // 998 empty block columns
H->setZero();
auto S = H->toScalar(); // heap-buffer-overflow READ in the counting lambda

SUMMARY: AddressSanitizer: heap-buffer-overflow SparseMatrices.hh:968 in
MeshFEM::CSCMatrix<...>::InOrderBuilder::InOrderBuilder<... expandSparsityPattern<3ul, false>() ...>

Fix: guard the detection on numBlocks > 0 (the filler loop below already iterates the empty range safely).

A possibly-related observation while reading BlockCSCHessian::toScalar: the uniformBlockSize() fast path computes its result and then falls through without returning it —

if (uniformBlockSize()) {
    CSCMat result = this->template expandSparsityPattern<VarStructure::MaxBlockDim>();
    copyValuesInto(result);
} // result discarded; the generic path below runs anyway

so the fast path's work is discarded (it was also how the OOB above got executed even though its result was unused). Possibly a missing return result; — left out of this PR since "fixing" it changes which code path serves all uniform-block-size callers, and you're best placed to judge intent.

🤖 Generated with Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[zfergus]

Pushed a second commit to this PR: the same assumption shows up again in visitDiagonalScalarEntries, and I hit it while exercising the fixed toScalar path.

for (_Index bj = 0; bj < n; ++bj) {
    auto cs = columnScanner(bj);
    _Index loc = cs.diagBlockScalarLoc();   // <-- for an empty column?

diagBlockScalarLoc() is scalarOffsetForColumn(bj + 1) - diagBlockSize(), which for an empty block column is N²·Ap[bj] - N². That points into the preceding column's storage, or below the start of Ax when the first block column is empty.

The visible symptom was trace(). On a contact Hessian I get 1951.93 where the equivalent dense trace is 447.82, because empty columns re-count the previous column's diagonal. apply() on the same matrix matches Eigen exactly, which is what narrowed it down: the values are fine, only the diagonal offsets are wrong. addDiag() and setDiag() write to those same bad locations.

The commit skips empty columns while still advancing the scalar column index, so trace() ignores their structurally-zero diagonals. Skipping is not a fix for the mutating operations, since there is no stored entry to write at all, so those now check and throw rather than corrupting a neighboring column.

One thing I deliberately did not change: missingRequiredDiagonalBlocks() excludes the StoreFullDiagonalBlocks case, so assertSupportsAssembly() never fires for these matrices even though diagBlockScalarLoc() is equally invalid there. Widening it would make setZero() throw for legitimately sparse patterns, and assembly itself is fine, since ElementHessianContribAssembler only touches columns that a stencil references and those always contain their diagonal block. Happy to widen it if you would rather the precondition be strict, but it seemed like it would break working uses.

Both fixes are validated through the IPC Toolkit's block-CSC assembly backend (ipc-sim/ipc-toolkit#246): trace() now agrees with the dense trace to 1e-12, addDiag() throws as expected, and the apply() round-trip against Eigen is unchanged.

🤖 Addressed by Claude Code

[jpanetta]

Thanks! Matrices with missing diagonal blocks are explicitly unsupported because random access to the noncontiguous storage layout cannot be implemented efficiently without assuming the presence of all diagonal blocks.

However, in the contiguous layout, omitting the diagonal blocks should work for most operations (apart from, e.g., the issues you found!). These patches here look reasonable and shouldn't degrade performance, so I will merge this with the warning that this use case is untested.