diff --git a/Project.toml b/Project.toml index 270a9f9..2adacb4 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" @@ -20,13 +20,14 @@ 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" [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 78854ad..b5a3e00 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 @@ -53,10 +53,14 @@ if CUSTOMIZABLE_PIVOT && isdefined(LinearAlgebra, :_ipiv_rows!) end end if CUSTOMIZABLE_PIVOT + # 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} - tri = @inbounds square_view(A.factors, size(A.factors, 1)) ldiv!(UpperTriangular(A.factors), ldiv!(UnitLowerTriangular(A.factors), B)) + return B end end @@ -100,7 +104,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..f02ba0a 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) @@ -65,6 +67,66 @@ 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 + +@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 + # 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 (Vector{T}, MT) + + 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 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) + 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 @@ -78,12 +140,67 @@ 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 + +@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