SupernodalLU: hand-rolled panel solve kernels below a size cutoff (−27 to −42% solve)#1112
Merged
ChrisRackauckas merged 1 commit intoJul 25, 2026
Conversation
…7 to -42%) The solve sweeps issued BLAS trsv/gemv per supernode panel, but the panels are tiny: on a 2D Poisson factorization the pivot block width np has mean ~12 and median 6, so a 12x12 unit-triangular solve (~72 flops) measured ~260 ns through BLAS - almost entirely call overhead. BLAS thread count is irrelevant at these sizes (123.1 us at 1 thread vs 122.8 us at 64). Below PANEL_BLAS_CUTOFF (64) the four panel operations now use column-oriented @inbounds/@simd kernels written out in-package (_unit_lower_solve!, _upper_solve!, _panel_gemv!, _panel_gemv_sub!); at or above the cutoff BLAS is used unchanged. This follows PureKLU's SciML#38/SciML#39/SciML#40, where hand-written per-column kernels beat generic library dispatch at these sizes. Measured single-RHS solve, BLAS pinned to 1 thread: poisson2d_20 n=400 27.92 -> 16.18 us (-42%) poisson2d_40 n=1600 123.15 -> 75.92 us (-38%) poisson2d_60 n=3600 286.33 -> 187.43 us (-35%) poisson2d_100 n=10000 843.54 -> 568.36 us (-33%) poisson3d_16 n=4096 491.79 -> 358.23 us (-27%) That moves solve from ~1.8x behind UMFPACK to ahead of it at scale (568 vs 703 us at n=10000; 187 vs 200 us at n=3600). Accuracy is unchanged (relative error vs backslash 3.8e-16 at n=400 to 4.9e-14 at n=10000, same as the BLAS path). The new testset checks single- and multi-RHS agreement and asserts that a 3D factorization actually has panels wider than the cutoff, so both branches stay exercised. GROUP=Core passes. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This was referenced Jul 25, 2026
ChrisRackauckas
marked this pull request as ready for review
July 25, 2026 13:39
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
The largest finding from auditing the vendored
src/SupernodalLUsolve path against PureKLU.jl's optimization history. Independent of #1109 and #1111; touchessolve.jland its test only.The problem
The solve sweeps issued BLAS
trsv/gemvper supernode panel — but the panels are tiny. On a 2D Poisson factorization the pivot block widthnphas mean ~12, median 6, so a 12×12 unit-triangular solve (~72 flops) measured ~260 ns through BLAS: almost entirely call overhead. Confirming it's dispatch and not arithmetic, BLAS thread count makes no difference at these sizes (123.1 µs at 1 thread vs 122.8 µs at 64)._solve_panels!was 96 % of solve time and its four BLAS calls were 96 % of that.The change
Below
PANEL_BLAS_CUTOFF(64) the four panel operations use column-oriented@inbounds/@simdkernels written out in-package (_unit_lower_solve!,_upper_solve!,_panel_gemv!,_panel_gemv_sub!); at or above the cutoff BLAS is used unchanged. Same approach as PureKLU's #38/#39/#40.Measured (single-RHS solve, BLAS pinned to 1 thread)
That moves solve from ~1.8× behind UMFPACK to ahead of it at scale: 568 µs vs UMFPACK's 703 µs at n=10 000, and 187 vs 200 µs at n=3 600.
Accuracy and coverage
Unchanged — relative error vs
\is 3.8e-16 (n=400) to 4.9e-14 (n=10 000), matching the BLAS path. The new testset checks single- and multi-RHS agreement against per-column solves, and asserts that a 3D factorization actually has panels wider than the cutoff so both branches stay exercised rather than one silently rotting.A note on method, from PureKLU's #23: isolated microbenchmarks of these loops are misleading (an
@simd ivdepvariant looked 1.8× better in isolation but was −0.5 % noise end-to-end). Every number above is an end-to-end A/B of the same binary with only the cutoff changed.GROUP=Corepasses.🤖 Generated with Claude Code
https://claude.ai/code/session_017rr42T1vFRRT8D7DMAmJzf