Skip to content

In-front Bunch-Kaufman pivoting in the fused pipeline - #139

Draft
sshin23 wants to merge 11 commits into
mainfrom
fused-bk-factor
Draft

sshin23 wants to merge 11 commits into
mainfrom
fused-bk-factor

Conversation

@sshin23

@sshin23 sshin23 commented Oct 14, 2026 •

Copy link
Copy Markdown
Contributor

Draft: in-front Bunch-Kaufman (BK) pivoting for the fused pipeline, so indefinite KKT systems get real pivoting at fused speed. Supersedes #135 (its commit is the first one here); #135 is left open pending the owner's call.

What

commit change
5947a93 the fused solve applies the within-front pivot order and a 2x2-aware D pass (it assumed identity order, so a pivoted factor solved silently wrong)
ebae400 build_fused_plan(; bk = true): the stock pivot machinery inside the fused kernels (pivoted subtree launches + the BK step machine in the mega kernel's PANEL roles)
09ee43f wide fronts batched per segment; BK contribution blocks via 2x2-aware TILE roles
71ed444 vendor-sytrf wide fronts, and the hybrid (BK narrow + unpivoted wides) mode
ff418fd MadNLP extension: sds_pivoting = :bk / :bk_hybrid, sds_pivot_threshold, sds_pivot_pairs; inertia from d/pivot_kind with 2x2 blocks
52d7b75 hybrid on merged in-kernel wides; subtree pivot search widened to the full workgroup
3765cde an optimistic certified fast path, implemented, measured and rejected (see below)

Measured (78k condensed ACOPF KKT, GV100, per refactorization)

path time
no-pivot fused LDL (the floor) 0.054 s
:bk_hybrid (BK subtrees + narrow fronts, unpivoted wides) 0.093 s
:bk (full BK, vendor-sytrf wides) 0.257 s
stock (non-fused) BK, what the ext routed to before 0.450 s
cuDSS LDL, same GPU and matrix 0.030 s

End-to-end on the INFORMS KKT benchmark (66-case ACOPF, L40S): condensed+SDS goes 60 -> 64 solved with :bk_hybrid, which equals condensed+cuDSS (64) with the same two failures (9241_pegase, 1888_rte). A 5x repeat of :bk vs :bk_hybrid on the two decisive cases came out identical, so wide-front pivoting is not load-bearing on this workload; :bk_hybrid + sds_pivot_pairs = "none" is the recommended configuration.

Two negative results worth recording

  • Wide-front in-front BK does not honor its own threshold bound. Pivoting is restricted to the fully-summed block, so when no candidate passes the test against the full column the fallback accepts a sub-threshold pivot. Measured on the 78k: the stock BK factor's wide-front multipliers reach 4.6e4, larger than the unpivoted path's 3.9e4, with 53 of 58 wide fronts over the nominal bound. An optimistic fast path certified by |l| <= 1/u was therefore reverted in 3765cde (the measurement kernel is kept as a diagnostic).
  • cuSOLVER's sytrf stores a different factor form than LAPACK: the same ipiv values, but the product form prod_k (P_k L_k) -- previously computed L columns are not row-swapped. LAPACK-style reconstruction fails at 0.97 relative error where the product form reconstructs to 3e-16; the conversion kernel handles it.

Validation

  • Narrow-front and subtree BK produce a pivot sequence byte-identical to the stock path (two gate fixtures: a grid-kron matrix with 332 2x2 pivots, and a saddle matrix with a 1195-wide root front); sytrf wides take LAPACK's sequence by design and are gated on exact inertia plus solve accuracy.
  • Full suite (CPU backend + CUDA, all four element types) 63913/63913; new test_solve.jl testset for the pivot-aware fused solve.
  • Extension smoke on a GV100: both modes info = 0, inertia equal to the eigenvalue counts, solve 6e-15.

Optimization round (commit 009707d)

Measured interleaved in one process against a same-shape no-pivot control (the earlier
1.71x figure in this thread compared against a differently shaped plan; like-for-like the
merged single-launch no-pivot floor is 0.0475 s, so the penalty was 1.91x):

time vs no-pivot
no-pivot fused LDL, merged plan 0.0475 s 1.00x
:bk_hybrid before this round 0.0910 s 1.91x
:bk_hybrid after 0.0874 s 1.84x

Two wins: the pivoted subtree kernel takes the next column's maxima out of the update
itself (the TRACK path front_ldlt_kernel! already used, so the separate per-step scan and
its barrier are gone), and when the chosen pivot already sits at k only the pivot block's
diagonal changes, which one work item writes while the rest of the workgroup skips the
interchange; plus BK panels take a 256-thread workgroup, since a pivoted front is factored
by ONE workgroup over all its rows where no-pivot spreads the rows below the block across
TRSM/TILE roles (0.0874 s at 256 vs 0.0906 at 128 and 0.0902 at 512).

A barrier may not be made runtime-conditional: @synchronize in one arm of an if/else
deadlocks the KA kernels, because threads taking the other arm never arrive. Caught as a
hung test_numeric_ldlt on CUDA; the work is conditional, the barrier is not.

Parameter levers that did not pay, recorded so they are not re-tried: widening the
subtree pivot search (neutral), 512-thread panels (no better than 256), and narrowing the
subtree workgroup to 64 or 32 (worse -- 128 is already optimal).

Attribution of the remaining +0.040 s: mega-kernel BK panels +0.027 s, pivoted subtree
launches +0.017 s. What is left is the per-pivot-step serial chain (cooperative search ->
single-threaded decide -> update, w times per front), which is inherent to threshold
pivoting rather than a tuning problem.

🤖 Generated with Claude Code

sshin23 and others added 11 commits October 13, 2026 21:03
The fused solve (algo2_solve! and its kernels fwd1/bwd1, the fwd2/bwd2 counter
sweep, the chain walkers, and the dense wide path) assumed identity pivot order
within every front, so a factor with in-front pivoting solved silently wrong.
Now: a compile-time Val switch (one dynamic dispatch per solve) folds the local
pivot order piv into the forward reads (z = P'y on the way into inv(L11)) and
the backward writes (x scattered back through piv), exactly _sv_piv_apply! of
the sweep solve; the dense wide fronts go through the sweep machinery with the
true solve kind (which applies piv itself); and the D pass goes through
diagonal_sweep! (1x1 and 2x2 pivot blocks) instead of the elementwise divide.
The no-pivot and Cholesky instantiations are unchanged (Val(false): _ppos is
the identity and piv is never read); pivoted mode engages only for an LDL
structure with pivot_type != 'N'.

New testset (test_solve.jl): on a stock Bunch-Kaufman factor of a kkt_matrix
with forced 2x2 pivots (asserted nonvacuous: n2x2 > 0, nothing perturbed) and
of a symindef with a front wider than MW (dense path), every fused solve
configuration matches the stock sweep solve. GPU-only: the fused solve kernels
keep workgroup-uniform front state in locals across @synchronize, which the KA
CPU barrier transform does not carry over (pre-existing, by design: the counter
sweep needs a resident grid).

Gated on a GV100: the new testset green in test_solve (508/508), full suite
(CPU + CUDA, all four element types) 63913/63913 with the kernels changed;
oracle script relerrs 1.8e-11 / 2.4e-11 / 9.0e-12 (level/counter/chain) vs
stock 2.1e-11 on n=3400 with 210 2x2 pivots and a 1195-wide dense front.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…= 64)

build_fused_plan(s; ldl = true, bk = true): the sorted subtree launches run the
stock pivoted subtree_ldlt_kernel! (same lists, LDLT local reserve), and the
mega kernel's PANEL roles run the stock regime-B pivot machinery verbatim on
the split front (F11 staged in sh, rows below in the factor panel): cooperative
search, interchanges, D/piv/pivot_kind/statistics, and the stock _lt_cb_*
contribution-block update into the private CB, which handles 2x2 pivot blocks.
No TRSM/TILE roles and no write-through for BK fronts (their trailing updates
assume a scalar d per column); the PANEL workgroup factors the whole front and
signals cb_done/arrived itself. Pivot parameters resolve per refactorization
(_ldlt_device_params + abs_max!), so setparam! takes effect without a plan
rebuild. Fronts wider than 64 and uniform batches raise NotSupportedError (the
blocked-path port is the follow-up); Cholesky and no-pivot LDL instantiations
are unchanged (Val(false) placeholders; the BK localmem collapses to 1 entry).

Gated on a GV100 (grid kron [0 1; 1 0] fixture, 332 2x2 pivots, 152 mega-kernel
fronts + 13 subtrees): pivot sequence and kinds EXACTLY equal to the stock BK
path, factor/d max rel diff 1.2e-15/2.2e-16, counter-sweep fused solve 3.6e-15
vs stock 3.0e-15. Full suite (CPU + CUDA, all four element types) 63913/63913.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…ILE roles

Wide fronts of a segment now go through one panel_ldlt_kernel! launch per
32-column block over a COMPACTED front list chunked by workspace at plan time
(the kernel's _ltc_slice walks the same list, so host offsets and kernel slices
agree), instead of one serial launch per front per block. BK PANEL fronts with
m > SYRK_IN hand their contribution block to TILE roles that fold W = L D with
1x1 and 2x2 pivot blocks at staging time (the 2x2 partner column read straight
from the panel, block-straddling pairs included) and spin on panel_done; small
CBs stay on the direct _lt_cb_update! in the PANEL workgroup. No TRSM roles in
BK mode (the pivot steps update every row below w already).

SDS_BK_TIMING=1 prints per-stage refactorization timing (dev instrumentation).

78k condensed on a GV100: 0.450s stock -> 0.328s fused BK; the decomposition is
subtrees 0.030s, mega kernel 0.064s, host wides 0.231s. The wide cost is the
inherently serial pivot stepping on the root chain (one workgroup per front,
~6-10 barriers per step); both gate fixtures stay pivot-sequence-identical to
the stock path (GATE OK), full fixture suite re-run green.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…wides) mode

