Replace the vector-RHS BLAS deferral with a rank-4 pure-Julia sweep (never-BLAS vectors) by ChrisRackauckas-Claude · Pull Request #48 · JuliaSIMD/TriangularSolve.jl

ChrisRackauckas-Claude · 2026-08-08T20:16:28Z

⚠️ Draft — please ignore until reviewed by @ChrisRackauckas.

What changed and why

Vector right-hand sides now run a pure-Julia rank-4 column sweep (_naive_vsolve_fwd!/_naive_vsolve_bwd!, plus inner-product forms for row-contiguous strided parents) instead of being routed into the matrix kernels below n = 128 and deferred to LinearAlgebra/BLAS trsv above it. The VECTOR_RHS_CUTOFF = 128 deferral is deleted: the vector path is now never-BLAS at every size, for all four wrapper types, Float32/Float64, 2- and 3-arg, Val(true/false). Matrix paths are untouched. Version bumped to 0.2.5 (ldiv! vector behavior above n = 128 changes from BLAS-backed to native; results differ only in floating-point rounding).

The swap logic this PR was scoped to add ("naive vs existing kernels, pick per regime") turned out to be unnecessary: after benchmarking, the naive kernel wins every measured cell against both the existing kernels and BLAS, so per the "data decides" instruction the swap was deleted and all vectors route to the naive sweep.

Decision table

Regime	Winner	Rule set
vector RHS, all n in [4, 2000], all 4 wrappers, F32+F64	rank-4 naive sweep, by 1.35–3.43x over the best existing TS path and 1.14–7.42x over trsv (medians of 3 pinned interleaved runs; all 80 cells)	all vector `ldiv!` → naive sweep; no swap, no BLAS
vector RHS, row-contiguous strided parent	(no Base strided type has this layout; not benchmarkable end-to-end)	inner-product form selected by `abs(stride(A,1)) <= abs(stride(A,2))`
matrix RHS, tiny (n ≤ 16, nrhs ∈ 2,4,8)	existing blocked SIMD kernels in 17 of 18 cells (up to 6.6x at n=16, nrhs=8); naive per-column won only (F32, n=4, nrhs=2): 70 vs 100 ns, replicated twice	matrix paths untouched
matrix RHS, large	not contested	matrix paths untouched

Kernel-form findings that set the final shape (AMD EPYC 7502, AVX2, taskset pinned, 1 BLAS thread, interleaved min-times):

@simd ivdep on the inner axpy beats plain @inbounds (measurable at Float32) and beats an explicit @turbo inner loop (e.g. F64 n=256 forward: 6.0 vs 8.9 µs).
Rank-4 outer unroll (retire 4 columns per pass over x) beats rank-2 by 10–35% and rank-1 by ~1.5–2x. It also eliminated a bimodal per-process code-placement slowdown that plagued the rank-1 form (F32 n=64 forward: 830 ns fast-mode vs 1300–1600 ns slow-mode across fresh processes, data-placement ruled out by within-process reallocation tests; rank-4: 526 ns in every process). Rank-8 not attempted (≤ ~10% theoretical headroom, 32-cycle dependency chain).
Unit-diagonal x[j] store discipline is direction-specific and measured: elide forward (storing cost 1.7x at F64 n=16), store backward (eliding cost 1.7x at F32 n=64).
The known "scalar sweep loses 2.3x at n=512" reproduced exactly (blocked single-thread driver: 83.7 µs vs trsv 29.8 µs); the rank-4 kernel at that size is 20.4 µs — 1.45x faster than trsv.
Before the rank-4 form, the large-n choice between the rank-1 sweep and the unblocked matrix driver was a statistical tie (ratio flipping 0.76–1.28 across processes at n ≥ 768), and the blocked driver's apparent win at exactly n = 2000 did not replicate (it collapses 1.3–2.4x at 512–768 and 3000–4000). No robust regime for the old kernels existed even then.

Benchmarks

Median over 3 fresh pinned processes (taskset -c 5, julia -t1 -O3, BLAS.set_num_threads(1)), interleaved min-of-31-passes per candidate, ns per solve. new = this PR's TriangularSolve.ldiv!; ts_cur = the pre-PR sub-cutoff path (matrix drivers via div_dispatch!/div_dispatch_L! with static(1), i.e. what ran at n ≤ 128) extended to all n; ts_ub = same drivers with large-N blocking skipped; trsv = LinearAlgebra.ldiv! (OpenBLAS). Pre-PR public behavior was ts_cur for n ≤ 128 and trsv above.

Float64

