TriangularSolve-only routing for NotIPIV backsolves; fix pivot-free lu! with user ipiv; factorization dispatch audit

[ChrisRackauckas-Claude]

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

Policy

Per the routing directive for RecursiveFactorization: both the backsolves consuming an RF factorization and the lu! factorization itself should run on TriangularSolve's native kernels / RF's own @turbo kernels — vector and matrix right-hand sides, with and without pivoting — and never silently fall back to a BLAS/LAPACK kernel where a native path exists (Float32/Float64, strided). This PR enforces that for everything RecursiveFactorization itself controls, audits the factorization path end to end, and fixes an undefined-memory crash found during the audit.

What changed

1. lu!(A, ipiv, Val(false), ...) now fills a user-supplied ipiv with the identity (src/lu.jl). Since ca26d78 (2023-01), the pivot-free algorithm on Julia ≥ 1.8 returned the caller's ipiv inside the LU without ever writing it. Any consumer of F.ipiv — LAPACK.getrs! via the stdlib ldiv!, or LinearAlgebra._ipiv_rows! — then read undefined memory. This is how LinearSolve calls lu!, and it segfaulted in dlaswp (companion LinearSolve PR below). NotIPIV (the path RF's own lu!(A, Val(false)) takes) is unaffected.
2. Vector backsolves stay on TriangularSolve's native vector kernels at every size, and the compat floor moves to TriangularSolve 0.2.5. TS ≤ 0.2.3 sent vectors to its LinearAlgebra/BLAS catch-all; 0.2.4 was native only to n=128 and deferred to BLAS trsv above. TriangularSolve 0.2.5 (Replace the vector-RHS BLAS deferral with a rank-4 pure-Julia sweep (never-BLAS vectors) JuliaSIMD/TriangularSolve.jl#48, merged and registered while this PR was in flight) replaced that deferral with a rank-4 pure-Julia sweep — never-BLAS and faster than trsv at all sizes — so the NotIPIV ldiv! now simply calls the vector entry directly (an interim revision of this branch routed vectors as n×1 reshapes onto the matrix kernels; that workaround is gone — measured ~2× slower than 0.2.5's vector kernels). Also removes a dead square_view line in that method and returns B explicitly.
3. Factorization-path audit + enforcement: removed the vestigial BLAS import from src/lu.jl (never used — the module now imports no BLAS entry point at all), and added a dispatch-audit testset for the factorization's panel solves (details below). No behavior change in the factorization; items 1–2 are the only behavior changes in the PR.
4. Butterfly testset made seeded, with a backward-stable residual bound (test/runtests.jl). The first CI round failed only the julia-pre (1.13.0-rc1) ubuntu lane, in the pre-existing 🦋 testset: norm(A*out .- b) = 7.6e-10 vs an absolute 1e-10 bound (https://github.com/JuliaLinearAlgebra/RecursiveFactorization.jl/actions/runs/31273560107/job/93143409318). Investigation: the same suite passes on 1.13.0-rc1 locally; a 550-draw sweep (n = 790–810, 50 seeds each) puts the worst residual at 2.2e-11 for the new backsolve route vs 1.1e-11 for the old one — same order, both far under the bound — while the backward-stable band for these n≈800 solves is ~1e-9 absolute, i.e. the old 1e-10 sat inside legitimate rounding. The test's inputs are also hardware-dependent (VectorizedRNG streams vary with SIMD width) and its b draws depended on suite-order RNG position, which the new testsets shifted. Fix: seed the testset locally and bound the relative residual at 1e-8 (≥20× above any observed rounding value, ≥3 orders below genuine breakage), with the justification in a comment.

Backsolve truth table for ldiv! on RF factorizations (direct use)

case	before (RF 0.2.28 + TS 0.2.4)	after
pivoted, vector	LinearAlgebra getrs! (BLAS)	unchanged — see "Not in scope"
pivoted, matrix	LinearAlgebra getrs! (BLAS)	unchanged — see "Not in scope"
NoPivot (NotIPIV), vector	TS vector kernel n ≤ 128, BLAS trsv above	TS native vector kernels, all n (TS ≥ 0.2.5)
NoPivot (NotIPIV), matrix	TS native matrix kernels	unchanged

Evidence (live which chains on Julia 1.12.4, RecursiveFactorization 0.2.28→this branch, TriangularSolve 0.2.5): the NoPivot entries resolve to RecursiveFactorization src/lu.jl and then TriangularSolve's native kernels — matrix legs at TriangularSolve.jl:691/818, vector legs at the new vector methods (TriangularSolve.jl:1083, both 2- and 3-arg) — never the LinearAlgebra.ldiv! catch-all. An IR walk of the 0.2.5 vector kernel reaches no BLAS/LAPACK symbol, and a solve is allocation-free.

Not in scope: ldiv! on a pivoted RF lu! result dispatches to LinearAlgebra's own method (LU{T, <:StridedMatrix, Vector{BlasInt}} → getrs!); RF cannot own that method without type piracy (every type in the signature belongs to Base/LinearAlgebra). LinearSolve routes the pivoted case through TriangularSolve on its side (companion PR below).

Factorization truth table for lu! (audited, dispatch-enforced; unchanged behavior)

eltype	size (mnmin)	route	kernels	BLAS/LAPACK?
Float64/Float32, strided	> threshold (40; 48 AVX512)	recursive	@turbo _generic_lufact! leaves + @(t)turbo schur_complement! + native TS panel ldiv! (TriangularSolve.jl:691 threaded / :705 unthreaded for the exact PtrArray view types) + pure-Julia apply_permutation!	none
Float64/Float32, strided	≤ threshold	RF _generic_lufact!	@turbo right-looking LU	none
any, pivoted, 3-arg lu!/lu with minmn < 10	—	stdlib LinearAlgebra.generic_lufact!	pure Julia	none (IR-walked). Not reachable via LinearSolve (it calls the 4-arg lu!)
ComplexF64/ComplexF32, strided	> threshold	recursive	leaves + Schur are LoopVectorization fallback loops (LV rejects complex — still no BLAS); panel ldiv! → TS catch-all (:987) → LinearAlgebra.ldiv! → LAPACK.trtrs!	yes — the panel solves only
ComplexF64/ComplexF32	≤ threshold	_generic_lufact! (fallback loop)	pure Julia	none
non-BlasFloat (BigFloat, Matrix{Real}, …) or non-strided	any	recursive-generic or _generic_lufact!	generic Julia (panels land on TS catch-all → LinearAlgebra generic triangular solve)	none

Evidence, all on Julia 1.12.4 (audit scripts: fact_probe.jl):

• Static: grep -rn "getrf\|LAPACK\|BLAS\." src/ matches nothing but comments; the BLAS import was vestigial (removed here). which resolves the exact recursion panel types (SquarePtrMatrix/PtrArray views) to native TS kernel methods for Float32/Float64, and to the catch-all for complex.
• Dynamic: running lu! over Float64/Float32 × {pivot, NoPivot} × {threaded, not} × {n=30, 300} plus a fat 200×300 and the tiny minmn<10 branch creates zero TriangularSolve catch-all specializations with real eltype — none exist in the whole session (Base.specializations sweep). The complex runs create exactly the complex panel-signature specializations.
• IR walk: walking the :invoke edges of the optimized IR (depth 6) for TS.ldiv! on the panel types, schur_complement!, RF _generic_lufact!, and stdlib generic_lufact! finds no gemm/gemv/getrf/trsm/trsv/trtrs/syrk/ger anywhere. The same walk on LinearAlgebra.ldiv! at the complex panel types hits LAPACK.trtrs! at depth 1 — that is the one BLAS/LAPACK dependence left in the factorization, and it is a TriangularSolve capability gap (no complex kernels), not an RF routing choice. Left unchanged and characterized in the tests; see proposal item 4.

Crash fix: failing before / passing after

Reproduction (rf_ipiv_repro.jl): poison an undef ipiv, factor with Val(false), consume F.ipiv via the stdlib ldiv!:

Before (registered RecursiveFactorization v0.2.28):

F.ipiv == 1:n  ->  false
Test Failed: Expression: F.ipiv == 1:n
   Evaluated: [9223372036854775800, …] == 1:30

(and from LinearSolve, where the ipiv is genuinely undef rather than poisoned, solve(prob, RFLUFactorization(pivot=Val(false))) dies with signal 11 (Segmentation fault) … dlaswp_plus … dgetrs_N_single.)

After (this branch):

F.ipiv == 1:n  ->  true
residual: 5.113875363395961e-16
PASS

The new "NoPivot lu! with a user-supplied ipiv leaves valid pivots" testset encodes exactly this.

Never-fallback enforcement

test/runtests.jl now contains three dispatch audits:

• "NotIPIV backsolves stay on TriangularSolve's native kernels": resolves, via which, every signature the NotIPIV ldiv! hands to TriangularSolve.ldiv! and asserts the resolved method is a TriangularSolve native kernel and not identical to the catch-all method object. A correctness testset also pins vector/matrix results at n = 8…300, i.e. across TriangularSolve's n=128 vector cutoff.
• "Factorization panel solves stay on TriangularSolve" (runs last): (a) the exact PtrArray panel-view types the Float32/Float64 recursion constructs must resolve to native TS kernels, threaded and not; (b) a whole-suite specialization sweep — after every factorization and backsolve in the test run, no TriangularSolve catch-all specialization may exist whose triangular argument has Float32/Float64 eltype, so any real-eltype solve anywhere that silently fell back to LinearAlgebra/BLAS fails the suite; (c) a characterization of the complex gap: complex panel signatures resolve to the catch-all (flips, and should be updated, if TriangularSolve ever gains complex kernels).

If a future TriangularSolve or RF restructuring reintroduces a silent BLAS fallback on real eltypes, these fail.

Performance

The factorization commit is behavior-neutral (import removal + tests only). For the backsolve: with TriangularSolve 0.2.5 the vector legs are faster than BLAS at every size — the earlier judgment call (accepting 2.3–2.6× at n = 512–1000 for never-BLAS purity, measured against the reshape workaround) no longer exists. Single-vector backsolve, interleaved min-times, 25×10 reps, 1 thread, Julia 1.12.4, Zen 2 (bench_ts5.jl; getrs! = pivoted stdlib baseline, reshape = this branch's interim route):

n      getrs!       TSvec(0.2.5)  reshape      vec/getrs  vec/reshape
8     291.0 ns     95.0 ns       189.0 ns     0.33       0.5
16    597.0 ns     211.9 ns      365.0 ns     0.35       0.58
32    1.11 us      456.9 ns      805.0 ns     0.41       0.57
64    2.4 us       1.0 us        1.8 us       0.42       0.56
100   4.42 us      1.89 us       3.3 us       0.43       0.57
128   6.03 us      2.65 us       4.73 us      0.44       0.56
200   12.16 us     6.51 us       11.57 us     0.54       0.56
256   17.64 us     9.61 us       19.72 us     0.54       0.49
512   60.14 us     42.32 us      180.98 us    0.7        0.23
1000  204.91 us    161.57 us     268.57 us    0.79       0.6

0 bytes allocated per solve (n=300 spot check). Never-BLAS and fastest option at every measured size.

Butterfly LU: BLAS census (findings only — nothing changed here)

src/butterflylu.jl still uses BLAS per solve, verified by the same IR walk:

• mul!(tmp, U', b) (line 50) and mul!(b, V, tmp) (line 52) in 🦋solve! — dense BLAS gemv applications of the U/V butterfly transforms, once each per solve. (Mirrored in LinearSolve's ButterflyFactorization solve!.)
• materializeUV's mul!(U, Bu2, Bu1) / mul!(V, Bv2, Bv1) (lines 176–177) — SparseBandedMatrix operands are not strided, so these go through stdlib generic matmul: not BLAS, and setup-only (once per workspace).
• The triangular legs ldiv!(F, tmp, thread) route TS-catch-all → LinearAlgebra.ldiv! → RF's NotIPIV method → TS native kernels (two dispatch hops, but native; covered by the NotIPIV audit).

Replacing the per-solve gemvs would need either a TriangularSolve/LoopVectorization dense-gemv kernel or restructuring the transform application; that is a design decision, not a routing fix, so it is left for review.

Verification

All with registered TriangularSolve 0.2.5 (fresh resolves), full suite julia --project=. -e 'using Pkg; Pkg.test()':

• Julia 1.12.4: exit 0 — Test LU factorization | 3120 pass, NoPivot lu! with a user-supplied ipiv | 6 pass, NotIPIV backsolves stay on TriangularSolve | 38 pass, NotIPIV ldiv! correctness | 24 pass, 🦋 | 21 pass, Factorization panel solves stay on TriangularSolve | 31 pass (3240 total).
• Julia 1.13.0-rc1: exit 0 (same testsets).
• Julia 1.10.11: exit 0 (same testsets).
• JULIA_NUM_THREADS=3, Julia 1.12.4 (threaded @batch/@tturbo paths): exit 0.
• Earlier revisions of this branch were additionally validated on 1.10/1.12/1.13 with TS 0.2.4 (reshape route) before the 0.2.5 rework.
• typos on the diff: clean (the siz hits in runtests.jl are pre-existing lines outside this diff).
• Not verified: Windows/macOS, x86, AVX512 hardware (threshold-48 branch taken structurally, not on this Zen 2 machine), downstream consumers other than LinearSolve (LinearSolve validated against this branch in the companion PR).
• CI note: the evaluate (Invalidations) job is red on this PR and on master since 2025-10-08 — an upstream SnoopCompile crash on Julia 1.12 (Core.Binding in report_invalidations JuliaDebug/SnoopCompile.jl#465), RF-independent. Separate mechanical fix: Pin the Invalidations workflow to Julia 1.11 (SnoopCompile #465 crashes on 1.12) #115.

Judgment calls a reviewer should check

1. The TriangularSolve compat floor moves 0.2.2 → 0.2.5 (the vector entry below 0.2.5 silently deferred large vectors to BLAS, which this PR's policy forbids). Consumers pinned to older TriangularSolve will hold back RecursiveFactorization 0.2.29.
2. Complex factorization cell left on LAPACK trtrs! — TriangularSolve has no complex kernels, so "TriangularSolve's path" does not exist there; rerouting complex to _generic_lufact! would be never-BLAS but a large regression at big n. Characterized in tests instead of changed; proposal item 4 is the structural fix.
3. The whole-suite specialization sweep in the new testset asserts a global invariant (no real-eltype TS catch-all specialization after the full run) — deliberate, so any new code path that silently falls back fails the suite, but it means unrelated future tests that legitimately hand TS a non-strided Float64 triangular would need the testset updated.
4. Version bump 0.2.28 → 0.2.29 (bugfix + additive behavior change on the NotIPIV vector route).
5. The butterfly tolerance change (absolute 1e-10 → relative 1e-8 + local seed) — loosening a bound is normally suspect; the distribution evidence and the backward-stability arithmetic above are the case that 1e-10 absolute was wrong-in-principle for n≈800. Push back if a tighter, still-principled bound is preferred (e.g. c·n·eps·‖b‖ with measured c).

TriangularSolve follow-up proposal (not implemented here, deliberately not a PR)

1. Strict entry point (ldiv_native! or similar) with no catch-all, so never-fallback callers get MethodError instead of silent BLAS.
2. Public has_native_kernel(::Type, ::Type) trait maintained next to the kernels, replacing downstream which-based audits.
3. Blocked single-RHS kernel to close the n ≥ 512 gap — delivered upstream while this PR was in flight: Replace the vector-RHS BLAS deferral with a rank-4 pure-Julia sweep (never-BLAS vectors) JuliaSIMD/TriangularSolve.jl#48 (registered as 0.2.5) removed VECTOR_RHS_CUTOFF and the BLAS deferral outright; this PR now consumes it.
4. Complex kernels (split-complex over the existing real kernels, or hand-SIMD via VectorizationBase — LoopVectorization will not do complex), which would close the last BLAS/LAPACK cell in RF's factorization table.

Companion LinearSolve PR: SciML/LinearSolve.jl#1171

🤖 Generated with Claude Code