perf: split regime-A subtrees by flops; experiment 0 phase profile

[michel2323]

Not a TASKS.md task. This is experiment 0 of PERFORMANCE.md, run on the owner's request.

Task: PERFORMANCE.md, experiment 0 (baseline split), plus the fix it pointed to.

What was built:

• bench/profile_phases.jl profiles one warm refactorization and solve per harness matrix on CUDA, using CUPTI. It writes bench/profile/phase_split.{md,csv}. The state before the fix is in bench/profile/phase_split_main_9d280d0.*.
• subtree_parallelism is a new analysis tuning knob in Options, default 4096. It is not a parameter string, so the cuDSS API is unchanged. A regime-A subtree may do at most 1/subtree_parallelism of the factorization flops, measured by front_flops in src/symbolic/schedule.jl. 0 restores the previous behaviour.
• bench/comparison is regenerated. performance.md is renamed to PERFORMANCE.md, like PLAN.md and TASKS.md. It has a new "Experiment 0 results" section and a "Tracked issues" table of the issues labelled performance.

Why: on every pglib KKT dump, the whole factorization and the forward solve ran in one 128-thread workgroup. Regime A took any subtree whose stack fit the local-memory budget, with no parallelism criterion. GPU busy time equalled wall time on every row, so the gap was not launch overhead.

SDS time divided by cuDSS time, geometric mean over the harness:

feature	refactorization	solve
Cholesky	9.06× → 3.89×	8.93× → 4.36×
LDLᵀ	29.8× → 9.90×	13.4× → 3.96×
LDLᵀ + 2 IR steps, K2	36.0× → 5.45×	33.3× → 7.57×

On case1354 K2 with LDLᵀ, refactorization goes from 146 ms to 7.2 ms and the solve from 78.7 ms to 2.6 ms. SuiteSparse and Laplacian matrices are unchanged within noise.

Tests: run with julia --project=. -e 'using Pkg; Pkg.test()'.

backend	pass	fail	broken
KA CPU	60801	0	1
CUDA, RTX 4080	all but test_api in the full run	0	1

The broken test is the @test_broken already present in test_refinement. On CUDA, test_api was rerun after the fix below and passes, with 3908 passes on CPU and CUDA together.

Test changes, for review:

• Regime-A kernel tests in test_numeric_cholesky_a.jl and the schedule test's "small example" now pin subtree_parallelism = 0. They exercise the unsplit subtrees exactly as before. A new testset checks the flop limit, the maximality of the subtrees, and the knob's validation.
• Three test_api testsets compare solve layouts bitwise. They now set deterministic_mode = 1. The small test matrices now reach the atomic regime-B forward sweep on CUDA, which is not bitwise repeatable. That is the documented contract of deterministic_mode = 0, but these matrices never hit the atomic path before.

Deviations from PLAN.md / TASKS.md: PLAN §2.3 step 5 defines the regime-A subtrees by the local-memory budget only, and this PR adds a flop limit. PLAN.md is not edited. The owner decides whether the plan text changes.

Follow-up issues opened: #82, KKT refactorization and solve are now bound by the number of tree levels. Next steps are ranked in PERFORMANCE.md: LDLᵀ regime B/C kernels (#75), then launch minimization and the solve path for the KKT systems, which are now bound by the number of tree levels.

Possible conflict with T17 (#80) in src/options.jl and src/symbolic/schedule.jl.

🤖 Generated with Claude Code

https://claude.ai/code/session_01TJrmq7e9jiW4VxBMEVup2E

[claude]

nit (non-blocking): the two descriptions in the docstring disagree by one. "The elimination of column k updates the trailing (f-k)×(f-k) block" for k = 1..w gives Σ_{j=f-w}^{f-1} j², while the formula and the code (sq(f) - sq(f-w)) compute Σ_{j=f-w+1}^{f} j². As a relative work measure it does not matter, and the tests pin the implemented sum, so just make the prose match, e.g. "column k (1-based) scales and updates the trailing (f-k+1)×(f-k+1) block".

[claude]

nit (non-blocking): when total < subtree_parallelism (any matrix with fewer than 4096 multiply-adds, e.g. the 6×6 arrow matrix in the schedule test) maxwork is 0, so regime A is switched off entirely and the tiny tree runs as tree_height regime-B launches instead of one. The PR body says a work floor did not help on the harness, but the harness has no such tiny matrices. Consider maxwork = max(total ÷ par, floor) with a small constant floor (a few thousand flops, i.e. a single workgroup's worth), or at least note the behaviour in the Options table. Not blocking: absolute cost is microseconds and the test documents the current behaviour.

[claude]

VERDICT: APPROVE (head 7b3c851)

Not a TASKS.md task: owner-requested experiment 0 of performance.md plus the fix it pointed to. Reviewed against AGENTS.md and PLAN.md §2.3 step 5.

Checked

• Scope: PLAN.md, TASKS.md, Manifest.toml, Project.toml untouched; no new dependencies (the profile script uses CUDA.jl's built-in CUDA.@profile). The deviation from PLAN §2.3 step 5 (a flop limit on regime-A subtrees in addition to the local-memory budget) is stated in the PR body and left to the owner. No cuDSS parameter or phase strings change; subtree_parallelism is a tuning knob in TUNING_OPTIONS, validated through _parse_tuning (InvalidValueError on -1 and 1.5), copied by Base.copy, documented in the Options table.
• Correctness of the split: subtree flops are accumulated child-to-parent in Int64, which is valid because the supernode numbering already puts children before parents (the existing peak/eligibility loop relies on the same order). eligible[s] is monotone (children eligible, work within limit), so the roots picked by the existing (p == 0 || !eligible[p]) test are still the maximal eligible subtrees; 0 restores typemax(Int64) and the previous behaviour. No kernel, layout or device-side change, so the KA/genericity rules are unaffected.
• Tests: the new subtree_parallelism testset checks the knob, front_flops values, the flop limit, regime shrinkage and subtree maximality on laplacian2d(60,60) for 0/64/4096. Existing regime-A kernel tests pin subtree_parallelism = 0 and therefore exercise exactly the subtrees they did before; nothing was weakened. The three test_api testsets that compare layouts bitwise now set deterministic_mode = 1, which is correct: _solve_phase! reads opts.deterministic_mode at solve time and the default atomic regime-B forward sweep is not bitwise repeatable on a GPU; the matrices only reach that path now that regime A is split. setparam! on a DirectSolver accepts config parameters after factorization.
• Bench: bench/profile_phases.jl uses the existing internal fields (numeric.plan.sub_first/group_first, schedule.regime/level/subtree_ptr); documented in bench/README.md; before/after profiles and regenerated comparison data committed.
• Report: claims match the diff. I could not run the suite in this review session (no execution approval); CPU and CUDA CI are pending on this head and are handled by the pipeline.

Blocking findings: none.

Non-blocking nits (inline): front_flops docstring prose is off by one from its formula; with total < subtree_parallelism the flop cap is 0 and regime A is disabled entirely for tiny matrices (a small floor might be worth considering).