Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

converted sparse/umfpack to testsets #20956

Closed
wants to merge 4 commits into from
Closed
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
271 changes: 142 additions & 129 deletions test/sparse/umfpack.jl
Original file line number Diff line number Diff line change
@@ -1,156 +1,169 @@
# This file is a part of Julia. License is MIT: http://julialang.org/license

se33 = speye(3)
do33 = ones(3)
@test isequal(se33 \ do33, do33)
@testset "isequal" begin
se33 = speye(3)
do33 = ones(3)
@test isequal(se33 \ do33, do33)
end

# 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 "umfpack tests" begin
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))
end
end

for Tv in (Float64, Complex128)
Ac0 = complex.(A0,A0)
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]
Ac = convert(SparseMatrixCSC{Complex128,Ti}, Ac0)
lua = lufact(Ac)
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))
@test (Diagonal(Rs) * Ac)[p,q] ≈ L * U
end
end

Ac0 = complex.(A0,A0)
for Ti in Base.uniontypes(Base.SparseArrays.UMFPACK.UMFITypes)
Ac = convert(SparseMatrixCSC{Complex128,Ti}, Ac0)
lua = lufact(Ac)
L,U,p,q,Rs = lua[:(:)]
@test (Diagonal(Rs) * Ac)[p,q] ≈ L * U
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
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)
@testset "issue #4523 complex sparse \\" begin
x = speye(2) + im * speye(2)
@test (x*(lufact(x) \ ones(2))) ≈ ones(2)
end

@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])
@testset "UMFPACK_ERROR_n_nonpositive" begin
@test_throws ArgumentError lufact(sparse(Int[], Int[], Float64[], 5, 0))
end

for T in (BigFloat, Complex{BigFloat})
@test_throws ArgumentError lufact(sparse(ones(T, 1, 1)))
@testset "issue #15099" begin
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
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)
@testset "BigFloat, Complex{BigFloat} throws error in lufact" begin
for T in (BigFloat, Complex{BigFloat})
@test_throws ArgumentError lufact(sparse(ones(T, 1, 1)))
end
end

@testset "size(::UmfpackLU)" begin
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
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)
@testset "aliased input and output array throws ArgumentError" begin
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
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]
@testset "issues #18246, #18244-lufact sparse pivot" begin
let A = speye(4)
A[1:2,1:2] = [-.01 -200; 200 .001]
F = lufact(A)
@test F[:p] == [3 ; 4 ; 2 ; 1]
end
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)
@testset "Test that A[c|t]_ldiv_B!{T<:Complex}(X::StridedMatrix{T}, lu::UmfpackLU{Float64}, B::StridedMatrix{T}) works as expected." begin
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
end