wrap	n	new	ts_cur	ts_ub	trsv	best-old/new	trsv/new
L	4	23	33	33	130	1.41	5.57
L	8	48	72	72	191	1.50	3.97
L	16	118	168	166	315	1.41	2.67
L	32	278	377	375	590	1.35	2.12
L	64	603	857	855	1206	1.42	2.00
L	128	1563	2517	2517	3034	1.61	1.94
L	256	4790	9736	8487	8558	1.77	1.79
L	512	20436	86468	33876	29684	1.66	1.45
L	1000	79349	133359	125659	96600	1.58	1.22
L	2000	547158	796626	937485	663587	1.46	1.21
UnitL	4	16	22	22	109	1.41	7.04
UnitL	8	26	43	42	146	1.60	5.56
UnitL	16	58	109	108	235	1.88	4.09
UnitL	32	137	314	311	476	2.26	3.47
UnitL	64	309	771	767	1075	2.48	3.48
UnitL	128	1198	2364	2370	2805	1.97	2.34
UnitL	256	4642	9696	8719	8394	1.88	1.81
UnitL	512	23007	83803	33790	30353	1.47	1.32
UnitL	1000	78870	132819	126470	96219	1.60	1.22
UnitL	2000	564557	807546	942045	695396	1.43	1.23
U	4	23	39	53	117	1.66	5.03
U	8	48	89	104	184	1.84	3.80
U	16	118	206	220	351	1.75	2.99
U	32	290	445	457	664	1.54	2.30
U	64	665	989	1000	1375	1.49	2.07
U	128	1709	2559	2572	3329	1.50	1.95
U	256	4967	9729	8219	9323	1.65	1.88
U	512	21457	86242	29127	30970	1.36	1.44
U	1000	79390	129320	110569	97830	1.39	1.23
U	2000	597907	862685	925915	709846	1.44	1.19
UnitU	4	16	28	36	99	1.73	6.01
UnitU	8	27	54	59	144	1.96	5.26
UnitU	16	60	134	138	272	2.23	4.52
UnitU	32	147	288	293	521	1.96	3.55
UnitU	64	364	679	684	1111	1.87	3.05
UnitU	128	1290	2035	2038	2772	1.58	2.15
UnitU	256	4692	9171	7332	8362	1.56	1.78
UnitU	512	21044	84617	28524	29436	1.36	1.40
UnitU	1000	77760	128379	110100	95699	1.42	1.23
UnitU	2000	567297	873245	874166	662896	1.54	1.17

Float32

wrap	n	new	ts_cur	ts_ub	trsv	best-old/new	trsv/new
L	4	21	30	30	136	1.42	6.40
L	8	42	64	63	188	1.52	4.53
L	16	107	150	145	313	1.36	2.93
L	32	250	355	362	611	1.42	2.45
L	64	536	824	821	1238	1.53	2.31
L	128	1190	2260	2256	2709	1.90	2.28
L	256	3208	6709	6156	6420	1.92	2.00
L	512	11402	23418	22058	18269	1.93	1.60
L	1000	45960	79410	76674	57074	1.67	1.24
L	2000	174469	329198	294548	200889	1.69	1.15
UnitL	4	15	24	22	112	1.49	7.42
UnitL	8	26	45	44	153	1.67	5.83
UnitL	16	60	109	109	243	1.83	4.08
UnitL	32	140	300	300	510	2.14	3.64
UnitL	64	308	905	904	1167	2.94	3.79
UnitL	128	817	2138	2136	2570	2.62	3.15
UnitL	256	2758	6299	5908	6180	2.14	2.24
UnitL	512	11416	23717	21834	18364	1.91	1.61
UnitL	1000	44794	79494	75944	57440	1.70	1.28
UnitL	2000	184539	336038	308068	214119	1.67	1.16
U	4	21	36	48	112	1.67	5.26
U	8	44	82	95	179	1.86	4.05
U	16	106	191	202	334	1.81	3.15
U	32	254	531	542	711	2.09	2.81
U	64	559	1144	1151	1414	2.05	2.53
U	128	1318	2513	2522	3102	1.91	2.35
U	256	3333	6946	6333	7270	1.90	2.18
U	512	11592	23856	21026	20230	1.81	1.75
U	1000	44754	76734	66410	59914	1.48	1.34
U	2000	176399	326489	252269	205419	1.43	1.16
UnitU	4	17	30	37	92	1.83	5.51
UnitU	8	29	59	66	140	2.02	4.80
UnitU	16	64	158	165	255	2.48	3.99
UnitU	32	147	503	510	564	3.43	3.84
UnitU	64	331	1073	1080	1174	3.24	3.55
UnitU	128	886	2465	2473	2610	2.78	2.95
UnitU	256	2733	6865	6070	6354	2.22	2.32
UnitU	512	10761	23356	19200	18531	1.78	1.72
UnitU	1000	44014	76180	67895	57505	1.54	1.31
UnitU	2000	176159	328748	265419	200929	1.51	1.14

