michel2323 · Oct 9, 2026 · Oct 9, 2026 · Oct 9, 2026
diff --git a/PERFORMANCE.md b/PERFORMANCE.md
@@ -9,15 +9,15 @@ Performance issues carry the GitHub label `performance` ([list](https://github.c
 | issue | what | experiment | status |
 | --- | --- | --- | --- |
 | #81 (PR) | regime-A subtrees ran a whole KKT tree on one workgroup; flop limit `subtree_parallelism` | 0 | merged |
-| #82 | KKT refactorization and solve are level-bound after #81: 54–59 launches per refactorization, 62–85 per solve | 5, 6 | open, triaged; both phases prototyped in PR #107 (see Experiments 5–6 prototypes), in-src adoption is T25 |
+| #82 | KKT refactorization and solve are level-bound after #81: 54–59 launches per refactorization, 62–85 per solve | 5, 6 | open, triaged; both phases prototyped in PR #107 (see Experiments 5–6 prototypes), in-src adoption is T24 (solve) and T25/T26 (factorization) |
 | #107 (PR) | prototypes in `bench/`: partitioned-inverse fused solve, split/fused factorization, ordering study, `subtree_max_fronts`, MadNLP end-to-end on CUDA and AMD, AMD CI leg | 0, 4, 5, 6, 7 | merged (bench only; one `Options` knob in `src/`) |
 | #108 | ordering chooser picks AMD on large KKT systems: the cost model scores column-etree depth, not supernodal schedule depth (65 vs 29 levels, 36% on the solve) | 5, 6, 10 | open (found-by-agent) |
 | #109 | update-stack placement assumes level-synchronous execution: out-of-level-order numeric phases need private contribution blocks (4× the stack) or DAG lifetimes | 5 | open (found-by-agent) |
 | #110 | regime-C Cholesky reportedly host-synchronizes per front through its `info` check; to verify | 5 | open (found-by-agent) |
 | #75 | device LDLᵀ: serial pivot search, one workgroup per regime-B/C front, no vendor `sytrf` | 1, 2 | open, triaged; steps 1–3 in PRs from `perf/exp1-pivot-search`, `perf/exp1-regime-b-local`, `perf/exp2-regime-c-blas`; criteria partly met (see Experiments 1–2 results) |
 | #86 | device LDLᵀ differs from `ref_ldlt!` on the K2 dumps (pivot sequence, `nperturbed`; rounding with max\|L\| 1e14–1e16), pre-existing | 1, 3 | open (found-by-agent) |
 | #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 |
+| #22, #24, #25, #26 | T22 ordering chooser, T24 fused solve, T25 split factorization, T26 segmented fused factorization (tasks; the former T25 performance pass) | 5, 6, 7, 10 | 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.
 
@@ -62,7 +62,7 @@ Related accuracy issues that gate the K2 results: #71 (max\|L\| 1e14–1e16 on t
 - **Evidence on weaknesses.** Pacaud and Shin (arXiv:2403.15913, CDC 2024) find that "cuDSS spends most of its time in the symbolic analysis" and that "ma27 is approximately twice as fast during the pre-processing". On their N = 50,000 instance, initialization took 15.4 s for ma27 against 27.9 s for Lifted-KKT and 29.7 s for HyKKT. On an A30 GPU, that same N = 50,000 instance (3,350,067 rows) took 20.187 s for SYM, 0.432 s for FAC and 0.165 s for SOLVE, i.e. "less than 0.5 seconds to recompute the factorization of the largest instance". Over the full solve, "Lifted-KKT and HyKKT solve the problem respectively 26x and 18x faster than HSL ma27 on the largest instance (N = 50, 000)" (totals of 33.8 s and 34.5 s against 125.5 s). QOCO-GPU (arXiv:2603.29197) reports that "for larger problems, up to 75% of qoco-gpu's runtime is spent in the setup phase, where the dominant cost is the reordering step in cuDSS's analysis phase."
 
 ### 3. Techniques most likely to matter, ranked for repeated IPM KKT solves
-1. **Parallel indefinite pivoting.** Use APTP (Duff–Hogg–Lopez; SSIDS): factor a block optimistically with BLAS-3, check the threshold a posteriori, and fail or delay only the bad columns. It keeps TPP-level robustness and was designed for GPU/multicore. [researchgate](https://www.researchgate.net/publication/339832172_A_New_Sparse_LDLT_Solver_Using_A_Posteriori_Threshold_Pivoting) [researchgate](https://www.researchgate.net/publication/335908712_Exploring_Benefits_of_Linear_Solver_Parallelism_on_Modern_Nonlinear_Optimization_Applications) For SDS this means replacing the serial `_choose_pivot` with a warp-parallel argmax BK search over the panel plus an a-posteriori check. Pivots that fail inside the front fall back to perturbation now and to delayed pivots in T27. This is the biggest lever, because it is what lets regime C use vendor trsm/gemm (or syrk on L·D) for most of the flops.
+1. **Parallel indefinite pivoting.** Use APTP (Duff–Hogg–Lopez; SSIDS): factor a block optimistically with BLAS-3, check the threshold a posteriori, and fail or delay only the bad columns. It keeps TPP-level robustness and was designed for GPU/multicore. [researchgate](https://www.researchgate.net/publication/339832172_A_New_Sparse_LDLT_Solver_Using_A_Posteriori_Threshold_Pivoting) [researchgate](https://www.researchgate.net/publication/335908712_Exploring_Benefits_of_Linear_Solver_Parallelism_on_Modern_Nonlinear_Optimization_Applications) For SDS this means replacing the serial `_choose_pivot` with a warp-parallel argmax BK search over the panel plus an a-posteriori check. Pivots that fail inside the front fall back to perturbation now and to delayed pivots in T28. This is the biggest lever, because it is what lets regime C use vendor trsm/gemm (or syrk on L·D) for most of the flops.
 2. **Matching and scaling for K2.** Use MC64-style symmetric scaling (max diag product, as in cuDSS algo5) or cheaper equilibration (Ruiz) computed once per pattern on the host, with the scaling recomputed on the device each iteration. cuDSS needs algo5 to be accurate on K2, so SDS must have it for parity, and per-iteration device-side equilibration is a differentiator. [github](https://github.com/exanauts/SparseDirectSolver.jl/pull/44) Recomputing the matching itself per iteration is a host cost, so cache the matching and refresh it only when IR stalls.
 3. **Device-resident, launch-minimal refactor+solve.** At pglib scale (n ≈ 1e4–1e6), elimination-tree depth times launches per level dominates. Merge levels: a persistent kernel walks a level window with grid-wide or atomic-counter sync, i.e. sync-free dependency counters on supernodes (Liu et al. 2016 style). [ssslab](https://www.ssslab.cn/assets/papers/2016-liu-sptrsv.pdf) Capture the whole refactorize+solve+IR sequence into one CUDA graph. cuDSS itself cannot capture analysis, and hybrid execute (its small-matrix mode) forces syncs, which is exactly SDS's opening. [nvidia](https://docs.nvidia.com/cuda/archive/13.0.1/cudss/general.html)
 4. **Triangular solve.** Use supernodal level-set SpTRSV with batched TRSV/GEMV at the bottom levels and streams at the top. Invert the diagonal blocks so TRSV becomes GEMV (the partitioned inverse of Alvarado–Pothen–Schreiber). Tacho uses exactly this split. Yamazaki, Rajamanickam and Ellingwood (ICPP '20), whose paper also covers "an algorithmic variant called the partitioned inverse", report that their Kokkos supernodal solver "can be 12.4× or 19.5× faster than the vendor optimized implementation in NVIDIA's CuSPARSE library" on V100/P100. In IPMs the solve is called 2–6× per factorization (IR, inertia correction retries, second-order correction), so solve latency matters as much as factorization.
@@ -238,7 +238,7 @@ A prototype campaign by a human collaborator with their own agent, recorded in `
 | refactorization + solve | 185 ms | ~35 ms | 28.1 ms |
 | solve, Radeon VII (gfx906, no vendor solver exists there) | 20.5 ms | **8.05 ms** | — |
 
-How the solve got there (`bench/solve_proto_*.jl`, `bench/solve_nd.jl`, `bench/final_sweep.jl`): per-front `L₁₁` inverses (the partitioned-inverse idea of plan row 6) turn every TRSV into a GEMV; one fused dependency-counter kernel per direction covers every front of width ≤ 256 and the vendor path only the root; the whole solve is two custom kernels, the root and two permutes. The ordering is a precondition: under the default METIS ND the schedule has 65 levels and the prototype reads 6.35 ms; under cuDSS's reordering (imported with `bench/get_cudss_perm.jl`) or SDS's own ND (`reordering_alg = "algo4"`, `bench/order_search.jl`) it has 29–30 levels and reads 3.5–3.6 ms. The default chooser picks AMD here because its cost model scores column-etree depth, which does not predict the supernodal schedule depth (#108). An amalgamation sweep only regresses the solve (8.7–28 ms), and on lap3d_40 and apache2 the same solve regresses, so T25 has to pick the strategy per schedule.
+How the solve got there (`bench/solve_proto_*.jl`, `bench/solve_nd.jl`, `bench/final_sweep.jl`): per-front `L₁₁` inverses (the partitioned-inverse idea of plan row 6) turn every TRSV into a GEMV; one fused dependency-counter kernel per direction covers every front of width ≤ 256 and the vendor path only the root; the whole solve is two custom kernels, the root and two permutes. The ordering is a precondition: under the default METIS ND the schedule has 65 levels and the prototype reads 6.35 ms; under cuDSS's reordering (imported with `bench/get_cudss_perm.jl`) or SDS's own ND (`reordering_alg = "algo4"`, `bench/order_search.jl`) it has 29–30 levels and reads 3.5–3.6 ms. The default chooser picks AMD here because its cost model scores column-etree depth, which does not predict the supernodal schedule depth (#108). An amalgamation sweep only regresses the solve (8.7–28 ms), and on lap3d_40 and apache2 the same solve regresses, so T24 has to pick the strategy per schedule.
 
 How the factorization got there (`bench/fact_split.jl`, `bench/fact_fused.jl`), each step on top of the previous one:
 
@@ -278,7 +278,7 @@ The loop arithmetic closes (193 factorizations × ~36 ms + 202 solves × ~4 ms +
 
 **Criteria.** Experiment 6 (solve ≤ cuDSS): met in the prototype on this matrix, not on all harness matrices, and not in `src/`. Experiment 5 (≥ 3× fewer launches, refactor+solve faster than cuDSS for n ≤ 1e5): launches reduced far beyond 3× in the prototype, the time criterion not met (1.31× on this matrix, and n is 6.7e5); graph replay, the lever the row names, is neutral. Experiment 4 (iterations within ±2, same objective): met for the condensed system; K2 with LDLᵀ not run.
 
-**Revised priority.** The ordering comes first: it is a third of the remaining solve gap, a precondition for the fused solve, and the fix (#108: score the supernodal schedule depth, or prefer ND on GPU backends) is host code. Then experiment 6 in `src/`, then experiment 5 with #109 solved by placement, then the factorization buckets above. The T25 owner note in `TASKS.md` orders the deliverables.
+**Revised priority.** The ordering comes first: it is a third of the remaining solve gap, a precondition for the fused solve, and the fix (#108: score the supernodal schedule depth, or prefer ND on GPU backends) is host code. Then experiment 6 in `src/`, then experiment 5 with #109 solved by placement, then the factorization buckets above. `TASKS.md` splits them into T22 (ordering), T24 (solve), T25 (split factorization) and T26 (segmented fused factorization, #109, #110).
 
 ## Recommendations: Prioritized Experiment Plan
 
@@ -289,13 +289,13 @@ The loop arithmetic closes (193 factorizations × ~36 ms + 202 solves × ~4 ms +
 | 2 | **Regime-C LDLᵀ via vendor BLAS**: KA pivoting of an F₁₁ panel (nb = 32–64), then cuBLAS trsm + gemm (L·D·Lᵀ update) for F₂₁/F₂₂; multi-workgroup per root front | Very high | Med | achieved TFLOP/s on fronts > 256; share of factor time in root fronts | ≥ 50% of cuBLAS DGEMM peak on fronts > 512; "C only" faster than reference on CPU and ≥ 5× on GPU |
 | 3 | **Matching + scaling (T21)**: MC64 max-product symmetric scaling cached per pattern; device Ruiz equilibration per iteration as a cheap variant | Very high (accuracy) | Med | max\|L\|, factor error, nperturbed, relres after 0/2/5 IR on case1354 K2 | relres ≤ cuDSS algo5 (≤ 4e-7) with ≤ 5 IR steps; nperturbed ≤ cuDSS's 0–14 [github](https://github.com/exanauts/SparseDirectSolver.jl/pull/73) |
 | 4 | **MadNLP end-to-end** (#28) on pglib via ExaModelsPower: K2 with SDS-LDLᵀ vs cuDSS-LDLᵀ, plus condensed/Lifted with SDS-Cholesky vs cuDSS-Cholesky. *Condensed part done in PR #107 (`bench/e2e/`, identical iterations on CUDA and AMD); K2 with LDLᵀ and the MadNLP task remain* | Essential | Low–Med | IPM iterations, inertia-correction count, time per iteration split (factor/solve/other) | iterations within ±2 of cuDSS, same objective to 1e-6 (the #28 criterion); [github](https://github.com/exanauts/SparseDirectSolver.jl/issues/28) fewer inertia corrections thanks to 2×2 pivots |
-| 5 | **Launch minimization**: level merging into persistent kernels with atomic dependency counters (regimes A/B); one graph for refactor+solve+IR, replayed via a cached exec (not `@captured`). *Prototyped in PR #107: segmented counters and split kernels reach 1.31× cuDSS; graph replay is neutral; needs #109; in-src adoption is T25* | High at small/medium n | Med | launches/refactor, CPU-side time, GPU idle gaps (Nsight timeline) | ≥ 3× fewer launches; refactor+solve faster than cuDSS for n ≤ 1e5 |
-| 6 | **Solve path**: partitioned-inverse diagonal blocks (`solve_alg="algo1"`), batched GEMV for small supernodes, sync-free forward sweep; fused permute/scale/IR residual kernels. *Prototyped in PR #107: criterion met on the 78k-bus condensed KKT (3.5 vs 3.8 ms), needs the ND ordering (#108); in-src adoption is T25* | High | Med | solve latency (1 RHS and 2–6 RHS), backward error | solve ≤ cuDSS solve on all harness matrices; backward error unchanged within 10× |
+| 5 | **Launch minimization**: level merging into persistent kernels with atomic dependency counters (regimes A/B); one graph for refactor+solve+IR, replayed via a cached exec (not `@captured`). *Prototyped in PR #107: segmented counters and split kernels reach 1.31× cuDSS; graph replay is neutral; needs #109; in-src adoption is T25 and T26* | High at small/medium n | Med | launches/refactor, CPU-side time, GPU idle gaps (Nsight timeline) | ≥ 3× fewer launches; refactor+solve faster than cuDSS for n ≤ 1e5 |
+| 6 | **Solve path**: partitioned-inverse diagonal blocks (`solve_alg="algo1"`), batched GEMV for small supernodes, sync-free forward sweep; fused permute/scale/IR residual kernels. *Prototyped in PR #107: criterion met on the 78k-bus condensed KKT (3.5 vs 3.8 ms), needs the ND ordering (#108, T22); in-src adoption is T24* | High | Med | solve latency (1 RHS and 2–6 RHS), backward error | solve ≤ cuDSS solve on all harness matrices; backward error unchanged within 10× |
 | 7 | **Amalgamation and bin tuning**: sweep (max_width, zero_fraction, min_width), regime thresholds, local-memory budgets; size-sorted regime-A subtrees | Med | Low | flops vs time Pareto; stack memory/nnz(L) | 10–30% factor-time gain without > 15% nnz(L) growth |
 | 8 | **Mixed precision**: FP32 factor + FGMRES-IR (T18 first); enable by IPM phase or by condensed SPD | Med | Med | IR iterations, time to 1e-10 relres, failures near convergence | ≥ 1.3× factor speedup with no change in IPM iteration count |
 | 9 | **Uniform batching** (T17) for multi-scenario/MPC: shared symbolic, batched fronts across instances | Med–High for batch workloads | Med | throughput vs cuDSS uniform batch | ≥ cuDSS batch throughput |
 | 10 | **Analysis caching / GPU symbolic**: serialize analysis per pattern; multithreaded AMD/ND; device etree/colcounts | Med (MPC/online) | Med–High | analysis time vs cuDSS (default and MT) | analysis ≤ cuDSS on the harness; zero cost when reusing a cached pattern |
-| 11 | **Delayed pivots / APTP fallback** (T27) and portability (AMDGPU on MI250/MI300) | Med | High | relres on the hard delayed-pivot case; AMD timings | T27 tests; portability as a unique selling point (cuDSS is NVIDIA-only) |
+| 11 | **Delayed pivots / APTP fallback** (T28) and portability (AMDGPU on MI250/MI300) | Med | High | relres on the hard delayed-pivot case; AMD timings | T28 tests; portability as a unique selling point (cuDSS is NVIDIA-only) |
 
 **Benchmark matrices.** MadNLP KKT dumps (K2 and condensed, early, middle and converged iterates) from pglib-opf via ExaModelsPower: case118, case1354_pegase, case2869_pegase, case9241_pegase, case13659_pegase, and ACTIVSg25k/70k if available. Add COPS and optimal-control instances from the MadNLP/HybridKKT benchmarks. From SuiteSparse: TSOPF_RS_*, rajat21, and an SPD set (e.g. thermal2, G3_circuit, audikw_1, Serena) for Cholesky parity.