perf: regime-B LDLᵀ with F₁₁ in local memory, no scale phase, tiled update (#75)

[michel2323]

Refs #75

Stacked on #87 (its four commits are at the bottom of this branch); merge #87 first. The step-2 commits: "stage F₁₁ in local memory", "regime-B LDLᵀ step restructured, tiled contribution-block update", "narrow regime-B bins keep the panel in global memory", "experiment 1 step 2 results".

Task: #75, experiment 1 of PERFORMANCE.md, step 2 of 3 — regime B with F₁₁ in local memory (T25 owner note, item 3).

What was built (src/numeric/ldlt.jl):

• front_ldlt_kernel! takes the width class W of its launch group (Val(W); 0 = panel in global memory) and keeps the packed W×W F₁₁ in @localmem (_SplitFront: rows 1:w local, rows below in the panel) from W = 32; classes 8 and 16 keep the panel in global memory (_LT_GLOBAL_MAX_W), the staging phases cost more than they save on KKT fronts.
• Staging alone gained nothing (lap2d_300 190 → 192 ms, the rows below F₁₁ and the contribution block dominate), so the step was restructured with the reference's arithmetic unchanged: pivot columns stay unscaled (L D) until the end of the front, where _lt_finalize! applies the same divisions and 2×2 transforms with the pivot recomputed from D (no scale phase); the update reduces the next column as it writes it (pass 1); work items own the rows below F₁₁ when there are at least WG of them (one division per row). Four barriers per step (five before).
• Contribution blocks with m ≥ 64 are updated in 16×16 tiles (_lt_cb_load!/_lt_cb_acc!/_lt_cb_store!) with W₂₁ = L₂₁D and L₂₁ staged 16 columns at a time in the F₁₁ buffer (sized at least _LT_TBUF); smaller blocks keep _lt_cb_update!. Regime A and regime C (still this kernel until step 3) use the same finalize/update.
• PERFORMANCE.md: step-2 section.

Tests: julia --project=. -e 'using Pkg; Pkg.test()' with CUDA in the local test environment (RTX 4080): 104914 pass / 0 fail / 2 broken (T16 @test_broken, both backends) on the step before "narrow bins"; after it, on CPU and CUDA, test_numeric_ldlt 3820 pass and test_api, test_ported, test_refinement, test_ubatch, test_solve, test_symbolic_schedule 36571 pass, 0 fail, 2 broken; rebased on T18 (#83, no overlap; the stacked step 3 passes test_numeric_ldlt/test_symbolic_schedule/test_refinement on CPU and CUDA after the rebase). No test changed. factorize! LDLᵀ still 64 B on the CPU backend.

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

matrix	main ms	step 1 ms	step 2 ms	step 2 / cuDSS	step 2 / SDS Cholesky
lap2d_300	228	192	116	22.5×	3.68×
lap3d_40	14,777*	14,494*	8,753*	150×	39.1×
HB/bcsstk17	206	167	115	34.2×	4.14×
Boeing/bcsstk38	137	106	72.9	22×	3.22×
GHS_psdef/apache2	92,576*	97,296*	50,736*	104×	49.3×
kkt_pglib_opf_case118_ieee_condensed_1	1.77	1.82	1.91	4.89×	2.09×
kkt_pglib_opf_case118_ieee_condensed_10	1.72	1.79	1.87	5.35×	2.12×
kkt_pglib_opf_case118_ieee_condensed_20	1.87	1.82	1.88	5.36×	2.15×
kkt_pglib_opf_case118_ieee_k2_1	1.69	1.69	1.8	4.97×
kkt_pglib_opf_case118_ieee_k2_10	1.62	1.68	1.76	4.97×
kkt_pglib_opf_case118_ieee_k2_20	1.6	1.66	1.74	4.36×
kkt_pglib_opf_case1354_pegase_condensed_1	10.9	7.2	7.51	11×	2.56×
kkt_pglib_opf_case1354_pegase_condensed_10	10.9	7.3	7.56	11.8×	2.59×
kkt_pglib_opf_case1354_pegase_condensed_20	11	7.32	7.56	10.7×	2.59×
kkt_pglib_opf_case1354_pegase_k2_1	8.31	6.6	6.95	8.51×
kkt_pglib_opf_case1354_pegase_k2_10	7.15	6.05	6.27	7.58×
kkt_pglib_opf_case1354_pegase_k2_20	9.1	6.71	6.84	8.48×
kkt_pglib_opf_case14_ieee_condensed_1	0.446	0.503	0.484	2.51×	1.33×
kkt_pglib_opf_case14_ieee_condensed_10	0.443	0.52	0.517	2.26×	1.52×
kkt_pglib_opf_case14_ieee_condensed_11	0.448	0.463	0.515	1.77×	1.57×
kkt_pglib_opf_case14_ieee_k2_1	0.669	0.697	0.728	2.6×
kkt_pglib_opf_case14_ieee_k2_10	0.722	0.708	0.727	3.06×
kkt_pglib_opf_case14_ieee_k2_15	0.97	1.03	1.07	3.33×

• K2 dumps (CUDA): inertia and nperturbed identical to main on all nine.
• kkt_matrix(Float64, 3000, 1000, 1e-8) (T15 number): CUDA 8.3 s (step 1) → 5.6 s; KA CPU backend 7.0 → 7.4 s, ref_ldlt! 2.7 s.
• Experiment-1 criterion (within 1.5× of SDS Cholesky) not yet met: lap2d 3.7×, bcsstk17 4.1×, case1354 condensed 2.6×.

Deviations from PLAN.md / TASKS.md: F₁₁ staging only from width class 32 (measured, see above). The KKT dumps are 3–6% slower than step 1 (one more phase per front: finalize).

Follow-up issues opened: none for this step (#86 opened with step 1).

🤖 Generated with Claude Code

https://claude.ai/code/session_01Cqp7TUSyCEdjPtLDEftJf5

[claude]

nit (non-blocking): the docstring (and the header comment at the top of the file, plus the factorize_ldlt! docstring: "F₁₁ in local memory for regime-B width classes") still equate W = 0 with regime C. After the "narrow bins" commit the regime-B classes 8 and 16 also run with W = 0 (_LT_GLOBAL_MAX_W), so only classes 32 and 64 stage F₁₁. Suggest: "W = 0: regime C and the regime-B classes ≤ _LT_GLOBAL_MAX_W, in the panel; W ≥ 32: F₁₁ staged in @localmem".

[claude]

nit (non-blocking): _lt_front(fa::Tuple, ctl, ::Val{W}) is ambiguous with _lt_front(fa, ctl, ::Val{true}) for a call _lt_front((L11, factor), ctl, Val(true)) (neither method is more specific). It is never reached (the tuple form is only paired with an integer Val(W), and Aqua runs with ambiguities = false), so nothing breaks; constraining the third argument to ::Val{W} where {W <: Integer}, or dispatching the tuple form on _lt_pk's result, would remove the latent ambiguity.

[claude]

nit (non-blocking, performance): _lt_cb_tile is an O(nt) search, and every work item evaluates it in every _lt_cb_load! call (once per tile and chunk) and again in _lt_cb_store!. For regime-C blocks with nt in the hundreds (apache2, lap3d_40: m of several thousand rows) that is comparable to the chunk's own work per work item (4 loads + 32 FMAs). Work item 1 could store (I, J) of the current tile in two ctl slots before the chunk loop (there is already a barrier after the first _lt_cb_load!), or use the closed-form inverse of the column-major triangular numbering.

[claude]

VERDICT: APPROVE (head cd3776b)

What was checked (diff origin/main...HEAD, 7 files; the four step-2 commits touch only src/numeric/ldlt.jl and PERFORMANCE.md, the other files are the stacked #87 commits):

• Scope: experiment 1 step 2 of PERFORMANCE.md (T25 owner note, item 3) on the owner-triaged performance issue #75, same shape as #81/#87. PLAN.md, TASKS.md and Manifest.toml untouched; no public name, parameter or phase string changed; no new dependency. Refs #75 rather than Closes is right, since T25 closes #75.
• Kernel correctness, read against the reference _ldlt_front!: the deferred scaling is exact (_lt_finalize! recomputes e11..e22 from d[g], d[n+g], d[g+1] with the same expressions _lt_select! used, and divides 1x1 columns by the stored, possibly perturbed, d[g]); row interchanges of later steps commute with the per-column scaling, so the unscaled panel gives the same L; row k+1 of a 2x2 block keeps the zero _lt_swap! wrote. The rank-1/rank-2 update matches the reference (l1 = x1 e11 + x2 e21, l2 = x1 e12 + x2 e22, conj mirrors only for Hermitian). The fused pass-1 reduction in _lt_update! (TRACK) covers every entry of the next column exactly once, visits each lane in increasing row order, keeps lambda block-only and the column maximum over all rows, and _lt_nused bounds every lane that can hold data (block lanes <= nt, row lanes <= m, or all WG lanes when m >= WG), so the arg-max equals the serial first maximum. _SplitFront routes only rows <= w to local memory and every reader goes through _fget/_fset!; the staging loops cover exactly the lower triangle and _packed(i, j, W) with w <= W stays in bounds. Tiled F22 update: W21 = L21 D for 1x1, 2x2-first and 2x2-second columns matches _lt_cb_update!, tile numbering is consistent between load/acc/store, the buffer regions used by store(tt) and load(tt+1) are disjoint with a barrier before acc, lower-only and real-diagonal (Hermitian) stores. Regime C and classes 8/16 use W = 0 with the panel in global memory; regime A reuses the same update/finalize with TRACK = false and its reduction slots are not written by the update. No atomics, no allocation or host sync inside the phase, @localmem sizes from Val/constants, 1-D workgroups, errors via InvalidValueError. Local memory at the largest class (W = 64, ComplexF64, Int64 maps) is about 39 KB, within the CUDA limit and the same order as the Cholesky regime-B kernel.
• Tests: no test changed in step 2 (the one-line test_symbolic_schedule.jl change belongs to #87 and only passes the structure-dependent reserve to subtree_capacity; no weakening). test_numeric_ldlt compares piv, pivot_kind, stats, D and the panel with the reference over BACKENDS x ELTYPES and the three regime mixes, which exercises W = 0, the staged class 32 and the tiled update (KKT and random fronts with m >= 64). The claimed counts are plausible for the suite.
• PR body / PERFORMANCE.md: every "what was built" claim matches the code (_LT_GLOBAL_MAX_W = 16, _LT_TILED_M = 64, 16x16 tiles, L11 sized max(W(W+1)/2, _LT_TBUF), four barriers per step, finalize at the end of the front); the deviation (staging only from class 32) is stated with its measurement; the two tables agree.
• CI (gh pr checks 88): build green; test-github-cpuonly, test-gpu (cuda) and the review workflow still pending at the time of this review. The body reports CPU and CUDA green locally and nothing in CI contradicts it yet; the pipeline handles the pending checks. I could not run Julia in this session.

Non-blocking inline nits (3): docstrings still describe W = 0 as regime C only; a latent method ambiguity of _lt_front(::Tuple, ctl, ::Val{true}) (never reached, Aqua ambiguities off); _lt_cb_tile is an O(nt) loop evaluated per work item, per tile and chunk.

Note for the owner: this PR is stacked on #87 (its four commits are included in this diff; #87 has no verdict yet). Merging #88 first squashes the #87 content in and leaves #87 empty; merging #87 first reduces this diff to the step-2 commits. Either order is fine for the code reviewed here, but the pipeline should not merge both.

Blocking findings: none.

[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 main after #87's squash merge (the step-1 commits are dropped; the tree is identical to the approved head cd3776b, git diff cd3776b ef5c6f5 is empty). The conflict came from the squash, not from a content change; removed needs-owner.

[claude]

nit (non-blocking): stale section header. Regime B now stages F₁₁ in local memory for width classes 32 and 64; only W = 0 (regime C and the 8/16 bins) keeps the panel in global memory. Same for the pass-1 comment at line 190 ("computed in the scale phase of the previous step"): the regime-B/C kernel now reduces pass 1 in _lt_update! (TRACK), the scale phase is gone; only regime A still calls _lt_pass1! after the update.

[claude]

nit (non-blocking): _lt_cb_tile is an O(nt) scan and every work item re-runs it in _lt_cb_load! for each of the cld(w, 16) chunks of every tile, and again in _lt_cb_store!. On a root front with m in the thousands (nt in the hundreds) that is a few hundred scalar iterations per chunk per work item, all redundant. Work item 1 could write (I, J) into two control slots once per tile (one more barrier per tile, or reuse the existing one before the first load), or the kernel could iterate J, I directly instead of a linear tile index. Also worth noting for the Report: the tiled path accumulates the two terms of a 2×2 block as two separate adds (acc += W[i,k]·conj(L[j,k]); acc += W[i,k+1]·conj(L[j,k+1])), where _lt_cb_update! adds their sum once, so the two paths round differently for 2×2 pivots (still deterministic; d_error/panel_error are tolerance tests, so fine).

[claude]

VERDICT: APPROVE (head ef5c6f5)

Perf PR for issue #75, experiment 1 step 2 (T25 owner note, item 3); like #87 it is not a TNN task, so the Report lives in the PR body. #87 is merged, so the diff against main is exactly the four step-2 commits: src/numeric/ldlt.jl and PERFORMANCE.md. PLAN.md, TASKS.md, tests and Manifest untouched; no public name, parameter or phase string changed.

Checked

• _SplitFront / _lt_stage!: packed W×W indexing with W ≥ w (width classes 32, 64 from _with_width_class; 8, 16 and regime C with group_width = 0 take Val(0)). Local memory at W = 64, ComplexF64: about 33 KB for F₁₁ plus about 5 KB of reductions, within the CUDA static limit; the _LT_TBUF floor covers the tiles when W is small.
• Unscaled pivot columns: the update already used l = F[i,k]/d on unscaled columns, so the arithmetic is unchanged; row interchanges of already-pivoted columns commute with the per-column scaling (1×1 and 2×2); _lt_finalize! recomputes det, c2/det, -b/det, -up/det, a/det from D in the same order as _lt_select! formed pv[1:4], so L is bitwise what the old scale phase produced; the sub-diagonal zero of a 2×2 block is skipped (j + 1 > jmax). Regime A calls finalize before _lt_cb_update! and the panel write.
• Pass 1 fused into _lt_update! (TRACK): every entry of the next column below the diagonal is written exactly once; each lane visits its rows in increasing order with strict >, so first-maximum tie-breaking is preserved; the lanes used are within 1:min(WG, f - k + 1), so the _lt_nused bounds of _lt_stage1!/_lt_final! still hold. Four barriers per step as claimed.
• Tiled F₂₂ update: W₂₁ = L₂₁D for 1×1, 2×2-first and 2×2-second columns matches _lt_cb_update!; tile to (I, J) mapping, chunk bounds, Hermitian real diagonal, and the local-buffer hazards (load/acc separated by barriers, store reads only the work item's own accumulators, the next tile's load touches a disjoint region) are sound. No atomics, deterministic.
• Kernel argument count stays below the 32-argument specialization limit (29); _with_width_class keeps the launch allocation-free, consistent with the "still 64 B" claim.
• Report vs diff: all claims match (_LT_GLOBAL_MAX_W = 16, _LT_TILED_M = 64, finalize/update shared by all regimes, PERFORMANCE.md step-2 section and tracked-issue row). Deviations (staging from class 32 only, KKT dumps 3–6 % slower) are stated.

Not verified here: the review sandbox refused to run Pkg.test, so the CPU and CUDA counts rest on the Report and on CI (build green, test jobs pending at review time; the pipeline handles them).

Blocking findings: none. Two non-blocking nits inline (stale comments; _lt_cb_tile recomputed per chunk per work item, and a note that the tiled 2×2 accumulation rounds differently from _lt_cb_update!).