diff --git a/src/lu.jl b/src/lu.jl index be1a9da..f75ee0e 100644 --- a/src/lu.jl +++ b/src/lu.jl @@ -1,7 +1,7 @@ using LoopVectorization using TriangularSolve: ldiv! using LinearAlgebra: BlasInt, BlasFloat, LU, UnitLowerTriangular, checknonsingular, BLAS, - LinearAlgebra, Adjoint, Transpose, UpperTriangular + LinearAlgebra, Adjoint, Transpose, UpperTriangular, AbstractVecOrMat using StrideArraysCore using Polyester: @batch @@ -46,7 +46,8 @@ if CUSTOMIZABLE_PIVOT && isdefined(LinearAlgebra, :_ipiv_cols!) end end if CUSTOMIZABLE_PIVOT && isdefined(LinearAlgebra, :_ipiv_rows!) - function LinearAlgebra._ipiv_rows!(::LU{<:Any, <:Any, NotIPIV}, ::OrdinalRange, + function LinearAlgebra._ipiv_rows!(::(LU{T, <:AbstractMatrix{T}, NotIPIV} where {T}), + ::OrdinalRange, B::StridedVecOrMat) return B end diff --git a/test/runtests.jl b/test/runtests.jl index f0d7c37..8b9467f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,7 @@ using Test import RecursiveFactorization import LinearAlgebra -using LinearAlgebra: norm, Adjoint, Transpose +using LinearAlgebra: norm, Adjoint, Transpose, ldiv! using Random Random.seed!(12) @@ -9,12 +9,24 @@ Random.seed!(12) const baselu = LinearAlgebra.lu const mylu = RecursiveFactorization.lu -function testlu(A, MF, BF) +function testlu(A, MF, BF, p) @test MF.info == BF.info - @test norm(MF.L * MF.U - A[MF.p, :], Inf) < 200sqrt(eps(real(one(float(first(A)))))) + if !iszero(MF.info) + return nothing + end + E = 20size(A, 1) * eps(real(one(float(first(A))))) + @test norm(MF.L * MF.U - A[MF.p, :], Inf) < (p ? E : 10sqrt(E)) + if ==(size(A)...) + b = ldiv!(MF, A[:, end]) + if all(isfinite, b) + n = size(A, 2) + rhs = [i == n for i in 1:n] + @test b≈rhs atol=p ? 100E : 100sqrt(E) + end + end nothing end -testlu(A::Union{Transpose, Adjoint}, MF, BF) = testlu(parent(A), parent(MF), BF) +testlu(A::Union{Transpose, Adjoint}, MF, BF, p) = testlu(parent(A), parent(MF), BF, p) @testset "Test LU factorization" begin for _p in (true, false), @@ -24,29 +36,31 @@ testlu(A::Union{Transpose, Adjoint}, MF, BF) = testlu(parent(A), parent(MF), BF) p = Val(_p) for (i, s) in enumerate([1:10; 50:80:200; 300]) iseven(i) && (p = RecursiveFactorization.to_stdlib_pivot(p)) - siz = (s, s + 2) - @info("size: $(siz[1]) × $(siz[2]), T = $T, p = $_p") - if isconcretetype(T) - A = rand(T, siz...) - else - _A = rand(siz...) - A = Matrix{T}(undef, siz...) - copyto!(A, _A) + for m in (s, s + 2) + siz = (s, m) + @info("size: $(siz[1]) × $(siz[2]), T = $T, p = $_p") + if isconcretetype(T) + A = rand(T, siz...) + else + _A = rand(siz...) + A = Matrix{T}(undef, siz...) + copyto!(A, _A) + end + MF = mylu(A, p) + BF = baselu(A, p) + testlu(A, MF, BF, _p) + testlu(A, mylu(A, p, Val(false)), BF, false) + A′ = permutedims(A) + MF′ = mylu(A′', p) + testlu(A′', MF′, BF, _p) + testlu(A′', mylu(A′', p, Val(false)), BF, false) + i = rand(1:s) # test `MF.info` + A[:, i] .= 0 + MF = mylu(A, p, check = false) + BF = baselu(A, p, check = false) + testlu(A, MF, BF, _p) + testlu(A, mylu(A, p, Val(false), check = false), BF, false) end - MF = mylu(A, p) - BF = baselu(A, p) - testlu(A, MF, BF) - testlu(A, mylu(A, p, Val(false)), BF) - A′ = permutedims(A) - MF′ = mylu(A′', p) - testlu(A′', MF′, BF) - testlu(A′', mylu(A′', p, Val(false)), BF) - i = rand(1:s) # test `MF.info` - A[:, i] .= 0 - MF = mylu(A, p, check = false) - BF = baselu(A, p, check = false) - testlu(A, MF, BF) - testlu(A, mylu(A, p, Val(false), check = false), BF) end end end