Two ways out of the wide-front pivot-step wall (the serial blocked panel costs
0.20s of the 0.29s refactorization on the 78k condensed KKT):

- wides_sytrf (default on a vendor_sytrf-capable backend): cuSOLVER/LAPACK
  Bunch-Kaufman on the dense w x w fully-summed block, then a conversion to the
  stock conventions. MEASURED: cuSOLVER emits LAPACK's ipiv VALUES but stores
  the PRODUCT form prod_k (P_k L_k) -- it does not row-swap previously computed
  L columns (LAPACK-style reconstruction fails at 0.97 rel err, product form
  reconstructs to 3e-16), so the conversion kernel walks ipiv applying each
  interchange to the prior columns cooperatively before extracting D/piv/kinds
  and scrubbing to unit-L. L21 is gathered into pivot order, TRSM'd (W = L21 D),
  divided by the D blocks race-free from the W copy, and the CB is one GEMM.
  Wide pivot sequences are LAPACK's (alpha rule), not the stock kernels': gated
  on inertia (exact) and solve accuracy instead, as wides_potrf trades for chol.
- bk + wides_potrf = the HYBRID: BK on subtrees and narrow fronts, the verbatim
  potrf-derived path on wides (a sign flip there fails -> MadNLP regularizes,
  the no-pivot behavior); identity piv/1x1 kinds keep the pivot-aware solve and
  2x2 D pass consistent.

