Wavefront-parallel numeric LDL factorization (pixiu series 4/4)

[batterseapower]

What: Parallelizes the numeric LDL factorization (42% of solve time) across rows of equal elimination-tree height ("wavefronts"), keeping the serial code verbatim for small problems (<20k dims) and logical factorization.

Safety argument (documented in the code): row k's pattern consists only of its etree descendants; equal-height rows have disjoint descendant sets; each column's entries are appended by its ancestor chain in strictly increasing height = strictly increasing row order, so CSC ordering is preserved and at most one row per wavefront touches any column. Shared outputs are mutated through documented SendPtr raw pointers; each worker has private y workspace.

Determinism: factors are bitwise identical to the serial code (verified max|Δx| = 0 over full solves on two reference problems) — every entry is computed from the same inputs in the same inner-loop order, independent of scheduling.

Measured: ~5–10% full-solve. On our KKT, 92% of Lnz-weighted work sits in parallelizable levels, but the kernel is bound by short dependent random accesses (avg L column ~7 entries) and saturates at ~4 threads. Also tested and rejected on the same grounds (see branch pixiu/99-replay-experiment, bitwise-verified, zero gain): target-cpu=native, static-structure replay refactorization with u32 indices, software prefetch of scatter targets.

Series result: 2.17s → ~1.8s per solve plus zero fallback re-solves. Remaining gap to mosek (0.75s) is structural; next investigation is a condensed-KKT (Woodbury) backend exploiting the diagonal + low-rank structure of our constraint matrix.

[batterseapower]

Closing on measured evidence: I implemented this same approach independently and it produces no speedup at all on this problem class, for a structural reason that will not go away with tuning.

I built the identical scheme — group columns by elimination-tree height, replay each wavefront concurrently — on top of the refactorization schedule from #6/#9, using std::thread::scope (no new dependency), with per-thread accumulators. Correctness was fine: a test asserted refactorization is bit-identical across 1/2/4/8 threads, which the disjoint-descendant-set argument guarantees, and it passed.

Timings on the largest portfolio problem, max_threads = 1/2/4/8:

1.299s / 1.304s / 1.307s / 1.299s

Flat. The wavefront diagnostic says why:

	n=61679 (full)	n=32069 (screened)
wavefronts	159	222
total update work	21.1M elements	11.2M
work in wavefronts ≥1e6 elements	0 (0.0%)	0 (0.0%)
heaviest 5 wavefronts	518k, 506k, 450k, 447k, 446k	272k, 265k, 247k, 247k, 246k
columns in those wavefronts	1, 1, 1, 1, 1	1, 1, 1, 1, 1

The elimination tree is a chain through the dense border. The heaviest wavefronts contain exactly one column each, so height-based parallelism finds no concurrency where the work actually is. Intra-column parallelism is also unavailable: a column's update steps are sequentially dependent through the sparse accumulator (step cidx writes the y_vals entries that later steps read) — it is a triangular solve.

This is a property of KKT systems with a dense-ish row border, i.e. exactly the factor-model structure these problems have; it is not specific to my implementation. The approach could still pay on problems with bushy elimination trees, so the idea isn't wrong in general — just ineffective here.

What did work on the same 40% of runtime was vectorising rather than parallelising it: #9 replays the factor updates over consecutive-row runs (95.7% of flop-weighted entries lie in runs ≥ 8, from that same dense border), which is bit-identical and −12.7%. Multi-core for the factor would need a supernodal/multifrontal front where the dense border becomes a dense block, parallelised inside the BLAS3 kernel — sidestepping the etree chain entirely.

Also note this PR is based on pixiu/03-best-iterate, whose PR (#3) is closed in favour of #7.