Skip to content

converted sparse/umfpack to testsets #20956

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

Closed
wants to merge 4 commits into from
Closed
Changes from 2 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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove the # characters from the description string

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

remove the # characters from the description string

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I did that

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

there are 2 of them and you only got rid of one

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

where is the second one?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I also do not see a second one – there are some testset strings with "issue #123" in them, but that's fine.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

lu::UmfpackLU{Float64}, # B::StridedMatrix{T})

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

its gone now

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