Make sparsity fills output-driven in block left-multiply and csr@csc matmul

[Transurgeon]

What changed

block_left_multiply_fill_sparsity and csr_csc_matmul_alloc build the output
sparsity pattern of C = A @ B. Previously they were input-driven: every
candidate output entry (i, j) ran a has_overlap merge-scan of row i of A
against column j of B — m x p scans regardless of how sparse the operands
are. They are now output-driven (the symbolic phase of Gustavson's sparse
matmul): an output row's pattern is the union of the operand index lists its
entries select, so we only walk index lists that exist. Candidates are
deduplicated with a generation-stamped workspace and sorted, so the produced
pattern is byte-identical to the old code's. The shared gather lives in one
helper, sorted_pattern_union. Each function builds one transposed view of the
constant operand up front; both functions run once per problem at derivative-
structure initialization, so that cost is amortized. The per-iteration
*_fill_values paths are untouched.

Complexity per output row goes from O(m) (resp. O(p)) overlap scans to
O(flops for that row + sort of its pattern).

Correctness

• The gather reproduces the old scan's output exactly (same dedup semantics,
  same ascending order). Unit tests include hand-built cases where a missing
  dedup or missing sort would each fail deterministically, plus randomized
  cross-checks against the original has_overlap implementation (kept verbatim
  in the test file as a reference): identical p/i arrays asserted over
  random trials for both functions.
• Full engine test suite passes (401 tests).
• End-to-end A/B: 17 problems mirroring the cvxpy nlp_tests suite and DNLP
  doc examples were run through the cvxpy -> sparsediffpy path against two
  builds identical except for this file. Objective, gradient, constraint
  values, and the full Jacobian and Hessian sparsity patterns and values
  are exactly equal (np.array_equal) on every problem, in two independent
  runs. (Comparison is NaN-aware: two division-by-quad_form formulations
  evaluate to NaN at the harness's arbitrary evaluation point — identically,
  in the same positions, in both builds.)

Benchmarks (DNLP suite, end-to-end through cvxpy)

"kernel" = wall time inside the two rewritten functions, measured with
CLOCK_MONOTONIC accumulators compiled identically into both builds
(benchmark-only instrumentation, not part of this PR). "deriv. init" =
problem_init_jacobian_coo + problem_init_hessian_coo_lower_triangular.
Times are min over 2 independent runs x 3 in-process repetitions.

problem	kernel calls (block / csr@csc)	kernel old (ms)	kernel new (ms)	kernel speedup	deriv. init old (ms)	deriv. init new (ms)	init speedup
A @ x == b, sparse A 5000x5000 (0.2%)	1 / 0	973.25	1.29	752x	975.01	2.55	381.84x
signal recovery, sparse A 2000x6000 (1%)	1 / 0	1,848.06	2.70	683x	1,896.65	50.49	37.57x
convolve(k, x), kernel 100, signal 2000	0 / 1	939.28	3.08	305x	1,814.09	882.11	2.06x
logistic regression, sparse A 10000x2000 (0.5%)	1 / 0	745.03	1.93	386x	1,550.35	809.82	1.91x
logistic regression, sparse A 2000x785 (1%)	1 / 0	45.00	0.33	136x	98.26	54.04	1.82x
A @ x == b, sparse A 1000x1000 (1%)	1 / 0	33.83	0.23	147x	34.09	0.46	73.69x
AC power flow, 9-bus (test suite)	0 / 2	<0.1	<0.1	--	0.53	0.52	1.02x
min |X @ Y - C|_F^2, 40x20x30	0 / 0	~0	~0	--	218.33	216.98	1.01x
t == sum(X @ Y) (test suite, 5x7x11)	0 / 0	~0	~0	--	0.07	0.07	1.01x
risk parity, n=8 (test suite)	0 / 0	~0	~0	--	0.06	0.06	1.01x
Sharpe ratio, n=100 (test suite)	0 / 0	~0	~0	--	0.18	0.19	0.99x
signal recovery, dense A 200x600	0 / 0	~0	~0	--	0.92	0.93	0.99x
max entropy sum(entr(A @ q)), dense A 100x100	0 / 0	~0	~0	--	0.03	0.02	1.05x
chemical equilibrium (doc example, verbatim)	0 / 0	~0	~0	--	0.03	0.03	1.01x
trace(exp(A @ X)), dense A 10x10 (stress suite)	0 / 0	~0	~0	--	0.11	0.11	1.02x
sum(exp(X @ A)), dense A 10x10 (stress suite)	0 / 0	~0	~0	--	0.13	0.13	1.03x
path planning w/ obstacles, n=50 (doc example)	0 / 0	~0	~0	--	0.26	0.26	0.99x
total		4,584.48	9.59	478x	6,589.09	2,018.76	3.26x