new wins all 80 cells against both comparators (min best-old/new = 1.35, min trsv/new = 1.14).

Matrix-path regression spot-check (untouched code, must be unchanged)

Upper-ldiv!, nrhs = 8, Float64, ns per solve, main → branch: n=64: 1500 → 1510; n=256: 20430 → 20010; n=1000: 290169 → 290179. Within noise, as expected — the diff does not touch matrix kernels.

Tiny matrix RHS (naive rank-4 per column vs existing kernels, ns, 2 replicated processes)

T	n	nrhs	existing	naive per-column
F64	4	2	59–60	79
F64	4	8	80	189
F64	8	4	99	219
F64	16	8	229–239	949
F32	4	2	99–100	70
F32	8	4	79	199–200
F32	16	8	139	879–919

Existing kernels win 17 of 18 cells (up to 6.6x); naive wins only (F32, 4, 2) by 30 ns. Matrix swapping rejected; matrix paths untouched.

Behavioral notes for review

Never-BLAS above n = 128 is a results change: bitwise output differs from trsv (different summation order). Accuracy is the same order; tests compare against Matrix(U) \ b references at rtol = sqrt(eps(T)) * n.
The Val(thread) argument is now ignored for vector RHS (previously it only toggled internal driver configuration; threading never engaged at M = 1). Both Val(true) and Val(false) run the same single-threaded sweep; both are tested.
2-arg/3-arg ldiv!(Lower*/UnitLower*, b) with mismatched sizes previously did not throw below the cutoff (div_dispatch! has no size checks; it silently read out of bounds). All four wrappers now throw DimensionMismatch in both arities, matching div_dispatch_L! and LinearAlgebra semantics. Tests added.
muladd sweeps contract to FMA; diagonal division uses plain / (no fastmath approximations).

Verification

All on this machine (AMD EPYC 7502, AVX2, Linux), at commit 30d74ad (the branch tip adds only a README sentence on top of it).

Full suite, Julia 1.12.4 — julia --project=. -e 'using Pkg; Pkg.test()':

Test Summary:      |   Pass   Total     Time
TriangularSolve.jl | 149626  149626  2m44.4s
     Testing TriangularSolve tests passed

Full suite, Julia 1.10.11 — temp env, Pkg.develop(path=...); Pkg.test("TriangularSolve"):

Test Summary:      |   Pass   Total      Time
TriangularSolve.jl | 149626  149626  11m18.2s
     Testing TriangularSolve tests passed

Both runs include the Aqua QA + ambiguity testsets. New tests: vector solves at n ∈ {1, 2, 5, 8, 16, 33, 64, 127, 128, 129, 200, 500, 1201} × 4 wrappers × F32/F64 × Val(true/false) × {2-arg, 3-arg, 3-arg-aliased}, packed lu!-factor vector solves at n ∈ {8..512}, non-contiguous strided views (step-2 parent and vector, 2- and 3-arg), direct tests of the inner-product kernels, vector DimensionMismatch throws for both triangles and arities, BigFloat and Bidiagonal-parent LinearAlgebra fallbacks, and warm-in-loop zero-allocation assertions for 2- and 3-arg vector solves. Beyond the suite: an exhaustive n = 1:40 sweep (every rank-4 remainder case) × 4 wrappers × F32/F64 × both arities passed, plus n ∈ {127, 128, 129, 500, 1000, 2003}.

typos over the diff: clean (one false positive on a git hash). JuliaFormatter v1 with the repo's .JuliaFormatter.toml was run; its churn on untouched code was reverted so the diff stays confined to the vector path (formatting-only sweep can be a separate PR). Added-line comment share: 17/254 = 6.7%.

Not verified: AVX-512 and aarch64/Apple Silicon (all crossover measurements are AVX2-EPYC; the all-naive rule has no tunable constants to mis-transfer, but win margins there are unmeasured), MKL as BLAS backend, Windows/macOS, multi-threaded BLAS comparisons (deliberately pinned to 1 thread), and downstream packages (RecursiveFactorization.jl, LinearSolve.jl) against this branch.

Downstream

