From 9028a940d67edc24365ffe4e4c3e90cdd58f1237 Mon Sep 17 00:00:00 2001 From: Michel Schanen Date: Mon, 5 Oct 2026 09:26:43 -0500 Subject: [PATCH 1/4] =?UTF-8?q?perf:=20stage=20F=E2=82=81=E2=82=81=20in=20?= =?UTF-8?q?local=20memory=20in=20the=20regime-B=20LDL=E1=B5=80=20kernel=20?= =?UTF-8?q?(#75)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Regime-B launch groups run front_ldlt_kernel! with their width class W: the packed W×W lower triangle of F₁₁ lives in @localmem for the pivot search, interchanges and updates, the rows below it stay in the panel. Regime C keeps the panel in global memory (W = 0). Co-Authored-By: Claude Opus 5.5 --- src/numeric/ldlt.jl | 122 ++++++++++++++++++++++++++++++++++---------- 1 file changed, 94 insertions(+), 28 deletions(-) diff --git a/src/numeric/ldlt.jl b/src/numeric/ldlt.jl index 7b96d5b..1abf6b3 100644 --- a/src/numeric/ldlt.jl +++ b/src/numeric/ldlt.jl @@ -10,8 +10,9 @@ # Cholesky kernel of `src/numeric/subtree.jl`). # * Regimes B and C: `front_ldlt_kernel!`, one workgroup per front of a launch # group, fused assembly (zero, scatter A, owner-pull extend-add), then the -# pivoted factorization of the panel in global memory and the update of the -# front's packed contribution block on the update stack. +# pivoted factorization of the panel (regime B: `F₁₁` staged in `@localmem`, +# the rows below it in global memory; regime C: all in global memory) and the +# update of the front's packed contribution block on the update stack. # # Per pivot step (a loop over the `w` columns of the front with a uniform trip # count; a 2×2 pivot leaves the last iterations idle): the workgroup chooses the @@ -117,9 +118,59 @@ end # F[i, j] from the lower triangle (upper entries: conj for Hermitian, plain for complex symmetric) @inline _fsym(F, i, j, h::Val) = i >= j ? _fget(F, i, j) : _cj(_fget(F, j, i), h) -# the front of the current node from the control words: offset `_ST_LF` (1-based) into `fa` +# regime B: the fully-summed block F₁₁ (rows and columns `1:w`) packed `W×W` lower triangle in local +# memory `l`, the rows `w+1:f` in the `f×w` panel in global memory +struct _SplitFront{L, A} + l::L + a::A + off::Int + f::Int + w::Int + W::Int +end + +@inline _fget(F::_SplitFront, i, j) = @inbounds i <= F.w ? F.l[_packed(i, j, F.W)] : F.a[F.off + (j - 1) * F.f + i] +@inline function _fset!(F::_SplitFront, i, j, v) + @inbounds if i <= F.w + F.l[_packed(i, j, F.W)] = v + else + F.a[F.off + (j - 1) * F.f + i] = v + end + return nothing +end + +# the front of the current node from the control words: offset `_ST_LF` (1-based) into `fa`; `Val(W)` +# with `fa = (L11, factor)`: F₁₁ staged in the local `L11` (width class `W`) @inline _lt_front(fa, ctl, ::Val{true}) = @inbounds _PackedFront(fa, Int(ctl[_ST_LF]) - 1, Int(ctl[_ST_F])) @inline _lt_front(fa, ctl, ::Val{false}) = @inbounds _PanelFront(fa, Int(ctl[_ST_LF]) - 1, Int(ctl[_ST_F])) +@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 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) +@inline _lt_pk(::Val{0}) = Val(false) +@inline _lt_pk(w::Val) = w + +# copy F₁₁ between the panel and local memory (`load`), or write it back with a unit diagonal +@inline function _lt_stage!(L11, factor, ctl, li, ::Val{W}, ::Val{WG}, load::Bool) where {W, WG} + @inbounds if W > 0 + f = Int(ctl[_ST_F]) + w = Int(ctl[_ST_W]) + p0 = 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 || continue + if load + L11[_packed(i, j, W)] = factor[p0 + (j - 1) * f + i] + else + factor[p0 + (j - 1) * f + i] = i == j ? one(eltype(factor)) : L11[_packed(i, j, W)] + end + end + end + return nothing +end # --------------------------------------------------------------------------- # cooperative pivot search (the reference's `_choose_pivot` on the lower triangle, pivot for pivot) @@ -320,7 +371,7 @@ end @inbounds begin n = length(perm) IT = eltype(ctl) - T = eltype(fa) + T = eltype(pv) k = Int(ctl[_LT_K]) c0 = Int(ctl[_LT_C0]) g = c0 + k - 1 @@ -562,7 +613,7 @@ end end if i == k _fset!(F, k, k, pv[5]) - r == 0 || _fset!(F, k + 1, k, zero(eltype(fa))) # the 2×2 block lives in D, L's is I + r == 0 || _fset!(F, k + 1, k, zero(eltype(pv))) # the 2×2 block lives in D, L's is I elseif r != 0 && i == k + 1 _fset!(F, k + 1, k + 1, pv[6]) end @@ -718,14 +769,16 @@ end front_ldlt_kernel!(backend, WG)(factor, stack, info, stats, d, piv, pivot_kind, psign, perm, aux, nzval, amap, amap_ptr, amap_src, nodes, bm, super_ptr, front_ptr, front_nrows, front_ncols, cb_ptr, child_ptr, child_list, relind_ptr, relind, maxchild, prm, - Val(WG); ndrange = WG * count * bm.nact) + Val(W), Val(WG); ndrange = WG * count * bm.nact) LDLᵀ/LDLᴴ kernel of regimes B and C: workgroup `G` takes front `s = nodes[bm.first + g - 1]` of batch member `k` (`g`, `k` from the [`BatchMap`](@ref) `bm`; the per-member arrays are the member's), zeroes and assembles it (A through the `amap`, then the `maxchild` children's packed contribution blocks in `child_list` order), -factors its `w` fully-summed columns in place in the panel with in-block -Bunch–Kaufman pivoting, threshold acceptance and perturbation (`prm`, a +factors its `w` fully-summed columns in place with in-block +Bunch–Kaufman pivoting (`W = 0`, regime C: in the panel; regime B, width +class `W ≥ w`: `F₁₁` staged as a packed `W×W` lower triangle in `@localmem`, +the rows below it in the panel; the pivot search is a workgroup reduction), threshold acceptance and perturbation (`prm`, a `_LDLTDevice{R, HERM}`; `HERM`: Hermitian or real symmetric, else complex symmetric), writes D, the local pivot order `piv`, the pivot kinds and the front's statistics, updates its packed contribution block on the update stack @@ -734,7 +787,7 @@ statistics, updates its packed contribution block on the update stack @kernel function front_ldlt_kernel!(factor, stack, info, stats, d, piv, pivot_kind, psign, perm, aux, nzval, amap, amap_ptr, amap_src, nodes, bm, super_ptr, front_ptr, front_nrows, front_ncols, cb_ptr, child_ptr, child_list, relind_ptr, relind, maxchild, - prm::_LDLTDevice{R, HERM}, ::Val{WG}) where {R, HERM, WG} + prm::_LDLTDevice{R, HERM}, ::Val{W}, ::Val{WG}) where {R, HERM, W, WG} @uniform TT = eltype(factor) @uniform IT = eltype(front_ptr) li = @index(Local, Linear) @@ -743,6 +796,7 @@ statistics, updates its packed contribution block on the update stack pv = @localmem TT (_LT_NPV,) red = @localmem R (3 * WG,) redi = @localmem IT (2 * WG,) + L11 = @localmem TT (max(W * (W + 1) ÷ 2, 1),) if li == 1 @inbounds s = nodes[bm.first + _bm_node(bm, G) - 1] _lt_front_setup!(ctl, s, super_ptr, member_panels(front_ptr, _bm_gmember(bm, G), bm.nbatch), front_nrows, @@ -768,7 +822,9 @@ statistics, updates its packed contribution block on the update stack front_nrows, front_ncols, cb_ptr, child_ptr, child_list, relind_ptr, relind, Val(WG)) @synchronize end - _lt_pass1!(factor, ctl, red, redi, li, Val(WG), prm, Val(false)) + _lt_stage!(L11, factor, ctl, li, Val(W), Val(WG), true) + @synchronize + _lt_pass1!(_lt_fa(L11, factor, Val(W)), ctl, red, redi, li, Val(WG), prm, _lt_pk(Val(W))) @synchronize for it in 1:ctl[_ST_W] if WG > _LT_S1 @@ -778,12 +834,12 @@ statistics, updates its packed contribution block on the update stack if li == 1 k = _bm_gmember(bm, G) nb = bm.nbatch - _lt_decide1!(factor, ctl, pv, red, redi, _mview(d, k, nb), _mview(pivot_kind, k, nb), _mview(piv, k, nb), - psign, perm, _mview(aux, k, nb), prm, Val(WG), Val(false), Val(HERM)) + _lt_decide1!(_lt_fa(L11, factor, Val(W)), ctl, pv, red, redi, _mview(d, k, nb), _mview(pivot_kind, k, nb), _mview(piv, k, nb), + psign, perm, _mview(aux, k, nb), prm, Val(WG), _lt_pk(Val(W)), Val(HERM)) end @synchronize if ctl[_LT_PHASE] == 2 - _lt_pass2!(factor, ctl, red, redi, li, Val(WG), prm, Val(false)) + _lt_pass2!(_lt_fa(L11, factor, Val(W)), ctl, red, redi, li, Val(WG), prm, _lt_pk(Val(W))) @synchronize if WG > _LT_S1 _lt_stage1!(red, redi, ctl, li, Val(WG), Val(2)) @@ -792,14 +848,14 @@ statistics, updates its packed contribution block on the update stack if li == 1 k = _bm_gmember(bm, G) nb = bm.nbatch - _lt_decide2!(factor, ctl, pv, red, redi, _mview(d, k, nb), _mview(pivot_kind, k, nb), - _mview(piv, k, nb), psign, perm, _mview(aux, k, nb), prm, Val(WG), Val(false), Val(HERM)) + _lt_decide2!(_lt_fa(L11, factor, Val(W)), ctl, pv, red, redi, _mview(d, k, nb), _mview(pivot_kind, k, nb), + _mview(piv, k, nb), psign, perm, _mview(aux, k, nb), prm, Val(WG), _lt_pk(Val(W)), Val(HERM)) end @synchronize end if ctl[_LT_PHASE] == 3 - _lt_pass3!(factor, ctl, red, redi, li, Val(WG), _mview(aux, _bm_gmember(bm, G), bm.nbatch), prm, - Val(false)) + _lt_pass3!(_lt_fa(L11, factor, Val(W)), ctl, red, redi, li, Val(WG), _mview(aux, _bm_gmember(bm, G), bm.nbatch), prm, + _lt_pk(Val(W))) @synchronize if WG > _LT_S1 _lt_stage1!(red, redi, ctl, li, Val(WG), Val(3)) @@ -808,19 +864,21 @@ statistics, updates its packed contribution block on the update stack if li == 1 k = _bm_gmember(bm, G) nb = bm.nbatch - _lt_decide3!(factor, ctl, pv, red, redi, _mview(d, k, nb), _mview(pivot_kind, k, nb), - _mview(piv, k, nb), psign, perm, _mview(aux, k, nb), prm, Val(WG), Val(false), Val(HERM)) + _lt_decide3!(_lt_fa(L11, factor, Val(W)), ctl, pv, red, redi, _mview(d, k, nb), _mview(pivot_kind, k, nb), + _mview(piv, k, nb), psign, perm, _mview(aux, k, nb), prm, Val(WG), _lt_pk(Val(W)), Val(HERM)) end @synchronize end - _lt_swap!(factor, ctl, pv, _mview(piv, _bm_gmember(bm, G), bm.nbatch), li, Val(WG), Val(false), Val(HERM)) + _lt_swap!(_lt_fa(L11, factor, Val(W)), ctl, pv, _mview(piv, _bm_gmember(bm, G), bm.nbatch), li, Val(WG), _lt_pk(Val(W)), Val(HERM)) @synchronize - _lt_update!(factor, ctl, pv, li, Val(WG), Val(false), Val(HERM)) + _lt_update!(_lt_fa(L11, factor, Val(W)), ctl, pv, li, Val(WG), _lt_pk(Val(W)), Val(HERM)) @synchronize - _lt_scale!(factor, ctl, pv, li, Val(WG), Val(false)) - _lt_pass1!(factor, ctl, red, redi, li, Val(WG), prm, Val(false)) + _lt_scale!(_lt_fa(L11, factor, Val(W)), ctl, pv, li, Val(WG), _lt_pk(Val(W))) + _lt_pass1!(_lt_fa(L11, factor, Val(W)), ctl, red, redi, li, Val(WG), prm, _lt_pk(Val(W))) @synchronize end + _lt_stage!(L11, factor, ctl, li, Val(W), Val(WG), false) + @synchronize @inbounds s = nodes[bm.first + _bm_node(bm, G) - 1] @inbounds coff = Int(cb_ptr[s]) - 1 k = _bm_gmember(bm, G) @@ -1094,12 +1152,13 @@ end # --------------------------------------------------------------------------- # driver -function _launch_front_ldlt!(N::Numeric, S::Symbolic, nzval, first, count, maxchild, prm, ::Val{WG}) where {WG} +function _launch_front_ldlt!(N::Numeric, S::Symbolic, nzval, first, count, maxchild, prm, ::Val{W}, + ::Val{WG}) where {W, WG} bm = batch_map(N; first) kernel! = front_ldlt_kernel!(KernelAbstractions.get_backend(N.factor), WG) kernel!(N.factor, N.stack, N.info, N.stats, N.d, N.piv, N.pivot_kind, N.psign, S.perm, N.aux, nzval, S.amap, S.amap_ptr, S.amap_src, S.group_nodes, bm, S.super_ptr, S.front_ptr, S.front_nrows, S.front_ncols, - S.cb_ptr, S.child_ptr, S.child_list, S.relind_ptr, S.relind, Int(maxchild), prm, Val(WG); + S.cb_ptr, S.child_ptr, S.child_list, S.relind_ptr, S.relind, Int(maxchild), prm, Val(W), Val(WG); ndrange = WG * count * bm.nact) return nothing end @@ -1128,7 +1187,14 @@ function _factorize_ldlt_groups!(N::Numeric, S::Symbolic, nzval::AbstractVector, for k in eachindex(plan.group_first) _poll_interrupt(flag) a, b = plan.group_first[k], plan.group_last[k] - _launch_front_ldlt!(N, S, nzval, a, b - a + 1, plan.group_maxchild[k], prm, Val(LDLT_WORKGROUP)) + W = plan.group_width[k] + if W == 0 # regime C: 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 + _launch_front_ldlt!(N, S, nzval, a, b - a + 1, plan.group_maxchild[k], prm, w, Val(LDLT_WORKGROUP)) + end + end end return nothing end @@ -1144,8 +1210,8 @@ choice, same storage): `opts.pivot_type`, `pivot_threshold`, `pivot_epsilon`, other than `n` raises [`InvalidValueError`](@ref)). Regime-A groups are one [`subtree_ldlt_kernel!`](@ref) launch per budget class; every regime-B and regime-C launch group is one [`front_ldlt_kernel!`](@ref) launch (fused -assembly and factorization, one workgroup per front; no vendor calls, see the -T15 report); then [`reduce_stats!`](@ref). `opts.user_host_interrupt` is polled +assembly and factorization, one workgroup per front, `F₁₁` in local memory for +regime-B width classes; no vendor calls, see the T15 report); then [`reduce_stats!`](@ref). `opts.user_host_interrupt` is polled before every launch group ([`InterruptedError`](@ref)). With `pivot_epsilon_alg = "algo1"` [`abs_max!`](@ref) computes the scale first. Fills `numeric.factor` (unit-lower panels), `d`, `piv`, `pivot_kind`, `stats` and `totals`. The From d3c565e89e10e550fd42af113ec89c4f8517a05e Mon Sep 17 00:00:00 2001 From: Michel Schanen Date: Mon, 5 Oct 2026 11:24:42 -0500 Subject: [PATCH 2/4] =?UTF-8?q?perf:=20regime-B=20LDL=E1=B5=80=20step=20re?= =?UTF-8?q?structured,=20tiled=20contribution-block=20update=20(#75)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The pivot columns stay unscaled until the end of the front, where _lt_finalize! applies the reference's divisions (and 2×2 transforms) with the pivot recomputed from D, so a step needs no scale phase; the update reduces the next column as it writes it (four barriers per step). Rows below F₁₁ are owned by work items when there are at least WG of them. Contribution blocks with m ≥ 64 are updated in 16×16 tiles staged in the F₁₁ local buffer (W₂₁ = L₂₁D and L₂₁, 16 columns per chunk). Co-Authored-By: Claude Opus 5.5 --- src/numeric/ldlt.jl | 281 ++++++++++++++++++++++++++++++++++++++------ 1 file changed, 243 insertions(+), 38 deletions(-) diff --git a/src/numeric/ldlt.jl b/src/numeric/ldlt.jl index 1abf6b3..b73b409 100644 --- a/src/numeric/ldlt.jl +++ b/src/numeric/ldlt.jl @@ -19,10 +19,13 @@ # pivot (the reference's `_choose_pivot`, its column maxima as workgroup # reductions with the reference's tie-breaking, in one to three passes; work # item 1 decides and stores D, perturbing a tiny 1×1 pivot, and counts the -# statistics), the workgroup swaps rows/columns, applies the rank-1/rank-2 -# update to the remaining fully-summed columns, scales the pivot columns into -# L and reduces the next column (pass 1 of the next step). Five barriers per -# step when the first pass settles the pivot. +# statistics), swaps rows/columns, and applies the rank-1/rank-2 update to the +# remaining fully-summed columns while reducing the next column (pass 1 of the +# next step). The pivot columns stay unscaled (`L D`) until the end of the +# front, when `_lt_finalize!` divides them by D exactly as the reference does +# after each step; then the contribution block is updated (regimes B/C: tiled +# through local memory for blocks of `_LT_TILED_M` rows or more). Four barriers +# per step when the first pass settles the pivot. # Only the lower triangle of the front is stored; the upper entries the # algorithm reads are `conj` (Hermitian) or plain (complex symmetric) mirrors. # Pivot decisions, control words and the per-front counts live in `@localmem`; @@ -47,7 +50,8 @@ const _LT_R = _ST_CTL + 5 # 2×2 partner column (0: 1×1 pivot) const _LT_STAT = _ST_CTL + 5 # _LT_STAT + q: statistic q ∈ 1:5 (STAT_NPOS … STAT_N2X2) of the front 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_CTL = _ST_CTL + 12 +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 # `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 @@ -629,53 +633,130 @@ end return nothing end -# update of the remaining fully-summed columns (lower part, rows to f) with the unscaled pivot columns: -# F[i, j] -= l_i F[k, j] (1×1) or l1_i F[k, j] + l2_i F[k+1, j] (2×2), as the reference -@inline function _lt_update!(fa, ctl, pv, li, ::Val{WG}, pk::Val, h::Val) where {WG} +# update of the remaining fully-summed columns with the unscaled pivot column(s), as the reference: +# F[i, j] -= l_i F[k, j] (1×1, `l_i = F[i, k] / d`) or l1_i F[k, j] + l2_i F[k+1, j] (2×2). The pivot +# columns stay unscaled (W = L D) until `_lt_finalize!`, so the step needs no scale phase. The trailing +# triangle of the block is shared over the workgroup; each work item owns rows `w+1:f` and forms their +# multipliers once. With `TRACK`, the pass-1 partials of the next column are reduced from the values +# written (each lane's entries of that column come in increasing row order, as in `_lt_pass1!`). +@inline function _lt_update!(fa, ctl, pv, red, redi, li, ::Val{WG}, ::Val{TRACK}, p::_LDLTDevice{R}, pk::Val, + h::Val) where {WG, TRACK, R} @inbounds begin k = Int(ctl[_LT_K]) w = Int(ctl[_ST_W]) + λ = zero(R) + r = 0 + cm = zero(R) if k <= w F = _lt_front(fa, ctl, pk) + f = Int(ctl[_ST_F]) step = Int(ctl[_LT_STEP]) kk = k + step - 1 - nr = Int(ctl[_ST_F]) - kk - nc = w - kk - for q in (li - 1):WG:(nr * nc - 1) - i = kk + 1 + q % nr - j = kk + 1 + q ÷ nr + j0 = kk + 1 # the next column + nt = w - kk + for q in (li - 1):WG:(nt * nt - 1) + i = kk + 1 + q % nt + j = kk + 1 + q ÷ nt i >= j || continue if step == 1 l = _fget(F, i, k) / pv[1] - _fset!(F, i, j, _fget(F, i, j) - l * _cj(_fget(F, j, k), h)) + v = _fget(F, i, j) - l * _cj(_fget(F, j, k), h) else x1 = _fget(F, i, k) x2 = _fget(F, i, k + 1) l1 = x1 * pv[1] + x2 * pv[3] l2 = x1 * pv[2] + x2 * pv[4] - _fset!(F, i, j, _fget(F, i, j) - (l1 * _cj(_fget(F, j, k), h) + l2 * _cj(_fget(F, j, k + 1), h))) + v = _fget(F, i, j) - (l1 * _cj(_fget(F, j, k), h) + l2 * _cj(_fget(F, j, k + 1), h)) + end + _fset!(F, i, j, v) + if TRACK && j == j0 && i > j0 + a = abs(v) + cm = max(cm, a) + if a > λ + λ = a + r = i + end + end + end + m = f - w + if m >= WG # one row per work item, its multipliers formed once + for i in (w + li):WG:f + if step == 1 + l = _fget(F, i, k) / pv[1] + for j in (kk + 1):w + v = _fget(F, i, j) - l * _cj(_fget(F, j, k), h) + _fset!(F, i, j, v) + TRACK && j == j0 && (cm = max(cm, abs(v))) + end + else + x1 = _fget(F, i, k) + x2 = _fget(F, i, k + 1) + l1 = x1 * pv[1] + x2 * pv[3] + l2 = x1 * pv[2] + x2 * pv[4] + for j in (kk + 1):w + v = _fget(F, i, j) - (l1 * _cj(_fget(F, j, k), h) + l2 * _cj(_fget(F, j, k + 1), h)) + _fset!(F, i, j, v) + TRACK && j == j0 && (cm = max(cm, abs(v))) + end + end + end + else # few rows: one entry per work item + for q in (li - 1):WG:(m * nt - 1) + i = w + 1 + q % m + j = kk + 1 + q ÷ m + if step == 1 + l = _fget(F, i, k) / pv[1] + v = _fget(F, i, j) - l * _cj(_fget(F, j, k), h) + else + x1 = _fget(F, i, k) + x2 = _fget(F, i, k + 1) + l1 = x1 * pv[1] + x2 * pv[3] + l2 = x1 * pv[2] + x2 * pv[4] + v = _fget(F, i, j) - (l1 * _cj(_fget(F, j, k), h) + l2 * _cj(_fget(F, j, k + 1), h)) + end + _fset!(F, i, j, v) + TRACK && j == j0 && (cm = max(cm, abs(v))) end end end + if TRACK + red[li] = λ + redi[li] = r % eltype(redi) + red[WG + li] = cm + end end return nothing end -# scale the pivot columns below the pivot block into L -@inline function _lt_scale!(fa, ctl, pv, li, ::Val{WG}, pk::Val) where {WG} +# scale the pivot columns into L, row by row (rows `2:f`, columns `1:min(i - 1, w)`): `l = x / d` (1×1) or +# `(x1 e11 + x2 e21, x1 e12 + x2 e22)` with the inverse 2×2 block recomputed from D exactly as +# `_lt_select!` formed it; the zero below a 2×2 block's diagonal stays +@inline function _lt_finalize!(fa, ctl, d, pivot_kind, li, ::Val{WG}, pk::Val, h::Val) where {WG} @inbounds begin - k = Int(ctl[_LT_K]) - if k <= Int(ctl[_ST_W]) - F = _lt_front(fa, ctl, pk) - step = Int(ctl[_LT_STEP]) - for i in (k + step + li - 1):WG:Int(ctl[_ST_F]) - if step == 1 - _fset!(F, i, k, _fget(F, i, k) / pv[1]) + n = length(pivot_kind) + F = _lt_front(fa, ctl, pk) + f = Int(ctl[_ST_F]) + w = Int(ctl[_ST_W]) + c0 = Int(ctl[_LT_C0]) + for i in (1 + li):WG:f + jmax = min(i - 1, w) + j = 1 + while j <= jmax + g = c0 + j - 1 + if pivot_kind[g] == PIVOT_KIND_2X2_FIRST + if j + 1 <= jmax + a, b, c2 = d[g], d[n + g], d[g + 1] + up = _cj(b, h) + det = a * c2 - up * b + x1 = _fget(F, i, j) + x2 = _fget(F, i, j + 1) + _fset!(F, i, j, x1 * (c2 / det) + x2 * (-b / det)) + _fset!(F, i, j + 1, x1 * (-up / det) + x2 * (a / det)) + end + j += 2 else - x1 = _fget(F, i, k) - x2 = _fget(F, i, k + 1) - _fset!(F, i, k, x1 * pv[1] + x2 * pv[3]) - _fset!(F, i, k + 1, x1 * pv[2] + x2 * pv[4]) + _fset!(F, i, j, _fget(F, i, j) / d[g]) + j += 1 end end end @@ -725,6 +806,109 @@ end return nothing end +# tiled F₂₂ ← F₂₂ − (L₂₁ D) L₂₁ᴴ of the regime-B/C kernel (`L₂₁` in the panel, after `_lt_finalize!`): the +# packed m×m block is cut into `_LT_TB × _LT_TB` tiles (lower ones only); per tile and chunk of `_LT_TB` +# columns, `W₂₁ = L₂₁ D` (1×1 and 2×2 blocks, as `_lt_cb_update!`) of the tile's rows and `L₂₁` of its columns +# are staged in the local buffer `tb` (`2 _LT_TB²` entries, then the `_LT_TB²` accumulators), and every work +# item accumulates `_LT_TB² / WG` entries over the chunk in column order +const _LT_TB = 16 +const _LT_TBUF = 3 * _LT_TB * _LT_TB + +# smallest contribution block taking the tiled update (smaller ones: `_lt_cb_update!`, one entry per work item) +const _LT_TILED_M = 64 + +@inline _lt_cb_ntiles(f, w, coff) = coff < 0 || f - w < _LT_TILED_M ? 0 : (nt = cld(f - w, _LT_TB); nt * (nt + 1) ÷ 2) +@inline _lt_cb_nchunks(ctl) = @inbounds cld(Int(ctl[_ST_W]), _LT_TB) + +# tile `tt` (column-major over the lower tiles) → (I, J), I ≥ J +@inline function _lt_cb_tile(ctl, tt) + nt = @inbounds cld(Int(ctl[_ST_F]) - Int(ctl[_ST_W]), _LT_TB) + J = 1 + t = tt + while t > nt - J + 1 + t -= nt - J + 1 + J += 1 + end + return (J + t - 1, J) +end + +# stage chunk `kc` of tile `tt`: tb[1:TB²] = W₂₁ rows of block I, tb[TB² .+ (1:TB²)] = L₂₁ rows of block J +@inline function _lt_cb_load!(tb, factor, ctl, d, pivot_kind, tt, kc, li, ::Val{WG}, h::Val) where {WG} + @inbounds begin + n = length(pivot_kind) + T = eltype(tb) + F = _lt_front(factor, ctl, Val(false)) + w = Int(ctl[_ST_W]) + m = Int(ctl[_ST_F]) - w + c0 = Int(ctl[_LT_C0]) + I, J = _lt_cb_tile(ctl, tt) + TB2 = _LT_TB * _LT_TB + for e in (li - 1):WG:(2 * TB2 - 1) + isw = e < TB2 + e2 = isw ? e : e - TB2 + ri = e2 % _LT_TB + k = (kc - 1) * _LT_TB + e2 ÷ _LT_TB + 1 + i = ((isw ? I : J) - 1) * _LT_TB + ri + 1 + v = zero(T) + if i <= m && k <= w + x = _fget(F, w + i, k) + if !isw + v = x + else + g = c0 + k - 1 + kd = pivot_kind[g] + if kd == PIVOT_KIND_2X2_FIRST + v = x * d[g] + _fget(F, w + i, k + 1) * d[n + g] + elseif kd == PIVOT_KIND_2X2_SECOND + v = _fget(F, w + i, k - 1) * _cj(d[n + g - 1], h) + x * d[g] + else + v = x * d[g] + end + end + end + tb[e + 1] = v + end + end + return nothing +end + +# accumulate the staged chunk into the tile accumulators tb[2TB² .+ (1:TB²)] +@inline function _lt_cb_acc!(tb, ctl, kc, li, ::Val{WG}, h::Val) where {WG} + @inbounds begin + TB2 = _LT_TB * _LT_TB + nk = min(_LT_TB, Int(ctl[_ST_W]) - (kc - 1) * _LT_TB) + for e in (li - 1):WG:(TB2 - 1) + r = e % _LT_TB + c = e ÷ _LT_TB + acc = kc == 1 ? zero(eltype(tb)) : tb[2 * TB2 + e + 1] + for q in 0:(nk - 1) + acc += tb[r + q * _LT_TB + 1] * _cj(tb[TB2 + c + q * _LT_TB + 1], h) + end + tb[2 * TB2 + e + 1] = acc + end + end + return nothing +end + +# subtract the tile from the packed block (real diagonal for Hermitian) +@inline function _lt_cb_store!(tb, ctl, ca, coff, tt, li, ::Val{WG}, ::Val{H}) where {WG, H} + @inbounds begin + T = eltype(tb) + m = Int(ctl[_ST_F]) - Int(ctl[_ST_W]) + I, J = _lt_cb_tile(ctl, tt) + TB2 = _LT_TB * _LT_TB + for e in (li - 1):WG:(TB2 - 1) + i = (I - 1) * _LT_TB + e % _LT_TB + 1 + j = (J - 1) * _LT_TB + e ÷ _LT_TB + 1 + (i <= m && j <= m && i >= j) || continue + acc = tb[2 * TB2 + e + 1] + pos = coff + _packed(i, j, m) + ca[pos] = (H && i == j) ? T(real(ca[pos]) - real(acc)) : ca[pos] - acc + end + end + return nothing +end + # work item 1: per-front statistics and status @inline function _lt_stats!(stats, info, ctl, s) @inbounds begin @@ -741,7 +925,7 @@ end # --------------------------------------------------------------------------- # regimes B/C: one workgroup per front, panel in global memory -@inline function _lt_front_setup!(ctl, s, super_ptr, front_ptr, front_nrows, front_ncols) +@inline function _lt_front_setup!(ctl, s, super_ptr, front_ptr, front_nrows, front_ncols, cb_ptr) @inbounds begin IT = eltype(ctl) ctl[_ST_NODE] = s % IT @@ -749,6 +933,7 @@ end ctl[_ST_W] = front_ncols[s] % IT ctl[_ST_LF] = front_ptr[s] % IT _lt_reset!(ctl, super_ptr[s]) + ctl[_LT_NTILE] = _lt_cb_ntiles(Int(front_nrows[s]), Int(front_ncols[s]), Int(cb_ptr[s]) - 1) % IT end return nothing end @@ -796,11 +981,11 @@ statistics, updates its packed contribution block on the update stack pv = @localmem TT (_LT_NPV,) red = @localmem R (3 * WG,) redi = @localmem IT (2 * WG,) - L11 = @localmem TT (max(W * (W + 1) ÷ 2, 1),) + L11 = @localmem TT (max(W * (W + 1) ÷ 2, _LT_TBUF),) # F₁₁, then the tiles of the F₂₂ update if li == 1 @inbounds s = nodes[bm.first + _bm_node(bm, G) - 1] _lt_front_setup!(ctl, s, super_ptr, member_panels(front_ptr, _bm_gmember(bm, G), bm.nbatch), front_nrows, - front_ncols) + front_ncols, cb_ptr) end @synchronize @inbounds s = nodes[bm.first + _bm_node(bm, G) - 1] @@ -871,16 +1056,33 @@ statistics, updates its packed contribution block on the update stack end _lt_swap!(_lt_fa(L11, factor, Val(W)), ctl, pv, _mview(piv, _bm_gmember(bm, G), bm.nbatch), li, Val(WG), _lt_pk(Val(W)), Val(HERM)) @synchronize - _lt_update!(_lt_fa(L11, factor, Val(W)), ctl, pv, li, Val(WG), _lt_pk(Val(W)), Val(HERM)) - @synchronize - _lt_scale!(_lt_fa(L11, factor, Val(W)), ctl, pv, li, Val(WG), _lt_pk(Val(W))) - _lt_pass1!(_lt_fa(L11, factor, Val(W)), ctl, red, redi, li, Val(WG), prm, _lt_pk(Val(W))) + _lt_update!(_lt_fa(L11, factor, Val(W)), ctl, pv, red, redi, li, Val(WG), Val(true), prm, _lt_pk(Val(W)), + Val(HERM)) @synchronize end + @inbounds s = nodes[bm.first + _bm_node(bm, G) - 1] + k = _bm_gmember(bm, G) + nb = bm.nbatch + _lt_finalize!(_lt_fa(L11, factor, Val(W)), ctl, _mview(d, k, nb), _mview(pivot_kind, k, nb), li, Val(WG), + _lt_pk(Val(W)), Val(HERM)) + @synchronize _lt_stage!(L11, factor, ctl, li, Val(W), Val(WG), false) @synchronize + for tt in 1:ctl[_LT_NTILE] + for kc in 1:_lt_cb_nchunks(ctl) + k = _bm_gmember(bm, G) + _lt_cb_load!(L11, factor, ctl, _mview(d, k, bm.nbatch), _mview(pivot_kind, k, bm.nbatch), tt, kc, li, + Val(WG), Val(HERM)) + @synchronize + _lt_cb_acc!(L11, ctl, kc, li, Val(WG), Val(HERM)) + @synchronize + end + @inbounds s = nodes[bm.first + _bm_node(bm, G) - 1] + @inbounds coff = Int(cb_ptr[s]) - 1 + _lt_cb_store!(L11, ctl, _mview(stack, _bm_gmember(bm, G), bm.nbatch), coff, tt, li, Val(WG), Val(HERM)) + end @inbounds s = nodes[bm.first + _bm_node(bm, G) - 1] - @inbounds coff = Int(cb_ptr[s]) - 1 + @inbounds coff = ctl[_LT_NTILE] == 0 ? Int(cb_ptr[s]) - 1 : -1 k = _bm_gmember(bm, G) nb = bm.nbatch _lt_cb_update!(factor, ctl, _mview(stack, k, nb), coff, _mview(d, k, nb), _mview(pivot_kind, k, nb), li, Val(WG), @@ -1014,13 +1216,16 @@ stack. end _lt_swap!(buf, ctl, pv, _mview(piv, _bm_gmember(bm, G), bm.nbatch), li, Val(WG), Val(true), Val(HERM)) @synchronize - _lt_update!(buf, ctl, pv, li, Val(WG), Val(true), Val(HERM)) + _lt_update!(buf, ctl, pv, red, redi, li, Val(WG), Val(false), prm, Val(true), Val(HERM)) @synchronize - _lt_scale!(buf, ctl, pv, li, Val(WG), Val(true)) _lt_pass1!(buf, ctl, red, redi, li, Val(_LT_S1), prm, Val(true)) @synchronize end mb = _bm_gmember(bm, G) + _lt_finalize!(buf, ctl, _mview(d, mb, bm.nbatch), _mview(pivot_kind, mb, bm.nbatch), li, Val(WG), Val(true), + Val(HERM)) + @synchronize + mb = _bm_gmember(bm, G) _lt_cb_update!(buf, ctl, buf, _lt_local_cb_offset(ctl), _mview(d, mb, bm.nbatch), _mview(pivot_kind, mb, bm.nbatch), li, Val(WG), Val(true), Val(HERM)) @synchronize From 886c9d898f827e45ce8f344ea6ee641fdde39ab6 Mon Sep 17 00:00:00 2001 From: Michel Schanen Date: Mon, 5 Oct 2026 14:24:30 -0500 Subject: [PATCH 3/4] perf: narrow regime-B bins keep the panel in global memory (#75) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Width classes 8 and 16 run front_ldlt_kernel! without staging F₁₁: the load and write-back phases cost more than they save on KKT fronts (case14 condensed 0.55 → 0.51 ms, case1354 condensed 7.6 → 7.4 ms). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01Cqp7TUSyCEdjPtLDEftJf5 --- src/numeric/ldlt.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/numeric/ldlt.jl b/src/numeric/ldlt.jl index b73b409..3602afd 100644 --- a/src/numeric/ldlt.jl +++ b/src/numeric/ldlt.jl @@ -814,6 +814,10 @@ end const _LT_TB = 16 const _LT_TBUF = 3 * _LT_TB * _LT_TB +# widest regime-B width class that keeps its panel in global memory (staging F₁₁ costs more than it saves on +# narrow fronts: two more phases per front) +const _LT_GLOBAL_MAX_W = 16 + # smallest contribution block taking the tiled update (smaller ones: `_lt_cb_update!`, one entry per work item) const _LT_TILED_M = 64 @@ -1393,7 +1397,7 @@ 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 == 0 # regime C: the panel in global memory + if W <= _LT_GLOBAL_MAX_W # regime C and narrow 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 From ef5c6f545506498a5cefe883378251042edb7d91 Mon Sep 17 00:00:00 2001 From: Michel Schanen Date: Mon, 5 Oct 2026 18:06:57 -0500 Subject: [PATCH 4/4] perf: experiment 1 step 2 results in PERFORMANCE.md (#75) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01Cqp7TUSyCEdjPtLDEftJf5 --- PERFORMANCE.md | 36 +++++++++++++++++++++++++++++++++++- 1 file changed, 35 insertions(+), 1 deletion(-) diff --git a/PERFORMANCE.md b/PERFORMANCE.md index 54df883..90508af 100644 --- a/PERFORMANCE.md +++ b/PERFORMANCE.md @@ -10,7 +10,7 @@ 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) in a PR from `perf/exp1-pivot-search` | +| #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` | | #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 | @@ -147,6 +147,40 @@ LDLᵀ refactorization, `main` → step 1: On the K2 dumps the inertia and `nperturbed` are those of `main` on all nine. On the T15 matrix `kkt_matrix(3000, 1000, 1e-8)` (fallback scans on most steps of its two large root fronts) CUDA goes from 33.9 s to 8.3 s; the KA CPU backend stays at 7.0 s against 2.6 s for `ref_ldlt!` (it runs the work items of a workgroup one after another, so a parallel search saves nothing there). The medium matrices gain 16–34%; the smallest KKT dumps (case14, w ≤ 16) lose up to 17% to the reductions' barriers; the two matrices with very large regime-C fronts (apache2, lap3d_40) are within ±5%, still one workgroup per front. +### Step 2: regime B in local memory + +Staging only F₁₁ in `@localmem` (as `front_cholesky_kernel!`) gained nothing (lap2d_300 190 → 192 ms): the rows below F₁₁ are updated in global memory at every step and the contribution-block update was two thirds of regime B. The step became: `front_ldlt_kernel!` takes the width class `W` of its launch group and keeps the packed `W×W` F₁₁ in local memory from `W = 32` (classes 8 and 16 keep the panel in global memory: the staging phases cost more than they save on KKT fronts); the pivot columns stay unscaled (`L D`) until the end of the front, where `_lt_finalize!` applies the reference's divisions and 2×2 transforms with the pivot recomputed from D, so a step has no scale phase; the update reduces the next column as it writes it (pass 1), and work items own the rows below F₁₁ when there are at least `WG` of them (one division per row); four barriers per step. Contribution blocks with `m ≥ 64` are updated in 16×16 tiles with `W₂₁ = L₂₁D` and `L₂₁` staged 16 columns at a time in the F₁₁ buffer. Regime-C fronts run the same kernel with the panel in global memory and gain from the same changes. + +LDLᵀ refactorization, `main` → step 1 → step 2: + +| 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: 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. + ## Recommendations: Prioritized Experiment Plan | # | Experiment | Payoff | Effort | Key measurement | Success criterion |