diff --git a/PERFORMANCE.md b/PERFORMANCE.md index 90508af..3527141 100644 --- a/PERFORMANCE.md +++ b/PERFORMANCE.md @@ -10,7 +10,8 @@ Performance issues carry the GitHub label `performance` ([list](https://github.c | --- | --- | --- | --- | | #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 | -| #75 | device LDLᵀ: serial pivot search, one workgroup per regime-B/C front, no vendor `sytrf` | 1, 2 | open, triaged; step 1 (cooperative pivot search) and step 2 (regime B in local memory) in PRs from `perf/exp1-pivot-search`, `perf/exp1-regime-b-local` | +| #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 | @@ -181,6 +182,47 @@ LDLᵀ refactorization, `main` → step 1 → step 2: K2 dumps: inertia and `nperturbed` unchanged on all nine. T15 matrix on CUDA: 8.3 → 5.6 s. Experiment-1 criterion still not met (lap2d 3.7×, bcsstk17 4.1×, case1354 condensed 2.6× the Cholesky time); KKT dumps within +3–6% of step 1. +### Step 3: regime C through GEMMs + +Every front on the regime-C path, and every regime-B bin of width class ≥ 32 and row class ≥ 256 (`ldlt_blocked_path`), now runs its pivot steps in static blocks of `nb = 32` columns (`src/numeric/ldlt_c.jl`). Per block, `panel_ldlt_kernel!` (one workgroup per front, the fronts of a launch group concurrently in chunks whose workspace is at most twice the largest front's) runs the reference's pivot search on a lazily updated panel: the block's pivots are kept as columns of `Lb` (multipliers) and `Wb` (unscaled pivot columns) in the workspace, and an entry of a column that is not a pivot column yet is read as its stored value minus the pending pivots, in pivot order with 2×2 pivots as one term (the arithmetic of the right-looking reference). After the block, one GEMM updates the trailing fully-summed columns and one accumulates the contribution block, through `src/dense/interface.jl` (`:vendor` on CUDA, `:generic`/`:ka` elsewhere); `pack_add!` adds it to the update stack. A 2×2 pivot that starts at a block's last column takes the next block's first column, which is saved across the GEMM and restored by the next block. The fallback scan of this kernel tests (candidate column, 64 rows) tasks over the whole workgroup. Launch shapes depend on the analysis only (no host synchronization); `"algo2"` with empty `subtree_budgets` still forces every front onto this path, and a new test compares panels, D, `piv` and `pivot_kind` with the reference there for every element type, block sizes that make 2×2 pivots straddle blocks, and every dense implementation. + +LDLᵀ refactorization, `main` → step 3: + +| matrix | main ms | step 1 ms | step 2 ms | step 3 ms | step 3 / cuDSS | step 3 / SDS Cholesky | +| --- | --- | --- | --- | --- | --- | --- | +| lap2d_300 | 228 | 192 | 116 | 41.2 | 7.99× | 1.31× | +| lap3d_40 | 14,777* | 14,494* | 8,753* | 424 | 7.27× | 1.89× | +| HB/bcsstk17 | 206 | 167 | 115 | 48.1 | 14.3× | 1.73× | +| Boeing/bcsstk38 | 137 | 106 | 72.9 | 39.2 | 11.8× | 1.73× | +| GHS_psdef/apache2 | 92,576* | 97,296* | 50,736* | 2,120* | 4.33× | 2.06× | +| kkt_pglib_opf_case118_ieee_condensed_1 | 1.77 | 1.82 | 1.91 | 1.94 | 4.96× | 2.12× | +| kkt_pglib_opf_case118_ieee_condensed_10 | 1.72 | 1.79 | 1.87 | 1.83 | 5.25× | 2.08× | +| kkt_pglib_opf_case118_ieee_condensed_20 | 1.87 | 1.82 | 1.88 | 1.87 | 5.34× | 2.14× | +| kkt_pglib_opf_case118_ieee_k2_1 | 1.69 | 1.69 | 1.8 | 1.81 | 4.98× | | +| kkt_pglib_opf_case118_ieee_k2_10 | 1.62 | 1.68 | 1.76 | 1.75 | 4.93× | | +| kkt_pglib_opf_case118_ieee_k2_20 | 1.6 | 1.66 | 1.74 | 1.73 | 4.33× | | +| kkt_pglib_opf_case1354_pegase_condensed_1 | 10.9 | 7.2 | 7.51 | 7.52 | 11× | 2.57× | +| kkt_pglib_opf_case1354_pegase_condensed_10 | 10.9 | 7.3 | 7.56 | 7.56 | 11.8× | 2.59× | +| kkt_pglib_opf_case1354_pegase_condensed_20 | 11 | 7.32 | 7.56 | 7.57 | 10.7× | 2.59× | +| kkt_pglib_opf_case1354_pegase_k2_1 | 8.31 | 6.6 | 6.95 | 6.91 | 8.46× | | +| kkt_pglib_opf_case1354_pegase_k2_10 | 7.15 | 6.05 | 6.27 | 6.32 | 7.64× | | +| kkt_pglib_opf_case1354_pegase_k2_20 | 9.1 | 6.71 | 6.84 | 6.95 | 8.62× | | +| kkt_pglib_opf_case14_ieee_condensed_1 | 0.446 | 0.503 | 0.484 | 0.492 | 2.56× | 1.35× | +| kkt_pglib_opf_case14_ieee_condensed_10 | 0.443 | 0.52 | 0.517 | 0.51 | 2.22× | 1.5× | +| kkt_pglib_opf_case14_ieee_condensed_11 | 0.448 | 0.463 | 0.515 | 0.544 | 1.87× | 1.66× | +| kkt_pglib_opf_case14_ieee_k2_1 | 0.669 | 0.697 | 0.728 | 0.726 | 2.59× | | +| kkt_pglib_opf_case14_ieee_k2_10 | 0.722 | 0.708 | 0.727 | 0.722 | 3.04× | | +| kkt_pglib_opf_case14_ieee_k2_15 | 0.97 | 1.03 | 1.07 | 1.08 | 3.37× | | + +- GEMMs on fronts with `f > 512` (the calls of lap3d_40 and apache2 replayed through `_gemm_impl!(:vendor, …)`): 0.41 TFLOP/s, 56–57% of the 0.73 TFLOP/s cuBLAS DGEMM measures on a 4096² product. +- K2 dumps: inertia and `nperturbed` unchanged on all nine (and equal to `main`). +- `kkt_matrix(Float64, 3000, 1000, 1e-8)`, default analysis and "C only" (`"algo2"`): CUDA 33.9 s (`main`) → 2.5 s; KA CPU backend 7.0 s → 3.3–3.4 s; `ref_ldlt!` 2.5–2.7 s. Half of its pivot steps end in a fallback scan where no column is acceptable (1444 scans, about 76 candidates each), so every candidate must be rejected exactly, and each entry read by the scan subtracts up to 32 pending pivots: with `nb = 8` the same matrix takes 0.86 s. A tiled variant of the scan (pending rows staged in local memory) was slower (4.6 s) because of its barriers; the remaining lever is to keep the scan's entries current instead of lazy, or a smaller `nb` for fronts where the fallback dominates. +- Workspace: lap3d_40 101 MB (Cholesky 62 MB), apache2 223 MB (109 MB), next to a 128/578 MB update stack and 202 MB/1.7 GB of factor. + +**Criteria.** Experiment 1 (LDLᵀ refactorization within 1.5× of SDS Cholesky on the same pattern): met on lap2d_300 (1.31×) and case14 condensed_1 (1.35×), not yet on bcsstk17/bcsstk38 (1.73×), lap3d_40 (1.89×), apache2 (2.06×) and the other KKT dumps (1.5–2.6×); inertia unchanged on all K2 dumps (met). Experiment 2: ≥ 50% of cuBLAS DGEMM peak on fronts > 512 (met, 56%); "C only" faster than the reference on the CPU backend and ≥ 5× faster on the GPU: not met on the fallback-heavy T15 matrix (CPU backend 0.8×, CUDA 1.1× the reference; 13.7× faster than `main` on CUDA). What remains: the fallback scan on such fronts, the latency of the small KKT fronts in regime B (KKT dumps 1.9–12× cuDSS, level- and launch-bound, experiment 5), and `_extend_add_kernel` on lap3d_40 (92 of 424 ms, as for Cholesky). #75 stays open for these. + +`bench/comparison/comparison.{md,png}` are regenerated with the step-3 LDLᵀ rows: SDS/cuDSS geometric mean of "LDLᵀ, static pivoting" (23 matrices) factorization 6.37× → 4.41×, refactorization 9.90× → 6.27×. The "LDLᵀ + 2 refinement steps" row (9 K2 dumps) reads 5.45× → 7.19×, but it is dominated by the three case1354 K2 dumps, whose `compare.jl` medians vary 2–3× between reruns of the same code (8.8–27.6 ms for `main`, 7.5–16.5 ms for step 3, fresh solver per sample); the profile above (6.3–7.0 ms after, 7.2–9.1 ms before) is the reliable number there. + ## Recommendations: Prioritized Experiment Plan | # | Experiment | Payoff | Effort | Key measurement | Success criterion | diff --git a/TASKS.md b/TASKS.md index a2584f5..a6f759b 100644 --- a/TASKS.md +++ b/TASKS.md @@ -2563,6 +2563,7 @@ 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×. Partitioned-inverse solve (`solve_alg = "algo1"`), CUDA sync-free forward sweep behind a capability check, CUDA graph capture of refactorize+solve, diff --git a/bench/comparison/comparison.md b/bench/comparison/comparison.md index 7b90e23..06a7ac4 100644 --- a/bench/comparison/comparison.md +++ b/bench/comparison/comparison.md @@ -12,9 +12,9 @@ Geometric mean of the SDS/cuDSS ratio over the matrices both solvers ran. | --- | --- | --- | --- | --- | --- | --- | --- | | Cholesky, Float64 | T13 | done | 14 | 0.62× | 2.47× | 3.89× | 4.36× | | Cholesky, 16 right-hand sides | T12 | done | 14 | 0.64× | 2.74× | 3.89× | 3.87× | -| LDLᵀ, static pivoting | T15 | done | 23 | 1.15× | 6.37× | 9.90× | 3.96× | -| LDLᵀ + 2 refinement steps, K2 dumps | T16 | done | 9 | 2.69× | 3.18× | 5.45× | 7.57× | -| Uniform batch of 8, Cholesky | T17 | pending | | | | | | +| LDLᵀ, static pivoting | T15 | done | 23 | 1.18× | 4.41× | 6.27× | 4.22× | +| LDLᵀ + 2 refinement steps, K2 dumps | T16 | done | 9 | 2.77× | 5.18× | 7.19× | 7.91× | +| Uniform batch of 8, Cholesky | T17 | done | | | | | | | LU, unsymmetric | T19 | pending | | | | | | | Schur complement, Cholesky | T20 | pending | | | | | | | LDLᵀ + matching (algo5), K2 dumps | T21 | pending | | | | | | @@ -71,31 +71,29 @@ Geometric mean of the SDS/cuDSS ratio over the matrices both solvers ran. | matrix | n | analysis cuDSS ms | SDS ms | ratio | factorization cuDSS ms | SDS ms | ratio | refactorization cuDSS ms | SDS ms | ratio | solve cuDSS ms | SDS ms | ratio | nnz(L) cuDSS | SDS | relres cuDSS | SDS | | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | -| Boeing/bcsstk38 | 8032 | 35.5 | 60.6 | 1.71× | 4.55 | 131 | 28.90× | 3.31 | 131 | 39.63× | 0.476 | 2.03 | 4.26× | 7.99e+05 | 8.05e+05 | 6.5e-11 | 1.42e-10 | -| GHS_psdef/apache2 | 715176 | 4.45e+03 | 7.11e+03 | 1.60× | 518 | 9.08e+04 | 175.17× | 489 | 9.12e+04 | 186.38× | 6.8 | 89.2 | 13.11× | 1.45e+08 | 1.29e+08 | 1.58e-10 | 1.63e-10 | -| HB/bcsstk17 | 10974 | 42.2 | 72 | 1.70× | 5.45 | 198 | 36.32× | 3.37 | 198 | 58.64× | 0.5 | 3.69 | 7.37× | 1.08e+06 | 1.04e+06 | 3.87e-11 | 3.69e-11 | -| kkt_pglib_opf_case118_ieee_condensed_1 | 1088 | 12.9 | 5.75 | 0.45× | 0.905 | 1.93 | 2.13× | 0.39 | 1.7 | 4.36× | 0.383 | 0.741 | 1.93× | 1.52e+04 | 1.16e+04 | 2.99e-05 | 4.3e-05 | -| kkt_pglib_opf_case118_ieee_condensed_10 | 1088 | 12 | 5.19 | 0.43× | 0.921 | 2.21 | 2.40× | 0.349 | 2.79 | 7.98× | 0.249 | 1.75 | 7.02× | 1.52e+04 | 1.16e+04 | 9.52e-06 | 1.14e-05 | -| kkt_pglib_opf_case118_ieee_condensed_20 | 1088 | 12.9 | 5.26 | 0.41× | 0.899 | 1.66 | 1.85× | 0.351 | 2.5 | 7.12× | 0.282 | 0.871 | 3.09× | 1.52e+04 | 1.16e+04 | 0.108 | 0.0833 | -| kkt_pglib_opf_case118_ieee_k2_1 | 3150 | 16.2 | 55.3 | 3.41× | 0.845 | 1.5 | 1.78× | 0.363 | 1.98 | 5.45× | 0.229 | 0.897 | 3.92× | 2.02e+04 | 2.1e+04 | 2.4 | 0.0906 | -| kkt_pglib_opf_case118_ieee_k2_10 | 3150 | 16.2 | 55.6 | 3.44× | 0.868 | 1.4 | 1.61× | 0.354 | 1.35 | 3.83× | 0.284 | 1.49 | 5.26× | 2.02e+04 | 2.05e+04 | 10.1 | 0.0192 | -| kkt_pglib_opf_case118_ieee_k2_20 | 3150 | 16.2 | 54.2 | 3.35× | 0.83 | 1.33 | 1.60× | 0.4 | 1.5 | 3.76× | 0.272 | 1.59 | 5.84× | 2.02e+04 | 2.05e+04 | 6.03 | 0.0149 | -| kkt_pglib_opf_case1354_pegase_condensed_1 | 11192 | 53.3 | 54.8 | 1.03× | 1.39 | 9.82 | 7.09× | 0.684 | 9.6 | 14.03× | 0.355 | 1.22 | 3.45× | 1.55e+05 | 1.21e+05 | 0.000296 | 0.000275 | -| kkt_pglib_opf_case1354_pegase_condensed_10 | 11192 | 52.4 | 55.4 | 1.06× | 1.39 | 9.6 | 6.89× | 0.643 | 9.62 | 14.96× | 0.279 | 2.19 | 7.87× | 1.55e+05 | 1.21e+05 | 0.000122 | 9.49e-05 | -| kkt_pglib_opf_case1354_pegase_condensed_20 | 11192 | 53.1 | 54 | 1.02× | 1.65 | 11.1 | 6.72× | 0.706 | 9.73 | 13.79× | 0.611 | 1.46 | 2.38× | 1.55e+05 | 1.21e+05 | 7.55e-05 | 5.02e-05 | -| kkt_pglib_opf_case1354_pegase_k2_1 | 33811 | 97.4 | 1e+03 | 10.31× | 1.53 | 15.6 | 10.19× | 0.816 | 7.24 | 8.87× | 0.421 | 2.64 | 6.26× | 2.02e+05 | 2.16e+05 | 108 | 13 | -| kkt_pglib_opf_case1354_pegase_k2_10 | 33811 | 96.4 | 1.07e+03 | 11.06× | 1.55 | 7.88 | 5.09× | 0.827 | 7.24 | 8.75× | 0.448 | 2.6 | 5.81× | 2.02e+05 | 2.2e+05 | 52.2 | 52.8 | -| kkt_pglib_opf_case1354_pegase_k2_20 | 33811 | 100 | 1.08e+03 | 10.74× | 1.6 | 13 | 8.16× | 0.806 | 11 | 13.65× | 0.412 | 1.69 | 4.09× | 2.02e+05 | 2.21e+05 | 140 | 48 | -| kkt_pglib_opf_case14_ieee_condensed_1 | 118 | 6.65 | 1.13 | 0.17× | 0.226 | 0.386 | 1.71× | 0.192 | 0.437 | 2.27× | 0.321 | 0.539 | 1.68× | 1.27e+03 | 1.03e+03 | 9.81e-07 | 1.17e-06 | -| kkt_pglib_opf_case14_ieee_condensed_10 | 118 | 6.69 | 1.35 | 0.20× | 0.194 | 0.463 | 2.39× | 0.229 | 0.472 | 2.06× | 0.265 | 0.466 | 1.76× | 1.27e+03 | 1.03e+03 | 0.000129 | 0.000139 | -| kkt_pglib_opf_case14_ieee_condensed_11 | 118 | 6.74 | 1.28 | 0.19× | 0.212 | 0.509 | 2.40× | 0.291 | 0.71 | 2.44× | 0.299 | 0.488 | 1.63× | 1.27e+03 | 1.03e+03 | 1.44e-05 | 1.13e-05 | -| kkt_pglib_opf_case14_ieee_k2_1 | 344 | 9.21 | 4.26 | 0.46× | 0.216 | 0.839 | 3.88× | 0.28 | 0.669 | 2.39× | 0.265 | 0.576 | 2.17× | 2.02e+03 | 2.14e+03 | 0.182 | 7.68e-09 | -| kkt_pglib_opf_case14_ieee_k2_10 | 344 | 7.53 | 5.01 | 0.67× | 0.283 | 0.561 | 1.98× | 0.238 | 0.681 | 2.86× | 0.262 | 0.743 | 2.83× | 2.02e+03 | 2.15e+03 | 0.295 | 5.88e-06 | -| kkt_pglib_opf_case14_ieee_k2_15 | 344 | 7.63 | 3.01 | 0.39× | 0.217 | 1.15 | 5.28× | 0.322 | 1.17 | 3.63× | 0.435 | 0.565 | 1.30× | 2.02e+03 | 3.27e+03 | 0.221 | 0.48 | -| lap2d_300 | 90000 | 257 | 352 | 1.37× | 7.13 | 221 | 31.01× | 5.15 | 221 | 42.94× | 0.712 | 3.75 | 5.26× | 2.43e+06 | 2.47e+06 | 1.6e-12 | 1.68e-12 | -| lap3d_40 | 64000 | 332 | 428 | 1.29× | 62.9 | 1.46e+04 | 232.32× | 58.4 | 1.46e+04 | 250.39× | 1.41 | 19 | 13.42× | 1.75e+07 | 1.44e+07 | 7.68e-14 | 9.06e-14 | - -Single run instead of a BenchmarkTools trial (factorization above the `--single-run-above` limit): SDS on GHS_psdef/apache2, SDS on lap3d_40. +| Boeing/bcsstk38 | 8032 | 35.5 | 63 | 1.77× | 4.55 | 39.4 | 8.66× | 3.31 | 39.5 | 11.91× | 0.476 | 2.08 | 4.36× | 7.99e+05 | 8.05e+05 | 6.5e-11 | 1.74e-10 | +| GHS_psdef/apache2 | 715176 | 4.45e+03 | 6.68e+03 | 1.50× | 518 | 2.11e+03 | 4.07× | 489 | 2.12e+03 | 4.33× | 6.8 | 95.1 | 13.99× | 1.45e+08 | 1.29e+08 | 1.58e-10 | 1.42e-10 | +| HB/bcsstk17 | 10974 | 42.2 | 73.1 | 1.73× | 5.45 | 48.7 | 8.93× | 3.37 | 48.4 | 14.33× | 0.5 | 3.77 | 7.55× | 1.08e+06 | 1.04e+06 | 3.87e-11 | 4.42e-11 | +| kkt_pglib_opf_case118_ieee_condensed_1 | 1088 | 12.9 | 5.23 | 0.41× | 0.905 | 2.13 | 2.35× | 0.39 | 2.11 | 5.40× | 0.383 | 0.753 | 1.97× | 1.52e+04 | 1.16e+04 | 2.99e-05 | 4.3e-05 | +| kkt_pglib_opf_case118_ieee_condensed_10 | 1088 | 12 | 5.52 | 0.46× | 0.921 | 1.85 | 2.00× | 0.349 | 1.92 | 5.50× | 0.249 | 0.764 | 3.06× | 1.52e+04 | 1.16e+04 | 9.52e-06 | 1.14e-05 | +| kkt_pglib_opf_case118_ieee_condensed_20 | 1088 | 12.9 | 5.24 | 0.41× | 0.899 | 2.01 | 2.24× | 0.351 | 1.92 | 5.47× | 0.282 | 0.978 | 3.47× | 1.52e+04 | 1.16e+04 | 0.108 | 0.0833 | +| kkt_pglib_opf_case118_ieee_k2_1 | 3150 | 16.2 | 56.6 | 3.49× | 0.845 | 1.84 | 2.17× | 0.363 | 2.03 | 5.60× | 0.229 | 1.17 | 5.13× | 2.02e+04 | 2.1e+04 | 2.4 | 0.256 | +| kkt_pglib_opf_case118_ieee_k2_10 | 3150 | 16.2 | 58.3 | 3.60× | 0.868 | 1.77 | 2.04× | 0.354 | 1.76 | 4.98× | 0.284 | 1.19 | 4.20× | 2.02e+04 | 2.05e+04 | 10.1 | 0.0164 | +| kkt_pglib_opf_case118_ieee_k2_20 | 3150 | 16.2 | 55.2 | 3.42× | 0.83 | 1.76 | 2.12× | 0.4 | 1.76 | 4.39× | 0.272 | 0.903 | 3.32× | 2.02e+04 | 2.05e+04 | 6.03 | 0.0293 | +| kkt_pglib_opf_case1354_pegase_condensed_1 | 11192 | 53.3 | 56.3 | 1.05× | 1.39 | 7.5 | 5.41× | 0.684 | 7.48 | 10.92× | 0.355 | 1.39 | 3.92× | 1.55e+05 | 1.21e+05 | 0.000296 | 0.000275 | +| kkt_pglib_opf_case1354_pegase_condensed_10 | 11192 | 52.4 | 55.8 | 1.07× | 1.39 | 7.56 | 5.43× | 0.643 | 7.58 | 11.80× | 0.279 | 1.6 | 5.74× | 1.55e+05 | 1.21e+05 | 0.000122 | 0.000105 | +| kkt_pglib_opf_case1354_pegase_condensed_20 | 11192 | 53.1 | 55.7 | 1.05× | 1.65 | 7.58 | 4.60× | 0.706 | 7.58 | 10.75× | 0.611 | 1.37 | 2.25× | 1.55e+05 | 1.21e+05 | 7.55e-05 | 4.78e-05 | +| kkt_pglib_opf_case1354_pegase_k2_1 | 33811 | 97.4 | 1.03e+03 | 10.57× | 1.53 | 11.8 | 7.68× | 0.816 | 9.58 | 11.74× | 0.421 | 2.83 | 6.72× | 2.02e+05 | 2.16e+05 | 108 | 16.8 | +| kkt_pglib_opf_case1354_pegase_k2_10 | 33811 | 96.4 | 1.09e+03 | 11.30× | 1.55 | 31 | 20.02× | 0.827 | 13.3 | 16.02× | 0.448 | 2.37 | 5.30× | 2.02e+05 | 2.2e+05 | 52.2 | 55.9 | +| kkt_pglib_opf_case1354_pegase_k2_20 | 33811 | 100 | 1.12e+03 | 11.20× | 1.6 | 42 | 26.33× | 0.806 | 15.2 | 18.92× | 0.412 | 3.53 | 8.57× | 2.02e+05 | 2.21e+05 | 140 | 30.6 | +| kkt_pglib_opf_case14_ieee_condensed_1 | 118 | 6.65 | 1.17 | 0.18× | 0.226 | 0.532 | 2.36× | 0.192 | 0.555 | 2.88× | 0.321 | 1.4 | 4.36× | 1.27e+03 | 1.03e+03 | 9.81e-07 | 1.17e-06 | +| kkt_pglib_opf_case14_ieee_condensed_10 | 118 | 6.69 | 1.35 | 0.20× | 0.194 | 0.56 | 2.89× | 0.229 | 0.649 | 2.83× | 0.265 | 0.515 | 1.95× | 1.27e+03 | 1.03e+03 | 0.000129 | 0.000139 | +| kkt_pglib_opf_case14_ieee_condensed_11 | 118 | 6.74 | 1.31 | 0.19× | 0.212 | 0.559 | 2.63× | 0.291 | 0.549 | 1.89× | 0.299 | 0.431 | 1.44× | 1.27e+03 | 1.03e+03 | 1.44e-05 | 1.13e-05 | +| kkt_pglib_opf_case14_ieee_k2_1 | 344 | 9.21 | 4.68 | 0.51× | 0.216 | 0.794 | 3.67× | 0.28 | 0.807 | 2.88× | 0.265 | 1.89 | 7.12× | 2.02e+03 | 2.14e+03 | 0.182 | 7.68e-09 | +| kkt_pglib_opf_case14_ieee_k2_10 | 344 | 7.53 | 5.24 | 0.70× | 0.283 | 0.814 | 2.87× | 0.238 | 0.788 | 3.31× | 0.262 | 0.597 | 2.28× | 2.02e+03 | 2.15e+03 | 0.295 | 1.18e-05 | +| kkt_pglib_opf_case14_ieee_k2_15 | 344 | 7.63 | 2.9 | 0.38× | 0.217 | 1.17 | 5.37× | 0.322 | 1.11 | 3.44× | 0.435 | 0.59 | 1.35× | 2.02e+03 | 3.27e+03 | 0.221 | 0.359 | +| lap2d_300 | 90000 | 257 | 402 | 1.57× | 7.13 | 41.3 | 5.79× | 5.15 | 41.3 | 8.01× | 0.712 | 4.12 | 5.78× | 2.43e+06 | 2.47e+06 | 1.6e-12 | 1.59e-12 | +| lap3d_40 | 64000 | 332 | 399 | 1.20× | 62.9 | 424 | 6.75× | 58.4 | 424 | 7.27× | 1.41 | 20.4 | 14.39× | 1.75e+07 | 1.44e+07 | 7.68e-14 | 6.26e-14 | ## LDLᵀ + 2 refinement steps, K2 dumps @@ -103,19 +101,19 @@ Single run instead of a BenchmarkTools trial (factorization above the `--single- | matrix | n | analysis cuDSS ms | SDS ms | ratio | factorization cuDSS ms | SDS ms | ratio | refactorization cuDSS ms | SDS ms | ratio | solve cuDSS ms | SDS ms | ratio | nnz(L) cuDSS | SDS | relres cuDSS | SDS | | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | -| kkt_pglib_opf_case118_ieee_k2_1 | 3150 | 15.8 | 55.2 | 3.49× | 0.827 | 1.46 | 1.77× | 0.371 | 1.67 | 4.50× | 0.6 | 3.91 | 6.52× | 2.02e+04 | 2.1e+04 | 0.00442 | 3.23e-09 | -| kkt_pglib_opf_case118_ieee_k2_10 | 3150 | 16 | 55.4 | 3.47× | 0.812 | 1.37 | 1.69× | 0.503 | 1.36 | 2.71× | 0.55 | 4.29 | 7.80× | 2.02e+04 | 2.05e+04 | 0.0164 | 5.59e-09 | -| kkt_pglib_opf_case118_ieee_k2_20 | 3150 | 15.5 | 54.4 | 3.51× | 0.828 | 1.35 | 1.63× | 0.339 | 1.67 | 4.93× | 0.596 | 4.15 | 6.97× | 2.02e+04 | 2.05e+04 | 0.189 | 1.16e-08 | -| kkt_pglib_opf_case1354_pegase_k2_1 | 33811 | 97.7 | 1.01e+03 | 10.34× | 1.56 | 7.27 | 4.67× | 0.931 | 8.99 | 9.65× | 0.814 | 22.3 | 27.41× | 2.02e+05 | 2.16e+05 | 5.11 | 0.0348 | -| kkt_pglib_opf_case1354_pegase_k2_10 | 33811 | 97.4 | 1.07e+03 | 10.96× | 1.53 | 10.3 | 6.74× | 0.813 | 9.07 | 11.16× | 0.978 | 22.8 | 23.35× | 2.02e+05 | 2.2e+05 | 1.91 | 7.94 | -| kkt_pglib_opf_case1354_pegase_k2_20 | 33811 | 97.4 | 1.07e+03 | 11.03× | 1.55 | 11.2 | 7.18× | 0.825 | 11.2 | 13.63× | 0.958 | 22.8 | 23.82× | 2.02e+05 | 2.21e+05 | 13.9 | 0.216 | -| kkt_pglib_opf_case14_ieee_k2_1 | 344 | 7.76 | 4.49 | 0.58× | 0.222 | 0.713 | 3.21× | 0.215 | 0.882 | 4.10× | 0.719 | 2.16 | 3.01× | 2.02e+03 | 2.14e+03 | 3.23e-07 | 6.15e-13 | -| kkt_pglib_opf_case14_ieee_k2_10 | 344 | 7.47 | 5.04 | 0.68× | 0.295 | 0.594 | 2.01× | 0.233 | 0.755 | 3.24× | 1.26 | 1.79 | 1.42× | 2.02e+03 | 2.15e+03 | 9.14e-06 | 2.41e-14 | -| kkt_pglib_opf_case14_ieee_k2_15 | 344 | 7.46 | 2.68 | 0.36× | 0.206 | 0.966 | 4.70× | 0.235 | 0.858 | 3.65× | 0.65 | 2.31 | 3.55× | 2.02e+03 | 3.27e+03 | 2.4e-08 | 5.69e-09 | +| kkt_pglib_opf_case118_ieee_k2_1 | 3150 | 15.8 | 58.4 | 3.69× | 0.827 | 1.86 | 2.25× | 0.371 | 1.86 | 5.02× | 0.6 | 4.72 | 7.87× | 2.02e+04 | 2.1e+04 | 0.00442 | 5.06e-09 | +| kkt_pglib_opf_case118_ieee_k2_10 | 3150 | 16 | 57.2 | 3.59× | 0.812 | 1.79 | 2.20× | 0.503 | 2.22 | 4.42× | 0.55 | 4.17 | 7.58× | 2.02e+04 | 2.05e+04 | 0.0164 | 1.14e-08 | +| kkt_pglib_opf_case118_ieee_k2_20 | 3150 | 15.5 | 57.2 | 3.69× | 0.828 | 2.03 | 2.45× | 0.339 | 1.78 | 5.24× | 0.596 | 5.12 | 8.59× | 2.02e+04 | 2.05e+04 | 0.189 | 7.84e-09 | +| kkt_pglib_opf_case1354_pegase_k2_1 | 33811 | 97.7 | 1.06e+03 | 10.86× | 1.56 | 22.2 | 14.26× | 0.931 | 17.2 | 18.42× | 0.814 | 28.8 | 35.36× | 2.02e+05 | 2.16e+05 | 5.11 | 0.0498 | +| kkt_pglib_opf_case1354_pegase_k2_10 | 33811 | 97.4 | 1.1e+03 | 11.28× | 1.53 | 33.9 | 22.08× | 0.813 | 26.5 | 32.65× | 0.978 | 24.9 | 25.47× | 2.02e+05 | 2.2e+05 | 1.91 | 14 | +| kkt_pglib_opf_case1354_pegase_k2_20 | 33811 | 97.4 | 1.09e+03 | 11.24× | 1.55 | 25.1 | 16.14× | 0.825 | 12.2 | 14.77× | 0.958 | 25.8 | 26.98× | 2.02e+05 | 2.21e+05 | 13.9 | 0.145 | +| kkt_pglib_opf_case14_ieee_k2_1 | 344 | 7.76 | 4.4 | 0.57× | 0.222 | 0.747 | 3.36× | 0.215 | 0.743 | 3.45× | 0.719 | 1.82 | 2.53× | 2.02e+03 | 2.14e+03 | 3.23e-07 | 6.15e-13 | +| kkt_pglib_opf_case14_ieee_k2_10 | 344 | 7.47 | 5.09 | 0.68× | 0.295 | 0.737 | 2.50× | 0.233 | 0.729 | 3.13× | 1.26 | 1.7 | 1.35× | 2.02e+03 | 2.15e+03 | 9.14e-06 | 2.64e-14 | +| kkt_pglib_opf_case14_ieee_k2_15 | 344 | 7.46 | 2.74 | 0.37× | 0.206 | 1.07 | 5.21× | 0.235 | 1.08 | 4.61× | 0.65 | 1.85 | 2.85× | 2.02e+03 | 3.27e+03 | 2.4e-08 | 6.19e-09 | ## Uniform batch of 8, Cholesky -`ubatch8`, T17 (pending), structure `SPD`, Float64, nrhs = 1, uniform batch of 8. +`ubatch8`, T17 (done), structure `SPD`, Float64, nrhs = 1, uniform batch of 8. | matrix | n | analysis cuDSS ms | SDS ms | ratio | factorization cuDSS ms | SDS ms | ratio | refactorization cuDSS ms | SDS ms | ratio | solve cuDSS ms | SDS ms | ratio | nnz(L) cuDSS | SDS | relres cuDSS | SDS | | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | @@ -220,3 +218,4 @@ Not run yet. * cuDSS 0.8.0 on NVIDIA GeForce RTX 4080, repository 57e2ea9+dirty, 2026-10-02. * SparseDirectSolver.jl 0.1.0 on NVIDIA GeForce RTX 4080, repository 9d280d0+dirty, 2026-10-02. +* SparseDirectSolver.jl 0.1.0 on NVIDIA GeForce RTX 4080, repository f2e3dfe+dirty, 2026-10-05. diff --git a/bench/comparison/comparison.png b/bench/comparison/comparison.png index 5c35c73..e6cff65 100644 Binary files a/bench/comparison/comparison.png and b/bench/comparison/comparison.png differ diff --git a/bench/comparison/sds.csv b/bench/comparison/sds.csv index 1141b16..1ed9118 100644 --- a/bench/comparison/sds.csv +++ b/bench/comparison/sds.csv @@ -27,35 +27,35 @@ sds,cholesky_nrhs16,T12,kkt_pglib_opf_case14_ieee_condensed_10,SPD,118,1282,Floa sds,cholesky_nrhs16,T12,kkt_pglib_opf_case14_ieee_condensed_11,SPD,118,1282,Float64,16,1,0.001020132,0.000425807,0.000357847,0.000449716,0.000952953,0.000349898,0.000330647,0.000428507,1027,15193.0,33,1.4705124707443054e-5,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 sds,cholesky_nrhs16,T12,lap2d_300,SPD,90000,448800,Float64,16,1,0.325324887,0.03152837,0.031424426,0.008513792,0.314441856,0.031439825,0.031318001,0.008498163,2465905,3.73611053e8,20495,1.7851598547374033e-12,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 sds,cholesky_nrhs16,T12,lap3d_40,SPD,64000,438400,Float64,16,1,0.386902999,0.223042271,0.22320726,0.049068244,0.385093288,0.222839828,0.222611719,0.048659465,14387160,1.6401480134e10,11066,8.448190094604811e-14,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,Boeing/bcsstk38,S,8032,355460,Float64,1,1,0.060624225,0.131407119,0.131362702,0.002029925,0.059439843,0.13137589,0.131347536,0.001910636,804613,1.40793119e8,604,1.4180843867122588e-10,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,GHS_psdef/apache2,S,715176,4817870,Float64,1,1,7.107202684,90.755114345,91.203945714,0.08917829,7.107202684,90.755114345,91.203945714,0.08917829,129200941,1.66491403722e11,123689,1.629045645261273e-10,1,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,HB/bcsstk17,S,10974,428650,Float64,1,1,0.071961975,0.197887792,0.197856718,0.003685172,0.070547262,0.197867827,0.197776184,0.003669032,1043601,1.6795197e8,1300,3.693578455267767e-11,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case118_ieee_condensed_1,S,1088,12860,Float64,1,1,0.005745238,0.001926796,0.001701419,0.000740874,0.005228572,0.001491029,0.001513493,0.000723474,11634,219042.0,303,4.302669682539129e-5,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case118_ieee_condensed_10,S,1088,12860,Float64,1,1,0.005194696,0.002214169,0.002787787,0.001750367,0.005080131,0.001637618,0.001799307,0.001066752,11634,219042.0,303,1.1399610969542334e-5,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case118_ieee_condensed_20,S,1088,12860,Float64,1,1,0.00525929,0.001663352,0.002498752,0.000871243,0.00520416,0.001584528,0.001724887,0.000754464,11634,219042.0,303,0.08325458299212221,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case118_ieee_k2_1,S,3150,17714,Float64,1,1,0.055343651,0.00150123,0.001978034,0.000896683,0.054610868,0.00139382,0.00142681,0.000823674,21032,308794.0,1917,0.09056867773597634,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case118_ieee_k2_10,S,3150,17714,Float64,1,1,0.0555984,0.00139845,0.001353419,0.001494139,0.054707287,0.001327759,0.001330961,0.000849923,20542,292847.0,1906,0.01922668190063813,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case118_ieee_k2_20,S,3150,17714,Float64,1,1,0.05421622,0.00133248,0.001501393,0.001588253,0.053774374,0.001320149,0.001313525,0.000957439,20507,291377.0,1889,0.014854294291476057,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case1354_pegase_condensed_1,S,11192,136724,Float64,1,1,0.054796466,0.009824687,0.009603797,0.001224495,0.054552087,0.009600263,0.009591012,0.001200021,120598,2.787739e6,3134,0.00027545234473510356,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case1354_pegase_condensed_10,S,11192,136724,Float64,1,1,0.055422452,0.009602167,0.009617228,0.002194214,0.05428184,0.009593329,0.009571867,0.00123898,120598,2.787739e6,3134,9.488278667035474e-5,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case1354_pegase_condensed_20,S,11192,136724,Float64,1,1,0.054016472,0.011066087,0.009729919,0.00145556,0.053790238,0.010603346,0.009716458,0.00121572,120598,2.787739e6,3134,5.0192353762337224e-5,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case1354_pegase_k2_1,S,33811,188063,Float64,1,1,1.004016075,0.015597837,0.007240037,0.002638421,0.99742274,0.007233705,0.007170497,0.001572179,215953,3.404739e6,19624,13.01856461229925,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case1354_pegase_k2_10,S,33811,188063,Float64,1,1,1.065767901,0.007882551,0.007240921,0.00260473,1.062880945,0.007190032,0.00628589,0.00135963,220146,3.494638e6,19496,52.848192429671556,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case1354_pegase_k2_20,S,33811,188063,Float64,1,1,1.077886349,0.013019252,0.011003308,0.001685644,1.073363761,0.008562279,0.007917806,0.001496138,221292,3.599399e6,19600,47.992929768911,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case14_ieee_condensed_1,S,118,1282,Float64,1,1,0.001127281,0.000385627,0.000436776,0.000539256,0.001057122,0.000364507,0.000371267,0.000429547,1027,15193.0,33,1.1657999541762487e-6,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case14_ieee_condensed_10,S,118,1282,Float64,1,1,0.001352675,0.000463447,0.000472476,0.000465566,0.001049342,0.000369728,0.000429867,0.000428587,1027,15193.0,33,0.00013880687786021974,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case14_ieee_condensed_11,S,118,1282,Float64,1,1,0.001277951,0.000508636,0.000710394,0.000487586,0.001036812,0.000442206,0.000481592,0.000427987,1027,15193.0,33,1.1319105246087817e-5,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case14_ieee_k2_1,S,344,1924,Float64,1,1,0.004259567,0.000838714,0.000668709,0.000575875,0.004160759,0.000555685,0.000589436,0.000542126,2139,26663.0,186,7.679955110816424e-9,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case14_ieee_k2_10,S,344,1924,Float64,1,1,0.005011122,0.000560776,0.000680865,0.000743114,0.004824065,0.000556056,0.000606406,0.000520316,2154,27106.0,185,5.877775182937206e-6,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,kkt_pglib_opf_case14_ieee_k2_15,S,344,1924,Float64,1,1,0.003013006,0.001146422,0.001168541,0.000565445,0.00272051,0.000994923,0.000940113,0.000559966,3273,56950.0,85,0.480402287562271,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,lap2d_300,S,90000,448800,Float64,1,1,0.352224026,0.221263451,0.221137994,0.003748761,0.347788493,0.220387852,0.221107635,0.003679562,2465905,3.73611053e8,20495,1.6779389216462592e-12,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt,T15,lap3d_40,S,64000,438400,Float64,1,1,0.427970632,14.612296974,14.613361251,0.018982602,0.427970632,14.612296974,14.613361251,0.018982602,14387160,1.6401480134e10,11066,9.063326423318782e-14,1,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt_ir2,T16,kkt_pglib_opf_case118_ieee_k2_1,S,3150,17714,Float64,1,1,0.055196048,0.001461409,0.001667507,0.003909787,0.054259484,0.00141077,0.001403119,0.003774639,21032,308794.0,1917,3.2260005214076745e-9,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt_ir2,T16,kkt_pglib_opf_case118_ieee_k2_10,S,3150,17714,Float64,1,1,0.055369095,0.00136816,0.001362058,0.004290508,0.054318588,0.00133382,0.001339428,0.003783444,20542,292847.0,1906,5.5897389139096295e-9,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt_ir2,T16,kkt_pglib_opf_case118_ieee_k2_20,S,3150,17714,Float64,1,1,0.054366693,0.001350819,0.001671586,0.004152631,0.053275971,0.001335069,0.001420728,0.0037175,20507,291377.0,1889,1.160357102037044e-8,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt_ir2,T16,kkt_pglib_opf_case1354_pegase_k2_1,S,33811,188063,Float64,1,1,1.010063189,0.007270851,0.008991956,0.022317976,1.003646642,0.007231236,0.007204166,0.021328615,215953,3.404739e6,19624,0.034762418896692626,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt_ir2,T16,kkt_pglib_opf_case1354_pegase_k2_10,S,33811,188063,Float64,1,1,1.067244555,0.01034738,0.00906991,0.022827412,1.061806312,0.006174879,0.006321956,0.021533663,220146,3.494638e6,19496,7.936310503696266,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt_ir2,T16,kkt_pglib_opf_case1354_pegase_k2_20,S,33811,188063,Float64,1,1,1.07477958,0.011161858,0.011246247,0.022811512,1.063724657,0.009126394,0.007992944,0.021728675,221292,3.599399e6,19600,0.2161981921511238,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt_ir2,T16,kkt_pglib_opf_case14_ieee_k2_1,S,344,1924,Float64,1,1,0.004491031,0.000713474,0.000881953,0.002160144,0.004418033,0.000551436,0.000564095,0.001663596,2139,26663.0,186,6.147932736718633e-13,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt_ir2,T16,kkt_pglib_opf_case14_ieee_k2_10,S,344,1924,Float64,1,1,0.005043358,0.000593745,0.000754834,0.001786805,0.004973018,0.000586476,0.000562155,0.001700796,2154,27106.0,185,2.4139806966654987e-14,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 -sds,ldlt_ir2,T16,kkt_pglib_opf_case14_ieee_k2_15,S,344,1924,Float64,1,1,0.002677508,0.000966181,0.000857869,0.002308462,0.002588568,0.000854268,0.000827763,0.00185908,3273,56950.0,85,5.6903688411491105e-9,5,ok,NVIDIA GeForce RTX 4080,0.1.0,9d280d0+dirty,2026-10-02 +sds,ldlt,T15,Boeing/bcsstk38,S,8032,355460,Float64,1,1,0.062968879,0.039384778,0.039483886,0.00207725,0.060976728,0.03916677,0.039337186,0.001828223,804613,1.40793119e8,604,1.743507164976252e-10,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,GHS_psdef/apache2,S,715176,4817870,Float64,1,1,6.677054091,2.109598742,2.119594559,0.095149225,6.493347834,2.106747765,2.118319932,0.089544662,129200941,1.66491403722e11,123689,1.4191728752798158e-10,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,HB/bcsstk17,S,10974,428650,Float64,1,1,0.073100648,0.048676182,0.048350411,0.003773646,0.071383332,0.048520903,0.048117996,0.003548618,1043601,1.6795197e8,1300,4.4243215645658814e-11,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case118_ieee_condensed_1,S,1088,12860,Float64,1,1,0.005229271,0.00212726,0.002107751,0.000753173,0.005134034,0.001901992,0.001928053,0.000724553,11634,219042.0,303,4.302669682539129e-5,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case118_ieee_condensed_10,S,1088,12860,Float64,1,1,0.00551601,0.001847123,0.001922434,0.000763553,0.005373031,0.001830963,0.001837252,0.000750693,11634,219042.0,303,1.1399610969542334e-5,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case118_ieee_condensed_20,S,1088,12860,Float64,1,1,0.005241282,0.002013351,0.001919132,0.000978281,0.005160173,0.001905302,0.001871533,0.000761723,11634,219042.0,303,0.08325458299212221,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case118_ieee_k2_1,S,3150,17714,Float64,1,1,0.056612879,0.001836794,0.002033351,0.001173729,0.055637387,0.001811383,0.001835314,0.001026621,21032,308794.0,1917,0.2560147626226783,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case118_ieee_k2_10,S,3150,17714,Float64,1,1,0.058311374,0.001774064,0.001762103,0.001193274,0.056139343,0.001748855,0.001749604,0.000876862,20542,292847.0,1906,0.016404740598452035,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case118_ieee_k2_20,S,3150,17714,Float64,1,1,0.055244091,0.001758255,0.001755505,0.000903001,0.054896424,0.001749994,0.001738914,0.000845742,20507,291377.0,1889,0.029292850294016723,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case1354_pegase_condensed_1,S,11192,136724,Float64,1,1,0.056254742,0.007498594,0.007475191,0.001393398,0.054932283,0.007463082,0.007445301,0.001226739,120598,2.787739e6,3134,0.00027487128532320196,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case1354_pegase_condensed_10,S,11192,136724,Float64,1,1,0.055807575,0.007559461,0.00758465,0.001598726,0.055074512,0.00754294,0.007550822,0.001456526,120598,2.787739e6,3134,0.00010523920564427243,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case1354_pegase_condensed_20,S,11192,136724,Float64,1,1,0.055650012,0.007578071,0.007584281,0.001371909,0.055168273,0.007559061,0.00755278,0.001187149,120598,2.787739e6,3134,4.784066937614418e-5,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case1354_pegase_k2_1,S,33811,188063,Float64,1,1,1.02918183,0.011758811,0.009582122,0.002831046,1.019129191,0.007035615,0.008557481,0.001579776,215953,3.404739e6,19624,16.778159435039974,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case1354_pegase_k2_10,S,33811,188063,Float64,1,1,1.089173672,0.031010106,0.013253068,0.002372669,1.082955646,0.018037445,0.012203497,0.001662275,220146,3.494638e6,19496,55.93851084192943,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case1354_pegase_k2_20,S,33811,188063,Float64,1,1,1.124377108,0.042017733,0.01524833,0.003531017,1.091390706,0.018712367,0.008053956,0.002838308,221292,3.599399e6,19600,30.612483225129548,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case14_ieee_condensed_1,S,118,1282,Float64,1,1,0.001165129,0.000532405,0.000554984,0.001399987,0.00107508,0.00052,0.000512615,0.000415187,1027,15193.0,33,1.1657999541762487e-6,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case14_ieee_condensed_10,S,118,1282,Float64,1,1,0.001353587,0.00055983,0.000648634,0.000515296,0.001195658,0.000515655,0.000565065,0.000421146,1027,15193.0,33,0.00013880687786021974,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case14_ieee_condensed_11,S,118,1282,Float64,1,1,0.001313368,0.000558585,0.000548745,0.000430796,0.001093181,0.000543525,0.000526245,0.000403406,1027,15193.0,33,1.1319105246087817e-5,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case14_ieee_k2_1,S,344,1924,Float64,1,1,0.004679736,0.000794003,0.000806823,0.001886611,0.004339701,0.000757934,0.000758823,0.000531426,2139,26663.0,186,7.679955110816424e-9,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case14_ieee_k2_10,S,344,1924,Float64,1,1,0.005237051,0.000814233,0.000788332,0.000596924,0.005145581,0.000760023,0.000731943,0.000574444,2154,27106.0,185,1.1828974137803157e-5,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,kkt_pglib_opf_case14_ieee_k2_15,S,344,1924,Float64,1,1,0.002898651,0.001166131,0.001105079,0.000589565,0.002784254,0.0010763,0.00107035,0.000549265,3273,56950.0,85,0.35857303761291065,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,lap2d_300,S,90000,448800,Float64,1,1,0.401634232,0.04132111,0.041272569,0.004115072,0.395508976,0.040985408,0.040994424,0.003752336,2465905,3.73611053e8,20495,1.5882180912338264e-12,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt,T15,lap3d_40,S,64000,438400,Float64,1,1,0.398787987,0.424378052,0.424261673,0.020353114,0.39314832,0.423520435,0.422388961,0.018779947,14387160,1.6401480134e10,11066,6.258065839656417e-14,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt_ir2,T16,kkt_pglib_opf_case118_ieee_k2_1,S,3150,17714,Float64,1,1,0.05844427,0.001860962,0.001861943,0.004724067,0.056737158,0.001823284,0.001824653,0.004438679,21032,308794.0,1917,5.058367418552032e-9,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt_ir2,T16,kkt_pglib_opf_case118_ieee_k2_10,S,3150,17714,Float64,1,1,0.057230783,0.001789773,0.00222486,0.004168751,0.055983854,0.001765304,0.001784545,0.003939333,20542,292847.0,1906,1.1410435220446759e-8,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt_ir2,T16,kkt_pglib_opf_case118_ieee_k2_20,S,3150,17714,Float64,1,1,0.057170223,0.002032782,0.001776044,0.005117653,0.056759876,0.001766973,0.001759813,0.003802636,20507,291377.0,1889,7.84380738484693e-9,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt_ir2,T16,kkt_pglib_opf_case1354_pegase_k2_1,S,33811,188063,Float64,1,1,1.061296408,0.022171565,0.017156981,0.028790385,1.049282509,0.020592188,0.007943356,0.022750059,215953,3.404739e6,19624,0.04977202418560073,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt_ir2,T16,kkt_pglib_opf_case1354_pegase_k2_10,S,33811,188063,Float64,1,1,1.098027935,0.033896698,0.026534946,0.024909151,1.086466121,0.015682056,0.015764034,0.022255283,220146,3.494638e6,19496,13.953438350645612,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt_ir2,T16,kkt_pglib_opf_case1354_pegase_k2_20,S,33811,188063,Float64,1,1,1.094831637,0.02509138,0.012189538,0.025835582,1.08671851,0.007782348,0.007765178,0.021902928,221292,3.599399e6,19600,0.14454895243417556,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt_ir2,T16,kkt_pglib_opf_case14_ieee_k2_1,S,344,1924,Float64,1,1,0.00439831,0.000747493,0.000742755,0.001822133,0.00431144,0.000725404,0.000732633,0.001635245,2139,26663.0,186,6.147932736718633e-13,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt_ir2,T16,kkt_pglib_opf_case14_ieee_k2_10,S,344,1924,Float64,1,1,0.005088854,0.000736693,0.000729404,0.001701404,0.004994784,0.000717063,0.000718253,0.001690665,2154,27106.0,185,2.6360800359959953e-14,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 +sds,ldlt_ir2,T16,kkt_pglib_opf_case14_ieee_k2_15,S,344,1924,Float64,1,1,0.002735516,0.001071925,0.00108229,0.001853272,0.002639876,0.001066192,0.00105709,0.001770813,3273,56950.0,85,6.185155238243628e-9,5,ok,NVIDIA GeForce RTX 4080,0.1.0,f2e3dfe+dirty,2026-10-05 diff --git a/bench/profile_phases.jl b/bench/profile_phases.jl index 8350299..ad98116 100644 --- a/bench/profile_phases.jl +++ b/bench/profile_phases.jl @@ -12,8 +12,8 @@ # The SDS phases are run with `asynchronous = false`, as in bench/compare.jl. The # wall time is the best of 5 warm runs, each after 0.2 s of device load so that # the GPU is at full clock; the profiled run is a sixth one. -# LDLᵀ is skipped on matrices whose LDLᵀ factorization took more than 10 s in -# bench/comparison/sds.csv (GHS_psdef/apache2, lap3d_40): issue #75. +# LDLᵀ is skipped on GHS_psdef/apache2 (93 s per refactorization before issue #75, +# 2.1 s after; lap3d_40 is profiled again). using LinearAlgebra using SparseArrays @@ -41,7 +41,7 @@ function parse_args(args) end const OPTS = parse_args(ARGS) -const SLOW_LDLT = ("GHS_psdef/apache2", "lap3d_40") +const SLOW_LDLT = ("GHS_psdef/apache2",) const OUTDIR = joinpath(@__DIR__, "profile") # kernel name without the KernelAbstractions/CUDA template arguments diff --git a/src/SparseDirectSolver.jl b/src/SparseDirectSolver.jl index ff7906c..bacea50 100644 --- a/src/SparseDirectSolver.jl +++ b/src/SparseDirectSolver.jl @@ -63,6 +63,7 @@ include("numeric/front.jl") include("numeric/subtree.jl") include("numeric/factorize.jl") include("numeric/ldlt.jl") +include("numeric/ldlt_c.jl") include("numeric/extract.jl") # solve phase on the device (PLAN §2.5) diff --git a/src/numeric/factorize.jl b/src/numeric/factorize.jl index 7c58d8a..3e1332e 100644 --- a/src/numeric/factorize.jl +++ b/src/numeric/factorize.jl @@ -230,7 +230,7 @@ allocated on the device; the panels are bitwise reproducible for a fixed """ function factorize!(N::Numeric{T}, S::Symbolic, nzval::AbstractVector; impl::Symbol = :auto, opts::Options = Options()) where {T} - _is_ldlt_structure(S.structure) && return factorize_ldlt!(N, S, nzval; opts) + _is_ldlt_structure(S.structure) && return factorize_ldlt!(N, S, nzval; impl, opts) _check_numeric(N, S, nzval) p = _front_impls(N, S, impl) plan = N.plan diff --git a/src/numeric/ldlt.jl b/src/numeric/ldlt.jl index 3602afd..d9b9a03 100644 --- a/src/numeric/ldlt.jl +++ b/src/numeric/ldlt.jl @@ -51,7 +51,8 @@ const _LT_STAT = _ST_CTL + 5 # _LT_STAT + q: statistic q ∈ 1:5 (STAT_NPOS const _LT_PHASE = _ST_CTL + 11 # pivot search state: 0 chosen, 2 Bunch–Kaufman pass 2 due, 3 fallback scan due const _LT_RBK = _ST_CTL + 12 # Bunch–Kaufman candidate row r of pass 1 (0: none) const _LT_NTILE = _ST_CTL + 13 # tiles of the F₂₂ update of the regime-B/C kernel (0: no contribution block) -const _LT_CTL = _ST_CTL + 13 +const _LT_KEND = _ST_CTL + 14 # last column a pivot step may start at (w; the block's last column in regime C) +const _LT_CTL = _ST_CTL + 14 # `pv` slots: 1–4 the pivot (d, or the inverse 2×2 block e11, e12, e21, e22), 5–6 the new diagonal # entries of the pivot block (d, or a and c), 7–8 λ and the column maximum of pass 1 @@ -150,6 +151,9 @@ end @inline _lt_front(fa::Tuple, ctl, ::Val{W}) where {W} = @inbounds _SplitFront(fa[1], fa[2], Int(ctl[_ST_LF]) - 1, Int(ctl[_ST_F]), Int(ctl[_ST_W]), W) +# the front as stored, for the interchanges (the regime-C panel kernel moves stored values; see `_LazyFront`) +@inline _lt_raw_front(fa, ctl, pk::Val) = _lt_front(fa, ctl, pk) + # the buffers and the front mode of the regime-B/C kernel: `W = 0` panel in global memory, else F₁₁ staged @inline _lt_fa(L11, factor, ::Val{0}) = factor @inline _lt_fa(L11, factor, ::Val{W}) where {W} = (L11, factor) @@ -290,7 +294,7 @@ end λ = zero(R) r = 0 cm = zero(R) - if k <= w && p.ptype != _LT_PIVOT_NONE + if k <= Int(ctl[_LT_KEND]) && p.ptype != _LT_PIVOT_NONE F = _lt_front(fa, ctl, pk) for i in (k + li):NL:Int(ctl[_ST_F]) a = _fabs(F, i, k) @@ -457,7 +461,7 @@ end ctl[_LT_K] = k % IT ctl[_LT_STEP] = zero(IT) ctl[_LT_PHASE] = zero(IT) - k <= Int(ctl[_ST_W]) || return nothing + k <= Int(ctl[_LT_KEND]) || return nothing if p.ptype == _LT_PIVOT_NONE _lt_select!(fa, ctl, pv, d, pivot_kind, piv, psign, perm, aux, p, k, 0, pk, h) return nothing @@ -565,6 +569,7 @@ end ctl[_LT_R] = zero(IT) ctl[_LT_PHASE] = zero(IT) ctl[_LT_RBK] = zero(IT) + ctl[_LT_KEND] = ctl[_ST_W] for q in 1:5 ctl[_LT_STAT + q] = zero(IT) end @@ -588,11 +593,11 @@ end @inline function _lt_swap!(fa, ctl, pv, piv, li, ::Val{WG}, pk::Val, h::Val) where {WG} @inbounds begin k = Int(ctl[_LT_K]) - if k <= Int(ctl[_ST_W]) + if k <= Int(ctl[_LT_KEND]) r = Int(ctl[_LT_R]) p = r == 0 ? k : k + 1 q = r == 0 ? Int(ctl[_LT_C]) : r - F = _lt_front(fa, ctl, pk) + F = _lt_raw_front(fa, ctl, pk) for i in li:WG:Int(ctl[_ST_F]) if p != q if i < p @@ -647,7 +652,7 @@ end λ = zero(R) r = 0 cm = zero(R) - if k <= w + if k <= Int(ctl[_LT_KEND]) F = _lt_front(fa, ctl, pk) f = Int(ctl[_ST_F]) step = Int(ctl[_LT_STEP]) @@ -1384,7 +1389,8 @@ function _launch_subtrees_ldlt!(N::Numeric{T}, S::Symbolic, nzval, first, count, return nothing end -function _factorize_ldlt_groups!(N::Numeric, S::Symbolic, nzval::AbstractVector, prm, flag) +function _factorize_ldlt_groups!(N::Numeric, S::Symbolic, nzval::AbstractVector, prm, flag, gimpl::Symbol, + nbv::Val, herm::Val) plan = N.plan for k in eachindex(plan.sub_first) _poll_interrupt(flag) @@ -1397,7 +1403,9 @@ function _factorize_ldlt_groups!(N::Numeric, S::Symbolic, nzval::AbstractVector, _poll_interrupt(flag) a, b = plan.group_first[k], plan.group_last[k] W = plan.group_width[k] - if W <= _LT_GLOBAL_MAX_W # regime C and narrow bins: the panel in global memory + if W == 0 || ldlt_blocked_path(S.schedule, S.schedule.group_nodes[a]) # blocked steps and GEMMs (ldlt_c.jl) + _factorize_ldlt_c_group!(N, S, nzval, a, b, plan.group_maxchild[k], prm, gimpl, nbv, herm) + elseif W <= _LT_GLOBAL_MAX_W # narrow regime-B bins: the panel in global memory _launch_front_ldlt!(N, S, nzval, a, b - a + 1, plan.group_maxchild[k], prm, Val(0), Val(LDLT_WORKGROUP)) else # regime B: F₁₁ in local memory _with_width_class(W) do w @@ -1408,6 +1416,16 @@ function _factorize_ldlt_groups!(N::Numeric, S::Symbolic, nzval::AbstractVector, return nothing end +function _factorize_ldlt_herm!(N::Numeric{T}, S::Symbolic, nzval, opts, prm, gimpl, nbv::Val, ::Val{H}) where {T, H} + if H + _factorize_ldlt_groups!(N, S, nzval, prm, opts.user_host_interrupt, gimpl, nbv, Val(true)) + else + _factorize_ldlt_groups!(N, S, nzval, _ldlt_device_params(S, T, opts, Val(false)), opts.user_host_interrupt, + gimpl, nbv, Val(false)) + end + return nothing +end + """ factorize_ldlt!(numeric, symbolic, nzval; opts = Options()) -> 0 @@ -1427,8 +1445,11 @@ before every launch group ([`InterruptedError`](@ref)). With `pivot_epsilon_alg factorization always completes (`info = 0`); the phase allocates nothing on the device, never synchronizes with the host, and is deterministic. """ -function factorize_ldlt!(N::Numeric{T}, S::Symbolic, nzval::AbstractVector; opts::Options = Options()) where {T} +function factorize_ldlt!(N::Numeric{T}, S::Symbolic, nzval::AbstractVector; impl::Symbol = :auto, + opts::Options = Options(), nb::Integer = LDLT_C_NB) where {T} _check_numeric(N, S, nzval) + 1 <= nb <= LDLT_C_NB || throw(InvalidValueError("nb = $nb: the regime-C block size must be in 1:$LDLT_C_NB")) + gimpl = select_impl(:gemm, N.factor, impl === :auto && !S.schedule.vendor_c ? :ka : impl) herm = _ldlt_herm(S, T) prm = _ldlt_device_params(S, T, opts, Val(true)) # validates the options ps = opts.pivot_sign @@ -1440,10 +1461,14 @@ function factorize_ldlt!(N::Numeric{T}, S::Symbolic, nzval::AbstractVector; opts copyto!(N.psign, ps) end prm.scaled && abs_max!(N.aux, nzval, batch_map(N)) - if herm - _factorize_ldlt_groups!(N, S, nzval, prm, opts.user_host_interrupt) - else - _factorize_ldlt_groups!(N, S, nzval, _ldlt_device_params(S, T, opts, Val(false)), opts.user_host_interrupt) + if nb == LDLT_C_NB + if herm + _factorize_ldlt_herm!(N, S, nzval, opts, prm, gimpl, Val(LDLT_C_NB), Val(true)) + else + _factorize_ldlt_herm!(N, S, nzval, opts, prm, gimpl, Val(LDLT_C_NB), Val(false)) + end + else # other block sizes: tests only (dynamic dispatch) + _factorize_ldlt_herm!(N, S, nzval, opts, prm, gimpl, Val(Int(nb)), Val(herm)) end reduce_stats!(N, S) return 0 diff --git a/src/numeric/ldlt_c.jl b/src/numeric/ldlt_c.jl new file mode 100644 index 0000000..12f3105 --- /dev/null +++ b/src/numeric/ldlt_c.jl @@ -0,0 +1,566 @@ +# Regime C of the device LDLᵀ/LDLᴴ (PLAN §2.4, §3.3; experiment 2 of PERFORMANCE.md, issue #75): the +# in-front pivoting of the reference on the fully-summed columns, blocked so that most of the flops go to +# vendor (`:vendor`) or KA (`:ka`) GEMMs through `src/dense/interface.jl`. +# +# The `w` pivot steps of a front run in static blocks of `nb` columns (block b starts at column +# `k0 = (b - 1) nb + 1`). Per block, one `panel_ldlt_kernel!` launch (one workgroup) runs the pivot steps +# of `front_ldlt_kernel!` (the cooperative search, interchanges, D, statistics) on a lazily updated panel: +# the pivots of the block are kept as columns of `Lb` (multipliers) and `Wb` (unscaled pivot columns, +# `W = L D`) in the workspace, and an entry of a column that is not a pivot column yet is read as its +# stored value minus the block's pending pivots, `F[i, j] − Σₜ Lb[i, t] Wb[j, t]ᴴ` (2×2 pivots as one +# term, in pivot order: the arithmetic of the reference's right-looking steps). A pivot column is +# materialized once, when it is taken, and stored unscaled. After the block, one GEMM applies its pivots +# to the trailing fully-summed columns (`F[k0':f, k0':w] −= Lb Wbᴴ`, `k0' = k0 + nb`) and one accumulates +# them into the contribution block (`C −= Lb₂₁ Wb₂₁ᴴ`, `m×m` workspace). `Lb`/`Wb` have `nb + 1` slots: a +# 2×2 pivot that starts at the block's last column takes the first column of the next block too, which +# is then saved before the GEMM (that updates it as a trailing column) and restored by the next block's +# kernel, which starts one column later. Unused slots are zero. A last kernel scales the pivot columns +# into L (`_lt_finalize!`), clears the upper triangle of F₁₁ (the GEMMs write it) and sets the unit +# diagonal; `pack_add!` adds the contribution block to the update stack. Shapes and launches depend on +# the analysis only: no host synchronization. + +"Workgroup size of the regime-C LDLᵀ/LDLᴴ panel kernel." +const LDLT_C_WORKGROUP = 256 + +# workspace offsets (0-based) of Lb, Wb and the saved column of a front +@inline _ltc_offsets(f, m, nb, cb::Bool) = (lo = (cb ? m * m : 0); (lo, lo + f * (nb + 1), lo + 2 * f * (nb + 1))) + +# control words of the panel kernel, after those of the LDLᵀ kernels +const _LTC_K0 = _LT_CTL + 1 # first column of the block +const _LTC_T0 = _LT_CTL + 2 # first slot of a pivot of this block (2 when it starts one column late) +const _LTC_LO = _LT_CTL + 3 # workspace offset (1-based) of Lb +const _LTC_WO = _LT_CTL + 4 # workspace offset (1-based) of Wb +const _LTC_SO = _LT_CTL + 5 # workspace offset (1-based) of the saved column +const _LTC_IDLE = _LT_CTL + 6 # 1: the front has fewer blocks (nothing to do in this launch) +const _LTC_CTL = _LT_CTL + 6 + +# workspace offset (0-based) of the slice of the front at position q of a chunk `nodes[qa:…]`: the slices of +# the fronts before it (`ldlt_c_work_len`, as `ldlt_c_chunk_end` lays them out) +@inline function _ltc_slice(nodes, qa, q, front_nrows, front_ncols, cb_ptr) + base = 0 + @inbounds for q2 in qa:(q - 1) + s2 = nodes[q2] + base += ldlt_c_work_len(front_nrows[s2], front_ncols[s2], cb_ptr[s2] > 0) + end + return base +end + +# front mode of the panel kernel (`Val(_Lazy{H})`) +struct _Lazy{H} end + +# the panel as seen by the pivot search: columns `≥ kcur` with the block's pending pivots (slots t0:t1) +# subtracted; `H`: conj (Hermitian) or not (complex symmetric) +struct _LazyFront{A, B, K, H} + a::A + off::Int + f::Int + lb::B + lo::Int + wo::Int + kind::K + g0::Int # pivot_kind index of slot t: g0 + t + t0::Int + t1::Int + kcur::Int +end + +@inline function _fget(F::_LazyFront{A, B, K, H}, i, j) where {A, B, K, H} + @inbounds begin + x = F.a[F.off + (j - 1) * F.f + i] + j < F.kcur && return x + h = Val(H) + f = F.f + t = F.t0 + while t <= F.t1 + if F.kind[F.g0 + t] == PIVOT_KIND_2X2_FIRST + x -= (F.lb[F.lo + (t - 1) * f + i] * _cj(F.lb[F.wo + (t - 1) * f + j], h) + + F.lb[F.lo + t * f + i] * _cj(F.lb[F.wo + t * f + j], h)) + t += 2 + else + x -= F.lb[F.lo + (t - 1) * f + i] * _cj(F.lb[F.wo + (t - 1) * f + j], h) + t += 1 + end + end + return x + end +end + +@inline function _fset!(F::_LazyFront, i, j, v) + @inbounds F.a[F.off + (j - 1) * F.f + i] = v + return nothing +end + +# the lazy front from the control words: pending slots up to the last pivot taken (column K + STEP - 1), +# or (`cur = false`) without the pivot of the current step +@inline function _ltc_front(fa, ctl, ::Val{H}, cur::Bool = true) where {H} + @inbounds begin + factor, lw, kind = fa + k0 = Int(ctl[_LTC_K0]) + kc = Int(ctl[_LT_K]) + (cur ? Int(ctl[_LT_STEP]) : 0) + return _LazyFront{typeof(factor), typeof(lw), typeof(kind), H}(factor, Int(ctl[_ST_LF]) - 1, Int(ctl[_ST_F]), + lw, Int(ctl[_LTC_LO]) - 1, Int(ctl[_LTC_WO]) - 1, + kind, Int(ctl[_LT_C0]) + k0 - 2, + Int(ctl[_LTC_T0]), kc - k0, kc) + end +end + +@inline _lt_front(fa::Tuple{A, B, K}, ctl, ::Val{_Lazy{H}}) where {A, B, K, H} = _ltc_front(fa, ctl, Val(H)) +@inline _lt_raw_front(fa::Tuple{A, B, K}, ctl, ::Val{_Lazy{H}}) where {A, B, K, H} = + @inbounds _PanelFront(fa[1], Int(ctl[_ST_LF]) - 1, Int(ctl[_ST_F])) + +# block setup (work item 1): the front at position q of the chunk, its control words, the block's columns, +# the workspace offsets of its slice +@inline function _ltc_setup!(ctl, nodes, qa, q, b, nb, super_ptr, front_ptr, front_nrows, front_ncols, cb_ptr, + pivot_kind) + @inbounds begin + IT = eltype(ctl) + s = Int(nodes[q]) + f = Int(front_nrows[s]) + w = Int(front_ncols[s]) + ctl[_ST_NODE] = s % IT + ctl[_ST_F] = f % IT + ctl[_ST_W] = w % IT + ctl[_ST_LF] = front_ptr[s] % IT + _lt_reset!(ctl, super_ptr[s]) + c0 = Int(super_ptr[s]) + k0 = (b - 1) * nb + 1 + idle = k0 > w + late = !idle && k0 > 1 && pivot_kind[c0 + k0 - 2] == PIVOT_KIND_2X2_FIRST # k0 ended the last block + ctl[_LTC_IDLE] = idle % IT + ctl[_LT_K] = (late ? k0 + 1 : k0) % IT + ctl[_LT_KEND] = (idle ? 0 : min(k0 + nb - 1, w)) % IT + ctl[_LTC_K0] = k0 % IT + ctl[_LTC_T0] = (late ? 2 : 1) % IT + lo, wo, so = _ltc_offsets(f, f - w, nb, cb_ptr[s] > 0) .+ _ltc_slice(nodes, qa, q, front_nrows, front_ncols, + cb_ptr) + ctl[_LTC_LO] = (lo + 1) % IT + ctl[_LTC_WO] = (wo + 1) % IT + ctl[_LTC_SO] = (so + 1) % IT + end + return nothing +end + +# zero Lb and Wb, restore the column saved by the last block, identity pivot order (first block) +@inline function _ltc_prepare!(factor, lw, piv, ctl, b, nb, li, ::Val{WG}) where {WG} + @inbounds if ctl[_LTC_IDLE] == 0 + T = eltype(lw) + f = Int(ctl[_ST_F]) + lo = Int(ctl[_LTC_LO]) - 1 + for q in li:WG:(2 * f * (nb + 1)) + lw[lo + q] = zero(T) + end + k0 = Int(ctl[_LTC_K0]) + if Int(ctl[_LTC_T0]) == 2 + so = Int(ctl[_LTC_SO]) - 1 + p0 = Int(ctl[_ST_LF]) - 1 + for i in (k0 - 1 + li):WG:f + factor[p0 + (k0 - 1) * f + i] = lw[so + i] + end + end + b == 1 && _lt_init_piv!(piv, ctl, li, Val(WG)) + end + return nothing +end + +# pass 1 of `_lt_pass1!` on the lazy panel, keeping the materialized column (rows k+1:f) after the saved +# column in the workspace for `_ltc_commit!` +@inline function _ltc_pass1!(fa, lw, ctl, red, redi, li, ::Val{NL}, p::_LDLTDevice{R}, h::Val) where {NL, R} + @inbounds if li <= NL + k = Int(ctl[_LT_K]) + Int(ctl[_LT_STEP]) + w = Int(ctl[_ST_W]) + f = Int(ctl[_ST_F]) + λ = zero(R) + r = 0 + cm = zero(R) + if k <= Int(ctl[_LT_KEND]) # also for 'N': the commit reads the column + F = _ltc_front(fa, ctl, h) + co = Int(ctl[_LTC_SO]) - 1 + f + for i in (k + li):NL:f + x = _fget(F, i, k) + lw[co + i] = x + a = abs(x) + cm = max(cm, a) + if i <= w && a > λ + λ = a + r = i + end + end + end + red[li] = λ + redi[li] = r % eltype(redi) + red[NL + li] = cm + end + return nothing +end + +# rows per task of the fallback scan of the panel kernel +const _LTC_CHUNK = 64 + +# the fallback scan (`_lt_pass3!`) of the panel kernel in two steps over the whole workgroup. Flags (after the +# pass-1 column in the workspace; 1: rejected) are cleared for the candidate columns k:w, then every +# (candidate column, chunk of `_LTC_CHUNK` rows) task tests its entries with the per-entry form of the +# threshold test (`d ≥ u |F[i, j]|`, see `_lt_threshold_ok`; a NaN fails it) and flags the column on a +# violation; a column whose diagonal is not acceptable (`|a_jj| < ε` or NaN) is flagged too. A whole column +# per work item (as `_lt_pass3!`) is far slower here: every entry reads the block's pending pivots. +@inline function _ltc_flags_clear!(lw, ctl, li, ::Val{WG}) where {WG} + @inbounds begin + f = Int(ctl[_ST_F]) + fo = Int(ctl[_LTC_SO]) - 1 + 2 * f + for j in (Int(ctl[_LT_K]) + li - 1):WG:Int(ctl[_ST_W]) + lw[fo + j] = zero(eltype(lw)) + end + end + return nothing +end + +@inline function _ltc_scan!(fa, lw, ctl, li, ::Val{WG}, aux, p::_LDLTDevice, h::Val) where {WG} + @inbounds begin + k = Int(ctl[_LT_K]) + w = Int(ctl[_ST_W]) + f = Int(ctl[_ST_F]) + fo = Int(ctl[_LTC_SO]) - 1 + 2 * f + ε = _lt_eps(p, aux) + u = p.u + F = _ltc_front(fa, ctl, h) + nc = w - k + 1 + nch = cld(f - k + 1, _LTC_CHUNK) + for q in (li - 1):WG:(nc * nch - 1) + j = k + q % nc + iszero(lw[fo + j]) || continue # already rejected + i0 = k + (q ÷ nc) * _LTC_CHUNK + a = _fabs(F, j, j) + ok = a >= ε + if ok + for i in i0:min(i0 + _LTC_CHUNK - 1, f) + i == j && continue + if !(a >= u * _fabs(F, i, j)) + ok = false + break + end + end + end + ok || (lw[fo + j] = one(eltype(lw))) + end + end + return nothing +end + +# the lane partials of pass 3 from the flags: the unflagged column with the largest |a_jj| (first on a tie) +# and, for pivot type 'D', the largest |a_jj| +@inline function _ltc_pick!(fa, lw, ctl, red, redi, li, ::Val{NL}, p::_LDLTDevice{R}, h::Val) where {NL, R} + @inbounds if li <= NL + k = Int(ctl[_LT_K]) + f = Int(ctl[_ST_F]) + fo = Int(ctl[_LTC_SO]) - 1 + 2 * f + F = _ltc_front(fa, ctl, h) + bv = zero(R) + bi = 0 + gv = zero(R) + gi = 0 + for j in (k + li - 1):NL:Int(ctl[_ST_W]) + a = _fabs(F, j, j) + if p.ptype == _LT_PIVOT_DIAGONAL && !isnan(a) && (gi == 0 || a > gv) + gv = a + gi = j + end + if iszero(lw[fo + j]) && (bi == 0 || a > bv) + bv = a + bi = j + end + end + IT = eltype(redi) + red[li] = bv + redi[li] = bi % IT + red[NL + li] = gv + redi[NL + li] = gi % IT + end + return nothing +end + +# swap rows p and q of Lb and Wb along with the interchange of the step (`_lt_swap!` moves the stored panel) +@inline function _ltc_swap_rows!(lw, ctl, nb, li, ::Val{WG}) where {WG} + @inbounds begin + k = Int(ctl[_LT_K]) + if k <= Int(ctl[_LT_KEND]) + r = Int(ctl[_LT_R]) + p = r == 0 ? k : k + 1 + q = r == 0 ? Int(ctl[_LT_C]) : r + if p != q + f = Int(ctl[_ST_F]) + for e in (li - 1):WG:(2 * (nb + 1) - 1) + base = Int(ctl[_LTC_LO]) - 1 + e * f # Lb and Wb are adjacent: 2(nb + 1) columns + x = lw[base + p] + lw[base + p] = lw[base + q] + lw[base + q] = x + end + end + end + end + return nothing +end + +# materialize the pivot column(s) of the step below the pivot block: the stored panel gets the unscaled +# values, Wb the same, Lb the multipliers (as `_lt_update!` forms them) +@inline function _ltc_commit!(fa, lw, ctl, pv, li, ::Val{WG}, h::Val{H}) where {WG, H} + @inbounds begin + k = Int(ctl[_LT_K]) + if k <= Int(ctl[_LT_KEND]) + F = _ltc_front(fa, ctl, h, false) + f = Int(ctl[_ST_F]) + step = Int(ctl[_LT_STEP]) + t = k - Int(ctl[_LTC_K0]) + 1 + lo = Int(ctl[_LTC_LO]) - 1 + wo = Int(ctl[_LTC_WO]) - 1 + # pass 1 left column k materialized when the pivot is the 1×1 one at k (no interchange) + kept = step == 1 && Int(ctl[_LT_C]) == k + co = Int(ctl[_LTC_SO]) - 1 + f + for i in (k + step - 1 + li):WG:f + if step == 1 + x = kept ? lw[co + i] : _fget(F, i, k) + _fset!(F, i, k, x) + lw[wo + (t - 1) * f + i] = x + lw[lo + (t - 1) * f + i] = x / pv[1] + else + x1 = _fget(F, i, k) + x2 = _fget(F, i, k + 1) + _fset!(F, i, k, x1) + _fset!(F, i, k + 1, x2) + lw[wo + (t - 1) * f + i] = x1 + lw[wo + t * f + i] = x2 + lw[lo + (t - 1) * f + i] = x1 * pv[1] + x2 * pv[3] + lw[lo + t * f + i] = x1 * pv[2] + x2 * pv[4] + end + end + end + end + return nothing +end + +# after the block: save the first column of the next block when a 2×2 pivot took it; statistics +@inline function _ltc_block_end!(factor, lw, stats, ctl, pivot_kind, b, li, ::Val{WG}) where {WG} + @inbounds if ctl[_LTC_IDLE] == 0 + kend = Int(ctl[_LT_KEND]) + f = Int(ctl[_ST_F]) + c0 = Int(ctl[_LT_C0]) + if kend < Int(ctl[_ST_W]) && pivot_kind[c0 + kend - 1] == PIVOT_KIND_2X2_FIRST + so = Int(ctl[_LTC_SO]) - 1 + p0 = Int(ctl[_ST_LF]) - 1 + for i in (kend + li):WG:f + lw[so + i] = factor[p0 + kend * f + i] + end + end + if li == 1 + base = (Int(ctl[_ST_NODE]) - 1) * FRONT_STATS_FIELDS + for q in 1:5 + x = Int64(ctl[_LT_STAT + q]) + stats[base + q] = b == 1 ? x : stats[base + q] + x + end + stats[base + STAT_INFO] = 0 + end + end + return nothing +end + +""" + panel_ldlt_kernel!(backend, WG)(factor, lw, d, piv, pivot_kind, psign, perm, aux, stats, nodes, qa, b, nb, + super_ptr, front_ptr, k, nbatch, front_nrows, front_ncols, cb_ptr, prm, Val(WG); + ndrange = WG * count) + +Block `b` (columns `(b - 1) nb + 1 : min(b nb, w)`) of the regime-C LDLᵀ/LDLᴴ +pivot steps of the `count` fronts `nodes[qa:(qa + count - 1)]`, one workgroup +each (panel of batch member `k` of `nbatch` at `member_panels(front_ptr, k, nbatch)[s]`, workspace slice in `lw` after those of the +fronts before it: `Lb`, `Wb`, the saved column, the pass-1 column): the +reference's pivot sequence on the lazily updated panel (see the top of +`src/numeric/ldlt_c.jl`), D, `piv`, the pivot kinds and the front's statistics +(added to those of the earlier blocks). Fronts with fewer blocks do nothing. +""" +@kernel function panel_ldlt_kernel!(factor, lw, d, piv, pivot_kind, psign, perm, aux, stats, nodes, qa, b, nb, + super_ptr, front_ptr, mk, nbatch, front_nrows, front_ncols, cb_ptr, + prm::_LDLTDevice{R, HERM}, ::Val{WG}) where {R, HERM, WG} + @uniform TT = eltype(factor) + @uniform IT = eltype(front_ncols) + li = @index(Local, Linear) + G = @index(Group, Linear) + ctl = @localmem IT (_LTC_CTL,) + pv = @localmem TT (_LT_NPV,) + red = @localmem R (3 * WG,) + redi = @localmem IT (2 * WG,) + if li == 1 + _ltc_setup!(ctl, nodes, qa, qa + G - 1, b, nb, super_ptr, member_panels(front_ptr, mk, nbatch), front_nrows, + front_ncols, cb_ptr, pivot_kind) + end + @synchronize + _ltc_prepare!(factor, lw, piv, ctl, b, nb, li, Val(WG)) + @synchronize + _ltc_pass1!((factor, lw, pivot_kind), lw, ctl, red, redi, li, Val(WG), prm, Val(HERM)) + @synchronize + for it in 1:nb + _lt_stage1!(red, redi, ctl, li, Val(WG), Val(1)) + @synchronize + if li == 1 + _lt_decide1!((factor, lw, pivot_kind), ctl, pv, red, redi, d, pivot_kind, piv, psign, perm, aux, prm, + Val(WG), Val(_Lazy{HERM}), Val(HERM)) + end + @synchronize + if ctl[_LT_PHASE] == 2 + _lt_pass2!((factor, lw, pivot_kind), ctl, red, redi, li, Val(WG), prm, Val(_Lazy{HERM})) + @synchronize + _lt_stage1!(red, redi, ctl, li, Val(WG), Val(2)) + @synchronize + if li == 1 + _lt_decide2!((factor, lw, pivot_kind), ctl, pv, red, redi, d, pivot_kind, piv, psign, perm, aux, prm, + Val(WG), Val(_Lazy{HERM}), Val(HERM)) + end + @synchronize + end + if ctl[_LT_PHASE] == 3 + _ltc_flags_clear!(lw, ctl, li, Val(WG)) + @synchronize + _ltc_scan!((factor, lw, pivot_kind), lw, ctl, li, Val(WG), aux, prm, Val(HERM)) + @synchronize + _ltc_pick!((factor, lw, pivot_kind), lw, ctl, red, redi, li, Val(WG), prm, Val(HERM)) + @synchronize + _lt_stage1!(red, redi, ctl, li, Val(WG), Val(3)) + @synchronize + if li == 1 + _lt_decide3!((factor, lw, pivot_kind), ctl, pv, red, redi, d, pivot_kind, piv, psign, perm, aux, prm, + Val(WG), Val(_Lazy{HERM}), Val(HERM)) + end + @synchronize + end + _lt_swap!((factor, lw, pivot_kind), ctl, pv, piv, li, Val(WG), Val(_Lazy{HERM}), Val(HERM)) + _ltc_swap_rows!(lw, ctl, nb, li, Val(WG)) + @synchronize + _ltc_commit!((factor, lw, pivot_kind), lw, ctl, pv, li, Val(WG), Val(HERM)) + @synchronize + _ltc_pass1!((factor, lw, pivot_kind), lw, ctl, red, redi, li, Val(WG), prm, Val(HERM)) + @synchronize + end + _ltc_block_end!(factor, lw, stats, ctl, pivot_kind, b, li, Val(WG)) +end + +""" + finish_ldlt_kernel!(backend, WG)(factor, d, pivot_kind, info, nodes, qa, super_ptr, front_ptr, k, nbatch, + front_nrows, front_ncols, Val(HERM), Val(WG); ndrange = WG * count) + +After the blocks of the regime-C fronts `nodes[qa:(qa + count - 1)]`, one +workgroup each: scale the pivot columns into L (`_lt_finalize!`), clear the +upper triangle of `F₁₁`, set the unit diagonal and the front's status. +""" +@kernel function finish_ldlt_kernel!(factor, d, pivot_kind, info, nodes, qa, super_ptr, front_ptr, mk, nbatch, + front_nrows, front_ncols, ::Val{HERM}, ::Val{WG}) where {HERM, WG} + @uniform IT = eltype(front_ncols) + li = @index(Local, Linear) + G = @index(Group, Linear) + ctl = @localmem IT (_LT_CTL,) + if li == 1 + @inbounds begin + s = nodes[qa + G - 1] + ctl[_ST_NODE] = s % IT + ctl[_ST_F] = front_nrows[s] % IT + ctl[_ST_W] = front_ncols[s] % IT + ctl[_ST_LF] = member_panels(front_ptr, mk, nbatch)[s] % IT + ctl[_LT_C0] = super_ptr[s] % IT + end + end + @synchronize + _lt_finalize!(factor, ctl, d, pivot_kind, li, Val(WG), Val(false), Val(HERM)) + @inbounds begin + f = Int(ctl[_ST_F]) + w = Int(ctl[_ST_W]) + p = Int(ctl[_ST_LF]) - 1 + for q in (li - 1):WG:(w * w - 1) + j = q ÷ w + 1 + i = q - (j - 1) * w + 1 + i <= j && (factor[p + (j - 1) * f + i] = i == j ? one(eltype(factor)) : zero(eltype(factor))) + end + li == 1 && (info[ctl[_ST_NODE]] = Int32(0)) + end +end + +# regime C of the fronts `nodes[qa:qe]` of batch member `k`, concurrently: per block one panel launch, then the +# GEMMs of every front on its trailing columns and its contribution block; the last scaling; the contribution +# blocks onto the update stack +function _factor_chunk_ldlt_c!(N::Numeric{T}, S::Symbolic, qa::Int, qe::Int, k::Int, prm, gimpl::Symbol, + ::Val{NB}, ::Val{HERM}) where {T, NB, HERM} + L, sc = S.layout, S.schedule + nodes = sc.group_nodes + nb = N.nbatch + count = qe - qa + 1 + lw = view(N.work, 1:L.work_len) + tB = HERM ? 'C' : 'T' + d, piv, kind = _mview(N.d, k, nb), _mview(N.piv, k, nb), _mview(N.pivot_kind, k, nb) + backend = KernelAbstractions.get_backend(N.factor) + panel! = panel_ldlt_kernel!(backend, LDLT_C_WORKGROUP) + nblocks = 0 + for q in qa:qe + nblocks = max(nblocks, cld(sc.width[nodes[q]], NB)) + end + for b in 1:nblocks + panel!(N.factor, lw, d, piv, kind, N.psign, S.perm, _mview(N.aux, k, nb), _mview(N.stats, k, nb), + S.group_nodes, qa, b, NB, S.super_ptr, S.front_ptr, k, nb, S.front_nrows, S.front_ncols, S.cb_ptr, prm, + Val(LDLT_C_WORKGROUP); ndrange = LDLT_C_WORKGROUP * count) + base = 0 + for q in qa:qe + s = nodes[q] + f, w = sc.rows[s], sc.width[s] + m = f - w + cb = m > 0 && L.cb_ptr[s] > 0 + if b <= cld(w, NB) + p0 = panel_offset(L.panel_ptr, s, k, nb) + P = reshape(view(N.factor, p0:(p0 + f * w - 1)), f, w) + lo, wo, _ = _ltc_offsets(f, m, NB, cb) + Lb = reshape(view(lw, (base + lo + 1):(base + lo + f * (NB + 1))), f, NB + 1) + Wb = reshape(view(lw, (base + wo + 1):(base + wo + f * (NB + 1))), f, NB + 1) + k1 = b * NB + 1 # first column of the next block + if k1 <= w + _gemm_impl!(gimpl, 'N', tB, -one(T), view(Lb, k1:f, :), view(Wb, k1:w, :), one(T), + view(P, k1:f, k1:w)) + end + if cb + C = reshape(view(lw, (base + 1):(base + m * m)), m, m) + _gemm_impl!(gimpl, 'N', tB, -one(T), view(Lb, (w + 1):f, :), view(Wb, (w + 1):f, :), + b == 1 ? zero(T) : one(T), C) + end + end + base += ldlt_c_work_len(f, w, cb) + end + end + finish_ldlt_kernel!(backend, LDLT_C_WORKGROUP)(N.factor, d, kind, _iview(N.info, k, nb), S.group_nodes, qa, + S.super_ptr, S.front_ptr, k, nb, S.front_nrows, S.front_ncols, + Val(HERM), + Val(LDLT_C_WORKGROUP); ndrange = LDLT_C_WORKGROUP * count) + base = 0 + for q in qa:qe + s = nodes[q] + f, w = sc.rows[s], sc.width[s] + m = f - w + cb = m > 0 && L.cb_ptr[s] > 0 + cb && pack_add!(_mview(N.stack, k, nb), L.cb_ptr[s], view(lw, (base + 1):(base + m * m)), m) + base += ldlt_c_work_len(f, w, cb) + end + return nothing +end + +# a regime-C launch group of the LDLᵀ/LDLᴴ phase: assembly launches, then its fronts in chunks +# (`ldlt_c_chunk_end`) for every active member +function _factorize_ldlt_c_group!(N::Numeric, S::Symbolic, nzval, a::Int, b::Int, maxchild::Int, prm, gimpl::Symbol, + nbv::Val, herm::Val) + zero_fronts!(N, S, a, b - a + 1) + scatter_A!(N, S, nzval, a, b - a + 1) + extend_add!(N, S, a, b - a + 1, maxchild) + plan = N.plan + q = a + while q <= b + e = ldlt_c_chunk_end(S.schedule, S.layout.cb_len, q, b, S.layout.work_len) + if N.nbatch == 1 + _factor_chunk_ldlt_c!(N, S, q, e, 1, prm, gimpl, nbv, herm) + else + for j in 1:plan.nact[] + _factor_chunk_ldlt_c!(N, S, q, e, Int(plan.members_host[j]), prm, gimpl, nbv, herm) + end + end + q = e + 1 + end + return nothing +end diff --git a/src/solver.jl b/src/solver.jl index 3635f79..38cd5a7 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -377,7 +377,7 @@ function _symbolic!(solver::DirectSolver{T, INT}) where {T, INT} view = _stored_view(solver), index = A.index) sp = supernode_partition(factor_pattern(P, ord), ord.perm, opts) sc = build_schedule(sp, opts, T; reserve = subtree_local_reserve(solver.structure)) - layout = build_layout(sp, sc) + layout = build_layout(sp, sc; ldlt = _is_ldlt_structure(solver.structure)) Sh = Symbolic(sp, sc, layout, solver.host_rowptr, solver.host_colval, A.nrows, solver.structure; view = _stored_view(solver), index = A.index) nrhs = solver.workspace === nothing ? solver.nbatch : max_rhs(solver.workspace) diff --git a/src/symbolic/layout.jl b/src/symbolic/layout.jl index 70b1a85..325550c 100644 --- a/src/symbolic/layout.jl +++ b/src/symbolic/layout.jl @@ -44,7 +44,10 @@ offsets, 1-based: * `step_top[t + 1]`: last update-stack entry in use during step `t` (`t = 0:nsteps`); `stack_len = maximum(step_top)` is the high-water mark; * `work_len`: entries of the regime-C `syrk` workspace (largest `m^2` of a - front on the regime-C path, [`takes_c_path`](@ref), with a block on the stack); + front on the regime-C path, [`takes_c_path`](@ref), with a block on the stack; + LDLᵀ/LDLᴴ: the largest sum of [`ldlt_c_work_len`](@ref) over a chunk of + concurrent fronts on the blocked path, [`ldlt_blocked_path`](@ref) and + [`ldlt_c_chunk_end`](@ref)); * `local_front`, `local_cb` (regime A, `0` elsewhere): local-memory offset of the packed front of `s` and of its contribution block after the move (`0` for a subtree root, whose block goes to the update stack); `local_len[t]`: entries @@ -164,8 +167,70 @@ end _high_water(cb_ptr, cb_len, ids) = maximum((cb_ptr[s] + cb_len[s] - 1 for s in ids); init = 0) +"Pivot columns per block of the regime-C LDLᵀ/LDLᴴ path (`panel_ldlt_kernel!`)." +const LDLT_C_NB = 32 + +""" + ldlt_c_work_len(f, w, cb::Bool, nb = LDLT_C_NB) -> Int + +Workspace entries of a regime-C LDLᵀ/LDLᴴ front with `f` rows and `w` +fully-summed columns: the `m×m` contribution block (`m = f - w`, when `cb`), +`Lb` and `Wb` (`f × (nb + 1)` each), the saved column, the next pivot +column of pass 1 and the rejection flags of the fallback scan (`f` each). +""" +ldlt_c_work_len(f::Integer, w::Integer, cb::Bool, nb::Integer = LDLT_C_NB) = + (cb ? (Int(f) - Int(w))^2 : 0) + 2 * Int(f) * (nb + 1) + 3 * Int(f) + +""" +Smallest width class and row class of a regime-B bin that the LDLᵀ/LDLᴴ factorization runs on the blocked +regime-C path ([`ldlt_blocked_path`](@ref)) instead of the fused front kernel. """ - build_layout(sp::SupernodePartition, schedule::Schedule) -> Layout +const LDLT_BLOCKED_MIN_WCLASS = 32 +const LDLT_BLOCKED_MIN_FCLASS = 256 + +""" + ldlt_blocked_path(sc, s) -> Bool + +Front `s` takes the blocked LDLᵀ/LDLᴴ path (`panel_ldlt_kernel!` and GEMMs, +`src/numeric/ldlt_c.jl`): it is on the regime-C path ([`takes_c_path`](@ref)), +or in a regime-B bin of width class `≥ LDLT_BLOCKED_MIN_WCLASS` and row class +`≥ LDLT_BLOCKED_MIN_FCLASS`, where the GEMMs beat the fused kernel's in-kernel +updates. +""" +function ldlt_blocked_path(sc, s::Integer) + takes_c_path(sc, s) && return true + sc.regime[s] == REGIME_B || return false + nf = length(sc.fclasses) + return sc.wclasses[(sc.bin[s] - 1) ÷ nf + 1] >= LDLT_BLOCKED_MIN_WCLASS && + sc.fclasses[(sc.bin[s] - 1) % nf + 1] >= LDLT_BLOCKED_MIN_FCLASS +end + +""" + ldlt_c_chunk_end(sc, cb_len, q, last, cap) -> Int + +The fronts `sc.group_nodes[q:e]` of a regime-C launch group (`e ≤ last`) that +the blocked LDLᵀ/LDLᴴ path factors concurrently, one workgroup each: the +longest run from `q` whose [`ldlt_c_work_len`](@ref) slices, laid out one after +another, fit `cap` entries (at least one front). `cb_len[s] > 0` when front `s` +has a contribution block on the update stack. +""" +function ldlt_c_chunk_end(sc, cb_len::AbstractVector, q::Integer, last::Integer, cap::Integer) + nodes = sc.group_nodes + s = nodes[q] + acc = ldlt_c_work_len(sc.rows[s], sc.width[s], cb_len[s] > 0) + e = Int(q) + while e < last + s = nodes[e + 1] + len = ldlt_c_work_len(sc.rows[s], sc.width[s], cb_len[s] > 0) + acc + len <= cap || break + acc += len + e += 1 + end + return e +end + +""" + build_layout(sp::SupernodePartition, schedule::Schedule; ldlt = false) -> Layout Panel offsets in supernode order, D offsets, and the update-stack offsets of the contribution blocks that leave their front through global memory (B/C @@ -174,10 +239,11 @@ Blocks whose lifetimes `[step(s), step(parent)]` overlap never share entries; the offsets are the placement with the lowest high-water mark among a step-by-step first fit and two offline placements (lowest free offset, largest blocks first and largest size × lifetime first; issue #48). Also the regime-C -workspace and the local-memory offsets of the regime-A subtrees -([`subtree_local_layout`](@ref)). +workspace (with `ldlt`, structures `"S"`/`"H"`, the larger one of the blocked +LDLᵀ/LDLᴴ path, [`ldlt_c_work_len`](@ref)) and the local-memory offsets of the +regime-A subtrees ([`subtree_local_layout`](@ref)). """ -function build_layout(sp::SupernodePartition, sc::Schedule) +function build_layout(sp::SupernodePartition, sc::Schedule; ldlt::Bool = false) ns = nsupernodes(sp) n = sp.n panel_ptr = Vector{Int}(undef, ns + 1) @@ -220,8 +286,25 @@ function build_layout(sp::SupernodePartition, sc::Schedule) for s in ids, t in cb_first[s]:cb_last[s] step_top[t + 1] = max(step_top[t + 1], cb_ptr[s] + cb_len[s] - 1) end - work_len = maximum((cb_len[s] > 0 && takes_c_path(sc, s) ? (sc.rows[s] - sc.width[s])^2 : 0 for s in 1:ns); - init = 0) + work_len = if ldlt + # the fronts of a regime-C group run concurrently in chunks of at most twice the largest workspace + single = maximum((ldlt_blocked_path(sc, s) ? ldlt_c_work_len(sc.rows[s], sc.width[s], cb_len[s] > 0) : 0 + for s in 1:ns); init = 0) + len = single + for g in sc.groups + (g.regime == REGIME_A || !ldlt_blocked_path(sc, sc.group_nodes[g.first])) && continue + q = g.first + while q <= g.last + e = ldlt_c_chunk_end(sc, cb_len, q, g.last, 2 * single) + len = max(len, sum(s -> ldlt_c_work_len(sc.rows[s], sc.width[s], cb_len[s] > 0), + view(sc.group_nodes, q:e))) + q = e + 1 + end + end + len + else + maximum((cb_len[s] > 0 && takes_c_path(sc, s) ? (sc.rows[s] - sc.width[s])^2 : 0 for s in 1:ns); init = 0) + end local_front, local_cb, local_len = subtree_local_layout(sp, sc) return Layout(panel_ptr, panel_ptr[end] - 1, d_ptr, 2n, cb_ptr, cb_len, cb_first, cb_last, step_top, maximum(step_top; init = 0), work_len, local_front, local_cb, local_len) diff --git a/src/symbolic/maps.jl b/src/symbolic/maps.jl index 6a68282..c79e71f 100644 --- a/src/symbolic/maps.jl +++ b/src/symbolic/maps.jl @@ -246,7 +246,7 @@ function symbolic_analysis(A::CSR, structure, view = VIEW_FULL; opts::Options = ord = compute_ordering(P, opts; T, pp.pairs, pp.candidates) sp = supernode_partition(factor_pattern(P, ord), ord.perm, opts) sc = build_schedule(sp, opts, T; reserve = subtree_local_reserve(structure)) - layout = build_layout(sp, sc) + layout = build_layout(sp, sc; ldlt = _is_ldlt_structure(_structure(structure))) return Symbolic(sp, sc, layout, rowptr, colval, A.nrows, structure; view, index = A.index) end diff --git a/test/test_numeric_ldlt.jl b/test/test_numeric_ldlt.jl index 02cfd1a..40bd1b4 100644 --- a/test/test_numeric_ldlt.jl +++ b/test/test_numeric_ldlt.jl @@ -52,8 +52,10 @@ end for nrhs in (1, 5), det in (false, true) b = nrhs == 1 ? rand(T, size(A, 1)) : rand(T, size(A, 1), nrhs) x = device_solve(backend, ws, Sd, Nd, b; deterministic = det) - if name == "kkt(300,100,1e-8)" && T == ComplexF32 - # T14 / issue #66: this draw needs one refinement step in ComplexF32 (the reference too) + if name == "kkt(300,100,1e-8)" && T in (Float32, ComplexF32) + # T14 / issue #66: this draw needs one refinement step in ComplexF32 (the reference too); in + # Float32 the reference factor itself exceeds tol(T) on 1 of 20 right-hand sides (5.3e-4, + # #75), and the regime-C GEMMs (rounding only) hit such a draw x = x + device_solve(backend, ws, Sd, Nd, b - A * x; deterministic = det) end @test relres(A, x, b) <= tol(T) @@ -135,6 +137,54 @@ end InvalidValueError end +@testset "regime C: blocked pivot steps and GEMMs ($(backend_name(backend)), $T)" for backend in BACKENDS, T in ELTYPES + # every front on the regime-C path (`panel_ldlt_kernel!` per block of nb columns, GEMMs on the trailing + # columns and the contribution block); block sizes below the root's width, so that 2×2 pivots straddle + # block boundaries; every dense implementation of the GEMMs + Random.seed!(666) + conly = (factorization_alg = "algo2", subtree_budgets = Int[]) + cases = (("kkt 1e-3 H, δ = 0, interleaved (2×2)", + kkt_matrix(T, 200, 100, 0.0; hessian = :indefinite, hessian_scale = 1.0e-3), + (user_perm = kkt_interleaved_perm(200, 100),)), + ("random_symindef(200,0.5)", random_symindef(T, 200, 0.5), (;))) + for (name, A, kw) in cases + opts = Options(; kw..., conly...) + S, Nr, Sd, Nd, nz = ldlt_setup(backend, A; opts) + sp = S.partition + roots = findall(==(0), sp.snparent) + s = roots[argmax([SDS.snwidth(sp, r) for r in roots])] + w = SDS.snwidth(sp, s) + # a small block size at which a 2×2 pivot starts at the last column of a block (every front is on + # the regime-C path) + straddle(nb) = any(1:SDS.nsupernodes(sp)) do v + c0, wv = sp.super_ptr[v], SDS.snwidth(sp, v) + any(k -> k % nb == 0 && Nr.pivot_kind[c0 + k - 1] == SDS.PIVOT_KIND_2X2_FIRST, 1:(wv - 1)) + end + small = findfirst(straddle, 2:16) + name == first(cases[1]) && @test SDS.pivot_stats(Nr).n2x2 > 0 && small !== nothing + for nb in (something(small, 2) + 1, SDS.LDLT_C_NB) + @test SDS.takes_c_path(S.schedule, s) && w > nb + for impl in SDS.dense_impls(:gemm, backend, T) + @test SDS.factorize_ldlt!(Nd, Sd, nz; impl, opts, nb) == 0 + Nh = SDS.host_numeric(Nd) + @test Nh.piv == Nr.piv + @test Nh.pivot_kind == Nr.pivot_kind + @test d_error(Nh, Nr) <= panel_tol(T) + @test panel_error(Nh, Nr) <= panel_tol(T) + @test Nh.stats == Nr.stats + @test SDS.pivot_totals(Nd) == SDS.pivot_stats(Nr) + ws = SDS.allocate_solve(Sd, T, backend, 1) + b = rand(T, size(A, 1)) + @test relres(A, device_solve(backend, ws, Sd, Nd, b; deterministic = true), b) <= tol(T) + end + end + end + # the block size is checked + S, Nr, Sd, Nd, nz = ldlt_setup(backend, cases[2][2]; opts = Options(; conly...)) + @test thrown(() -> SDS.factorize_ldlt!(Nd, Sd, nz; nb = 0)) isa InvalidValueError + @test thrown(() -> SDS.factorize_ldlt!(Nd, Sd, nz; nb = SDS.LDLT_C_NB + 1)) isa InvalidValueError +end + @testset "pivot_type 'D' and 'N' on quasi-definite KKT ($(backend_name(backend)), $T)" for backend in BACKENDS, T in ELTYPES nh, nj = 300, 100