Reading the table honestly:

• On sparse problems at realistic scale the old kernels were the
  init-phase bottleneck (e.g. 97% of signal-recovery init) and are now
  negligible; where init speedup is below the kernel speedup, the remaining
  time is in other init work (e.g. Hessian assembly) that this PR does not
  touch.
• The call-count column shows which problems reach these kernels at all:
  only sparse constant matmuls (block fill) and convolve / power flow
  (csr@csc). Most suite problems at their natural unit-test sizes use dense
  data, which takes the permuted_dense path (0 calls) — on all of those,
  outputs are identical and timings differ only by noise (0.99x-1.05x).
• The bivariate X @ Y problems also made 0 calls: the matmul atom's
  csr_csc_matmul_alloc is on its composite-children Hessian path, and these
  tests use leaf variables.
• Four more suite problems (log_sum_exp x2, huber, sum_largest) could
  not run on the direct C-problem harness because those atoms are lowered by
  the dnlp2smooth reduction before reaching the engine in real solves; their
  lowered forms use the same dense-path atoms as the rows above.
• Per-iteration eval times are unaffected: an interleaved A/B/A/B measurement
  of the largest eval (bivariate X @ Y Hessian) shows the new build within
  noise (marginally faster: 157-159 ms vs 162-164 ms).
• The engine's own PROFILE_ONLY suite agrees: profile_left_matmul Jacobian
  init 0.238 s -> 0.028 s; profile_memory derivative-structure phase
  2.22 s -> 1.91 s; everything else unchanged.

Known trade-off

Kernel-level microbenchmarks (old vs new on random matrices, outputs verified
identical in every case) show the win grows as operands get sparser (230x at
0.1% density, 11-23x at 1%, size-independent), but when both operands are
dense-ish the old scan is faster: at 20% density old wins ~2x, and on fully
dense inputs ~70x (its has_overlap early-exits on the first match, O(1) per
pair, while the gather does full O(n) work per pair). This regime is not
reachable through cvxpy for constant matmuls — SPARSE_DENSITY_THRESHOLD = 0.05 routes denser constants to the dense (permuted_dense) path — and dense
operands should use that path anyway. If we ever want insurance, a cheap
per-row guard (sum the selected index-list lengths; fall back to the
has_overlap scan when it exceeds the output width) is possible as follow-up.

Environment

Intel Core i7-8850H (2.6 GHz), macOS 15.7.7, Apple clang 15, Release (-O3,
-DNDEBUG), Python 3.14 / cvxpy nlp branch. Benchmark harness: two sparsediffpy
builds identical except src/utils/linalg_sparse_matmuls.c (old = main,
new = this branch), swapped via PYTHONPATH.

[Transurgeon]

@dance858 another PR for you to review.
While trying to speedup a problem that used Kronecker products (as in the atom), Claude noticed that the sparse matmul with the Kronecker product would be much faster using Gustavson's algorithm (even with sorting).
I then stress-tested this new sparse matmul in all our DNLP problems and it seems to be faster throughout.
I think it could be nice to switch to this (especially since dense problems use permuted_dense now).

I've also added two correctness tests which 1. compare with the old code and 2. have both a unsorted and duplicated row gathers to ensure it is deduplicating and sorting the output.
I haven't looked too much in the implementation, but it seems quite concise and I at least understand it at a high level.
Let me know if you have any thoughts.. and sorry to give you more work ;).