More robust tests

Author: YingboMa

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

test/runtests.jl

@@ -1,20 +1,32 @@
 using Test
 import RecursiveFactorization
 import LinearAlgebra
-using LinearAlgebra: norm, Adjoint, Transpose
+using LinearAlgebra: norm, Adjoint, Transpose, ldiv!
 using Random
 
 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