78k condensed, GV100, per refactorization: stock BK 0.450s, fused BK blocked
0.293s, fused BK sytrf 0.257s (the cusolver call itself is 0.157s), HYBRID
0.1245s, no-pivot fused 0.054s, cuDSS LDL 0.030s. Gate fixtures: narrow path
still pivot-sequence-identical to stock; sytrf wides solve 2.7e-8 vs stock
3.6e-8 with exact inertia (2195, 1205).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The stock routing of :bk is gone (the fused pipeline is the main path): :bk
builds the fused plan with bk = true (stock pivot machinery in the subtree and
PANEL kernels, vendor-sytrf wide fronts) and :bk_hybrid adds wides_potrf (BK on
subtrees and narrow fronts, the potrf-derived path on wides: a sign flip there
fails -> regularization, the :none behavior). Both solve through the pivot-aware
fused solve; factorize! is refact_fused! for every mode. True inertia comes from
a device mapreduce over d and pivot_kind (2x2 blocks by det/trace sign), since
cholesky_stats! leaves the per-front stats Cholesky-flavored.

sds_pivot_pairs (new, default nothing = per-mode): analysis-time 2x2 pairs are
LOAD-BEARING for in-front BK on saddle structures -- they place zero-diagonal
columns adjacently in one front so the kernels can pair them; without them the
n=60 saddle perturbs 13 pivots and the solve degrades to 5e-2, stock and fused
alike. BK modes therefore default to "default" pairs; :none keeps "none"
(the detection costs ~42 s of the 78k analysis and no-pivot never consumes it).

78k condensed, GV100, per refactorization: :none 0.054s, :bk_hybrid 0.1245s,
:bk 0.257s, old stock routing 0.450s, cuDSS LDL 0.030s. GV100 ext smoke: both
BK modes info 0, inertia == eigenvalue counts, solve 6e-15; threshold knob live.
Fixture gates GATE OK; full suite (CPU + CUDA, all four element types)
63913/63913.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…arch

