michel2323 · Oct 5, 2026 · Oct 3, 2026 · Oct 5, 2026
diff --git a/performance.md → PERFORMANCE.md b/performance.md → PERFORMANCE.md
@@ -2,6 +2,20 @@
 
 **Bottom line:** For MadNLP-style repeated KKT solves, SparseDirectSolver.jl (SDS) will not beat cuDSS by out-BLASing it on large fronts. Its realistic edge comes from three places. First, robust indefinite numerics that cuDSS lacks: true 2×2 Bunch–Kaufman pivots, per-row pivot signs, and correct inertia. Second, a refactorize+solve path that is device-resident, launch-minimal and graph-captured, which removes the host synchronizations and launches that dominate at OPF scale. Third, a GPU-resident or cached analysis phase. The immediate blockers are not raw FLOP rate. They are (a) the serial, reference-faithful pivot search and the one-workgroup-per-front regime-C LDLᵀ (issue #75), and (b) missing matching and scaling, without which K2 systems produce max|L| of 1e14–1e16 (issue #71).
 
+## Tracked issues
+
+Performance issues carry the GitHub label `performance` ([list](https://github.com/exanauts/SparseDirectSolver.jl/issues?q=label%3Aperformance)). Keep this table in step with them: add a row when an issue is opened, and update the status when its PR merges or it closes.
+
+| issue | what | experiment | status |
+| --- | --- | --- | --- |
+| #81 (PR) | regime-A subtrees ran a whole KKT tree on one workgroup; flop limit `subtree_parallelism` | 0 | open PR |
+| #82 | KKT refactorization and solve are level-bound after #81: 54–59 launches per refactorization, 62–85 per solve | 5, 6 | open |
+| #75 | device LDLᵀ: serial pivot search, one workgroup per regime-B/C front, no vendor `sytrf` | 1, 2 | open, triaged |
+| #60 | regime-A follow-ups: CUDA timings (partly answered by experiment 0) and per-backend local-memory caps | 7 | open, triaged |
+| #25 | T25 performance pass (task) | 5, 6, 7 | open |
+
+Related accuracy issues that gate the K2 results: #71 (max\|L\| 1e14–1e16 on the K2 dumps, needs scaling) and #67 (pivot-pair matching), both experiment 3 / T21.
+
 ## TL;DR
 - **Current state:** SDS is a KernelAbstractions multifrontal solver. Symbolic analysis runs on the host (AMD or METIS ND, amalgamation, MA57-style pivot pairs). Numeric factorization runs in three regimes: fused subtree kernels, level-batched fused front kernels, and vendor potrf/trsm/syrk. SPD Cholesky and LDLᵀ with in-front Bunch–Kaufman and static perturbation both work on CUDA. LU, batching, matching and scaling, mixed precision and the non-CUDA backends do not exist yet. The LDLᵀ is deliberately slow: it reproduces the CPU reference pivot for pivot.
 - **cuDSS weak spots to exploit:** analysis (reordering) always runs on the host and is synchronous. It is the documented bottleneck in MadNLP/ExaModels: Pacaud, Shin, Montoison, Schanen and Anitescu (arXiv:2405.14236) report that "the analysis phase is four times slower for cuDSS compared to CHOLMOD". QOCO-GPU reports the same bottleneck. Symmetric-indefinite pivoting is diagonal-only within a supernode plus epsilon perturbation, with no LBLᵀ. The defaults (no matching) give relres 0.26–170 on the K2 dumps. Hybrid, MG and MGMN modes force synchronous phases.
@@ -59,6 +73,42 @@
 - **Compile latency.** The `Val`-specialized kernels make the first factorization cost minutes, which is bad for MadNLP users. Limit the specialization set (a few width classes), use PrecompileTools workloads on the CPU backend, and consider `@device_override`-free generic paths for rare sizes.
 - **Type stability and allocation.** Keep integer types uniform (Int32 on device). Keep the `pivot_stats` and inertia readback to a single host transfer per factorization. MadNLP needs inertia every iteration, and that one sync is unavoidable unless inertia correction moves to the device.
 
+## Experiment 0 results (2026-10-02, RTX 4080)
+
+Measured with `bench/profile_phases.jl` (CUPTI trace of one warm refactorization and solve per harness matrix; `bench/profile/phase_split.md`, the state before the fix in `bench/profile/phase_split_main_9d280d0.md`) and `bench/compare.jl` (`bench/comparison/comparison.md`).
+
+**The gap was not launch overhead.** GPU busy time (sum of kernel durations) equals wall time on every harness row, so the host never starves the device. The time goes into kernels that run on far too few thread blocks.
+
+**Root cause on every pglib KKT dump: one thread block.** The regime-A rule took a front into a fused subtree whenever its whole subtree's stack fit the local-memory budget, with no parallelism criterion. KKT trees of small fronts fit as a whole, so the entire factorization (3134 supernodes on case1354 condensed, 19.6k on case1354 K2) and the forward solve ran in one 128-thread workgroup on one of 76 SMs, at about 12 µs per front.
+
+**Fix (branch `perf/exp0-baseline-split`).** A new analysis tuning knob `subtree_parallelism` (default 4096): a regime-A subtree may do at most 1/4096 of the factorization flops. 4096 was best or near-best in a sweep over 128..4096 and a work floor never helped. Refactorization and solve, before and after:
+
+| matrix | refactorization | solve |
+| --- | --- | --- |
+| case118 condensed, Cholesky | 3.9 → 1.2 ms | 2.7 → 0.8 ms |
+| case118 K2, LDLᵀ | 13.8 → 2.0 ms | 7.7 → 0.9 ms |
+| case1354 condensed, Cholesky | 39.9 → 3.0 ms | 19.7 → 1.6 ms |
+| case1354 condensed, LDLᵀ | 56.3 → 9.6 ms | 22.6 → 1.2 ms |
+| case1354 K2, LDLᵀ | 146 → 7.2 ms | 78.7 → 2.6 ms |
+| SuiteSparse and Laplacian matrices | unchanged within noise (lap2d_300 36 → 31 ms) | unchanged |
+
+SDS/cuDSS geometric means over the harness (`comparison.md`), before → after:
+
+| feature | factorization | refactorization | solve |
+| --- | --- | --- | --- |
+| Cholesky, Float64 | 5.64× → 2.47× | 9.06× → 3.89× | 8.93× → 4.36× |
+| LDLᵀ, static pivoting | 20.4× → 6.37× | 29.8× → 9.90× | 13.4× → 3.96× |
+| LDLᵀ + 2 IR steps, K2 dumps | 23.4× → 3.18× | 36.0× → 5.45× | 33.3× → 7.57× |
+
+A side effect: the default solve (`deterministic_mode = 0`) now uses its atomic regime-B forward sweep on the small KKT and test matrices too, so two solves are no longer bitwise identical there. That was always the documented contract (bitwise reproducibility needs `deterministic_mode = 1`), but before the fix these matrices never reached the atomic path. MadNLP should set `deterministic_mode = 1` if it relies on repeatable solves.
+
+**What is left, by measurement:**
+1. **KKT refactor+solve is now level-bound.** 30–60 kernels per refactorization and 35–85 per solve, the longest kernel under 15% of busy time, regime-B fronts on one or two blocks per level. This is where the launch minimization of experiment 5 (level merging, graph capture) now pays, together with experiment 6 for the solve (issue #82). Remaining refactorization gap on case1354: 5× (Cholesky) and 9–14× (LDLᵀ).
+2. **LDLᵀ regime B/C is the single largest gap on everything else** (lap2d_300 43×, bcsstk17 59×, lap3d_40 250×, apache2 186× vs cuDSS): `front_ldlt_kernel` runs one front per workgroup on 1–2 blocks (issue #75). Experiments 1–2 unchanged in priority.
+3. **Large Cholesky** (lap3d_40 4.4×, apache2 2.4×, solve 11–13×): 1222/6608 kernels per refactorization, `extend_add` on 2 blocks. Level merging and wider extend-add grids.
+
+Revised order: experiments 1–2 (LDLᵀ kernels, #75) and 5–6 (launches and solve) now have the largest measured payoff; experiment 3 (matching and scaling) remains the accuracy blocker on K2 (relres unchanged by this fix).
+
 ## Recommendations: Prioritized Experiment Plan
 
 | # | Experiment | Payoff | Effort | Key measurement | Success criterion |

diff --git a/bench/README.md b/bench/README.md
@@ -15,6 +15,7 @@ Baselines every later milestone is measured against (PLAN.md §5 M0, §7). The
 | `dump_madnlp_kkt.jl` | dumps MadNLP K2 and condensed KKT matrices of pglib-opf cases |
 | `features.jl` | module `BenchFeatures`: the comparison features, one per planned capability in TASKS.md order (`cholesky_f64` … `mixed_precision`), with structure, element type, matrix selector and parameters; `task_status` reads the task's marker in TASKS.md |
 | `compare.jl` | cuDSS vs SparseDirectSolver.jl per feature × matrix with BenchmarkTools; one solver per run (`--solver=cudss` or `--solver=sds`), results merged into `bench/comparison/<solver>.csv` |
+| `profile_phases.jl` | where SDS spends refactorization and solve time on CUDA (PERFORMANCE.md, experiment 0): schedule per matrix (regime-A subtrees, B/C fronts, launch groups) and a CUPTI trace of one warm refactorization and solve (kernels, copies, busy time, longest kernel and its grid); writes `bench/profile/phase_split.{md,csv}`: `julia --project=bench bench/profile_phases.jl [--only=m1,m2] [--structures=SPD,S]` |
 | `compare_report.jl` | renders `bench/comparison/comparison.md` (overview + one table per feature) and `comparison.png` (SDS/cuDSS ratio per feature and phase); environment `bench/report/` (CairoMakie) |
 | `pivot_pairs.jl` | 2×2 pivot pairs of the `S` analysis (issue #66): `pivot_pairs` = `none`/`default`/`all` on the K2 dumps and the KKT generators, with nnz(L), zero/perturbed/2×2 pivots, max abs L and factor error of the CPU reference LDLᵀ; package environment: `julia --project=. bench/pivot_pairs.jl [--only=case118,...] [--generators=false]` |