perf: regime-C LDLᵀ through blocked pivot steps and GEMMs (#75)

[michel2323]

Refs #75

Stacked on #88 (and #87); merge those first. The step-3 commits are the top four of this branch. Not "Closes #75": two of the criteria are not met (below), #75 stays open for them.

Task: #75, experiment 2 of PERFORMANCE.md, step 3 of 3 — regime-C LDLᵀ through vendor BLAS with the reference pivot sequence (T25 owner note, item 2).

What was built:

• src/numeric/ldlt_c.jl (new): every front on the regime-C path, and every regime-B bin of width class ≥ 32 and row class ≥ 256 (ldlt_blocked_path, measured), runs its pivot steps in static blocks of LDLT_C_NB = 32 columns. panel_ldlt_kernel! (one workgroup per front; the fronts of a launch group concurrently, in chunks whose workspace fits twice the largest front's, ldlt_c_chunk_end) runs the pivot search, interchanges and D of front_ldlt_kernel! on a lazily updated panel: the block's pivots are kept as Lb (multipliers) and Wb (unscaled pivot columns) in the workspace, and an entry of a non-pivot column is read as its stored value minus the pending pivots in pivot order (2×2 pivots as one term: the reference's right-looking arithmetic). After each block one GEMM updates the trailing fully-summed columns and one accumulates the contribution block, through src/dense/interface.jl only (_gemm_impl!: :vendor on CUDA, :generic/:ka elsewhere; impl is passed from factorize!); finish_ldlt_kernel! scales into L and clears F₁₁'s upper triangle; pack_add! puts the block on the update stack (packed lower triangle, as the Cholesky regime-C path). A 2×2 pivot at a block's last column takes the next block's first column, saved across the GEMM and restored by the next block. Pass 1 keeps the materialized column for the commit; the fallback scan runs as (candidate column, 64 rows) tasks over the workgroup. Static launch shapes, no host synchronization, no atomics (one benign duplicate write of the same flag value).
• src/symbolic/layout.jl: ldlt_c_work_len, ldlt_c_chunk_end, ldlt_blocked_path, build_layout(…; ldlt) sizes the LDLᵀ workspace (lap3d_40 101 MB vs 62 MB Cholesky, apache2 223 vs 109 MB).
• factorize_ldlt!(…; impl, nb) (nb ≤ LDLT_C_NB for the tests, InvalidValueError otherwise).
• Tests (test/test_numeric_ldlt.jl): new testset "regime C: blocked pivot steps and GEMMs" for every backend and T ∈ ELTYPES: "C only" analysis, a root front wider than the block, a block size chosen from the reference's pivots so that 2×2 pivots straddle blocks, and LDLT_C_NB, every dense implementation of the GEMMs: piv, pivot_kind, stats, totals equal to ref_ldlt!, D and panels within panel_tol(T), relres ≤ tol(T); nb out of range raises.
• PERFORMANCE.md (step-3 results, criteria, tracked issues), bench/comparison/comparison.{md,png} and sds.csv regenerated, TASKS.md one line under the T25 owner note, bench/profile_phases.jl profiles lap3d_40 LDLᵀ again.

Tests: julia --project=. -e 'using Pkg; Pkg.test()' with CUDA in the local test environment (RTX 4080): 105738 pass / 0 fail / 2 broken (the T16 @test_broken, once per backend); after the rebase on T18, test_numeric_ldlt, test_symbolic_schedule, test_refinement on CPU and CUDA: 16871 pass, 0 fail, 2 broken. factorize! LDLᵀ allocates 64 B on the CPU backend (unchanged).
Test adaptation, please review: in "device factor equals the reference", the existing T14/#66 exception for kkt(300,100,1e-8) (one refinement step before the relres check, ComplexF32 only) now covers Float32 too: the GEMM rounding moved one right-hand side to relres 4.04e-4 > tol 3.45e-4 in "C only". Measured: on 20 draws the reference factor itself reaches 5.3e-4 (median 1.3e-4) and main's device factor 7.6e-4, so the bound depends on the draw for this matrix; the comment in the test says so. No other assertion changed.

Measurements (RTX 4080; refactorization best of 5 warm runs, bench/profile_phases.jl; *: bench/compare.jl median):

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 (lap3d_40 and apache2 calls replayed through _gemm_impl!(:vendor, …)): 0.41 TFLOP/s = 56–57% of cuBLAS DGEMM (0.73 TFLOP/s on 4096²).
• K2 dumps (CUDA): inertia and nperturbed identical to main on all nine.
• kkt_matrix(Float64, 3000, 1000, 1e-8), default and "C only": CUDA 33.9 s (main) → 2.5 s; KA CPU backend 7.0 → 3.3–3.4 s; ref_ldlt! 2.5–2.7 s.
• comparison.md overview, LDLᵀ static pivoting vs cuDSS: factorization 6.37× → 4.41×, refactorization 9.90× → 6.27×. The K2 + refinement row (5.45× → 7.19×) is dominated by the case1354 K2 dumps, whose harness medians vary 2–3× between reruns of the same code (8.8–27.6 ms for main; profile 7.2–9.1 → 6.3–7.0 ms).

Criteria: experiment 1 (LDLᵀ within 1.5× of SDS Cholesky): met on lap2d_300 (1.31×) and case14 condensed_1, not on bcsstk17/38 (1.73×), lap3d_40 (1.89×), apache2 (2.06×), other KKT dumps (1.5–2.6×); inertia unchanged on all K2 dumps: met. Experiment 2: ≥ 50% of DGEMM peak on fronts > 512: met (56%); "C only" faster than the reference on the CPU backend (0.8×) and ≥ 5× on the GPU (1.1×): not met on the T15 matrix, whose pivot steps end in a fallback scan half the time with no acceptable column (every candidate must be rejected exactly; each lazily read entry subtracts up to 32 pending pivots; nb = 8 gives 0.86 s; a local-memory tiled scan was slower, 4.6 s). Details and what remains in PERFORMANCE.md.

Deviations from PLAN.md / TASKS.md: regime C does not use sytrf (owner decision on #75); regime-B bins of width class ≥ 32 and row class ≥ 256 also take the blocked path (measured: bcsstk38 47 → 39 ms, lap2d 53 → 41 ms); the LDLᵀ workspace is larger (above). Suggested plan note (PLAN §2.4, not edited): regime C of "S"/"H" = KA pivot blocks + GEMMs for the trailing columns and the contribution block.

Follow-up issues opened: #86 (with step 1). The unmet criteria stay with #75.

🤖 Generated with Claude Code

https://claude.ai/code/session_01Cqp7TUSyCEdjPtLDEftJf5

[claude]

nit (non-blocking): the docstring above is now stale. Its signature line still reads factorize_ldlt!(numeric, symbolic, nzval; opts = Options()) -> 0 (no impl, no nb), and the text says every regime-B/C launch group is one front_ldlt_kernel! launch with "no vendor calls, see the T15 report". After this PR the regime-C fronts and the wide tall regime-B bins run panel_ldlt_kernel! blocks plus GEMMs through src/dense/interface.jl (:vendor on CUDA). Suggest: add impl = :auto (the GEMM implementation) and nb = LDLT_C_NB (block size, 1:LDLT_C_NB, tests only) to the signature, and replace the "no vendor calls" sentence with a pointer to ldlt_blocked_path / src/numeric/ldlt_c.jl.

[claude]

nit (non-blocking): impl is now forwarded to factorize_ldlt!, but the docstring of factorize! (line 197) still says "Structures "S"/"H": factorize_ldlt! with the pivoting options opts (impl is not used)". Update it to say impl selects the GEMM implementation of the blocked regime-C LDLᵀ path.

[claude]

nit (non-blocking): a docstring attaches only to the next binding, so this one documents LDLT_BLOCKED_MIN_WCLASS and LDLT_BLOCKED_MIN_FCLASS ends up undocumented. Either give each constant its own one-liner or move the shared text into a comment and keep short docstrings per constant.

[claude]

nit (non-blocking): the workspace offsets lo/wo/so (up to Layout.work_len, 2× the largest slice) are stored in ctl with % IT, i.e. truncated to INT. The symbolic maps get an overflow check in _index_vector, but work_len is not a map, so with INT = Int32 and a workspace above 2³¹ entries this would wrap silently instead of raising. Unreachable at today's sizes (apache2: 28 M entries), so just a suggestion: check L.work_len <= typemax(INT) once in factorize_ldlt! (or at adapt time) with an InvalidValueError, like the maps.

[claude]

VERDICT: APPROVE (head 452b273)

Reviewed the step-3 commits (cd3776b..452b273, the four commits stacked on #88/#87, which carry their own reviews) against the T25 owner note for #75, AGENTS.md and PLAN.md §2.4/§3.3.

Checked

• Scope: perf PR under T25/#75 (owner-sanctioned, labelled performance, milestone M11). PLAN.md untouched; TASKS.md changed by one line under the T25 owner note only; no Manifest.toml; no public name, parameter or phase string changed (factorize_ldlt! gains impl/nb keywords, factorize! forwards impl).
• Mathematics of src/numeric/ldlt_c.jl: the lazy panel (_LazyFront: stored value minus the block's pending pivots in pivot order, 2×2 as one term) reproduces the reference's right-looking update exactly up to rounding, Hermitian vs complex-symmetric via _cj/'C'/'T' consistently; the interchange swaps rows of Lb/Wb together with the stored panel; the straddling 2×2 pivot (slots nb, nb+1, column saved before the trailing GEMM and restored by the late-starting next block, slot 1 left zero) is consistent in both GEMMs and the next block's pending set; the fallback scan's per-entry threshold form is equivalent to _lt_threshold_ok (monotone u·max, NaN rejected) and its duplicate flag writes are benign; _LT_KEND keeps the fused kernels' behaviour (_lt_reset! sets it after _ST_W in all three kernels).
• Layout/chunking: ldlt_c_work_len matches _ltc_offsets and the device/host slice accumulation (cb_ptr > 0 ⇔ cb_len > 0); the runtime greedy with cap = work_len yields the same chunks the layout sized (work_len ≤ 2·single), so every slice fits; group_width == 0 ⇒ takes_c_path ⇒ ldlt_blocked_path, so the layout's and the numeric phase's choice of groups agree; batched members reuse the first slice sequentially (stream-ordered).
• Conventions: generic T/INT, no Float64/Int in kernels (only pre-existing Int32 for piv/info), 1-D workgroups, @localmem from Val, no atomics, no host sync in the phase, dense ops only through _gemm_impl! (:vendor on CUDA, :ka/:generic elsewhere, β = 0 on the first block is safe in all three), InvalidValueError for nb, docstrings on the new public names.
• Tests: new testset loops over BACKENDS × ELTYPES, Random.seed!(666), tol/panel_tol, shared ldlt_setup/kkt_matrix/kkt_interleaved_perm; checks piv, pivot_kind, D, panels, stats, totals and relres for every dense impl, at a straddling block size and at LDLT_C_NB; range check on nb. The "C only" regime of the existing suites (incl. complex symmetric, 'D'/'N', uniform batches in test_ubatch) exercises the new path too. The one adaptation (Float32 kkt(300,100,1e-8) gets the refinement step ComplexF32 already had) is documented in place with the measured reference bound; panel_error ≤ panel_tol for that case is still asserted, so the factor itself is unchanged to rounding. No @test_broken added.
• Report: matches the diff (files, workspace sizes, unmet criteria stated, #86 cited, deviations listed). Test counts plausible; "CUDA: local RTX 4080" — CI (test-github-cpuonly, test-gpu cuda) is still pending at review time and is handled by the pipeline. I could not run the CPU suite in this review session (sandbox permission), so test verification rests on CI.

Blocking findings: none.

Non-blocking nits (inline): stale factorize_ldlt! and factorize! docstrings (signature/impl, "no vendor calls"); the shared docstring before LDLT_BLOCKED_MIN_WCLASS/_FCLASS attaches to the first constant only; workspace offsets are truncated to INT without the overflow check the maps get.

[github-actions]

CI is green and the reviewer approved, but the merge failed (conflict with main or a branch rule). Labelled needs-owner.

[michel2323]

Rebased onto the rebased #88 (step-1 commits dropped after #87's squash merge; the tree is identical to the approved head 452b273). Removed needs-owner. Merge after #88.

[michel2323]

Draft until #88 merges: its squash would conflict with the step-2 commits still at the bottom of this branch. I will rebase onto main then and mark it ready (tree unchanged).

[claude]

nit (non-blocking): the regime-B branch of ldlt_blocked_path (width class ≥ 32, row class ≥ 256, so bins (32..64) × (256..512) under the default regime_c_width = 64 / regime_c_rows = 512) is the one place where the blocked path is reached with the group's assembly done by zero_fronts!/scatter_A!/extend_add! instead of the fused kernel. The new testset forces "C only", so it does not cover this branch; whether any default-analysis LDLᵀ test matrix (all have n ≤ 400) produces such a bin is not asserted anywhere. Suggest either an @test in the new testset that at least one front of some default-analysis case satisfies ldlt_blocked_path(sc, s) && !takes_c_path(sc, s) (and then compare piv/pivot_kind/D/panels with the reference as the "C only" case does), or a case with regime_c_width/regime_c_rows raised so that a 32-wide, 256-row front stays in regime B. The owner's bcsstk38/lap2d runs exercise it on CUDA, but the suite should too.

[claude]

VERDICT: APPROVE (head 5fbca4d)

Re-review after the owner's rebase onto the rebased #88. The tree of this head is bitwise identical to the previously approved head 452b273 (git rev-parse <sha>^{tree} = 7c542d7b… for both), so the earlier findings stand; this pass re-read the step-3 commits (diff ef5c6f5…5fbca4d, i.e. on top of #88's head) independently rather than relying on the earlier verdict.

Checked

• Scope: owner-sanctioned perf PR under T25/#75 (label performance, milestone M11, Refs #75, not Closes, with the unmet criteria stated). PLAN.md untouched; TASKS.md changed by one line under the T25 owner note; no Manifest.toml; no public name, parameter or phase string changed (factorize_ldlt! gains impl/nb keywords, factorize! forwards impl, build_layout gains ldlt).
• Mathematics of src/numeric/ldlt_c.jl, re-derived: the lazy read F[i,j] − Σₜ Lb[i,t]·cj(Wb[j,t]) (2×2 as one term, pivot order) is exactly the fused kernel's _lt_update! arithmetic (l = x/d, l1 = x1 e11 + x2 e21, l2 = x1 e12 + x2 e22); the commit materializes with cur = false (pending set = pivots before the current one), pass 1 keeps the column only when no interchange happened (C == k, 1×1); the interchange moves the stored panel (_lt_raw_front) and the rows of Lb/Wb together, which keeps the lazy value permutation-consistent (and Hermitian-consistent, since E is Hermitian and d real); the straddling 2×2 at the block's last column occupies slots nb, nb+1, is saved (rows kend+1:f) before the trailing GEMM corrupts column k1, restored by the next block (late ⇔ the save condition, t0 = 2, slot 1 zero), and the next block's K/KEND and slot numbering agree; g0 + t indexes pivot_kind of column k0 + t − 1. The fallback scan's per-entry form equals _lt_threshold_ok (u·max = max(u·|x|) exactly for u ≥ 0, NaN rejected), _ltc_pick! reproduces _lt_pass3!'s best/tie/'D' rules; duplicate flag writes are benign. finish_ldlt_kernel!: _lt_finalize! touches strictly-lower entries only, the clear loop the diagonal and upper triangle (the GEMMs write it), so no barrier is needed. _LT_KEND is set in _lt_reset! after _ST_W in all three setup paths.
• Layout/launch consistency: _ltc_offsets matches ldlt_c_work_len; cb_ptr > 0 ⇔ cb_len > 0 (build_layout), so the device, host and layout agree on each slice; the runtime greedy with cap = work_len yields exactly the chunks the layout sized with 2·single (every chunk sum ≤ work_len, and the layout's stop condition implies the runtime's); W == 0 ⇒ regime C ⇒ ldlt_blocked_path; both build_layout callers (solver.jl, maps.jl) pass ldlt, the only plain call is an SPD test. pack_add! adds the lower triangle exactly as the Cholesky regime-C path.
• Conventions: generic T/INT, no Float64/Int in kernels, 1-D workgroups, @localmem from Val, no atomics, no host sync in the phase, dense ops through _gemm_impl! only, InvalidValueError for nb, docstrings on new public names.
• Tests: new testset over BACKENDS × ELTYPES, Random.seed!(666), tol/panel_tol, shared helpers; piv, pivot_kind, D, panels, stats, totals, relres for every dense impl at a straddling nb and at LDLT_C_NB; range check on nb. The Float32 kkt(300,100,1e-8) adaptation is documented with the measured reference bound and keeps the panel_error assertion. No @test_broken added.
• Report (PR body): matches the diff; #86 exists (found-by-agent); numbers are from the owner's machine. CI: build green, test-github-cpuonly and test-gpu (cuda) still pending and handled by the pipeline. I could not run Julia in this review sandbox (commands denied), so test verification rests on CI.

Blocking findings: none.

Non-blocking: the four nits from the previous round (stale factorize_ldlt!/factorize! docstrings, shared docstring on LDLT_BLOCKED_MIN_WCLASS, work_len offsets truncated to INT without a check) still apply to this identical tree; one new inline nit on test coverage of the blocked regime-B branch of ldlt_blocked_path.