Skip to content

Commit

Permalink
Reduce allocation of dense arrays on-the-fly in linalg tests (#485)
Browse files Browse the repository at this point in the history
  • Loading branch information
@dkarrasch
dkarrasch authored Jan 7, 2024
1 parent c73d6e3 commit 63459e5
Show file tree
Hide file tree
Showing 2 changed files with 111 additions and 99 deletions.
197 changes: 104 additions & 93 deletions test/linalg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,13 +83,15 @@ end

@testset "Sparse promotion in sparse matmul" begin
A = SparseMatrixCSC{Float32, Int8}(2, 2, Int8[1, 2, 3], Int8[1, 2], Float32[1., 2.])
MA = Array(A)
B = SparseMatrixCSC{ComplexF32, Int32}(2, 2, Int32[1, 2, 3], Int32[1, 2], ComplexF32[1. + im, 2. - im])
@test A*transpose(B) ≈ Array(A) * transpose(Array(B))
@test A*adjoint(B) ≈ Array(A) * adjoint(Array(B))
@test transpose(A)*B ≈ transpose(Array(A)) * Array(B)
@test transpose(A)*transpose(B) ≈ transpose(Array(A)) * transpose(Array(B))
@test adjoint(B)*A ≈ adjoint(Array(B)) * Array(A)
@test adjoint(B)*adjoint(complex.(A)) ≈ adjoint(Array(B)) * adjoint(Array(complex.(A)))
MB = Array(B)
@test A*transpose(B) ≈ MA * transpose(MB)
@test A*adjoint(B) ≈ MA * adjoint(MB)
@test transpose(A)*B ≈ transpose(MA) * MB
@test transpose(A)*transpose(B) ≈ transpose(MA) * transpose(MB)
@test adjoint(B)*A ≈ adjoint(MB) * MA
@test adjoint(B)*adjoint(complex.(A)) ≈ adjoint(MB) * adjoint(Array(complex.(A)))
end

@testset "multiplication of triangular sparse and dense matrices" begin
Expand Down Expand Up @@ -226,10 +228,11 @@ end

@testset "UniformScaling" begin
local A = sprandn(10, 10, 0.5)
@test A + I == Array(A) + I
@test I + A == I + Array(A)
@test A - I == Array(A) - I
@test I - A == I - Array(A)
MA = Array(A)
@test A + I == MA + I
@test I + A == I + MA
@test A - I == MA - I
@test I - A == I - MA
end

@testset "unary minus for SparseMatrixCSC{Bool}" begin
Expand All @@ -240,25 +243,30 @@ end

@testset "sparse matrix norms" begin
Ac = sprandn(10,10,.1) + im* sprandn(10,10,.1)
MAc = Array(Ac)
Ar = sprandn(10,10,.1)
Ai = ceil.(Int,Ar*100)
@test opnorm(Ac,1) ≈ opnorm(Array(Ac),1)
@test opnorm(Ac,Inf) ≈ opnorm(Array(Ac),Inf)
@test norm(Ac) ≈ norm(Array(Ac))
@test opnorm(Ar,1) ≈ opnorm(Array(Ar),1)
@test opnorm(Ar,Inf) ≈ opnorm(Array(Ar),Inf)
@test norm(Ar) ≈ norm(Array(Ar))
@test opnorm(Ai,1) ≈ opnorm(Array(Ai),1)
@test opnorm(Ai,Inf) ≈ opnorm(Array(Ai),Inf)
@test norm(Ai) ≈ norm(Array(Ai))
MAr = Array(Ar)
Ai = ceil.(Int, Ar*100)
MAi = Array(Ai)
@test opnorm(Ac,1) ≈ opnorm(MAc,1)
@test opnorm(Ac,Inf) ≈ opnorm(MAc,Inf)
@test norm(Ac) ≈ norm(MAc)
@test opnorm(Ar,1) ≈ opnorm(MAr,1)
@test opnorm(Ar,Inf) ≈ opnorm(MAr,Inf)
@test norm(Ar) ≈ norm(MAr)
@test opnorm(Ai,1) ≈ opnorm(MAi,1)
@test opnorm(Ai,Inf) ≈ opnorm(MAi,Inf)
@test norm(Ai) ≈ norm(MAi)
Ai = trunc.(Int, Ar*100)
@test opnorm(Ai,1) ≈ opnorm(Array(Ai),1)
@test opnorm(Ai,Inf) ≈ opnorm(Array(Ai),Inf)
@test norm(Ai) ≈ norm(Array(Ai))
MAi = Array(Ai)
@test opnorm(Ai,1) ≈ opnorm(MAi,1)
@test opnorm(Ai,Inf) ≈ opnorm(MAi,Inf)
@test norm(Ai) ≈ norm(MAi)
Ai = round.(Int, Ar*100)
@test opnorm(Ai,1) ≈ opnorm(Array(Ai),1)
@test opnorm(Ai,Inf) ≈ opnorm(Array(Ai),Inf)
@test norm(Ai) ≈ norm(Array(Ai))
MAi = Array(Ai)
@test opnorm(Ai,1) ≈ opnorm(MAi,1)
@test opnorm(Ai,Inf) ≈ opnorm(MAi,Inf)
@test norm(Ai) ≈ norm(MAi)
# make certain entries in nzval beyond
# the range specified in colptr do not
# impact norm of a sparse matrix
Expand All @@ -269,16 +277,18 @@ end

# Test (m x 1) sparse matrix
colM = sprandn(10, 1, 0.6)
@test opnorm(colM, 1) ≈ opnorm(Array(colM), 1)
@test opnorm(colM) ≈ opnorm(Array(colM))
@test opnorm(colM, Inf) ≈ opnorm(Array(colM), Inf)
McolM = Array(colM)
@test opnorm(colM, 1) ≈ opnorm(McolM, 1)
@test opnorm(colM) ≈ opnorm(McolM)
@test opnorm(colM, Inf) ≈ opnorm(McolM, Inf)
@test_throws ArgumentError opnorm(colM, 3)

# Test (1 x n) sparse matrix
rowM = sprandn(1, 10, 0.6)
@test opnorm(rowM, 1) ≈ opnorm(Array(rowM), 1)
@test opnorm(rowM) ≈ opnorm(Array(rowM))
@test opnorm(rowM, Inf) ≈ opnorm(Array(rowM), Inf)
MrowM = Array(rowM)
@test opnorm(rowM, 1) ≈ opnorm(MrowM, 1)
@test opnorm(rowM) ≈ opnorm(MrowM)
@test opnorm(rowM, Inf) ≈ opnorm(MrowM, Inf)
@test_throws ArgumentError opnorm(rowM, 3)
end

Expand All @@ -293,20 +303,21 @@ end
if elty <: Complex
dd+=im*convert(Vector{elty}, randn(n))
end
D = Diagonal(dd)
b = rand(elty, n, n)
b = sparse(b)
@test ldiv!(D, copy(b)) ≈ Array(D)\Array(b)
D = Diagonal(dd); MD = Array(D)
bd = rand(elty, n, n)
b = sparse(bd)
@test ldiv!(D, copy(b)) ≈ MD\bd
@test_throws SingularException ldiv!(Diagonal(zeros(elty, n)), copy(b))
b = rand(elty, n+1, n+1)
b = sparse(b)
@test_throws DimensionMismatch ldiv!(D, copy(b))
b = view(rand(elty, n+1), Vector(1:n+1))
@test_throws DimensionMismatch ldiv!(D, b)
for b in (sparse(rand(elty,n,n)), sparse(rand(elty,n)))
@test lmul!(copy(D), copy(b)) ≈ Array(D)*Array(b)
@test lmul!(transpose(copy(D)), copy(b)) ≈ transpose(Array(D))*Array(b)
@test lmul!(adjoint(copy(D)), copy(b)) ≈ Array(D)'*Array(b)
bd = Array(b)
@test lmul!(copy(D), copy(b)) ≈ MD*bd
@test lmul!(transpose(copy(D)), copy(b)) ≈ transpose(MD)*bd
@test lmul!(adjoint(copy(D)), copy(b)) ≈ MD'*bd
end
end
end
Expand Down Expand Up @@ -662,29 +673,29 @@ end
(100, 0.01, 100, 0.01, 20),
(100, 0.1, 100, 0.2, 100),
)
a = sprand(m, n, p)
b = sprand(n, k, q)
a = sprand(m, n, p); ad = Array(a)
b = sprand(n, k, q); bd = Array(b)
as = sparse(a')
bs = sparse(b')
ab = a * b
aab = Array(a) * Array(b)
aab = ad * bd
@test maximum(abs.(ab - aab)) < 100*eps()
@test a*bs' == ab
@test as'*b == ab
@test as'*bs' == ab
f = Diagonal(rand(n))
@test Array(a*f) == Array(a)*f
@test Array(f*b) == f*Array(b)
@test Array(a*f) == ad*f
@test Array(f*b) == f*bd
A = rand(2n, 2n)
sA = view(A, 1:2:2n, 1:2:2n)
@test Array((sA*b)::Matrix) ≈ Array(sA)*Array(b)
@test Array((a*sA)::Matrix) ≈ Array(a)*Array(sA)
@test Array((sA'b)::Matrix) ≈ Array(sA')*Array(b)
c = sprandn(ComplexF32, n, n, q)
@test Array((sA*c')::Matrix) ≈ Array(sA)*Array(c)'
@test Array((c'*sA)::Matrix) ≈ Array(c)'*Array(sA)
@test Array((sA'c)::Matrix) ≈ Array(sA')*Array(c)
@test Array((sA'c')::Matrix) ≈ Array(sA')*Array(c)'
sA = view(A, 1:2:2n, 1:2:2n); dA = Array(sA)
@test (sA*b)::Matrix ≈ dA*bd
@test (a*sA)::Matrix ≈ ad*dA
@test (sA'b)::Matrix ≈ dA'*bd
c = sprandn(ComplexF32, n, n, q); cd = Array(c)
@test (sA*c')::Matrix ≈ dA*cd'
@test (c'*sA)::Matrix ≈ cd'*dA
@test (sA'c)::Matrix ≈ dA'*cd
@test (sA'c')::Matrix ≈ dA'*cd'
end
end

Expand Down Expand Up @@ -746,16 +757,16 @@ end
d = sparse(d_di); d_d = Array(d_di)
# mat ⊗ mat
for t in (identity, adjoint, transpose)
@test Array(kron(t(a), b)::SparseMatrixCSC) == kron(t(a_d), b_d)
@test Array(kron(a, t(b))::SparseMatrixCSC) == kron(a_d, t(b_d))
@test Array(kron(t(a), t(b))::SparseMatrixCSC) == kron(t(a_d), t(b_d))
@test Array(kron(t(a), b_d)::SparseMatrixCSC) == kron(t(a_d), b_d)
@test Array(kron(a_d, t(b))::SparseMatrixCSC) == kron(a_d, t(b_d))
@test Array(kron(t(a), c_di)::SparseMatrixCSC) == kron(t(a_d), c_d)
@test Array(kron(a, t(c_di))::SparseMatrixCSC) == kron(a_d, t(c_d))
@test Array(kron(t(a), t(c_di))::SparseMatrixCSC) == kron(t(a_d), t(c_d))
@test Array(kron(c_di, y)::SparseMatrixCSC) == kron(c_di, y_d)
@test Array(kron(x, d_di)::SparseMatrixCSC) == kron(x_d, d_di)
@test kron(t(a), b)::SparseMatrixCSC == kron(t(a_d), b_d)
@test kron(a, t(b))::SparseMatrixCSC == kron(a_d, t(b_d))
@test kron(t(a), t(b))::SparseMatrixCSC == kron(t(a_d), t(b_d))
@test kron(t(a), b_d)::SparseMatrixCSC == kron(t(a_d), b_d)
@test kron(a_d, t(b))::SparseMatrixCSC == kron(a_d, t(b_d))
@test kron(t(a), c_di)::SparseMatrixCSC == kron(t(a_d), c_d)
@test kron(a, t(c_di))::SparseMatrixCSC == kron(a_d, t(c_d))
@test kron(t(a), t(c_di))::SparseMatrixCSC == kron(t(a_d), t(c_d))
@test kron(c_di, y)::SparseMatrixCSC == kron(c_di, y_d)
@test kron(x, d_di)::SparseMatrixCSC == kron(x_d, d_di)
end
end
# vec ⊗ vec
Expand All @@ -764,22 +775,22 @@ end
@test Vector(kron(x, y_d)::SparseVector) == kron(x_d, y_d)
for t in (identity, adjoint, transpose)
# mat ⊗ vec
@test Array(kron(t(a), y)::SparseMatrixCSC) == kron(t(a_d), y_d)
@test Array(kron(t(a_d), y)::SparseMatrixCSC) == kron(t(a_d), y_d)
@test Array(kron(t(a), y_d)::SparseMatrixCSC) == kron(t(a_d), y_d)
@test kron(t(a), y)::SparseMatrixCSC == kron(t(a_d), y_d)
@test kron(t(a_d), y)::SparseMatrixCSC == kron(t(a_d), y_d)
@test kron(t(a), y_d)::SparseMatrixCSC == kron(t(a_d), y_d)
# vec ⊗ mat
@test Array(kron(x, t(b))::SparseMatrixCSC) == kron(x_d, t(b_d))
@test Array(kron(x_d, t(b))::SparseMatrixCSC) == kron(x_d, t(b_d))
@test Array(kron(x, t(b_d))::SparseMatrixCSC) == kron(x_d, t(b_d))
@test kron(x, t(b))::SparseMatrixCSC == kron(x_d, t(b_d))
@test kron(x_d, t(b))::SparseMatrixCSC == kron(x_d, t(b_d))
@test kron(x, t(b_d))::SparseMatrixCSC == kron(x_d, t(b_d))
end
# vec ⊗ vec'
@test Array(kron(v, y')::SparseMatrixCSC) == kron(v_d, y_d')
@test Array(kron(x, y')::SparseMatrixCSC) == kron(x_d, y_d')
@test kron(v, y')::SparseMatrixCSC == kron(v_d, y_d')
@test kron(x, y')::SparseMatrixCSC == kron(x_d, y_d')
# test different types
z = convert(SparseVector{Float16, Int8}, y); z_d = Vector(z)
@test Vector(kron(x, z)) == kron(x_d, z_d)
@test Array(kron(a, z)) == kron(a_d, z_d)
@test Array(kron(z, b)) == kron(z_d, b_d)
@test kron(a, z) == kron(a_d, z_d)
@test kron(z, b) == kron(z_d, b_d)
# test bounds checks
@test_throws DimensionMismatch kron!(copy(a), a, b)
@test_throws DimensionMismatch kron!(copy(x), x, y)
Expand All @@ -790,30 +801,30 @@ end
@testset "sparse Frobenius dot/inner product" begin
full_view = M -> view(M, :, :)
for i = 1:5
A = sprand(ComplexF64,10,15,0.4)
B = sprand(ComplexF64,10,15,0.5)
C = rand(10,15) .> 0.3
@test dot(A,B) ≈ dot(Matrix(A), Matrix(B))
@test dot(A,B) ≈ dot(A, Matrix(B))
@test dot(A,B) ≈ dot(Matrix(A), B)
@test dot(A,C) ≈ dot(Matrix(A), C)
@test dot(C,A) ≈ dot(C, Matrix(A))
A = sprand(ComplexF64,10,15,0.4); MA = Matrix(A)
B = sprand(ComplexF64,10,15,0.5); MB = Matrix(B)
C = rand(10,15) .> 0.3; MC = Matrix(C)
@test dot(A,B) ≈ dot(MA, MB)
@test dot(A,B) ≈ dot(A, MB)
@test dot(A,B) ≈ dot(MA, B)
@test dot(A,C) ≈ dot(MA, C)
@test dot(C,A) ≈ dot(C, MA)
# square matrices required by most linear algebra wrappers
SA = A * A'
SB = B * B'
SC = C * C'
SA = A * A'; MSA = Matrix(SA)
SB = B * B'; MSB = Matrix(SB)
SC = C * C'; MSC = Matrix(SC)
for W in (full_view, LowerTriangular, UpperTriangular, UpperHessenberg, Symmetric, Hermitian)
WA = W(Matrix(SA))
WB = W(Matrix(SB))
WC = W(Matrix(SC))
@test dot(WA,SB) ≈ dot(WA, Matrix(SB))
@test dot(SA,WB) ≈ dot(Matrix(SA), WB)
@test dot(SA,WC) ≈ dot(Matrix(SA), WC)
WA = W(MSA)
WB = W(MSB)
WC = W(MSC)
@test dot(WA,SB) ≈ dot(WA, MSB)
@test dot(SA,WB) ≈ dot(MSA, WB)
@test dot(SA,WC) ≈ dot(MSA, WC)
end
for W in (transpose, adjoint)
WA = W(Matrix(A))
WB = W(Matrix(B))
WC = W(Matrix(C))
WA = W(MA)
WB = W(MB)
WC = W(MC)
TA = copy(W(A))
TB = copy(W(B))
@test dot(WA,TB) ≈ dot(WA, Matrix(TB))
Expand Down
13 changes: 7 additions & 6 deletions test/sparsematrix_ops.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,12 +101,13 @@ do33 = fill(1.,3)
end
@testset "binary operations on sparse matrices with union eltype" begin
A = sparse([1,2,1], [1,1,2], Union{Int, Missing}[1, missing, 0])
MA = Array(A)
for fun in (+, -, *, min, max)
if fun in (+, -)
@test collect(skipmissing(Array(fun(A, A)))) == collect(skipmissing(Array(fun(Array(A), Array(A)))))
@test collect(skipmissing(Array(fun(A, A)))) == collect(skipmissing(Array(fun(MA, MA))))
end
@test collect(skipmissing(Array(map(fun, A, A)))) == collect(skipmissing(map(fun, Array(A), Array(A))))
@test collect(skipmissing(Array(broadcast(fun, A, A)))) == collect(skipmissing(broadcast(fun, Array(A), Array(A))))
@test collect(skipmissing(Array(map(fun, A, A)))) == collect(skipmissing(map(fun, MA, MA)))
@test collect(skipmissing(Array(broadcast(fun, A, A)))) == collect(skipmissing(broadcast(fun, MA, MA)))
end
b = convert(SparseMatrixCSC{Union{Float64, Missing}}, sprandn(Float64, 20, 10, 0.2)); b[rand(1:200, 3)] .= missing
C = convert(SparseMatrixCSC{Union{Float64, Missing}}, sprandn(Float64, 20, 10, 0.9)); C[rand(1:200, 3)] .= missing
Expand Down Expand Up @@ -263,10 +264,10 @@ end

@testset "findall" begin
# issue described in https://groups.google.com/d/msg/julia-users/Yq4dh8NOWBQ/GU57L90FZ3EJ
A = sparse(I, 5, 5)
@test findall(A) == findall(x -> x == true, A) == findall(Array(A))
A = sparse(I, 5, 5); MA = Array(A)
@test findall(A) == findall(x -> x == true, A) == findall(MA)
# Non-stored entries are true
@test findall(x -> x == false, A) == findall(x -> x == false, Array(A))
@test findall(x -> x == false, A) == findall(x -> x == false, MA)

# Not all stored entries are true
@test findall(sparse([true false])) == [CartesianIndex(1, 1)]
Expand Down

0 comments on commit 63459e5

Please sign in to comment.