From cc126db8c902fdf6546afd16de688299d92de718 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 8 Aug 2026 09:01:41 -0400 Subject: [PATCH 1/5] Initialize a user-supplied ipiv in pivot-free lu! Since ca26d78 (Jan 2023), lu!(A, ipiv, Val(false), ...) on Julia >= 1.8 returned the caller's ipiv inside the LU without ever writing it. Stdlib consumers of F.ipiv (LAPACK.getrs! via ldiv!, LinearAlgebra._ipiv_rows!) then read undefined memory: LinearSolve.jl's RFLUFactorization(pivot = Val(false)) segfaulted in dlaswp on every vector solve and threw BoundsError on matrix solves. Fill the supplied vector with the identity permutation; NotIPIV (RF's own pivot-free path) is unaffected. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01UGVN6qeNL2jGCtaYg1386X --- src/lu.jl | 6 +++++- test/runtests.jl | 16 ++++++++++++++++ 2 files changed, 21 insertions(+), 1 deletion(-) diff --git a/src/lu.jl b/src/lu.jl index 78854ad..74019f8 100644 --- a/src/lu.jl +++ b/src/lu.jl @@ -100,7 +100,11 @@ function lu!(A::AbstractMatrix{T}, ipiv::AbstractVector{<:Integer}, info = zero(BlasInt) m, n = size(A) mnmin = min(m, n) - if pivot === Val(false) && !CUSTOMIZABLE_PIVOT + # The pivot-free algorithm never writes `ipiv`, but a user-supplied vector + # is still returned inside the `LU`; fill it with the identity so consumers + # of `F.ipiv` (e.g. `LAPACK.getrs!`, `LinearAlgebra._ipiv_rows!`) see valid + # pivots instead of undefined memory. + if pivot === Val(false) && !(ipiv isa NotIPIV) copyto!(ipiv, 1:mnmin) end if recurse(A) && mnmin > threshold diff --git a/test/runtests.jl b/test/runtests.jl index 2454246..8f9e597 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -65,6 +65,22 @@ testlu(A::Union{Transpose, Adjoint}, MF, BF, p) = testlu(parent(A), parent(MF), end end +@testset "NoPivot lu! with a user-supplied ipiv leaves valid pivots" begin + # Stdlib consumers (`LAPACK.getrs!`, `_ipiv_rows!`) read `F.ipiv`; it must + # hold the identity, not undefined memory (which crashed in `dlaswp`). + for T in (Float64, Float32, ComplexF64) + n = 30 + A = rand(T, n, n) + T(10) * LinearAlgebra.I + b = rand(T, n) + ipiv = Vector{LinearAlgebra.BlasInt}(undef, n) + fill!(ipiv, typemax(LinearAlgebra.BlasInt) - 7) # poison undefined memory + F = RecursiveFactorization.lu!(copy(A), ipiv, Val(false), Val(false)) + @test F.ipiv == 1:n + x = LinearAlgebra.ldiv!(F, copy(b)) # stdlib LU path consumes F.ipiv + @test norm(A * x - b) < 1000 * n * eps(real(T)) + end +end + function wilkinson(N) A = zeros(N, N) A[1:(N+1):N*N] .= 1 From 19fd0147e3db3e560e65a0fde050b921bb189694 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 8 Aug 2026 09:01:52 -0400 Subject: [PATCH 2/5] Keep NotIPIV vector backsolves on TriangularSolve's native kernels ldiv!(::LU{T, <:StridedMatrix, <:NotIPIV}, b::StridedVector) handed the vector straight to TriangularSolve.ldiv!, whose vector entry point is native only up to n = 128 and defers to BLAS trsv above (and on TriangularSolve <= 0.2.3 always fell through the LinearAlgebra catch-all). Present a contiguous vector as an n-by-1 matrix (zero-copy reshape) so both triangular legs stay on TriangularSolve's native matrix kernels at every size; return B in its original shape. Also drop a dead square_view binding in that method. Add a which()-based dispatch audit that fails if any signature RF hands to TriangularSolve.ldiv! resolves to the LinearAlgebra catch-all again, plus correctness tests across the n = 128 cutoff, and bump to 0.2.29. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01UGVN6qeNL2jGCtaYg1386X --- Project.toml | 5 +++-- src/lu.jl | 16 ++++++++++++++-- test/runtests.jl | 50 +++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 66 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 270a9f9..02e34c9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "RecursiveFactorization" uuid = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" authors = ["Yingbo Ma "] -version = "0.2.28" +version = "0.2.29" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -27,6 +27,7 @@ julia = "1.10" [extras] Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +TriangularSolve = "d5829a12-d9aa-46ab-831f-fb7c9ab06edf" [targets] -test = ["Test", "Random"] +test = ["Test", "Random", "TriangularSolve"] diff --git a/src/lu.jl b/src/lu.jl index 74019f8..b7165ac 100644 --- a/src/lu.jl +++ b/src/lu.jl @@ -53,10 +53,22 @@ if CUSTOMIZABLE_PIVOT && isdefined(LinearAlgebra, :_ipiv_rows!) end end if CUSTOMIZABLE_PIVOT + # TriangularSolve's native kernels take strided *matrix* right-hand sides; + # its vector entry point defers to BLAS `trsv` above a size cutoff. A + # contiguous vector presented as an n×1 matrix stays on TriangularSolve's + # kernels at every size. view-then-reshape keeps the wrapper + # allocation-free (immutable, passed by value), unlike reshape(::Vector), + # which heap-allocates a Matrix header on every solve. + _ts_backsolve_rhs(B) = B + function _ts_backsolve_rhs(b::StridedVector{T}) where {T <: Union{Float32, Float64}} + R = reshape(view(b, :), length(b), 1) + return R isa StridedMatrix{T} ? R : b + end function LinearAlgebra.ldiv!(A::LU{T, <:StridedMatrix, <:NotIPIV}, B::StridedVecOrMat{T}) where {T <: BlasFloat} - tri = @inbounds square_view(A.factors, size(A.factors, 1)) - ldiv!(UpperTriangular(A.factors), ldiv!(UnitLowerTriangular(A.factors), B)) + B′ = _ts_backsolve_rhs(B) + ldiv!(UpperTriangular(A.factors), ldiv!(UnitLowerTriangular(A.factors), B′)) + return B end end diff --git a/test/runtests.jl b/test/runtests.jl index 8f9e597..565d224 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,9 @@ using Test import RecursiveFactorization import LinearAlgebra -using LinearAlgebra: norm, Adjoint, Transpose, ldiv! +import TriangularSolve +using LinearAlgebra: norm, Adjoint, Transpose, ldiv!, UnitLowerTriangular, + UpperTriangular using Random Random.seed!(12) @@ -81,6 +83,52 @@ end end end +@testset "NotIPIV backsolves stay on TriangularSolve's native kernels" begin + # The signatures `ldiv!(::LU{..., <:NotIPIV}, ...)` hands to + # TriangularSolve.ldiv! must keep resolving to native kernel methods, never + # to a catch-all that forwards to LinearAlgebra (= BLAS for BLAS types). + catchall2 = which(TriangularSolve.ldiv!, Tuple{Any, Any}) + catchall3 = which(TriangularSolve.ldiv!, Tuple{Any, Any, Val{true}}) + for T in (Float64, Float32) + MT = Matrix{T} + F = RecursiveFactorization.lu(rand(T, 40, 40) + T(10) * LinearAlgebra.I, + Val(false)) + @test F.ipiv isa RecursiveFactorization.NotIPIV + for BT in (Vector{T}, MT) + m = which(LinearAlgebra.ldiv!, Tuple{typeof(F), BT}) + @test m.module === RecursiveFactorization + end + # a contiguous vector reshapes, allocation-free, to a strided n×1 matrix + rhs = RecursiveFactorization._ts_backsolve_rhs(zeros(T, 4)) + @test rhs isa StridedMatrix{T} + RT = typeof(rhs) + for W in (UnitLowerTriangular{T, MT}, UpperTriangular{T, MT}), + BT in (MT, RT) + + m2 = which(TriangularSolve.ldiv!, Tuple{W, BT}) + m3 = which(TriangularSolve.ldiv!, Tuple{W, BT, Val{true}}) + @test m2 !== catchall2 + @test m3 !== catchall3 + @test m2.module === TriangularSolve + @test m3.module === TriangularSolve + end + end +end + +@testset "NotIPIV ldiv! correctness across the TriangularSolve size cutoff" begin + for T in (Float64, Float32), n in (8, 64, 200, 300) + A = rand(T, n, n) + T(10) * LinearAlgebra.I + b = rand(T, n) + B = rand(T, n, 3) + F = RecursiveFactorization.lu(A, Val(false)) + x = ldiv!(F, copy(b)) + @test x isa Vector{T} + @test norm(A * x - b) < 1000 * n * eps(T) + X = ldiv!(F, copy(B)) + @test norm(A * X - B) < 1000 * n * eps(T) + end +end + function wilkinson(N) A = zeros(N, N) A[1:(N+1):N*N] .= 1 From 074430bb385b68ed419f8159517ca1fd6115ffe3 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 8 Aug 2026 14:56:25 -0400 Subject: [PATCH 3/5] Dispatch-audit the factorization path; drop the vestigial BLAS import MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The recursive lu! never calls BLAS/LAPACK for Float32/Float64: leaf factorizations run RF's own @turbo _generic_lufact!, Schur complements run @(t)turbo schur_complement!, and panel solves dispatch to TriangularSolve's native kernels (verified by walking the optimized IR of all three for the exact PtrArray panel types — no gemm/getrf/trsm/trsv/trtrs/syrk/ger anywhere). The BLAS binding imported in lu.jl was never used; remove it, so the module no longer imports any BLAS entry point. New testset enforces the factorization side of the routing table: (a) the exact panel-view types the Float32/Float64 recursion constructs must resolve to native TriangularSolve kernel methods; (b) a whole-suite sweep asserts no TriangularSolve catch-all specialization exists with a Float32/Float64 triangular argument — any silent LinearAlgebra/BLAS fallback of a real-eltype solve anywhere in the test run fails it; (c) characterizes the known gap: complex panel solves resolve to the catch-all (LinearAlgebra -> LAPACK trtrs!), since TriangularSolve has no complex kernels. Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01UGVN6qeNL2jGCtaYg1386X --- src/lu.jl | 2 +- test/runtests.jl | 46 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 47 insertions(+), 1 deletion(-) diff --git a/src/lu.jl b/src/lu.jl index b7165ac..66511e2 100644 --- a/src/lu.jl +++ b/src/lu.jl @@ -1,6 +1,6 @@ using LoopVectorization using TriangularSolve: ldiv! -using LinearAlgebra: BlasInt, BlasFloat, LU, UnitLowerTriangular, checknonsingular, BLAS, +using LinearAlgebra: BlasInt, BlasFloat, LU, UnitLowerTriangular, checknonsingular, LinearAlgebra, Adjoint, Transpose, UpperTriangular, AbstractVecOrMat using StrideArraysCore using StrideArraysCore: square_view diff --git a/test/runtests.jl b/test/runtests.jl index 565d224..2b4fc27 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -151,3 +151,49 @@ end end end + +@testset "Factorization panel solves stay on TriangularSolve" begin + catchall2 = which(TriangularSolve.ldiv!, Tuple{Any, Any}) + catchall3 = which(TriangularSolve.ldiv!, Tuple{Any, Any, Val{true}}) + # (a) The exact panel types the recursive lu! constructs for Float32/Float64 + # (PtrArray-wrapped views) must resolve to native TriangularSolve kernel + # methods, threaded and not. + for T in (Float64, Float32) + A = rand(T, 100, 100) + B = view(RecursiveFactorization.PtrArray(A), axes(A)...) + n1 = RecursiveFactorization.nsplit(T, 100) + A11 = RecursiveFactorization.square_view(B, n1) + A12 = @view B[1:n1, (n1 + 1):100] + for V in (Val{false}, Val{true}) + m = which(TriangularSolve.ldiv!, + Tuple{UnitLowerTriangular{T, typeof(A11)}, typeof(A12), V}) + @test m !== catchall3 + @test m.module === TriangularSolve + end + end + # (b) Whole-suite sweep: every factorization and backsolve above has already + # run. A TriangularSolve catch-all specialization whose triangular argument + # has Float32/Float64 eltype would mean some real-eltype solve silently fell + # back to LinearAlgebra/BLAS; none may exist. Complex specializations are + # expected — TriangularSolve has no complex kernels, see (c). + flagged = Union{UpperTriangular{Float32}, UpperTriangular{Float64}, + UnitLowerTriangular{Float32}, UnitLowerTriangular{Float64}} + for ca in (catchall2, catchall3), mi in Base.specializations(ca) + sig = Base.unwrap_unionall(mi.specTypes) + TA = sig.parameters[2] + @test !(TA isa Type && TA <: flagged) + end + # (c) Characterization of the known gap (a finding, not policy): complex + # panel solves resolve to the catch-all, i.e. LinearAlgebra -> LAPACK + # trtrs!/trsm. If TriangularSolve gains native complex kernels these flip + # and the routing table should be updated. + for T in (ComplexF64, ComplexF32) + A = rand(T, 100, 100) + n1 = RecursiveFactorization.nsplit(T, 100) + A11 = RecursiveFactorization.square_view(A, n1) + A12 = @view A[1:n1, (n1 + 1):100] + m = which(TriangularSolve.ldiv!, + Tuple{UnitLowerTriangular{T, typeof(A11)}, typeof(A12), Val{false}}) + @test m === catchall3 + end +end From c0cd423e496fbda1fc45d377de23b4aabee4820a Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 8 Aug 2026 15:48:47 -0400 Subject: [PATCH 4/5] Make the butterfly testset seeded and its residual bound backward-stable MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The 1e-10 absolute bound sat inside the legitimate rounding band for the n≈800 Wilkinson solves (c*n*eps*normA*normx ≈ 1e-9) and depended on suite-order RNG position for its b draws, while the butterfly transforms themselves vary per machine (VectorizedRNG streams follow SIMD width): CI observed a spurious 7.6e-10 on one runner while 550 draws on Zen 2 stay below 2.3e-11 under both the old (TS vector entry / trsv) and new (n-by-1 reshape) backsolve routes, with same-order worst cases — rounding, not breakage. Seed the testset so earlier testsets cannot shift its draws, and bound the relative residual at 1e-8, far above rounding and far below any genuine routing/pivoting failure (>=1e-5). Co-Authored-By: Chris Rackauckas Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01UGVN6qeNL2jGCtaYg1386X --- test/runtests.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 2b4fc27..3e4e18d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -142,12 +142,21 @@ function wilkinson(N) end @testset "🦋" begin + # Seeded so earlier testsets can't shift which b's this draws. The bound is + # a relative residual: the old absolute 1e-10 sat inside the backward-stable + # rounding band for n≈800 (c·n·eps·‖A‖·‖x‖ ≈ 1e-9) and failed spuriously on + # CI hardware — the butterfly transforms themselves are hardware-dependent + # (VectorizedRNG streams vary with SIMD width), so residuals near 1e-9 + # absolute are legitimate rounding, not breakage (measured: ≤2.3e-11 over + # 550 draws on Zen 2, 7.6e-10 once on a GitHub runner; genuine + # routing/pivot breakage produces ≥1e-5). + Random.seed!(1234) for i in 790 : 810 A = wilkinson(i) b = rand(i) ws = RecursiveFactorization.🦋workspace(copy(A), copy(b)) out = RecursiveFactorization.🦋solve!(ws, Val(true)) - @test norm(A * out .- b) <= 1e-10 + @test norm(A * out .- b) <= 1e-8 * norm(b) end end From 1efa602a15d8e632355a1d0de9ac09e10b6e69de Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 8 Aug 2026 16:55:35 -0400 Subject: [PATCH 5/5] Route NotIPIV vector backsolves through TriangularSolve 0.2.5's vector 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 Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01UGVN6qeNL2jGCtaYg1386X --- Project.toml | 2 +- src/lu.jl | 18 +++++------------- test/runtests.jl | 10 ++++------ 3 files changed, 10 insertions(+), 20 deletions(-) diff --git a/Project.toml b/Project.toml index 02e34c9..2adacb4 100644 --- a/Project.toml +++ b/Project.toml @@ -20,7 +20,7 @@ Polyester = "0.3.2,0.4.1, 0.5, 0.6, 0.7" PrecompileTools = "1" SparseBandedMatrices = "1" StrideArraysCore = "0.5.5" -TriangularSolve = "0.2.2" +TriangularSolve = "0.2.5" VectorizedRNG = "0.2.25" julia = "1.10" diff --git a/src/lu.jl b/src/lu.jl index 66511e2..b5a3e00 100644 --- a/src/lu.jl +++ b/src/lu.jl @@ -53,21 +53,13 @@ if CUSTOMIZABLE_PIVOT && isdefined(LinearAlgebra, :_ipiv_rows!) end end if CUSTOMIZABLE_PIVOT - # TriangularSolve's native kernels take strided *matrix* right-hand sides; - # its vector entry point defers to BLAS `trsv` above a size cutoff. A - # contiguous vector presented as an n×1 matrix stays on TriangularSolve's - # kernels at every size. view-then-reshape keeps the wrapper - # allocation-free (immutable, passed by value), unlike reshape(::Vector), - # which heap-allocates a Matrix header on every solve. - _ts_backsolve_rhs(B) = B - function _ts_backsolve_rhs(b::StridedVector{T}) where {T <: Union{Float32, Float64}} - R = reshape(view(b, :), length(b), 1) - return R isa StridedMatrix{T} ? R : b - end + # Both triangular legs run on TriangularSolve's native kernels for matrix + # *and* vector right-hand sides; the vector entry requires + # TriangularSolve >= 0.2.5 (compat-enforced), where it is BLAS-free at + # every size — older versions deferred vectors to BLAS `trsv` above n=128. function LinearAlgebra.ldiv!(A::LU{T, <:StridedMatrix, <:NotIPIV}, B::StridedVecOrMat{T}) where {T <: BlasFloat} - B′ = _ts_backsolve_rhs(B) - ldiv!(UpperTriangular(A.factors), ldiv!(UnitLowerTriangular(A.factors), B′)) + ldiv!(UpperTriangular(A.factors), ldiv!(UnitLowerTriangular(A.factors), B)) return B end end diff --git a/test/runtests.jl b/test/runtests.jl index 3e4e18d..f02ba0a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -98,12 +98,10 @@ end m = which(LinearAlgebra.ldiv!, Tuple{typeof(F), BT}) @test m.module === RecursiveFactorization end - # a contiguous vector reshapes, allocation-free, to a strided n×1 matrix - rhs = RecursiveFactorization._ts_backsolve_rhs(zeros(T, 4)) - @test rhs isa StridedMatrix{T} - RT = typeof(rhs) + # vector and matrix right-hand sides both resolve to native kernels + # (the vector methods require TriangularSolve >= 0.2.5) for W in (UnitLowerTriangular{T, MT}, UpperTriangular{T, MT}), - BT in (MT, RT) + BT in (Vector{T}, MT) m2 = which(TriangularSolve.ldiv!, Tuple{W, BT}) m3 = which(TriangularSolve.ldiv!, Tuple{W, BT, Val{true}}) @@ -115,7 +113,7 @@ end end end -@testset "NotIPIV ldiv! correctness across the TriangularSolve size cutoff" begin +@testset "NotIPIV ldiv! correctness across small-to-large sizes" begin for T in (Float64, Float32), n in (8, 64, 200, 300) A = rand(T, n, n) + T(10) * LinearAlgebra.I b = rand(T, n)