TriangularSolve-only routing for NotIPIV backsolves; fix pivot-free lu! with user ipiv; factorization dispatch audit JuliaLinearAlgebra/RecursiveFactorization.jl#114 reshapes vectors to n×1 to stay on TriangularSolve kernels and pays the old 2.4x penalty at n = 512; once this registers it can pass vectors straight through and drop the reshape.
RFLU: pure-Julia naive back-solve for single right-hand sides SciML/LinearSolve.jl#1169 adds an LS-side pure-Julia back-solve for single RHS to work around the same gap; with this registered, RFLU's vector path can rely on TriangularSolve directly. (Context: GenericLUFactorization: pure-Julia naive back-solve instead of getrs! SciML/LinearSolve.jl#1151, Orientation-specialized naive LU back-solve for Adjoint/Transpose factors SciML/LinearSolve.jl#1162, Enforce the LU backsolve-selection policy: GenericLU stays generic, LUFactorization gets the size-aware back-solve SciML/LinearSolve.jl#1164.)
Possible follow-ups here (not this PR): fusing the rank-4 sweep across a few RHS columns might win the tiny-(n, nrhs) matrix corner; the same rank-4 restructuring applied to the matrix kernels' M = 1 remainder would be redundant once vectors bypass them entirely.

🤖 Generated with Claude Code

Vector ldiv! previously ran the matrix kernels' M == 1 scalar remainder up to n = 128 and deferred to LinearAlgebra (BLAS trsv) above. Both lose to a plain column-oriented substitution sweep that LLVM vectorizes: outer-unrolled rank-4, it beats the old kernel path by 1.35-3.4x and trsv by 1.14-7.4x at every measured size (n = 4..2000, both triangles, unit/non-unit, Float32/64, AVX2), so the vector path is now never-BLAS with no swap logic. Inner-product forms cover row-contiguous strided parents. Vector methods now also throw DimensionMismatch for all four wrappers (the Lower ones previously read out of bounds on mismatched sizes). Matrix paths are unchanged; a naive per-column sweep lost to the blocked SIMD kernels in 17 of 18 tiny-size cells, so no matrix swapping. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

codecov · 2026-08-08T20:21:51Z

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 95.22%. Comparing base (24abdb7) to head (4e4f8f6).
⚠️ Report is 7 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #48      +/-   ##
==========================================
+ Coverage   94.69%   95.22%   +0.52%     
==========================================
  Files           1        1              
  Lines         849      942      +93     
==========================================
+ Hits          804      897      +93     
  Misses         45       45

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

