Make block_left_multiply_fill_values linear in the multiply-adds

[Transurgeon]

Problem

block_left_multiply_fill_values computed each entry of C as a merge-based sparse dot of a whole row of A with the matching block of a column of J. The cost is O(nnz(C) · row length), which is quadratic for a long dense row.

c @ x with c dense of length n, on the CSR path (cvxpy, DIFFENGINE with constants forced to CSR), one first canonicalization on an Apple M2:

n	before	after	dense (permuted_dense)	cvxpy CPP / SCIPY / COO
1e5	3.2 s	0.02 s	0.02 s	0.05–0.07 s
4e5	49 s	0.065 s	0.05 s	—
1e7	≈ 8.5 h (extrapolated)	1.59 s	1.36 s	6.4–7.4 s

The symbolic counterpart (block_left_multiply_fill_sparsity) was already made output-driven in c549c64. The value fill was not.

Fix

The fill is now Gustavson's numeric phase through A's CSC mirror. For each block of a column of C, it scatters A[:, c] * J[c, j] into an accumulator over A's rows, then gathers into the block's entries of C. The cost is now the number of multiply-adds.

• block_left_multiply_fill_values_csc(A_csc, J, C, acc) is the new kernel. acc holds A_csc->m doubles and needs no initialization.
• block_left_multiply_fill_values(A, J, C) keeps its signature. It builds the CSnvenience for one-off callers.
• sparse_matrix's block_left_mult_values reuses the existing version-guarded csc_cache (through refresh_csc_values) and a lazily allocated accumulator (bl_acc). A constant A therefore converts once. Sparse left_matmul matrices are always constant, since sparse parameters are rejected, but the version guard would reconvert if that changed.

Each row's terms are still added in increasing column order of A, the same order avalues are bit-identical.

Tests

• Three new hand-computed cases in tests/utils/test_linalg_sparse_matmuls.h:
  • test_block_left_multiply_values_two_blocks: two blocks, with an entry of C that sums several products.
  • test_block_left_multiply_values_dense_row: the c @ x gradient that made the merge kernel quadratic.
  • test_sparse_matrix_block_left_mult_values_refill: the sparse_matrix path w refilled with new J values, then with new A values aftermatrix_values_changed.
• all_tests: 464/464 pass, also under -fsanitize=address,undefined. A 1e-6 perturbation of the kernel output makes the suite fail.
• End to end through cvxpy, on all 25 problems of the cvxpy benchmark suite: SHA-256 hashes of (P, c, A, b) match the previous engine (64f7432), both with default routing and with every constant forced to CSR. The exception is SimpleLPBenchmark on the CSR path, which previously never finished; it now matches the dense path's output exactly.
• clang-format -style=file: no replacements on the changed files.

🤖 Generated with Claude Code

[Transurgeon]

is this necessary to add to every sparse matrix? Also where is it lazily allocated?

[Transurgeon]

what do you mean by one-off convenience form?

[Transurgeon]

can't we just replace the original content of block_left_multiply_fill_values? instead of creating a new method?

[Transurgeon]

@dance858 I think we should focus on this PR first.
I encountered it while testing the "ablation" and turning off the matrix classes.
Our multiply fill values would have quadratic complexity.