perf: cooperative pivot search in the device LDLᵀ kernels (#75)

[michel2323]

Refs #75

Task: #75, experiment 1 of PERFORMANCE.md, step 1 of 3 — parallel pivot search in the device LDLᵀ/LDLᴴ kernels (T25 owner note, item 1).

What was built:

• src/numeric/ldlt.jl: the pivot search of front_ldlt_kernel! and subtree_ldlt_kernel! (was _lt_choose on work item 1) is a cooperative search: pass 1 (λ, r and the column maximum of column k), pass 2 (σ and the column-r maxima of Bunch–Kaufman) and pass 3 (the _best_1x1 fallback, one lane per candidate column) are workgroup reductions in @localmem (1-D, no warp intrinsics, two-stage fold 128 → 16 → 1, only over lanes that can hold data) with the reference's tie-breaking (first maximum), so every decision equals the serial one. Work item 1 takes the pivot from the values before the interchange and the interchange phase writes the new diagonal, which keeps five barriers per step; pass 1 of the next step runs in the scale phase.
• src/reference/ldlt.jl: _best_1x1 skips columns that cannot beat the best so far, and the threshold test exits at the first violation (_threshold_ok; u ≥ 0 is finite, so d ≥ u·max ⇔ d ≥ u·|x| for every entry, NaN included). Same results, the equality tests prove it.
• src/symbolic/schedule.jl: SUBTREE_LOCAL_RESERVE_LDLT = 1024 and subtree_local_reserve(structure): the regime-A LDLᵀ kernel keeps 1024 B (was 256) of its local memory for the control words and reduction slots; build_schedule(…; reserve) and subtree_capacity(…, reserve), passed by the analysis for "S"/"H". Cholesky schedules are unchanged.
• PERFORMANCE.md: tracked-issues table, "Experiments 1–2 results" with the step-1 table.

Tests: julia --project=. -e 'using Pkg; Pkg.test()' with CUDA in the local test environment (RTX 4080, not committed): 104914 pass / 0 fail / 2 broken (the T16 @test_broken, once per backend) before the last commit (lane-bounded folds); after it test_numeric_ldlt, test_symbolic_schedule, test_api, test_ported, test_refinement, test_ubatch, test_solve again on CPU and CUDA: 15020 + 25371 pass, 0 fail, 2 broken. Rebased on T18 (#83, no overlap). test/test_symbolic_schedule.jl: check_schedule uses the structure's reserve (required by the change, no assertion loosened). No allocation change (factorize! LDLᵀ 64 B on the CPU backend).

Measurements (RTX 4080; refactorization best of 5 warm runs, bench/profile_phases.jl; *: bench/compare.jl median where the profile skips the matrix):

matrix	main ms	step 1 ms	step 1 / cuDSS	step 1 / SDS Cholesky
lap2d_300	228	192	37.4×	6.12×
lap3d_40	14,777*	14,494*	248×	64.7×
HB/bcsstk17	206	167	49.5×	6×
Boeing/bcsstk38	137	106	32.1×	4.69×
GHS_psdef/apache2	92,576*	97,296*	199×	94.5×
kkt_pglib_opf_case118_ieee_condensed_1	1.77	1.82	4.66×	1.99×
kkt_pglib_opf_case118_ieee_condensed_10	1.72	1.79	5.13×	2.04×
kkt_pglib_opf_case118_ieee_condensed_20	1.87	1.82	5.18×	2.08×
kkt_pglib_opf_case118_ieee_k2_1	1.69	1.69	4.66×
kkt_pglib_opf_case118_ieee_k2_10	1.62	1.68	4.74×
kkt_pglib_opf_case118_ieee_k2_20	1.6	1.66	4.16×
kkt_pglib_opf_case1354_pegase_condensed_1	10.9	7.2	10.5×	2.46×
kkt_pglib_opf_case1354_pegase_condensed_10	10.9	7.3	11.4×	2.5×
kkt_pglib_opf_case1354_pegase_condensed_20	11	7.32	10.4×	2.51×
kkt_pglib_opf_case1354_pegase_k2_1	8.31	6.6	8.08×
kkt_pglib_opf_case1354_pegase_k2_10	7.15	6.05	7.32×
kkt_pglib_opf_case1354_pegase_k2_20	9.1	6.71	8.32×
kkt_pglib_opf_case14_ieee_condensed_1	0.446	0.503	2.61×	1.38×
kkt_pglib_opf_case14_ieee_condensed_10	0.443	0.52	2.26×	1.53×
kkt_pglib_opf_case14_ieee_condensed_11	0.448	0.463	1.59×	1.41×
kkt_pglib_opf_case14_ieee_k2_1	0.669	0.697	2.49×
kkt_pglib_opf_case14_ieee_k2_10	0.722	0.708	2.98×
kkt_pglib_opf_case14_ieee_k2_15	0.97	1.03	3.2×

• K2 dumps (CUDA): inertia and nperturbed identical to main on all nine.
• kkt_matrix(Float64, 3000, 1000, 1e-8), default analysis (the T15 number): CUDA 33.9 s → 8.3 s; KA CPU backend 7.0 s → 7.0 s against ref_ldlt! 2.6 s (work items of a workgroup run one after another there, so a parallel search cannot help).
• Experiment-1 criterion (LDLᵀ within 1.5× of SDS Cholesky): not met by this step (lap2d 6.1×, bcsstk17 6.0×, case1354 condensed 2.5×); the remaining cost is regime B/C (steps 2 and 3).

Deviations from PLAN.md / TASKS.md: the regime-A LDLᵀ local-memory reserve is an analysis change for "S"/"H" (1024 B instead of 256 B per budget class), needed for the reduction slots in the 48 KiB class. The smallest KKT dumps (case14, w ≤ 16) are 3–17% slower (reduction barriers); apache2/lap3d_40 within ±5% (one workgroup per huge front, step 3).

Follow-up issues opened: #86 (found-by-agent, pre-existing on main): the device LDLᵀ differs from ref_ldlt! on the K2 dumps in pivot sequence and nperturbed (rounding with max|L| 1e14–1e16, #71); inertia equal.

🤖 Generated with Claude Code

https://claude.ai/code/session_01Cqp7TUSyCEdjPtLDEftJf5

[claude]

nit (non-blocking): for pass 1 the lanes that can hold data are the rows k+1:f, i.e. f - k lanes (lane li covers row k + li), but _lt_nused returns f - k + 1 for P == 1 as for pass 2 (rows k:f). Harmless, since _lt_pass1! writes the neutral values (0, 0, 0) into every lane ≤ NL, so the extra lane folds as a no-op; the correctness of the fold never depends on the bound. If you want the comment and the bound to agree, use Int(ctl[_ST_F]) - k for pass 1.

[claude]

VERDICT: APPROVE (head f47a2d6)

Perf PR for issue #75, step 1 of experiment 1 (T25 owner note, item 1); not a TNN task, so there is no TASKS.md Report block and the PR body serves as the report.

Checked

• Scope: only the cooperative pivot search (src/numeric/ldlt.jl), the matching early-exit threshold test and _best_1x1 pruning in the reference, the regime-A LDLᵀ local-memory reserve (SUBTREE_LOCAL_RESERVE_LDLT, subtree_local_reserve, build_schedule(…; reserve), subtree_capacity(…, reserve), both analysis callers), and the PERFORMANCE.md tables. PLAN.md, TASKS.md and Manifest.toml untouched; no public name, parameter or phase string changed. The reserve is used consistently: the LDLᵀ launch uses LB - SUBTREE_LOCAL_RESERVE_LDLT, the schedule capacity for "S"/"H" subtracts the same 1024 B; Cholesky keeps 256 B. The 944 B figure in the docstring matches ctl (22 × 8) + pv (8 × 16) + red (48 × 8) + redi (32 × 8).
• Correctness of the search, pass by pass against _choose_pivot: pass 1 (λ, r = first maximum over rows k+1:w, column max over k+1:f), pass 2 (σ over k:w but r; column-r max over k:f but r, k; column-k max over k:f but k, r, so _accept_1x1(F, r, k) = max(mr, |F[k, r]|) and _accept_2x2 gets m1 = mk, m2 = mr), pass 3 (_best_1x1 with the first-maximum tie-break, the pivot_type = "D" largest-diagonal fallback with the NaN-at-k case). Per-lane strict > plus the fold rule (greater value, or equal value and smaller index; index 0 = none) reproduces the serial first maximum; max propagates NaN as _colmax did; _threshold_ok’s early exit is equivalent because the prefix max is monotone and u ≥ 0. _lt_select! reads the 2×2 block (F[k,k], F[r,k], F[r,r]) and psign[perm[piv[c0+c-1]]] before the interchange, which equals the post-swap values the old _lt_pivot! read; _lt_swap! writes the new diagonal/zero from the work item that owns row k (and k+1), after its own swap writes, so no cross-lane race. End-of-step pass 1 reads column K+STEP while _lt_scale! writes columns K..K+STEP−1 only. Phase branches read ctl[_LT_PHASE] after a barrier (uniform), so the nested @synchronize are safe; the reductions have a fixed order (deterministic).
• Conventions: 1-D workgroups, @localmem sized by Val/constants, no warp intrinsics, no atomics, no allocation or host sync in the kernels, generic T/INT (red::R, redi::IT), no @allowscalar.
• Tests: none weakened; check_schedule passes the structure’s reserve as the change requires. The existing equality suite (piv, pivot_kind, D, L, stats vs ref_ldlt! on BACKENDS × ELTYPES, three regime mixes, "D"/"N" pivot types, perturbation) exercises all three passes on both kernels. The PR body’s counts (104914 pass / 2 broken, the T16 @test_broken) are plausible; CPU and CUDA CI are still pending on this head and are handled by the pipeline.
• Follow-up: #86 exists, open, labelled found-by-agent (not yet triaged, which will hold the chain after merge, as designed).

Blocking findings: none.

Non-blocking: one inline nit on _lt_nused (pass 1 bound is one lane over; harmless because every lane holds a neutral value).