Use testsets and \approx in umfpack.jl

test/sparse/umfpack.jl

@@ -1,163 +1,169 @@
 # This file is a part of Julia. License is MIT: https://julialang.org/license
 
-se33 = speye(3)
-do33 = ones(3)
-@test isequal(se33 \ do33, do33)
-
-# based on deps/Suitesparse-4.0.2/UMFPACK/Demo/umfpack_di_demo.c
-
-using Base.SparseArrays.UMFPACK.increment!
-
-A0 = sparse(increment!([0,4,1,1,2,2,0,1,2,3,4,4]),
-           increment!([0,4,0,2,1,2,1,4,3,2,1,2]),
-           [2.,1.,3.,4.,-1.,-3.,3.,6.,2.,1.,4.,2.], 5, 5)
-
-for Tv in (Float64, Complex128)
-    for Ti in Base.uniontypes(Base.SparseArrays.UMFPACK.UMFITypes)
-        A = convert(SparseMatrixCSC{Tv,Ti}, A0)
-        lua = lufact(A)
-        @test nnz(lua) == 18
-        @test_throws KeyError lua[:Z]
-        L,U,p,q,Rs = lua[:(:)]
-        @test (Diagonal(Rs) * A)[p,q] ≈ L * U
-
-        det(lua) ≈ det(Array(A))
-
-        b = [8., 45., -3., 3., 19.]
-        x = lua\b
-        @test x ≈ float([1:5;])
-
-        @test norm(A*x-b,1) < eps(1e4)
-        z = complex.(b,zeros(b))
-        x = Base.SparseArrays.A_ldiv_B!(lua, z)
-        @test x ≈ float([1:5;])
-        @test z === x
-        y = similar(z)
-        A_ldiv_B!(y, lua, complex.(b,zeros(b)))
-        @test y ≈ x
-
-        @test norm(A*x-b,1) < eps(1e4)
-
-        b = [8., 20., 13., 6., 17.]
-        x = lua'\b
-        @test x ≈ float([1:5;])
-
-        @test norm(A'*x-b,1) < eps(1e4)
-        z = complex.(b,zeros(b))
-        x = Base.SparseArrays.Ac_ldiv_B!(lua, z)
-        @test x ≈ float([1:5;])
-        @test x === z
-        y = similar(x)
-        Base.SparseArrays.Ac_ldiv_B!(y, lua, complex.(b,zeros(b)))
-        @test y ≈ x
-
-        @test norm(A'*x-b,1) < eps(1e4)
-        x = lua.'\b
-        @test x ≈ float([1:5;])
-
-        @test norm(A.'*x-b,1) < eps(1e4)
-        x = Base.SparseArrays.At_ldiv_B!(lua,complex.(b,zeros(b)))
-        @test x ≈ float([1:5;])
-        y = similar(x)
-        Base.SparseArrays.At_ldiv_B!(y, lua,complex.(b,zeros(b)))
-        @test y ≈ x
-
-        @test norm(A.'*x-b,1) < eps(1e4)
-
-        # Element promotion and type inference
-        @inferred lua\ones(Int, size(A, 2))
+@testset "UMFPACK wrappers" begin
+    se33 = speye(3)
+    do33 = ones(3)
+    @test isequal(se33 \ do33, do33)
+
+    # based on deps/Suitesparse-4.0.2/UMFPACK/Demo/umfpack_di_demo.c
+
+    using Base.SparseArrays.UMFPACK.increment!
+
+    A0 = sparse(increment!([0,4,1,1,2,2,0,1,2,3,4,4]),
+                increment!([0,4,0,2,1,2,1,4,3,2,1,2]),
+                [2.,1.,3.,4.,-1.,-3.,3.,6.,2.,1.,4.,2.], 5, 5)
+
+    @testset "Core functionality for $Tv elements" for Tv in (Float64, Complex128)
+        # We might be able to support two index sizes one day
+        for Ti in Base.uniontypes(Base.SparseArrays.UMFPACK.UMFITypes)
+            A = convert(SparseMatrixCSC{Tv,Ti}, A0)
+            lua = lufact(A)
+            @test nnz(lua) == 18
+            @test_throws KeyError lua[:Z]
+            L,U,p,q,Rs = lua[:(:)]
+            @test (Diagonal(Rs) * A)[p,q] ≈ L * U
+
+            det(lua) ≈ det(Array(A))
+
+            b = [8., 45., -3., 3., 19.]
+            x = lua\b
+            @test x ≈ float([1:5;])
+
+            @test A*x ≈ b
+            z = complex.(b,zeros(b))
+            x = Base.SparseArrays.A_ldiv_B!(lua, z)
+            @test x ≈ float([1:5;])
+            @test z === x
+            y = similar(z)
+            A_ldiv_B!(y, lua, complex.(b,zeros(b)))
+            @test y ≈ x
+
+            @test A*x ≈ b
+
+            b = [8., 20., 13., 6., 17.]
+            x = lua'\b
+            @test x ≈ float([1:5;])
+
+            @test A'*x ≈ b
+            z = complex.(b,zeros(b))
+            x = Base.SparseArrays.Ac_ldiv_B!(lua, z)
+            @test x ≈ float([1:5;])
+            @test x === z
+            y = similar(x)
+            Base.SparseArrays.Ac_ldiv_B!(y, lua, complex.(b,zeros(b)))
+            @test y ≈ x
+
+            @test A'*x ≈ b
+            x = lua.'\b
+            @test x ≈ float([1:5;])
+
+            @test A.'*x ≈ b
+            x = Base.SparseArrays.At_ldiv_B!(lua,complex.(b,zeros(b)))
+            @test x ≈ float([1:5;])
+            y = similar(x)
+            Base.SparseArrays.At_ldiv_B!(y, lua,complex.(b,zeros(b)))
+            @test y ≈ x
+
+            @test A.'*x ≈ b
+
+            # Element promotion and type inference
+            @inferred lua\ones(Int, size(A, 2))
+        end
     end
-end
 
-Ac0 = complex.(A0,A0)
-for Ti in Base.uniontypes(Base.SparseArrays.UMFPACK.UMFITypes)
-    Ac = convert(SparseMatrixCSC{Complex128,Ti}, Ac0)
-    x  = complex.(ones(size(Ac, 1)), ones(size(Ac,1)))
-    lua = lufact(Ac)
-    L,U,p,q,Rs = lua[:(:)]
-    @test (Diagonal(Rs) * Ac)[p,q] ≈ L * U
-    b  = Ac*x
-    @test Ac\b ≈ x
-    b  = Ac'*x
-    @test Ac'\b ≈ x
-    b  = Ac.'*x
-    @test Ac.'\b ≈ x
-end
+    @testset "More tests for complex cases" begin
+        Ac0 = complex.(A0,A0)
+        for Ti in Base.uniontypes(Base.SparseArrays.UMFPACK.UMFITypes)
+            Ac = convert(SparseMatrixCSC{Complex128,Ti}, Ac0)
+            x  = complex.(ones(size(Ac, 1)), ones(size(Ac,1)))
+            lua = lufact(Ac)
+            L,U,p,q,Rs = lua[:(:)]
+            @test (Diagonal(Rs) * Ac)[p,q] ≈ L * U
+            b  = Ac*x
+            @test Ac\b ≈ x
+            b  = Ac'*x
+            @test Ac'\b ≈ x
+            b  = Ac.'*x
+            @test Ac.'\b ≈ x
+        end
+    end
 
-for elty in (Float64, Complex128)
-    for (m, n) in ((10,5), (5, 10))
-        A = sparse([1:min(m,n); rand(1:m, 10)], [1:min(m,n); rand(1:n, 10)], elty == Float64 ? randn(min(m, n) + 10) : complex.(randn(min(m, n) + 10), randn(min(m, n) + 10)))
-        F = lufact(A)
-        L, U, p, q, Rs = F[:(:)]
-        @test (Diagonal(Rs) * A)[p,q] ≈ L * U
+    @testset "Rectangular cases" for elty in (Float64, Complex128)
+        for (m, n) in ((10,5), (5, 10))
+            A = sparse([1:min(m,n); rand(1:m, 10)], [1:min(m,n); rand(1:n, 10)], elty == Float64 ? randn(min(m, n) + 10) : complex.(randn(min(m, n) + 10), randn(min(m, n) + 10)))
+            F = lufact(A)
+            L, U, p, q, Rs = F[:(:)]
+            @test (Diagonal(Rs) * A)[p,q] ≈ L * U
+        end
     end
-end
 
-#4523 - complex sparse \
-x = speye(2) + im * speye(2)
-@test (x*(lufact(x) \ ones(2))) ≈ ones(2)
-
-@test det(sparse([1,3,3,1], [1,1,3,3], [1,1,1,1])) == 0
-
-# UMFPACK_ERROR_n_nonpositive
-@test_throws ArgumentError lufact(sparse(Int[], Int[], Float64[], 5, 0))
-
-#15099
-for (Tin, Tout) in (
-        (Complex32, Complex128),
-        (Complex64, Complex128),
-        (Complex128, Complex128),
-        (Float16, Float64),
-        (Float32, Float64),
-        (Float64, Float64),
-        (Int, Float64),
-    )
-
-    F = lufact(sparse(ones(Tin, 1, 1)))
-    L = sparse(ones(Tout, 1, 1))
-    @test F[:p] == F[:q] == [1]
-    @test F[:Rs] == [1.0]
-    @test F[:L] == F[:U] == L
-    @test F[:(:)] == (L, L, [1], [1], [1.0])
-end
+    @testset "Issue #4523 - complex sparse \\" begin
+        x = speye(2) + im * speye(2)
+        @test (x*(lufact(x) \ ones(2))) ≈ ones(2)
 
-for T in (BigFloat, Complex{BigFloat})
-    @test_throws ArgumentError lufact(sparse(ones(T, 1, 1)))
-end
+        @test det(sparse([1,3,3,1], [1,1,3,3], [1,1,1,1])) == 0
+    end
 
-#size(::UmfpackLU)
-let
-    m = n = 1
-    F = lufact(sparse(ones(m, n)))
-    @test size(F) == (m, n)
-    @test size(F, 1) == m
-    @test size(F, 2) == n
-    @test size(F, 3) == 1
-    @test_throws ArgumentError size(F,-1)
-end
+    @testset "UMFPACK_ERROR_n_nonpositive" begin
+        @test_throws ArgumentError lufact(sparse(Int[], Int[], Float64[], 5, 0))
+    end
 
-let
-    a = rand(5)
-    @test_throws ArgumentError Base.SparseArrays.UMFPACK.solve!(a, lufact(speye(5,5)), a, Base.SparseArrays.UMFPACK.UMFPACK_A)
-    aa = complex(a)
-    @test_throws ArgumentError Base.SparseArrays.UMFPACK.solve!(aa, lufact(complex(speye(5,5))), aa, Base.SparseArrays.UMFPACK.UMFPACK_A)
-end
+    @testset "Issue #15099" for (Tin, Tout) in (
+            (Complex32, Complex128),
+            (Complex64, Complex128),
+            (Complex128, Complex128),
+            (Float16, Float64),
+            (Float32, Float64),
+            (Float64, Float64),
+            (Int, Float64),
+        )
+
+        F = lufact(sparse(ones(Tin, 1, 1)))
+        L = sparse(ones(Tout, 1, 1))
+        @test F[:p] == F[:q] == [1]
+        @test F[:Rs] == [1.0]
+        @test F[:L] == F[:U] == L
+        @test F[:(:)] == (L, L, [1], [1], [1.0])
+    end
 
-#18246,18244-lufact sparse pivot
-let A = speye(4)
-    A[1:2,1:2] = [-.01 -200; 200 .001]
-    F = lufact(A)
-    @test F[:p] == [3 ; 4 ; 2 ; 1]
-end
+    @testset "BigFloat not supported" for T in (BigFloat, Complex{BigFloat})
+        @test_throws ArgumentError lufact(sparse(ones(T, 1, 1)))
+    end
+
+    @testset "size(::UmfpackLU)" begin
+        m = n = 1
+        F = lufact(sparse(ones(m, n)))
+        @test size(F) == (m, n)
+        @test size(F, 1) == m
+        @test size(F, 2) == n
+        @test size(F, 3) == 1
+        @test_throws ArgumentError size(F,-1)
+    end
+
+    @testset "Test aliasing" begin
+        a = rand(5)
+        @test_throws ArgumentError Base.SparseArrays.UMFPACK.solve!(a, lufact(speye(5,5)), a, Base.SparseArrays.UMFPACK.UMFPACK_A)
+        aa = complex(a)
+        @test_throws ArgumentError Base.SparseArrays.UMFPACK.solve!(aa, lufact(complex(speye(5,5))), aa, Base.SparseArrays.UMFPACK.UMFPACK_A)
+    end
+
+    @testset "Issues #18246,18244 - lufact sparse pivot" begin
+        A = speye(4)
+        A[1:2,1:2] = [-.01 -200; 200 .001]
+        F = lufact(A)
+        @test F[:p] == [3 ; 4 ; 2 ; 1]
+    end
+
+    @testset "Test that A[c|t]_ldiv_B!{T<:Complex}(X::StridedMatrix{T}, lu::UmfpackLU{Float64}, B::StridedMatrix{T}) works as expected." begin
+        N = 10
+        p = 0.5
+        A = N*speye(N) + sprand(N, N, p)
+        X = zeros(Complex{Float64}, N, N)
+        B = complex.(rand(N, N), rand(N, N))
+        luA, lufA = lufact(A), lufact(Array(A))
+        @test A_ldiv_B!(copy(X), luA, B) ≈ A_ldiv_B!(copy(X), lufA, B)
+        @test At_ldiv_B!(copy(X), luA, B) ≈ At_ldiv_B!(copy(X), lufA, B)
+        @test Ac_ldiv_B!(copy(X), luA, B) ≈ Ac_ldiv_B!(copy(X), lufA, B)
+    end
 
-# Test that A[c|t]_ldiv_B!{T<:Complex}(X::StridedMatrix{T}, lu::UmfpackLU{Float64},
-# B::StridedMatrix{T}) works as expected.
-let N = 10, p = 0.5
-    A = N*speye(N) + sprand(N, N, p)
-    X = zeros(Complex{Float64}, N, N)
-    B = complex.(rand(N, N), rand(N, N))
-    luA, lufA = lufact(A), lufact(Array(A))
-    @test A_ldiv_B!(copy(X), luA, B) ≈ A_ldiv_B!(copy(X), lufA, B)
-    @test At_ldiv_B!(copy(X), luA, B) ≈ At_ldiv_B!(copy(X), lufA, B)
-    @test Ac_ldiv_B!(copy(X), luA, B) ≈ Ac_ldiv_B!(copy(X), lufA, B)
 end