bk + inkernel_wides (hybrid-v2) runs Bunch-Kaufman on subtrees and narrow
fronts with the in-kernel no-pivot LDL W-roles on wide fronts, in one merged
launch: 0.093s per refactorization on the 78k condensed KKT (potrf-wides hybrid
0.12s, no-pivot 0.054s, full BK 0.26s; GV100). In BK mode piv/pivot_kind are
initialized identity/1x1 for every column up front, so unpivoted wide fronts
keep solve-consistent state and the per-front init in the potrf branch is gone.

The stock subtree_ldlt_kernel! pivot search now scans with the full workgroup
(Val(WG) passes + the _lt_stage1! combine, exactly the regime-B kernel) instead
of 16 lanes; SUBTREE_LOCAL_RESERVE_LDLT grows to 4608 for the WG-sized red/redi.
The reduction keeps the reference tie-breaking (regime B is gated against the
reference with exact pivot sequences, and the LDLT suites pass unchanged:
2828/2828 targeted, 63913/63913 full). No effect on the tiny subtree fronts of
the 78k (their cost is the per-step barriers, not the scan); it pays on wide
subtree fronts.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…cted on evidence

:bk_hybrid now builds the merged in-kernel-wides plan (0.093s per 78k
refactorization, was 0.12s with potrf wides): wide fronts factor through sign
flips with true inertia from d instead of failing potrf into regularization.

An optimistic full-:bk (unpivoted wides + a post-hoc |l| <= 1/u certification,
retry with real pivoting on failure) was implemented, measured, and REVERTED:
the criterion is not an invariant of the real path. In-front BK can only pivot
within the fully-summed block, and on wide fronts where no candidate passes the
threshold against the full column the fallback accepts a sub-threshold pivot --
measured on the 78k, the stock BK factor's wide multipliers reach 4.6e4, LARGER
than the unpivoted path's 3.9e4, with 53 of 58 wide fronts over the nominal
bound. So no sound certification separates the paths, and the retry would fire
on virtually every real iteration (paying fast + slow). The measurement kernel
stays as a diagnostic (wide_multipliers_ok, _wide_lmax_kernel!, fp.wall/lmax).

GV100 ext smoke green for :bk and :bk_hybrid (6e-15, exact inertia); full suite
63913/63913.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…de BK panels

Two measured wins, 0.0910 -> 0.0874s per refactorization on the 78k condensed
ACOPF KKT (GV100, interleaved in-process against a same-shape no-pivot control):

- the pivoted subtree kernel now takes the next column's maxima out of the
  update itself (the TRACK path front_ldlt_kernel! already used), so the
  separate per-step pass-1 scan and its barrier are gone; and when the chosen
  pivot already sits at k, only the pivot block's diagonal changes, which work
  item 1 writes while the rest of the workgroup skips the interchange.
- BK panels take a 256-thread workgroup: a pivoted front is factored by ONE
  workgroup over all its rows, where no-pivot spreads the rows below the block
  across TRSM/TILE roles, so the tuning bundle's no-pivot width is too narrow
  here (0.0874s at 256 vs 0.0906 at 128, 0.0902 at 512).

A barrier may NOT be made runtime-conditional: @synchronize in one arm of an
if/else deadlocks the KA kernels (threads taking the other arm never arrive) --
measured as a hung test_numeric_ldlt on CUDA. The work is conditional, the
barrier is not; the stage-1 guard stays the original compile-time one.

Honest baseline correction: the penalty over no-pivot is 1.84x, not the 1.71x
reported earlier, which compared against a differently shaped plan (merged
single-launch no-pivot is 0.0475s, not 0.0532s). Attribution of the remainder:
mega-kernel BK panels +0.027s, pivoted subtree launches +0.017s.

Gates: reference-exact LDL suites 1868/1868 (pivot sequences unchanged), both
oracle fixtures GATE OK, extension smoke OK, full suite 63913/63913.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…diag)

Per column take the largest remaining diagonal of the fully-summed block and
permute it in symmetrically; never a 2x2. Because only 1x1 pivots occur, the
factor keeps the no-pivot SHAPE plus a within-front permutation, so every
no-pivot role applies unchanged (TRSM, TILE, in-kernel wides) and -- unlike
Bunch-Kaufman -- the selection reads nothing below the fully-summed block, so
those rows are never touched during the column loop. Per front that is O(f w) of
interchange traffic against BK's O(f w^2) of per-step updates.