…r kernels TriangularSolve 0.2.5 (JuliaSIMD/TriangularSolve.jl#48) replaced its vector-entry BLAS deferral (native <= 128, trsv above) with a rank-4 pure-Julia sweep: never-BLAS and faster than trsv/getrs! at every size (measured 0.33-0.79x of getrs! and ~2x faster than the n-by-1-reshape workaround this branch previously used, 1 thread, Zen 2; allocation-free). Drop the reshape helper, call the vector entry directly, and raise the TriangularSolve compat floor to 0.2.5 so the vector legs can never silently defer to BLAS on older TriangularSolve. The dispatch audit now asserts the vector signatures resolve to the native vector kernel methods. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01UGVN6qeNL2jGCtaYg1386X

…nels TriangularSolve 0.2.5 (JuliaSIMD/TriangularSolve.jl#48) made the vector ldiv! entry never-BLAS at every size and faster than getrs!/trsv (measured 0.33-0.79x of getrs!, 1 thread, Zen 2; allocation-free), so _rf_ldiv!'s vector method calls it directly instead of presenting the vector as an n-by-1 reshape onto the matrix kernels (~2x slower). Raise the TriangularSolve compat floor to 0.2.5 — below it the vector entry silently deferred large vectors to BLAS trsv, which the routing policy forbids. Supersedes the ext-local naive back-solve approach of SciML#1169. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01UGVN6qeNL2jGCtaYg1386X

…rash, factorization-cell audit (#1171) * RFLU backsolves: TriangularSolve-only routing, fix pivot=Val(false) crash The extension's backsolve consumed fact.ipiv even for pivot = Val(false), where RecursiveFactorization returns the caller-supplied ipiv unwritten (undefined memory): every RFLUFactorization(pivot = Val(false)) and RF32MixedLUFactorization(pivot = Val(false)) solve segfaulted in dlaswp (vector RHS, via LAPACK.getrs!) or threw BoundsError (matrix RHS, via _ipiv_rows!). _rf_ldiv! now takes the pivot flag from the algorithm type and never touches fact.ipiv when pivoting is off. Routing policy: wherever TriangularSolve has a native kernel (Float32/Float64, strided), both triangular legs run on TriangularSolve — never on a BLAS kernel. Vector right-hand sides and single-column matrices, which previously took LAPACK.getrs!, are presented to TriangularSolve as n-by-1 view-reshapes (allocation-free, preserving the zero-allocation re-solve contract); the n x 1 early-return in the matrix path is gone. Complex and non-strided types keep the stdlib path (TriangularSolve has no kernels there), with NoPivot solved through the triangular legs so no path reads unwritten pivots. RF32MixedLUFactorization shares _rf_ldiv! instead of plain ldiv!. Adds _ts_native_backsolve, a which()-based enforcement helper the test suite uses to fail if TriangularSolve dispatch ever lands these argument types on its LinearAlgebra catch-all again. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01UGVN6qeNL2jGCtaYg1386X * Test RFLU backsolve routing: both pivots, all RHS shapes, dispatch audit Correctness for {vector, matrix, n-by-1} x {pivoted, NoPivot} x n in {8, 40, 300} (300 spans TriangularSolve's vector-entry cutoff) plus ComplexF64, and a which()-based dispatch audit asserting the extension's TriangularSolve-routing methods are selected and both triangular legs — including the n-by-1 view-reshape type the vector path actually passes — resolve to native TriangularSolve kernels, never its LinearAlgebra catch-all. The NoPivot cases crashed (segfault/BoundsError) before the extension fix. Re-enables the previously commented-out RF32MixedLUFactorization(pivot = Val(false)) test. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01UGVN6qeNL2jGCtaYg1386X * Pin the RFLU factorization cells: ext-gated defaults, complex stays opt-in RFLU's factorization is RecursiveFactorization.lu!, dispatch-audited as BLAS-free for Float32/Float64 in RF's own suite; the complex panel solves are the one LAPACK cell (TriangularSolve has no complex kernels). Assert the two LinearSolve-side facts that keep that cell contained: the ext enables userecursivefactorization, and defaultalg never routes complex matrices to RFLU — so LAPACK panels stay reachable only by explicitly requesting the algorithm. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01UGVN6qeNL2jGCtaYg1386X * RFLU vector backsolves: use TriangularSolve 0.2.5's native vector kernels TriangularSolve 0.2.5 (JuliaSIMD/TriangularSolve.jl#48) made the vector ldiv! entry never-BLAS at every size and faster than getrs!/trsv (measured 0.33-0.79x of getrs!, 1 thread, Zen 2; allocation-free), so _rf_ldiv!'s vector method calls it directly instead of presenting the vector as an n-by-1 reshape onto the matrix kernels (~2x slower). Raise the TriangularSolve compat floor to 0.2.5 — below it the vector entry silently deferred large vectors to BLAS trsv, which the routing policy forbids. Supersedes the ext-local naive back-solve approach of #1169. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01UGVN6qeNL2jGCtaYg1386X --------- Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com> Co-authored-by: Claude Fable 5 <noreply@anthropic.com>

ChrisRackauckas added 2 commits August 8, 2026 16:00

Update README for the native vector right-hand-side path

4e4f8f6

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

ChrisRackauckas marked this pull request as ready for review August 8, 2026 20:30

ChrisRackauckas merged commit ee8dd65 into JuliaSIMD:main Aug 8, 2026
7 checks passed

ChrisRackauckas-Claude mentioned this pull request Aug 8, 2026

TriangularSolve-only routing for NotIPIV backsolves; fix pivot-free lu! with user ipiv; factorization dispatch audit JuliaLinearAlgebra/RecursiveFactorization.jl#114

Closed

Navigation Menu

Uh oh!

Replace the vector-RHS BLAS deferral with a rank-4 pure-Julia sweep (never-BLAS vectors) - #48

Replace the vector-RHS BLAS deferral with a rank-4 pure-Julia sweep (never-BLAS vectors)#48
ChrisRackauckas merged 2 commits into
JuliaSIMD:mainfrom
ChrisRackauckas-Claude:naive-vector-kernel

ChrisRackauckas-Claude commented Aug 8, 2026

Uh oh!

codecov Bot commented Aug 8, 2026 •

edited

Loading

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Conversation

ChrisRackauckas-Claude commented Aug 8, 2026

What changed and why

Decision table

Benchmarks

Float64

Float32

Matrix-path regression spot-check (untouched code, must be unchanged)

Tiny matrix RHS (naive rank-4 per column vs existing kernels, ns, 2 replicated processes)

Behavioral notes for review

Verification

Downstream

Uh oh!

codecov Bot commented Aug 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

codecov Bot commented Aug 8, 2026 •

edited

Loading