diff --git a/src/SparseArrays.jl b/src/SparseArrays.jl index 727abddb..7db213b6 100644 --- a/src/SparseArrays.jl +++ b/src/SparseArrays.jl @@ -38,21 +38,9 @@ include("linalg.jl") include("deprecated.jl") -# temporarily moved here and commented out from from base/linalg/diagonal.jl, base/linalg/tridiag.jl -# and base/linalg/bidiag.jl due to their usage of spzeros -similar(B::Bidiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...) -similar(D::Diagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...) -similar(S::SymTridiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...) -similar(M::Tridiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...) - zero(a::AbstractSparseArray) = spzeros(eltype(a), size(a)...) -const BiTriSym = Union{Bidiagonal,SymTridiagonal,Tridiagonal} -function *(A::BiTriSym, B::BiTriSym) - TS = promote_op(matprod, eltype(A), eltype(B)) - mul!(similar(A, TS, size(A)...), A, B) -end - -LinearAlgebra.diagzero(D::Diagonal{<:AbstractSparseMatrix{T}},i,j) where {T} = spzeros(T, size(D.diag[i], 1), size(D.diag[j], 2)) +LinearAlgebra.diagzero(D::Diagonal{<:AbstractSparseMatrix{T}},i,j) where {T} = + spzeros(T, size(D.diag[i], 1), size(D.diag[j], 2)) end diff --git a/src/sparsevector.jl b/src/sparsevector.jl index c8f4eaaa..edf2dd40 100644 --- a/src/sparsevector.jl +++ b/src/sparsevector.jl @@ -4,6 +4,7 @@ import Base: sort, findall, copy! import LinearAlgebra: promote_to_array_type, promote_to_arrays_ +using LinearAlgebra: _SpecialArrays, _DenseConcatGroup ### The SparseVector @@ -1062,30 +1063,16 @@ vcat(X::Union{Vector,SparseVector}...) = vcat(map(sparse, X)...) ### Concatenation of un/annotated sparse/special/dense vectors/matrices -# TODO: These methods and definitions should be moved to a more appropriate location, -# particularly some future equivalent of base/linalg/special.jl dedicated to interactions -# between a broader set of matrix types. - -# TODO: A definition similar to the third exists in base/linalg/bidiag.jl. These definitions -# should be consolidated in a more appropriate location, e.g. base/linalg/special.jl. const _SparseArrays = Union{SparseVector, AbstractSparseMatrixCSC, Adjoint{<:Any,<:SparseVector}, Transpose{<:Any,<:SparseVector}} -const _SpecialArrays = Union{Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal} const _SparseConcatArrays = Union{_SpecialArrays, _SparseArrays} const _Symmetric_SparseConcatArrays{T,A<:_SparseConcatArrays} = Symmetric{T,A} const _Hermitian_SparseConcatArrays{T,A<:_SparseConcatArrays} = Hermitian{T,A} const _Triangular_SparseConcatArrays{T,A<:_SparseConcatArrays} = LinearAlgebra.AbstractTriangular{T,A} const _Annotated_SparseConcatArrays = Union{_Triangular_SparseConcatArrays, _Symmetric_SparseConcatArrays, _Hermitian_SparseConcatArrays} - -const _Symmetric_DenseArrays{T,A<:Matrix} = Symmetric{T,A} -const _Hermitian_DenseArrays{T,A<:Matrix} = Hermitian{T,A} -const _Triangular_DenseArrays{T,A<:Matrix} = LinearAlgebra.AbstractTriangular{T,A} -const _Annotated_DenseArrays = Union{_Triangular_DenseArrays, _Symmetric_DenseArrays, _Hermitian_DenseArrays} -const _Annotated_Typed_DenseArrays{T} = Union{_Triangular_DenseArrays{T}, _Symmetric_DenseArrays{T}, _Hermitian_DenseArrays{T}} - -const _SparseConcatGroup = Union{Number, Vector, Adjoint{<:Any,<:Vector}, Transpose{<:Any,<:Vector}, Matrix, _SparseConcatArrays, _Annotated_SparseConcatArrays, _Annotated_DenseArrays} -const _DenseConcatGroup = Union{Number, Vector, Adjoint{<:Any,<:Vector}, Transpose{<:Any,<:Vector}, Matrix, _Annotated_DenseArrays} -const _TypedDenseConcatGroup{T} = Union{Vector{T}, Adjoint{T,Vector{T}}, Transpose{T,Vector{T}}, Matrix{T}, _Annotated_Typed_DenseArrays{T}} +# It's important that _SparseConcatGroup is a larger union than _DenseConcatGroup to make +# sparse cat-methods less specific and to kick in only if there is some sparse array present +const _SparseConcatGroup = Union{_DenseConcatGroup, _SparseConcatArrays, _Annotated_SparseConcatArrays} # Concatenations involving un/annotated sparse/special matrices/vectors should yield sparse arrays _makesparse(x::Number) = x @@ -1123,23 +1110,8 @@ _hvcat_rows(::Tuple{}, X::_SparseConcatGroup...) = () # make sure UniformScaling objects are converted to sparse matrices for concatenation promote_to_array_type(A::Tuple{Vararg{Union{_SparseConcatGroup,UniformScaling}}}) = SparseMatrixCSC -promote_to_array_type(A::Tuple{Vararg{Union{_DenseConcatGroup,UniformScaling}}}) = Matrix promote_to_arrays_(n::Int, ::Type{SparseMatrixCSC}, J::UniformScaling) = sparse(J, n, n) -# Concatenations strictly involving un/annotated dense matrices/vectors should yield dense arrays -Base._cat(dims, xs::_DenseConcatGroup...) = Base.cat_t(promote_eltype(xs...), xs...; dims=dims) -vcat(A::Vector...) = Base.typed_vcat(promote_eltype(A...), A...) -vcat(A::_DenseConcatGroup...) = Base.typed_vcat(promote_eltype(A...), A...) -hcat(A::Vector...) = Base.typed_hcat(promote_eltype(A...), A...) -hcat(A::_DenseConcatGroup...) = Base.typed_hcat(promote_eltype(A...), A...) -hvcat(rows::Tuple{Vararg{Int}}, xs::_DenseConcatGroup...) = Base.typed_hvcat(promote_eltype(xs...), rows, xs...) -# For performance, specially handle the case where the matrices/vectors have homogeneous eltype -Base._cat(dims, xs::_TypedDenseConcatGroup{T}...) where {T} = Base.cat_t(T, xs...; dims=dims) -vcat(A::_TypedDenseConcatGroup{T}...) where {T} = Base.typed_vcat(T, A...) -hcat(A::_TypedDenseConcatGroup{T}...) where {T} = Base.typed_hcat(T, A...) -hvcat(rows::Tuple{Vararg{Int}}, xs::_TypedDenseConcatGroup{T}...) where {T} = Base.typed_hvcat(T, rows, xs...) - - ### math functions ### Unary Map diff --git a/test/sparse.jl b/test/sparse.jl index f68b84ee..b7c946c3 100644 --- a/test/sparse.jl +++ b/test/sparse.jl @@ -132,6 +132,8 @@ end # constructor methods well-exercised by the immediately preceding testset @test sparse(2I, 3, 4)::SparseMatrixCSC{Int,Int} == Matrix(2I, 3, 4) @test sparse(2I, (3, 4))::SparseMatrixCSC{Int,Int} == Matrix(2I, 3, 4) + @test sparse(3I, 4, 5) == sparse(1:4, 1:4, 3, 4, 5) + @test sparse(3I, 5, 4) == sparse(1:4, 1:4, 3, 5, 4) end se33 = SparseMatrixCSC{Float64}(I, 3, 3) @@ -154,6 +156,33 @@ do33 = fill(1.,3) @test_throws DimensionMismatch map(|, sqrboolmat, colboolmat) @test_throws DimensionMismatch map(xor, sqrboolmat, colboolmat) end + + # ascertain inference friendliness, ref. https://github.com/JuliaLang/julia/pull/25083#issuecomment-353031641 + sparsevec = SparseVector([1.0, 2.0, 3.0]) + @test map(-, Adjoint(sparsevec), Adjoint(sparsevec)) isa Adjoint{Float64,SparseVector{Float64,Int}} + @test map(-, Transpose(sparsevec), Transpose(sparsevec)) isa Transpose{Float64,SparseVector{Float64,Int}} + @test broadcast(-, Adjoint(sparsevec), Adjoint(sparsevec)) isa Adjoint{Float64,SparseVector{Float64,Int}} + @test broadcast(-, Transpose(sparsevec), Transpose(sparsevec)) isa Transpose{Float64,SparseVector{Float64,Int}} + @test broadcast(+, Adjoint(sparsevec), 1.0, Adjoint(sparsevec)) isa Adjoint{Float64,SparseVector{Float64,Int}} + @test broadcast(+, Transpose(sparsevec), 1.0, Transpose(sparsevec)) isa Transpose{Float64,SparseVector{Float64,Int}} + + @testset "binary ops with matrices" begin + λ = complex(randn(),randn()) + J = UniformScaling(λ) + B = bitrand(2, 2) + @test B + I == B + Matrix(I, size(B)) + @test I + B == B + Matrix(I, size(B)) + AA = randn(2, 2) + for SS in (sprandn(3,3, 0.5), sparse(Int(1)I, 3, 3)) + for S in (SS, view(SS, 1:3, 1:3)) + @test @inferred(I*S) !== S # Don't alias + @test @inferred(S*I) !== S # Don't alias + + @test @inferred(S*J) == S*λ + @test @inferred(J*S) == S*λ + end + end + end end @testset "Issue #30006" begin @@ -194,6 +223,14 @@ end @test nnz(blockdiag()) == 0 end + @testset "Diagonal of sparse matrices" begin + s = sparse([1 2; 3 4]) + D = Diagonal([s, s]) + @test D[1, 1] == s + @test D[1, 2] == zero(s) + @test isa(D[2, 1], SparseMatrixCSC) + end + @testset "concatenation promotion" begin sz41_f32 = spzeros(Float32, 4, 1) se33_i32 = sparse(Int32(1)I, 3, 3) @@ -213,6 +250,141 @@ end @test [a[1:2,1:2] a[1:2,3:4]; a[3:5,1] [a[3:4,2:4]; a[5:5,2:4]]] == a end end + + # should all yield sparse arrays + @testset "concatenations of combinations of special and other matrix types" begin + N = 4 + diagmat = Diagonal(1:N) + bidiagmat = Bidiagonal(1:N, 1:(N-1), :U) + tridiagmat = Tridiagonal(1:(N-1), 1:N, 1:(N-1)) + symtridiagmat = SymTridiagonal(1:N, 1:(N-1)) + specialmats = (diagmat, bidiagmat, tridiagmat, symtridiagmat) + # Test concatenating pairwise combinations of special matrices with sparse matrices, + # dense matrices, or dense vectors + spmat = spdiagm(0 => fill(1., N)) + spvec = sparse(fill(1., N)) + for specialmat in specialmats + # --> Tests applicable only to pairs of matrices + @test issparse(vcat(specialmat, spmat)) + @test issparse(vcat(spmat, specialmat)) + # --> Tests applicable also to pairs including vectors + for specialmat in specialmats, othermatorvec in (spmat, spvec) + @test issparse(hcat(specialmat, othermatorvec)) + @test issparse(hcat(othermatorvec, specialmat)) + @test issparse(hvcat((2,), specialmat, othermatorvec)) + @test issparse(hvcat((2,), othermatorvec, specialmat)) + @test issparse(cat(specialmat, othermatorvec; dims=(1,2))) + @test issparse(cat(othermatorvec, specialmat; dims=(1,2))) + end + end + end + + # Test that concatenations of annotated sparse/special matrix types with other matrix + # types yield sparse arrays, and that the code which effects that does not make concatenations + # strictly involving un/annotated dense matrices yield sparse arrays + @testset "concatenations of annotated types" begin + N = 4 + # The tested annotation types + testfull = Bool(parse(Int,(get(ENV, "JULIA_TESTFULL", "0")))) + utriannotations = (UpperTriangular, UnitUpperTriangular) + ltriannotations = (LowerTriangular, UnitLowerTriangular) + triannotations = (utriannotations..., ltriannotations...) + symannotations = (Symmetric, Hermitian) + annotations = testfull ? (triannotations..., symannotations...) : (LowerTriangular, Symmetric) + # Concatenations involving these types, un/annotated, should yield sparse arrays + spvec = spzeros(N) + spmat = sparse(1.0I, N, N) + diagmat = Diagonal(1:N) + bidiagmat = Bidiagonal(1:N, 1:(N-1), :U) + tridiagmat = Tridiagonal(1:(N-1), 1:N, 1:(N-1)) + symtridiagmat = SymTridiagonal(1:N, 1:(N-1)) + sparseconcatmats = testfull ? (spmat, diagmat, bidiagmat, tridiagmat, symtridiagmat) : (spmat, diagmat) + # Concatenations involving strictly these types, un/annotated, should yield dense arrays + densevec = fill(1., N) + densemat = fill(1., N, N) + # Annotated collections + annodmats = [annot(densemat) for annot in annotations] + annospcmats = [annot(spmat) for annot in annotations] + # Test that concatenations of pairwise combinations of annotated sparse/special + # yield sparse matrices + for annospcmata in annospcmats, annospcmatb in annospcmats + @test issparse(vcat(annospcmata, annospcmatb)) + @test issparse(hcat(annospcmata, annospcmatb)) + @test issparse(hvcat((2,), annospcmata, annospcmatb)) + @test issparse(cat(annospcmata, annospcmatb; dims=(1,2))) + end + # Test that concatenations of pairwise combinations of annotated sparse/special + # matrices and other matrix/vector types yield sparse matrices + for annospcmat in annospcmats + # --> Tests applicable to pairs including only matrices + for othermat in (densemat, annodmats..., sparseconcatmats...) + @test issparse(vcat(annospcmat, othermat)) + @test issparse(vcat(othermat, annospcmat)) + end + # --> Tests applicable to pairs including other vectors or matrices + for other in (spvec, densevec, densemat, annodmats..., sparseconcatmats...) + @test issparse(hcat(annospcmat, other)) + @test issparse(hcat(other, annospcmat)) + @test issparse(hvcat((2,), annospcmat, other)) + @test issparse(hvcat((2,), other, annospcmat)) + @test issparse(cat(annospcmat, other; dims=(1,2))) + @test issparse(cat(other, annospcmat; dims=(1,2))) + end + end + # The preceding tests should cover multi-way combinations of those types, but for good + # measure test a few multi-way combinations involving those types + @test issparse(vcat(spmat, densemat, annospcmats[1], annodmats[2])) + @test issparse(vcat(densemat, spmat, annodmats[1], annospcmats[2])) + @test issparse(hcat(spvec, annodmats[1], annospcmats[1], densevec, diagmat)) + @test issparse(hcat(annodmats[2], annospcmats[2], spvec, densevec, diagmat)) + @test issparse(hvcat((5,), diagmat, densevec, spvec, annodmats[1], annospcmats[1])) + @test issparse(hvcat((5,), spvec, annodmats[2], diagmat, densevec, annospcmats[2])) + @test issparse(cat(annodmats[1], diagmat, annospcmats[2], densevec, spvec; dims=(1,2))) + @test issparse(cat(spvec, diagmat, densevec, annospcmats[1], annodmats[2]; dims=(1,2))) + end + + @testset "hcat and vcat involving UniformScaling" begin + @test_throws ArgumentError hcat(I) + @test_throws ArgumentError [I I] + @test_throws ArgumentError vcat(I) + @test_throws ArgumentError [I; I] + @test_throws ArgumentError [I I; I] + + A = SparseMatrixCSC(rand(3,4)) + B = SparseMatrixCSC(rand(3,3)) + C = SparseMatrixCSC(rand(0,3)) + D = SparseMatrixCSC(rand(2,0)) + E = SparseMatrixCSC(rand(1,3)) + F = SparseMatrixCSC(rand(3,1)) + α = rand() + @test (hcat(A, 2I))::SparseMatrixCSC == hcat(A, Matrix(2I, 3, 3)) + @test (hcat(E, α))::SparseMatrixCSC == hcat(E, [α]) + @test (hcat(E, α, 2I))::SparseMatrixCSC == hcat(E, [α], fill(2, 1, 1)) + @test (vcat(A, 2I))::SparseMatrixCSC == vcat(A, Matrix(2I, 4, 4)) + @test (vcat(F, α))::SparseMatrixCSC == vcat(F, [α]) + @test (vcat(F, α, 2I))::SparseMatrixCSC == vcat(F, [α], fill(2, 1, 1)) + @test (hcat(C, 2I))::SparseMatrixCSC == C + @test_throws DimensionMismatch hcat(C, α) + @test (vcat(D, 2I))::SparseMatrixCSC == D + @test_throws DimensionMismatch vcat(D, α) + @test (hcat(I, 3I, A, 2I))::SparseMatrixCSC == hcat(Matrix(I, 3, 3), Matrix(3I, 3, 3), A, Matrix(2I, 3, 3)) + @test (vcat(I, 3I, A, 2I))::SparseMatrixCSC == vcat(Matrix(I, 4, 4), Matrix(3I, 4, 4), A, Matrix(2I, 4, 4)) + @test (hvcat((2,1,2), B, 2I, I, 3I, 4I))::SparseMatrixCSC == + hvcat((2,1,2), B, Matrix(2I, 3, 3), Matrix(I, 6, 6), Matrix(3I, 3, 3), Matrix(4I, 3, 3)) + @test hvcat((3,1), C, C, I, 3I)::SparseMatrixCSC == hvcat((2,1), C, C, Matrix(3I, 6,6)) + @test hvcat((2,2,2), I, 2I, 3I, 4I, C, C)::SparseMatrixCSC == + hvcat((2,2,2), Matrix(I, 3, 3), Matrix(2I, 3,3 ), Matrix(3I, 3,3), Matrix(4I, 3,3), C, C) + @test hvcat((2,2,4), C, C, I, 2I, 3I, 4I, 5I, D)::SparseMatrixCSC == + hvcat((2,2,4), C, C, Matrix(I, 3, 3), Matrix(2I,3,3), + Matrix(3I, 2, 2), Matrix(4I, 2, 2), Matrix(5I,2,2), D) + @test (hvcat((2,3,2), B, 2I, C, C, I, 3I, 4I))::SparseMatrixCSC == + hvcat((2,2,2), B, Matrix(2I, 3, 3), C, C, Matrix(3I, 3, 3), Matrix(4I, 3, 3)) + @test hvcat((3,2,1), C, C, I, B ,3I, 2I)::SparseMatrixCSC == + hvcat((2,2,1), C, C, B, Matrix(3I,3,3), Matrix(2I,6,6)) + @test (hvcat((1,2), A, E, α))::SparseMatrixCSC == hvcat((1,2), A, E, [α]) == hvcat((1,2), A, E, α*I) + @test (hvcat((2,2), α, E, F, 3I))::SparseMatrixCSC == hvcat((2,2), [α], E, F, Matrix(3I, 3, 3)) + @test (hvcat((2,2), 3I, F, E, α))::SparseMatrixCSC == hvcat((2,2), Matrix(3I, 3, 3), F, E, [α]) + end end let @@ -2023,6 +2195,31 @@ end @test_throws LinearAlgebra.SingularException UpperTriangular(A)\b end +@testset "Diagonal linear solve" begin + n = 12 + for relty in (Float32, Float64), elty in (relty, Complex{relty}) + dd=convert(Vector{elty}, randn(n)) + 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) + @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) + end + end +end + @testset "issue described in https://groups.google.com/forum/#!topic/julia-dev/QT7qpIpgOaA" begin @test sparse([1,1], [1,1], [true, true]) == sparse([1,1], [1,1], [true, true], 1, 1) == fill(true, 1, 1) @test sparsevec([1,1], [true, true]) == sparsevec([1,1], [true, true], 1) == fill(true, 1) @@ -2211,6 +2408,13 @@ end @test_throws BoundsError setindex!(A, [4.0, 5.0, 6.0], 3, 4) end +@testset "issue #29644" begin + F = lu(Tridiagonal(sparse(1.0I, 3, 3))) + @test F.L == Matrix(I, 3, 3) + @test startswith(sprint(show, MIME("text/plain"), F), + "$(LinearAlgebra.LU){Float64, $(LinearAlgebra.Tridiagonal){Float64, SparseArrays.SparseVector") +end + @testset "isstored" begin m = 5 n = 4 @@ -2466,6 +2670,32 @@ end @test similar(A, Float32, Int8, 6) == similar(A, Float32, Int8, (6,)) end +@testset "similar should preserve underlying storage type and uplo flag" begin + m, n = 4, 3 + sparsemat = sprand(m, m, 0.5) + for SymType in (Symmetric, Hermitian) + symsparsemat = SymType(sparsemat) + @test isa(similar(symsparsemat), typeof(symsparsemat)) + @test similar(symsparsemat).uplo == symsparsemat.uplo + @test isa(similar(symsparsemat, Float32), SymType{Float32,<:SparseMatrixCSC{Float32}}) + @test similar(symsparsemat, Float32).uplo == symsparsemat.uplo + @test isa(similar(symsparsemat, (n, n)), typeof(sparsemat)) + @test isa(similar(symsparsemat, Float32, (n, n)), SparseMatrixCSC{Float32}) + end +end + +@testset "similar should preserve underlying storage type" begin + local m, n = 4, 3 + sparsemat = sprand(m, m, 0.5) + for TriType in (UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular) + trisparsemat = TriType(sparsemat) + @test isa(similar(trisparsemat), typeof(trisparsemat)) + @test isa(similar(trisparsemat, Float32), TriType{Float32,<:SparseMatrixCSC{Float32}}) + @test isa(similar(trisparsemat, (n, n)), typeof(sparsemat)) + @test isa(similar(trisparsemat, Float32, (n, n)), SparseMatrixCSC{Float32}) + end +end + @testset "count specializations" begin # count should throw for sparse arrays for which zero(eltype) does not exist @test_throws MethodError count(SparseMatrixCSC(2, 2, Int[1, 2, 3], Int[1, 2], Any[true, true]))