MEASURED (78k condensed ACOPF KKT, GV100, interleaved in-process, ratios only --
another process held the GPU, so the absolutes in that run are inflated ~1.5x
uniformly): no-pivot 1.00x, diag 1.04x, BK 1.84x. The 78k solve residual is
0.0213 against no-pivot's 0.0184 and BK's 0.0173, i.e. the same class.

The mode is weaker than BK by construction and this is measured, not assumed: a
front whose remaining diagonals are ALL zero has nothing to select, so on an
unregularized saddle (zero 2,2 block) diag fails at the same column as no-pivot
where BK's 2x2 succeeds. It therefore sits between :none and :bk in robustness,
and which one a workload wants is an end-to-end question, not a speed one.

Implementation note worth keeping: the no-pivot subtree update writes
piv[g] = g on every column, which silently erased the permutation the selection
had just recorded -- the factor was permuted while the solve assumed it was not
(78k residual 1.1e31 with info = 0). It is now skipped under DIAG, which
maintains the order itself from the identity initialization.

Gates: extension smoke OK for :diag / :bk / :bk_hybrid (diag's honest failure on
a zero-diagonal block asserted as the contract, not accepted as success), full
suite (CPU + CUDA, all four element types) 63913/63913.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Sign-aware static perturbation for the no-pivot and :diag paths (sds_stabilize,
default off): a pivot below pivot_epsilon takes +/-epsilon with the pivot_sign
policy instead of failing the front, which is what the :bk path's _lt_select!
has always done -- so the three pivoting modes can share one stabilization
policy and differ only in pivoting. nperturbed(fp) counts the perturbations, so
a cell cannot be reported as stabilized when the policy never fired.

Also upstreams the benchmark branch's n_primal / n_dual / sds_pivot_signs onto
the extension, so that block stops conflicting between checkouts.

Targeted gate (GV100, 450-dim saddle with a zero 2,2 block where both 1x1 paths
fail at column 423): :none and :diag flip to info = 0 with 146 and 126 pivots
perturbed. NOTE the standalone solve stays poor (relerr ~2-3): perturbing zero
diagonals to +/-1e-13 completes the factorization while leaving it very
ill-conditioned, so the value is that MadNLP receives info = 0 and an inertia to
regularize against, NOT standalone accuracy. Cells must be labelled that way.

NOT yet regression-gated: the full suite run was cut short when shared-node jobs
had to be killed, and ORCD is the sanctioned place to re-run it (login-node
compute is forbidden, submit only), which currently needs Duo re-auth.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The device edit (_lt_select!) and this reference are compared on the exact pivot
sequence and factor, so the CPU reference needs the same perturbation entry for
entry: when |det| < eps * max(|a|,|b|,|c2|), move whichever diagonal is SMALLER
(dividing by the larger) so det becomes -eps * M, keeping it negative.

Justified on INVERSE-ENTRY MAGNITUDE, and the comment records that two
plausible-looking justifications are false. A pair chosen by Bunch-Kaufman
satisfies akk*sigma < alpha*lambda^2 and |c2| < alpha*sigma, hence
|a*c2| <= alpha^2 lambda^2 < lambda^2 = b^2, so det < 0 STRICTLY with
|det| >= (1 - alpha^2) lambda^2: the determinant can neither vanish nor round to
the wrong sign for such a pair, and the (1 pos, 1 neg) inertia credit is never
wrong for it. Verified over 2.3e5 sampled BK choices (min |det|/lambda^2 = 0.72)
and independently from the inequalities themselves (min 0.5918 against the
0.589903 bound). What IS reachable is a tiny diagonal a with a partner c2 far above
lambda: healthy |det|/lambda^2, and still |c2/det| above 1/eps.

The test fixture is DERIVED by calling the reference's own predicates, not
guessed -- three hand-built attempts failed because reaching the fallback also
needs NO acceptable 1x1 anywhere in the block (_best_1x1 short-circuits before
the nonsingular guard), which forces |c2| < u*sigma rather than just < alpha*sigma.
The test asserts the path is taken (choice, acceptance, best_1x1, guard) so it
cannot pass vacuously, checks the bump restores M/|det| <= 1/eps while keeping
the pair indefinite, and checks the eps = 0 configuration leaves the block
untouched.

test_reference_ldlt 975/975 (CPU).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant