From cc126db8c902fdf6546afd16de688299d92de718 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Sat, 8 Aug 2026 09:01:41 -0400 Subject: [PATCH 1/4] 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/4] 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/4] 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/4] 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