From 7af9d99d61f1aeb5e4015b79faef7c5f1971d1c4 Mon Sep 17 00:00:00 2001 From: Michel Schanen Date: Fri, 9 Oct 2026 13:16:33 +0000 Subject: [PATCH 1/2] Docs: split T25 into T22/T24-T26 after PR #107, renumber the post-v1 tasks T22-T30 [skip ci] The former T25 (performance pass) carried six deliverables from the #107 prototypes, several sessions of work under the one-task-per-session rule. TASKS.md now has: T22 ordering chooser scoring the supernodal schedule depth (#108, host only, precondition for the fused solve); T23 backends (unchanged); T24 partitioned-inverse fused solve, solve_alg = "algo1" (#82 solve half); T25 split factorization (tiled SYRK, chunked TRSM, split regime-C fronts, longest-first regime A); T26 segmented fused factorization with dependency counters, with #109 and #110 inside it plus the #75 remainder and the #96 re-measurement; T27 hybrid memory (was T26); T28 robustness extras (was T27, with the FP32 finding of #107); T29 non-uniform batch (was T22); T30 ND partition-tree export and ordering cache (was T24); External unchanged. Each new task names the bench prototype it adopts, the measured numbers, the negative results not to re-run, the tests and the Report asks. PLAN.md, PERFORMANCE.md, bench/README.md: task-id references follow the renumbering; the stale "empty repository" status line and the M1 row's ordering sentence are updated. STATE.md: section 8 records the restructuring and the GitHub steps (issue renames, triage labels, chain resumes with T22). Refs #96, #108, #109, #110. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01Wzx8WtMEXnAfyzyujdySFh --- PERFORMANCE.md | 16 +-- PLAN.md | 22 ++-- STATE.md | 20 ++++ TASKS.md | 302 ++++++++++++++++++++++++++++++++++++------------ bench/README.md | 2 +- 5 files changed, 266 insertions(+), 96 deletions(-) diff --git a/PERFORMANCE.md b/PERFORMANCE.md index 481c9f9..8dd47fc 100644 --- a/PERFORMANCE.md +++ b/PERFORMANCE.md @@ -9,7 +9,7 @@ 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) | @@ -17,7 +17,7 @@ Performance issues carry the GitHub label `performance` ([list](https://github.c | #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. diff --git a/PLAN.md b/PLAN.md index 1a298bd..a4012d6 100644 --- a/PLAN.md +++ b/PLAN.md @@ -65,7 +65,7 @@ IPM KKT systems, 2026-10-01). What it changed in this plan: sync-free solves and graph capture live only in the CUDA/ROCm extensions. * **Benchmark harness and cuDSS baseline measurements** move into M0. -Status: empty repository apart from `PLAN.md` and `RESEARCH.md`. +Status: see §5 (M0–M9 done, M11 in progress through T22–T26). --- @@ -123,7 +123,7 @@ solve phases raise `NotSupportedError` (T20, §3.6). | `use_cuda_register_memory` | reinterpret (M12) | Pinned host memory through the backend extension. | | `hybrid_execute_mode` | defer (M12) | Small levels on the KA CPU backend. | | `host_nthreads` | reinterpret | Julia threads for the host symbolic phase / CPU backend. | -| `nd_nlevels`, `nd_ubfactor` | port (T05) | ND is `METIS_NodeND` through CliqueTrees; `nd_ubfactor` is passed through, `nd_nlevels` is read as cuDSS documents it, a *minimum* number of dissection levels, which NodeND's full recursion meets (a level-capped ND was 4–13× slower with more fill). It becomes meaningful with the partition-tree export (T24). | +| `nd_nlevels`, `nd_ubfactor` | port (T05) | ND is `METIS_NodeND` through CliqueTrees; `nd_ubfactor` is passed through, `nd_nlevels` is read as cuDSS documents it, a *minimum* number of dissection levels, which NodeND's full recursion meets (a level-capped ND was 4–13× slower with more fill). It becomes meaningful with the partition-tree export (T30). | | `ubatch_size`, `ubatch_index` | port (M6) | Uniform batch. | | `use_superpanels` | reinterpret | Supernode amalgamation on/off. | | `device_count`, `device_indices` | not planned | Raise "not supported". | @@ -210,7 +210,7 @@ user API: LinearAlgebra generics + handle-style layer (CUDSS.jl strings, C-wra │ regime A fused subtree-per-workgroup kernels (KA) │ regime B fused per-front size-binned level kernels (KA; vendor batched optional) │ regime C Cholesky: vendor potrf + trsm + syrk/herk (extensions; KA fallback) - │ LDLᵀ: KA blocked in-front pivoting + vendor trsm/gemm; LU: KA kernel (BLAS-3 in T25) + │ LDLᵀ: KA blocked in-front pivoting + vendor trsm/gemm; LU: KA kernel (BLAS-3 in T26) │ + KA kernels for assembly, extend-add, gather/scatter, permutation ├─ solve (device): permute/scale → subtree/level fwd → diag → bwd → unpermute → IR/FGMRES └─ extensions: CUDA, AMDGPU, oneAPI, Metal (dense + sparse adapters + graph capture @@ -253,7 +253,7 @@ tree levels, so the numeric phase is organized by front size: step. LDLᵀ: the fully-summed block is factored in blocks of 32 columns by a KA kernel that reproduces the reference pivot sequence, with one vendor GEMM per block for the trailing columns and the contribution block (#89). LU: the - fused KA kernel, one workgroup per front; BLAS-3 for `L₂₁`/`U₁₂` is a T25 + fused KA kernel, one workgroup per front; BLAS-3 for `L₂₁`/`U₁₂` is a T26 item. Vendor `sytrf`/`getrf` are not used (§3.3). This mirrors cuDSS keeping a distinct algorithm for very sparse factors, and @@ -268,7 +268,7 @@ runtime decisions. 2. **Ordering** via CliqueTrees.jl: `permutation(graph; alg)` with `AMD()` (AMD.jl, a hard dependency) /`MMD()` or `METIS_NodeND` (Metis extension, `nd_ubfactor` passed through, `nd_nlevels` a minimum); user permutation; - natural; ND-tree export/import in the cuDSS encoding (T24). **Automatic + natural; ND-tree export/import in the cuDSS encoding (T30). **Automatic choice** computes both AMD and ND candidates (cheap for KKT sizes) and picks by a cost model `flops × (1 + nlevels/n)` on the column etree: ND gives bushier, shallower trees; AMD often gives lower fill on power grids (and @@ -276,7 +276,7 @@ runtime decisions. large KKT systems (issue #108, PR #107): at n ≈ 7e5 it weighs depth at 0.2% and scores column-etree depth, which does not predict the supernodal schedule depth that governs GPU time, so it picks AMD where ND halves the - schedule depth (65 → 29) and is 36% faster on the solve; until T25 scores + schedule depth (65 → 29) and is 36% faster on the solve; until T22 scores the schedule depth, GPU users should set `reordering_alg = "algo4"`. Schur mode constrains the ordering (§3.6). Matching (T21) composes a column permutation for `"G"` and @@ -366,7 +366,7 @@ test contract of T15/T19 (the pivot search is cooperative across the workgroup, the tie-breaking is the reference's). LDLᵀ root fronts run a blocked KA panel factorization (32 columns per block, lazily updated) with one vendor GEMM per block for the trailing columns and the contribution block (#89); LU root -fronts run the fused KA kernel, BLAS-3 pending (T25, #75). On the pglib K2 +fronts run the fused KA kernel, BLAS-3 pending (T26, #75). On the pglib K2 dumps the device and the reference can still differ by rounding (max|L| up to 1e16 unscaled), which flips threshold and tie decisions: there the contract is equal inertia and `nperturbed` within a tolerance (#86). The whole phase is a @@ -403,7 +403,7 @@ GEMV, and one fused dependency-counter kernel per direction over every front of width ≤ 256 solves the 78k-bus condensed KKT in 3.5 ms against 3.8 ms for cuDSS and 18 ms for the level-batched sweeps, on CUDA and (8 ms) on a Radeon VII; the shallow ND ordering (#108) is a precondition. On lap3d_40 and apache2 -the same approach regresses, so T25 picks the strategy per schedule. +the same approach regresses, so T24 picks the strategy per schedule. Refinement: plain IR with a KA CSR SpMV residual (gather, no atomics; the workspace is allocated by the first refining solve, so the handle layer with @@ -448,7 +448,7 @@ defect above 128×128 on macOS 27): Extensions: `…CUDAExt` (weak dependencies `CUDACore`, `cuSPARSE`, `cuBLAS`, `cuSOLVER` of CUDA.jl 6: dense bindings, sparse adapters; pinned memory, graph -capture and the sync-free solve are T25/T26), `…AMDGPUExt`, `…OneAPIExt`, +capture and the sync-free solve are T24/T26), `…AMDGPUExt`, `…OneAPIExt`, `…MetalExt` (T23), `…MetisExt` (ND), `…KrylovExt` (FGMRES-IR). Core depends on KernelAbstractions 0.9, GPUArrays(Core), Adapt, Atomix, LinearAlgebra, SparseArrays, CliqueTrees, AMD, Metis (ordering only through the extension). @@ -615,7 +615,7 @@ front sizes, extend-add maps, layout, schedule. The numeric phase and the solve are pure "values in, factors/solution out" kernel sequences with no allocation and no host synchronization, so the CUDA and ROCm extensions can capture the refactorize+solve sequence in a graph and replay it per IPM iteration (graph -capture itself is T25, #82). This holds for the dense implementations `:auto`, +capture itself is T26, #82). This holds for the dense implementations `:auto`, `:vendor` and `:ka` (the latter still allocates until T23, #53); `:generic` (host LinearAlgebra) is the reference path and may allocate. Exceptions by design: the `ir_tol > 0` early-exit test and every FGMRES iteration read on the @@ -680,7 +680,7 @@ rough: S ≈ days, M ≈ 1–2 weeks, L ≈ 3–4 weeks, XL ≈ more. | # | Milestone | Size | Definition of done | | --- | --- | --- | --- | | M0 | Scaffolding, audit, baselines | M | Package, four backend extensions with CSR adapters, CI matrix (GitHub Actions CPU; Buildkite juliagpu queue for CUDA/AMDGPU/oneAPI/Metal; self-hosted `kkt`), `Options` with the full name tables, error types, Aqua. **Capability audit** script filling the §2.6 table per backend and eltype (kept as a test). **Benchmark harness**: NREL opf_matrices, pglib-opf KKT and condensed matrices dumped from MadNLP, a CUTEst subset; cuDSS analysis/factorization/solve times, `flops`, supernode statistics and MA57 refinement counts recorded as the baseline every later milestone is measured against. | -| M1 | Symbolic engine | L | Pattern, orderings with the AMD/ND cost model, etree, column counts, GPU-tuned amalgamation, subtree partition, size bins, level lists, static layout, device maps, `memory_estimates`, `flops`, `lu_nnz`, `nsuperpanels`, ND-tree I/O, `max_lu_nnz`. Validated against CHOLMOD's nnz(L) and etree; supernode/level statistics compared with the cuDSS baseline on the harness matrices. A CPU reference numeric factorization (plain Julia, same layout) for testing. | +| M1 | Symbolic engine | L | Pattern, orderings with the AMD/ND cost model (scored on the supernodal schedule depth from T22, issue #108), etree, column counts, GPU-tuned amalgamation, subtree partition, size bins, level lists, static layout, device maps, `memory_estimates`, `flops`, `lu_nnz`, `nsuperpanels`, ND-tree I/O, `max_lu_nnz`. Validated against CHOLMOD's nnz(L) and etree; supernode/level statistics compared with the cuDSS baseline on the harness matrices. A CPU reference numeric factorization (plain Julia, same layout) for testing. | | M2 | Numeric kernels | L | Regime A fused subtree kernels, regime B fused per-front kernels (Cholesky first, LDLᵀ/LU hooks), regime C vendor bindings in all four extensions with KA tiled fallbacks, assembly/extend-add kernels. Unit tests against the CPU reference on every backend; per-bin micro-benchmarks deciding where vendor batched calls replace regime B. | | M3 | SPD/HPD end-to-end (**v0.1**) | M | All phases, single and multi RHS, `info`, `cholesky`/`cholesky!`/`ldiv!`/`\`, `Hermitian` wrapper, async flag; subtree/level solve sweeps. Ported `test_cudss.jl` subsets pass on all backends. **Target**: within 1.5× of cuDSS Cholesky refactorization+solve on condensed pglib-opf systems on CUDA, running on AMD and Intel. | | M4 | Symmetric indefinite LDLᵀ/LDLᴴ (**MadNLP-ready**) | L | In-front Bunch–Kaufman, `pivot_threshold`, `pivot_epsilon(_alg)`, `pivot_sign`, `inertia`, `pivot_stats`, `npivots`, `diag`, `solve_diag`, `ldlt`/`ldlt!`; vendor `sytrf` with post-check on root fronts. MadNLPGPU gets a `SparseDirectSolver` option; validated on OPF/ExaModels KKTs (MadNLP K2/K2r, MadIPM, MadNCL settings) against the cuDSS path, including refinement counts against MA27/MA57. | diff --git a/STATE.md b/STATE.md index 772051a..aed8870 100644 --- a/STATE.md +++ b/STATE.md @@ -156,3 +156,23 @@ Out-of-chain PRs merged since this review, all by hand and on `main`: New issues from #107, all `found-by-agent` + `performance`, none `triaged` yet (they hold nothing while #22 carries `on-hold`, but they gate the chain once it resumes): #108 (ordering chooser picks AMD on large KKT: score the supernodal schedule depth), #109 (update-stack placement assumes level-synchronous execution), #110 (regime-C Cholesky reportedly host-synchronizes per front; verify). Owner notes for them are under T23, T24, T25 and the External task in `TASKS.md`. Still open from §5–§6: #96 and #108–#110 untriaged; #60 item 1 unmeasured; `CLAUDE_GH_PAT` unset; the two `lts` required checks in the ruleset (unless already removed); the External MadNLP task before T22–T24 is still the recommendation, now with `bench/e2e/MadNLPSDS.jl`, `reordering_alg = "algo4"` and the phase-coupled knobs as known inputs; T25 absorbs #82, #75's remainder and the #107 designs. + +## 8. Addendum (2026-10-09): task restructuring after PR #107 + +The former T25 held six deliverables from the #107 prototypes; it is split and +the post-v1 tasks are renumbered so the pipeline's `TNN → TNN+1` chaining +still holds. New order in `TASKS.md`: **T22** ordering chooser scoring the +supernodal schedule depth (#108, host only, precondition for the fused solve); +**T23** AMDGPU/oneAPI/Metal extensions (unchanged); **T24** partitioned-inverse +fused solve `solve_alg = "algo1"` (#82, solve half); **T25** split +factorization (tiled SYRK, chunked TRSM, split regime-C fronts, longest-first +regime A); **T26** segmented fused factorization with dependency counters, +with #109 (stack lifetimes) and #110 (verify the regime-C host sync) inside +it, plus the #75 remainder and the #96 re-measurement; **T27** hybrid memory +(was T26); **T28** robustness extras (was T27); **T29** non-uniform batch (was +T22); **T30** ND partition-tree export and ordering cache (was T24); External +unchanged and last (owner decision: integrate against finished kernels; the +§5 item 5 recommendation is withdrawn). GitHub: issues #22–#27 renamed to the +new titles, #22 loses `on-hold`, new task issues for T28–T30; #96 and +#108–#110 labelled `triaged` (#108 → T22, #109/#110 → T26). The chain resumes +with T22. diff --git a/TASKS.md b/TASKS.md index f80aa06..982a18a 100644 --- a/TASKS.md +++ b/TASKS.md @@ -2833,9 +2833,41 @@ inertia with matching enabled equals the eigenvalue count (the cuDSS defect). matching cycles give the 2×2 pairs; jobs 1–4 keep the default pairs". - PLAN §5 M13: keep a posteriori pivoting after the MadNLP integration (T21 measurement above). -### T22 — Non-uniform batch (`BatchedDirectSolver`) `[ ]` - -Block-diagonal packing, forest schedule; `test_nonuniform_batch_cudss.jl` ported. +### T22 — Ordering chooser: supernodal schedule depth (issue #108) `[ ]` + +**Why this is first (PR #107)**: the automatic ordering (`compute_ordering`, +`src/symbolic/ordering.jl`, `reordering_alg = "default"`) picks AMD on large +KKT systems where ND halves the supernodal schedule depth (65 → 29 on the +78k-bus pglib condensed KKT, 36% on the solve, 11–17% on the factorization, +GV100). Two causes: `ordering_cost(flops, nlevels, n) = flops × (1 + +nlevels/n)` weighs depth at ~0.2% for n ≈ 7e5, so the chooser is a flop +contest; and `nlevels` is the column-etree depth, which does not predict the +schedule depth (cuDSS's permutation: 1279 column levels → 30 schedule levels; +METIS: 1216 → 29; AMD: competitive column metrics → 65). The fused solve of +T24 needs the shallow ordering to pay off, so the chooser is fixed before it. + +Deliverables, host code only: (1) `evaluate_ordering` scores each candidate +by the depth of its supernodal schedule — fundamental supernode partition of +the candidate's etree and `tree_levels(snparent)` (the T06 Report already +suggested it; the amalgamation step is not needed for the ranking) — with a +cost model whose depth term matters at n ≈ 1e6; the task picks the form and +documents it in PLAN-style in the `ordering_cost` docstring. (2) If no model +is robust across the harness, the fallback is a backend-dependent default +(ND on GPU backends, the current model on the CPU backend), stated in the +options table. (3) `Ordering.stats` reports the schedule depth of every +candidate and the chosen one. (4) A METIS knob sweep is not part of this task +(seed, ufactor and nseps move the depth 28–41 within noise in +`bench/order_search.jl`; plain ND is enough). + +Tests: on the KKT generators of `test/matrices.jl` the chooser selects the +candidate with the smaller schedule depth, asserted through `Ordering.stats`; +explicit `reordering_alg = "algo1"`/`"algo4"` reproduce the current orderings +bit for bit; `test_symbolic_*`, `test_options` and every numeric suite +unchanged. Report: table of schedule depth, column-etree depth, nnz(L) and +flops for AMD, ND and the chosen ordering on every T04 harness matrix, and on +the 78k-bus dump if the owner provides it (`bench/order_search.jl`, +`bench/dump_madnlp_kkt.jl`). Closes #108. The PLAN §2.3 step 2 sentence +"until T25 scores the schedule depth" is updated by the owner after the merge. ### T23 — AMDGPU, oneAPI and Metal extensions `[!]` @@ -2958,86 +2990,203 @@ needed to load them on Linux); Metal only by review. Mark both as untested. - PLAN §2.3 step 5 / §2.7 hazard row: regime-A ladder 8–64 KiB, capped per backend by `max_local_bytes`. - TASKS.md: a task for the oneAPI/Metal remainder of T23 with the #53 owner note. -### T24 — ND partition-tree export/import and ordering cache `[ ]` - -`nd_partition_tree`/`user_nd_partition_tree` in the cuDSS binary-tree -encoding; test: importing the exported tree reproduces the same supernode -partition and `nnz_L`. - -**Owner note (PR #107)**: the ordering study found that SDS's own ND -(`reordering_alg = "algo4"`) matches the imported cuDSS ordering on the -78k-bus KKT (schedule depth 29 vs 30), so the tree import is a compatibility -feature, not a performance one. The ordering cache (analysis under a stored -permutation ran in 3.6 s against 12 s for a fresh METIS analysis) is the part -MadNLP benefits from. - -### T25 — Performance pass `[ ]` - -**Owner note (PR #107, issues #82, #108, #109, #110)**: experiments 5 and 6 -exist as bench-level prototypes measured on a Quadro GV100 against the 78k-bus -pglib condensed KKT (`bench/solve_proto_*.jl`, `bench/fact_split.jl`, -`bench/fact_fused.jl`, `bench/e2e/SDSProto.jl`): solve 18 → 3.5 ms (cuDSS -3.8), refactorization 168 → 31.9 ms (cuDSS 24.3), per IPM iteration 185 → 35 ms -(cuDSS 28), all cuDSS-free. The deliverables of this task follow from them, in -this order: (1) the ordering chooser scores supernodal schedule depth, or GPU -backends default to ND (#108; worth 2× schedule depth and 36% on the solve); -(2) `solve_alg = "algo1"` as the partitioned-inverse solve with per-front -`L₁₁` inverses and one fused dependency-counter kernel per direction for fronts -of width ≤ 256, chosen per schedule (it regresses on lap3d_40 and apache2); -(3) the split factorization: tiled 32×32 SYRK and chunked TRSM as separate -many-block kernels for fronts with `m > 64`, the regime-C fronts narrower than -64 through the same split kernels batched per level with the regular assembly -kernels (vendor `potrf` only above width 64), longest-first regime-A launches; -(4) the segmented fused factorization with dependency counters, segments split -at the groups holding wide fronts (full overlap deadlocks, #109), which needs -private contribution blocks or DAG-aware stack lifetimes (#109) and atomic -extend-add of concurrent children; (5) an asynchronous vendor `potrf` whose -`info` stays on the device if #110 is confirmed; (6) `subtree_budgets` default -reconsidered for occupancy (16 KiB beat 48 KiB there) and `regime_c_rows` tuned -per phase, since the knobs are phase-coupled (64 is right for the prototype -factorization, 256 for the stock solve). Negative results to keep: CUDA-graph -replay is neutral above n ≈ 5e4 (kernel-busy), tiled extend-add is at its -scatter-bandwidth floor, packed-triangle iteration inside element loops loses, -amalgamation changes regress, a 96-wide local-memory class loses to vendor -`potrf`. `bench/e2e/MadNLPSDS.jl` is the end-to-end check. - -**Owner note (issue #75)**: the device LDLᵀ of T15 follows the reference -pivot for pivot and pays for it in three places that are deliverables here: -(1) the pivot search runs on one work item (`_lt_choose`, `O(w·f)` per -column when the Bunch–Kaufman choice fails the threshold and `_best_1x1` -scans the block): parallel reductions for the column maxima, `λ`/`σ` and the -`_best_1x1` scan, and a cheaper fallback, in the reference too; (2) regime-C -fronts run the fused KA kernel with one workgroup each and no vendor call -(the task's `sytrf` path was dropped by owner decision, since `sytrf` picks -its own pivot order on `F₁₁` and the reference equality test must hold): KA -in-front pivoting of `F₁₁` keeping the reference sequence, then vendor -`trsm`/`gemm` for `L₂₁` and the contribution block; (3) regime B keeps `F₁₁` -in global memory: stage it in `@localmem` as the Cholesky kernel does. -Baseline: `kkt_matrix(Float64, 3000, 1000, 1e-8)`, default analysis, 13.5 s -device / 12.9 s reference on the KA CPU backend (T15 report); the Report -gives the same numbers after, plus CUDA. Closes #75. -Delivered for #75 (perf PRs from `perf/exp1-pivot-search`, `perf/exp1-regime-b-local`, `perf/exp2-regime-c-blas`; numbers in PERFORMANCE.md "Experiments 1–2 results"): (1) cooperative pivot search, cheaper exact fallback in the reference; (3) regime B with F₁₁ in local memory, no scale phase, tiled contribution-block update; (2) regime C (and wide tall regime-B bins) as blocked pivot steps plus GEMMs through the dense interface, concurrent per launch group. Open: the fallback-heavy kkt(3000,1000,1e-8) (CUDA 2.5 s, KA CPU 3.3 s vs reference 2.6 s) and LDLᵀ/Cholesky 1.3–2.6×. -#96: matching pairs cost 2× nnz(L); re-measure after the T25 level merging whether the extra fronts matter; otherwise accepted. - -Partitioned-inverse solve (`solve_alg = "algo1"`), CUDA sync-free forward -sweep behind a capability check, CUDA graph capture of refactorize+solve, -level merging, amalgamation/bin tuning; tests: all previous suites unchanged; -Report: timing table vs the T04 cuDSS baseline for every harness matrix. - -### T26 — Hybrid memory / hybrid execute `[ ]` +### T24 — Partitioned-inverse fused solve (`solve_alg = "algo1"`, issue #82 solve half) `[ ]` + +**Prototype adopted (PR #107)**: `bench/solve_proto_78k.jl`, +`bench/solve_proto_all.jl`, `bench/final_sweep.jl`, `bench/e2e/SDSProto.jl` +(solve part). On the 78k-bus condensed KKT under the ND ordering: 18.0 ms +(stock level-batched sweeps) → 3.5 ms, cuDSS 3.8 ms, GV100; 20.5 → 8.05 ms on +a Radeon VII. Two states of the prototype were faster and wrong (an +early-release counter bug, relres 5e+11), so every timing in this task is +gated on the residual. + +Deliverables: (1) `solve_alg = "algo1"`: per-front inverses of `L₁₁` for +fronts of width ≤ `W` (default 256, an `Options` keyword), computed at the end +of (re)factorization with the dense interface (`impl` keyword, no direct +vendor calls), stored in a `Numeric` buffer sized at analysis and reported in +`memory_estimates`; every TRSV of the sweeps becomes a GEMV. (2) One fused +dependency-counter kernel per direction covering every front of width ≤ `W`, +regime-A subtrees included; the vendor path only for the wider fronts +(the root on KKT systems). Counters are the PLAN §2.5 place for atomics and +must have an atomic-free variant (`deterministic = true` keeps the level +sweeps); backward counters dispatch parents first or they deadlock. (3) The +strategy is chosen per schedule at analysis: the prototype regresses on +lap3d_40 and apache2 (wide fronts dominate), so `"default"` keeps the level +sweeps and `"algo1"` is chosen automatically only when the schedule +statistics predict a win; the Report states the rule and the full SPD-harness +table that justifies it. (4) Cholesky and LDLᵀ (the inverse is of `L₁₁`, D is +applied by `solve_diag` as today), real and complex, `nrhs > 1`, both +permutation paths, `solve_mode` variants; LU may stay on the level sweeps with +a documented reason. (5) No host synchronization in the solve +(PLAN §3.9); the per-call synchronization measured in `bench/e2e/MadNLPSDS.jl` +belongs to the wrapper, not here. + +Tests: every solve suite (`test_solve_*`, `test_api`, `test_refinement`, +`test_batch_*`) unchanged under `"default"` and repeated under `"algo1"` +(`SDS_TEST_SOLVE_ALG` or a loop, the task decides), with a relative-residual +gate per matrix; `test_options` accepts `"algo1"` and rejects the unsupported +combinations with `NotSupportedError`. Report: solve table (1 RHS and 6 RHS) +vs the T04 cuDSS baseline for every harness matrix on CPU and CUDA, with the +chosen strategy per matrix. Refs #82 (solve half). + +### T25 — Split factorization: tiled SYRK, chunked TRSM, split regime-C fronts `[ ]` + +**Prototype adopted (PR #107)**: `bench/fact_split.jl`. Refactorization of +the 78k-bus condensed KKT: 167.7 ms (stock) → 138.8 (ND ordering) → 72.7 +(`regime_c_rows = 256`, `subtree_parallelism = 16384`) → 64.0 (split kernels) +→ 40.6 ms (regime-C fronts split, `regime_c_rows = 64`, longest-first +subtree launches); cuDSS 24.3 ms. Factor bit-identical to stock on B-only +configurations, ~3e-10 relative where vendor arithmetic was replaced. + +Deliverables, Cholesky first: (1) the regime-B fused kernel's one-workgroup +SYRK (~12 ms on this matrix) and TRSM are replaced, for fronts with `m > 64`, +by separate many-block kernels: tiled 32×32 SYRK with disjoint tile writes +(deterministic, no atomics, bit-identical) and chunked TRSM; the fused kernel +stays for the small fronts. (2) Regime C is the inverse of regime B — its +vendor dense path costs ~8 API launches per tall-narrow front while its +multi-block assembly kernels are the efficient part — so regime-C fronts +narrower than 64 run the regular assembly kernels plus the split factor +kernels batched per level, and vendor `potrf` is used only above width 64 +(the boundary is measured right: a 96-wide local-memory class loses to +cuSOLVER's blocked `potrf`). (3) Regime-A subtree launches ordered longest +walk first. (4) `regime_c_rows` default re-measured after (2): it matters as +much as the width, and 64 is right once the C fronts are cheap; the knob is +phase-coupled with the stock solve (300 fronts on its per-front vendor path, +215 ms per backsolve), so the default must be chosen with T24's solve, not +alone. (5) The LDLᵀ regime-C path already uses blocked pivot steps plus GEMMs +(#89) and must not regress; if the tiled SYRK applies to `L₂₁ D L₂₁ᴴ` with +the same determinism, say so in the Report, otherwise leave it. + +Negative results that must not be re-run: a tiled extend-add is neutral +(~10 ms of scatter bandwidth is the floor); CUDA-graph replay of the +~1600-launch sequence is bit-identical and timing-neutral (kernel-busy); +amalgamation sweeps under the split kernels all regress, `(32, 0.25, 8)` is +optimal; packed-triangle iteration inside element loops loses (33 → 42 ms); +Float32 factorization is 18% faster and unusable (κ ≈ 1e14, with or without +Jacobi scaling; δ-regularization fails even in FP64). + +Tests: every numeric suite unchanged (`test_numeric_*`, `test_reference_*`, +`test_api`), with the bitwise panel checks of well-conditioned generators +still passing for the split path; a test that exercises each of the three +front classes (fused B, split B, split C) on one matrix. Report: +refactorization table vs the T04 cuDSS baseline on every harness matrix, +Cholesky and LDLᵀ, with the per-class time split from `bench/profile_phases.jl`. + +### T26 — Segmented fused factorization with dependency counters (issues #82, #109, #110) `[ ]` + +**Prototype adopted (PR #107)**: `bench/fact_fused.jl`, `bench/e2e/SDSProto.jl` +(factorization part). On top of T25: 40.6 → 33.5 ms (segmented +dependency-counter kernel) → 33.1 (native ND) → 31.9 ms (`subtree_budgets = +[16384]`); cuDSS 24.3 ms. Remaining buckets at 31.9: mega-kernel 11.0, +subtrees 8.4, vendor dense ~10, statistics 0.6. + +Deliverables: (1) First, verify #110 with a CUDA profile of one +refactorization on a matrix with several regime-C fronts +(`bench/profile_phases.jl`): count the host synchronizations per +factorization. If `_factor_panel_c!`'s `info` check does synchronize, an +asynchronous vendor `potrf` whose `info` stays on the device with the +workspace size queried once at analysis, and the same for the LDLᵀ/LU +regime-C paths; the numeric phase is host-synchronization-free as PLAN §3.9 +states, or the Report says where it is not. Closes or documents #110. (2) +#109 is a prerequisite of (3): the update-stack placement (`build_layout`, +`src/symbolic/layout.jl`, PR #54) derives block lifetimes from schedule steps, +so any out-of-level-order execution reuses a live slot and fails as +mid-column pivot errors far from the cause. Either compute the lifetimes on +the execution order the fused kernel actually follows (interval colouring over +that order), or budget private slots for the non-subtree fronts and report +them in `memory_estimates` (0.14 GB vs 34 MB on the 78k-bus KKT). Closes #109. +(3) One kernel per segment of levels with four workgroup roles (extend-add, +panel, TRSM, SYRK tiles) and per-front dependency counters; segments split at +the launch groups holding wide (vendor) fronts, vendor fronts between +segments: full cross-stream overlap deadlocks structurally on this tree +(62 wide fronts from level ~3, 31.5k dependent blocks saturate the SMs ahead of +the vendor stream). Concurrent extend-add of a parent's children needs atomic +adds; keep the owner-pull level-synchronous path as the atomic-free variant +(`deterministic = true`). The write-through extend-add (a child with a single +consumer writes into a pre-assembled parent) is correct and worth ~0.4%; take +it only if free. (4) `subtree_budgets` default: the 48 KiB class caps the +subtree kernel at one block per SM; `[16384]` was 13.3 → 8.4 ms on that +bucket. Re-measure on the RTX 4080 and the CI runner before changing the +default; `subtree_max_fronts` stays off (neutral, ~3,300 subtrees hide the +188-front walk). (5) The #75 remainder on the LDLᵀ side: the fallback-heavy +`kkt(3000, 1000, 1e-8)` (CUDA 2.5 s, KA CPU 3.3 s vs reference 2.6 s) and the +extend-add width; and #96: re-measure whether the extra fronts of the +matching-based pairs matter once the levels are merged; otherwise accepted. + +**Owner note (issue #75, history)**: the device LDLᵀ of T15 follows the +reference pivot for pivot and paid for it in three places: (1) the pivot +search ran on one work item (`_lt_choose`, `O(w·f)` per column when the +Bunch–Kaufman choice fails the threshold and `_best_1x1` scans the block); +(2) regime-C fronts ran the fused KA kernel with one workgroup each and no +vendor call (the `sytrf` path was dropped by owner decision, since `sytrf` +picks its own pivot order on `F₁₁` and the reference equality test must +hold); (3) regime B kept `F₁₁` in global memory. Delivered for #75 (perf PRs +from `perf/exp1-pivot-search`, `perf/exp1-regime-b-local`, +`perf/exp2-regime-c-blas`; numbers in PERFORMANCE.md "Experiments 1–2 +results"): (1) cooperative pivot search, cheaper exact fallback in the +reference; (3) regime B with `F₁₁` in local memory, no scale phase, tiled +contribution-block update; (2) regime C (and wide tall regime-B bins) as +blocked pivot steps plus GEMMs through the dense interface, concurrent per +launch group. Open: the fallback-heavy case above and LDLᵀ/Cholesky +1.3–2.6×. + +Negative results that must not be re-run: subtree classes on concurrent +streams are neutral and break graph capture; a stream pool for independent +boundary wides is neutral until (1) is done; CUDA-graph replay is neutral +above n ≈ 5e4 and only pays below, so it is not a deliverable here (PLAN +§2.5 keeps it for the CUDA extension when a small-n case asks for it). + +Tests: every numeric suite unchanged under the default; the fused path on +every generator of `test/matrices.jl` with a zero-failing-fronts gate and the +factor compared to the level-synchronous path (`panel_tol`/`growth_tol`); +a test that the lifetimes or private slots of (2) are respected (a +deliberately reordered execution must not clobber a live block: assert +through the factor equality on a matrix with deep CB reuse). Report: +refactorization and refactorization+solve tables vs cuDSS on every harness +matrix, launches per refactorization (#82's criterion: ≥ 3× fewer), the +#110 synchronization count before and after, stack bytes with and without +private slots; `bench/e2e/MadNLPSDS.jl` on the 78k-bus ACOPF as the +end-to-end check. Closes #82. + +### T27 — Hybrid memory / hybrid execute `[ ]` Host-resident panels streamed per level with pinned memory; `hybrid_device_memory_min`; CPU-backend execution of regime A/B levels; test: a factorization under a device memory budget smaller than the factor size completes with the same solution. -### T27 — Robustness extras `[ ]` +### T28 — Robustness extras `[ ]` A posteriori threshold pivoting with optional delayed pivots (host re-analysis), `factor_precision = Float32` with Float64 refinement; tests: delayed-pivot case where static perturbation gives `relres > 1e-6` and APTP gives `≤ 1e-10`; mixed precision reaches Float64 accuracy with FGMRES-IR. +**Owner note (PR #107)**: on the condensed pglib KKT systems (κ ≈ 1e14) a +Float32 factor is unusable with or without Jacobi scaling, and FP64's +single-solve relres of 1.4e-2 is itself the κ·ε floor, consistent with cuDSS's +documented FP32 failures there. The mixed-precision mode is therefore for +better-conditioned systems and for Metal; it needs FGMRES-IR against the +Float32 factor (T18) and a documented failure path, not a precision trick at +the factorization level. + +### T29 — Non-uniform batch (`BatchedDirectSolver`) `[ ]` + +Block-diagonal packing, forest schedule; `test_nonuniform_batch_cudss.jl` ported. + +### T30 — ND partition-tree export/import and ordering cache `[ ]` + +`nd_partition_tree`/`user_nd_partition_tree` in the cuDSS binary-tree +encoding; test: importing the exported tree reproduces the same supernode +partition and `nnz_L`. + +**Owner note (PR #107)**: the ordering study found that SDS's own ND +(`reordering_alg = "algo4"`) matches the imported cuDSS ordering on the +78k-bus KKT (schedule depth 29 vs 30), so the tree import is a compatibility +feature, not a performance one. The ordering cache (analysis under a stored +permutation ran in 3.6 s against 12 s for a fresh METIS analysis) is the part +MadNLP benefits from. + ### External — MadNLPGPU integration (in the MadNLP repository) `[ ]` Add a `SparseDirectSolver`-backed `AbstractLinearSolver` next to @@ -3051,8 +3200,9 @@ upper triangle with view `'U'`, as `CUDSSSolver` does; values aliased with an `update!` fallback; Cholesky `info` mapped to inertia), backend-agnostic, run on the 78k-bus ACOPF with `SparseCondensedKKTSystem` on CUDA and AMD with the same 101 iterations and objective as cuDSS. Start from it. Known points: set -`reordering_alg = "algo4"` (#108); the tuning knobs are phase-coupled -(`regime_c_rows = 256` balances the stock factorization and solve); the 2×2 -pivot pairs of `"S"` are chosen from the values present at analysis, so -analyse after the first KKT assembly; the per-call synchronization of the -wrapper (`asynchronous = false`) is part of the measured gap. +`reordering_alg = "algo4"` until T22 merges (#108); the tuning knobs are +phase-coupled (`regime_c_rows = 256` balances the stock factorization and +solve; T25 retunes them); the 2×2 pivot pairs of `"S"` are chosen from the +values present at analysis, so analyse after the first KKT assembly; the +per-call synchronization of the wrapper (`asynchronous = false`) is part of +the measured gap. diff --git a/bench/README.md b/bench/README.md index 7fd76fd..351dd2c 100644 --- a/bench/README.md +++ b/bench/README.md @@ -96,7 +96,7 @@ julia --project=bench bench/compare.jl --solver=sds --features=ldlt solve (no residual). `nubatch` batches all selected condensed dumps into one row. The SDS calls for uniform batches assume `(n, nbatch)` right-hand sides and the non-uniform batch assumes `BatchedDirectSolver` (PLAN §3); adjust - `compare.jl` when T17/T22 fix the API. + `compare.jl` when T17/T29 fix the API. * **Adding a feature.** One `Feature(...)` entry in `features.jl` (task id, structure, matrix selector, parameters); `compare.jl` handles the kinds `:single`, `:ubatch`, `:nubatch` and `:schur`. From 9eb21453e269156a5bdad235b684dc70b7c2dc1c Mon Sep 17 00:00:00 2001 From: Michel Schanen Date: Fri, 9 Oct 2026 13:21:20 +0000 Subject: [PATCH 2/2] Docs: T31 for the oneAPI/Metal remainder of T23 (#112, #53) [skip ci] PR #113 delivers the AMDGPU extension and the local-memory cap alone and opened #112 for the rest of T23. The remainder becomes T31 at the end of the chain: oneAPI and Metal extensions, the #53 allocation-free :ka fallbacks and the select_impl order noted in #112. STATE.md section 8 updated. Refs #112, #53. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01Wzx8WtMEXnAfyzyujdySFh --- STATE.md | 14 ++++++++------ TASKS.md | 30 ++++++++++++++++++++++++++++++ 2 files changed, 38 insertions(+), 6 deletions(-) diff --git a/STATE.md b/STATE.md index aed8870..e8ab98e 100644 --- a/STATE.md +++ b/STATE.md @@ -170,9 +170,11 @@ regime A); **T26** segmented fused factorization with dependency counters, with #109 (stack lifetimes) and #110 (verify the regime-C host sync) inside it, plus the #75 remainder and the #96 re-measurement; **T27** hybrid memory (was T26); **T28** robustness extras (was T27); **T29** non-uniform batch (was -T22); **T30** ND partition-tree export and ordering cache (was T24); External -unchanged and last (owner decision: integrate against finished kernels; the -§5 item 5 recommendation is withdrawn). GitHub: issues #22–#27 renamed to the -new titles, #22 loses `on-hold`, new task issues for T28–T30; #96 and -#108–#110 labelled `triaged` (#108 → T22, #109/#110 → T26). The chain resumes -with T22. +T22); **T30** ND partition-tree export and ordering cache (was T24); **T31** +oneAPI and Metal extensions with the allocation-free `:ka` fallbacks (#112, +#53), the remainder of T23 after PR #113 delivered the AMDGPU extension alone; +External unchanged and last (owner decision: integrate against finished +kernels; the §5 item 5 recommendation is withdrawn). GitHub: issues #22–#27 +renamed to the new titles, #22 loses `on-hold`, new task issues for T28–T31; +#96, #108–#110 and #112 labelled `triaged` (#108 → T22, #109/#110 → T26, +#112 → T31). The chain resumes with T22 once PR #113 (T23) has merged. diff --git a/TASKS.md b/TASKS.md index 982a18a..957dc12 100644 --- a/TASKS.md +++ b/TASKS.md @@ -3187,6 +3187,36 @@ feature, not a performance one. The ordering cache (analysis under a stored permutation ran in 3.6 s against 12 s for a fresh METIS analysis) is the part MadNLP benefits from. +### T31 — oneAPI and Metal extensions, allocation-free KA dense fallbacks (issues #112, #53) `[ ]` + +The remainder of T23, split off by owner decision when PR #113 delivered the +AMDGPU extension and the per-backend local-memory cap alone (#112). Three +parts: (1) `ext/SparseDirectSolverOneAPIExt.jl`: `oneSparseMatrixCSR` adapters, +the T13 API on them, oneMKL dense bindings for the `vendor_*` table, +`max_local_bytes(::oneAPIBackend)` from the device; (2) +`ext/SparseDirectSolverMetalExt.jl`: MPS bindings where they exist or the KA +fallbacks only, `max_local_bytes(::MetalBackend) = 32768` (the cap already +clamps the default `subtree_budgets`), Float32/ComplexF32 only; (3) the owner +note of #53: the `impl = :ka` dense fallbacks (`ka_potrf!`, `ka_trsm!`, +`ka_gemm!` and the batched variants) are the only dense path on these backends +and must be allocation-free per call — the T23 survey found no allocation +inside `src/dense/fallback/*.jl`; the per-call cost is the launch +configuration (`Val`s built from runtime flags such as `Val(ul == 'L')` in +`potrf.jl`, `Val(tA)`/`Val(tB)` and the `Union`-typed tile of `_default_tile` +in `gemm.jl`, kernel objects built per call, keyword calls), so the fix is +static launch configurations resolved at analysis and preallocated workspace +in `Numeric`, asserted with `ka_cpu_alloc_budget` from `test/utils.jl`. Also +from #112: `select_impl` resolves `:auto` to the allocating `:generic` path +before `:ka` on a backend without `vendor_*` bindings whose `mul!`/`cholesky!` +probes pass, so on oneAPI/Metal `:auto` must prefer `:ka`. `:generic` stays +the allocating reference path (PLAN §3.9). + +Written by analogy with the CUDA and AMDGPU extensions. Test on the owner's +machine: a temporary environment that adds oneAPI and Metal and precompiles +both extensions (no hardware needed to load them on Linux; Metal by review +only). Mark both as untested in the Report; the `:ka` allocation budget and +the `select_impl` order are tested on the KA CPU backend. Closes #112 and #53. + ### External — MadNLPGPU integration (in the MadNLP repository) `[ ]` Add a `SparseDirectSolver`-backed `AbstractLinearSolver` next to