Make LinearAlgebra.jl independent of SparseArrays.jl

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

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

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]))