Merge pull request JuliaLang#43127 from JuliaLang/dk/freelinalg

Make LinearAlgebra.jl independent of SparseArrays.jl

NEWS.md

@@ -95,11 +95,24 @@ Standard library changes
 #### Package Manager
 
 #### LinearAlgebra
-* The BLAS submodule now supports the level-2 BLAS subroutine `spr!` ([#42830]).
 
+* The BLAS submodule now supports the level-2 BLAS subroutine `spr!` ([#42830]).
 * `cholesky[!]` now supports `LinearAlgebra.PivotingStrategy` (singleton type) values
  as its optional `pivot` argument: the default is `cholesky(A, NoPivot())` (vs.
  `cholesky(A, RowMaximum())`); the former `Val{true/false}`-based calls are deprecated. ([#41640])
+* The standard library `LinearAlgebra.jl` is now completely independent of `SparseArrays.jl`,
+ both in terms of the source code as well as unit testing ([#43127]). As a consequence,
+ sparse arrays are no longer (silently) returned by methods from `LinearAlgebra` applied
+ to `Base` or `LinearAlgebra` objects. Specifically, this results in the following breaking
+ changes:
+
+ * Concatenations involving special "sparse" matrices (`*diagonal`) now return dense matrices;
+ As a consequence, the `D1` and `D2` fields of `SVD` objects, constructed upon `getproperty`
+ calls are now dense matrices.
+ * 3-arg `similar(::SpecialSparseMatrix, ::Type, ::Dims)` returns a dense zero matrix.
+ As a consequence, products of bi-, tri- and symmetric tridiagonal matrices with each
+ other result in dense output. Moreover, constructing 3-arg similar matrices of special
+ "sparse" matrices of (nonstatic) matrices now fails for the lack of `zero(::Type{Matrix{T}})`.
 
 #### Markdown
 
@@ -128,6 +141,12 @@ Standard library changes
 
 #### SparseArrays
 
+* New sparse concatenation functions `sparse_hcat`, `sparse_vcat`, and `sparse_hvcat` return
+ `SparseMatrixCSC` output independent from the types of the input arguments. They make
+ concatenation behavior available, in which the presence of some special "sparse" matrix
+ argument resulted in sparse output by multiple dispatch. This is no longer possible after
+ making `LinearAlgebra.jl` independent from `SparseArrays.jl` ([#43127]).
+
 #### Dates
 
 #### Downloads

stdlib/LinearAlgebra/Project.toml

@@ -8,7 +8,6 @@ libblastrampoline_jll = "8e850b90-86db-534c-a0d3-1478176c7d93"
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [targets]
-test = ["Test", "Random", "SparseArrays"]
+test = ["Test", "Random"]

stdlib/LinearAlgebra/docs/src/index.md

@@ -1,7 +1,7 @@
 # [Linear Algebra](@id man-linalg)
 
 ```@meta
-DocTestSetup = :(using LinearAlgebra, SparseArrays, SuiteSparse)
+DocTestSetup = :(using LinearAlgebra)
 ```
 
 In addition to (and as part of) its support for multi-dimensional arrays, Julia provides native implementations
@@ -308,7 +308,9 @@ of the Linear Algebra documentation.
 
 ## Standard functions
 
-Linear algebra functions in Julia are largely implemented by calling functions from [LAPACK](https://www.netlib.org/lapack/). Sparse matrix factorizations call functions from [SuiteSparse](https://suitesparse.com). Other sparse solvers are available as Julia packages.
+Linear algebra functions in Julia are largely implemented by calling functions from [LAPACK](https://www.netlib.org/lapack/).
+Sparse matrix factorizations call functions from [SuiteSparse](https://suitesparse.com).
+Other sparse solvers are available as Julia packages.
 
 ```@docs
 Base.:*(::AbstractMatrix, ::AbstractMatrix)

stdlib/LinearAlgebra/src/LinearAlgebra.jl

@@ -372,15 +372,14 @@ algorithm.
 
 See also: `copy_similar`, `copy_to_array`.
 """
-copy_oftype(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A,T), A)
+copy_oftype(A::AbstractArray, ::Type{T}) where {T} = copyto!(similar(A, T), A)
 
 """
  copy_similar(A, T)
 
 Copy `A` to a mutable array with eltype `T` based on `similar(A, T, size(A))`.
 
-Compared to `copy_oftype`, the result can be more flexible. For example,
-supplying a tridiagonal matrix results in a sparse array. In general, the type
+Compared to `copy_oftype`, the result can be more flexible. In general, the type
 of the output corresponds to that of the three-argument method `similar(A, T, size(s))`.
 
 See also: `copy_oftype`, `copy_to_array`.

stdlib/LinearAlgebra/src/adjtrans.jl

@@ -227,7 +227,7 @@ _adjoint_hcat(avs::Union{Number,AdjointAbsVec}...) = adjoint(vcat(map(adjoint, a
 _transpose_hcat(tvs::Union{Number,TransposeAbsVec}...) = transpose(vcat(map(transpose, tvs)...))
 typed_hcat(::Type{T}, avs::Union{Number,AdjointAbsVec}...) where {T} = adjoint(typed_vcat(T, map(adjoint, avs)...))
 typed_hcat(::Type{T}, tvs::Union{Number,TransposeAbsVec}...) where {T} = transpose(typed_vcat(T, map(transpose, tvs)...))
-# otherwise-redundant definitions necessary to prevent hitting the concat methods in sparse/sparsevector.jl
+# otherwise-redundant definitions necessary to prevent hitting the concat methods in LinearAlgebra/special.jl
 hcat(avs::Adjoint{<:Any,<:Vector}...) = _adjoint_hcat(avs...)
 hcat(tvs::Transpose{<:Any,<:Vector}...) = _transpose_hcat(tvs...)
 hcat(avs::Adjoint{T,Vector{T}}...) where {T} = _adjoint_hcat(avs...)

stdlib/LinearAlgebra/src/bidiag.jl

@@ -204,12 +204,8 @@ AbstractMatrix{T}(A::Bidiagonal) where {T} = convert(Bidiagonal{T}, A)
 
 convert(T::Type{<:Bidiagonal}, m::AbstractMatrix) = m isa T ? m : T(m)
 
-# For B<:Bidiagonal, similar(B[, neweltype]) should yield a Bidiagonal matrix.
-# On the other hand, similar(B, [neweltype,] shape...) should yield a sparse matrix.
-# The first method below effects the former, and the second the latter.
 similar(B::Bidiagonal, ::Type{T}) where {T} = Bidiagonal(similar(B.dv, T), similar(B.ev, T), B.uplo)
-# The method below is moved to SparseArrays for now
-# similar(B::Bidiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)
+similar(B::Bidiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = zeros(T, dims...)
 
 
 ###################
@@ -706,6 +702,11 @@ function *(A::SymTridiagonal, B::Diagonal)
  A_mul_B_td!(Tridiagonal(zeros(TS, size(A, 1)-1), zeros(TS, size(A, 1)), zeros(TS, size(A, 1)-1)), A, B)
 end
 
+function *(A::BiTriSym, B::BiTriSym)
+ TS = promote_op(matprod, eltype(A), eltype(B))
+ mul!(similar(A, TS, size(A)...), A, B)
+end
+
 function dot(x::AbstractVector, B::Bidiagonal, y::AbstractVector)
  require_one_based_indexing(x, y)
  nx, ny = length(x), length(y)

stdlib/LinearAlgebra/src/diagonal.jl

@@ -87,12 +87,8 @@ Construct an uninitialized `Diagonal{T}` of length `n`. See `undef`.
 """
 Diagonal{T}(::UndefInitializer, n::Integer) where T = Diagonal(Vector{T}(undef, n))
 
-# For D<:Diagonal, similar(D[, neweltype]) should yield a Diagonal matrix.
-# On the other hand, similar(D, [neweltype,] shape...) should yield a sparse matrix.
-# The first method below effects the former, and the second the latter.
 similar(D::Diagonal, ::Type{T}) where {T} = Diagonal(similar(D.diag, T))
-# The method below is moved to SparseArrays for now
-# similar(D::Diagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)
+similar(::Diagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = zeros(T, dims...)
 
 copyto!(D1::Diagonal, D2::Diagonal) = (copyto!(D1.diag, D2.diag); D1)
 
@@ -114,8 +110,8 @@ end
  end
  r
 end
-diagzero(::Diagonal{T},i,j) where {T} = zero(T)
-diagzero(D::Diagonal{<:AbstractMatrix{T}},i,j) where {T} = zeros(T, size(D.diag[i], 1), size(D.diag[j], 2))
+diagzero(::Diagonal{T}, i, j) where {T} = zero(T)
+diagzero(D::Diagonal{<:AbstractMatrix{T}}, i, j) where {T} = zeros(T, size(D.diag[i], 1), size(D.diag[j], 2))
 
 function setindex!(D::Diagonal, v, i::Int, j::Int)
  @boundscheck checkbounds(D, i, j)

stdlib/LinearAlgebra/src/special.jl

@@ -339,3 +339,27 @@ end
 
 ==(A::Bidiagonal, B::SymTridiagonal) = iszero(_evview(B)) && iszero(A.ev) && A.dv == B.dv
 ==(B::SymTridiagonal, A::Bidiagonal) = A == B
+
+# concatenation
+const _SpecialArrays = Union{Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal}
+const _Symmetric_DenseArrays{T,A<:Matrix} = Symmetric{T,A}
+const _Hermitian_DenseArrays{T,A<:Matrix} = Hermitian{T,A}
+const _Triangular_DenseArrays{T,A<:Matrix} = AbstractTriangular{T,A}
+const _Annotated_DenseArrays = Union{_SpecialArrays, _Triangular_DenseArrays, _Symmetric_DenseArrays, _Hermitian_DenseArrays}
+const _Annotated_Typed_DenseArrays{T} = Union{_Triangular_DenseArrays{T}, _Symmetric_DenseArrays{T}, _Hermitian_DenseArrays{T}}
+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}}
+
+promote_to_array_type(::Tuple{Vararg{Union{_DenseConcatGroup,UniformScaling}}}) = Matrix
+
+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...)

stdlib/LinearAlgebra/src/svd.jl

@@ -333,13 +333,13 @@ Q factor:
  1.0 0.0
  0.0 1.0
 D1 factor:
-2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 2 stored entries:
- 0.707107  ⋅
-  ⋅  0.707107
+2×2 Matrix{Float64}:
+ 0.707107 0.0
+ 0.0 0.707107
 D2 factor:
-2×2 SparseArrays.SparseMatrixCSC{Float64, Int64} with 2 stored entries:
- 0.707107  ⋅
-  ⋅  0.707107
+2×2 Matrix{Float64}:
+ 0.707107 0.0
+ 0.0 0.707107
 R0 factor:
 2×2 Matrix{Float64}:
  1.41421 0.0
@@ -352,8 +352,8 @@ julia> F.U*F.D1*F.R0*F.Q'
 
 julia> F.V*F.D2*F.R0*F.Q'
 2×2 Matrix{Float64}:
- 0.0 1.0
- 1.0 0.0
+ -0.0 1.0
+  1.0 0.0
 ```
 """
 struct GeneralizedSVD{T,S} <: Factorization{T}

stdlib/LinearAlgebra/src/tridiag.jl

@@ -148,12 +148,8 @@ function size(A::SymTridiagonal, d::Integer)
  end
 end
 
-# For S<:SymTridiagonal, similar(S[, neweltype]) should yield a SymTridiagonal matrix.
-# On the other hand, similar(S, [neweltype,] shape...) should yield a sparse matrix.
-# The first method below effects the former, and the second the latter.
 similar(S::SymTridiagonal, ::Type{T}) where {T} = SymTridiagonal(similar(S.dv, T), similar(S.ev, T))
-# The method below is moved to SparseArrays for now
-# similar(S::SymTridiagonal, ::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} = zeros(T, dims...)
 
 copyto!(dest::SymTridiagonal, src::SymTridiagonal) =
  (copyto!(dest.dv, src.dv); copyto!(dest.ev, _evview(src)); dest)
@@ -584,12 +580,8 @@ end
 Matrix(M::Tridiagonal{T}) where {T} = Matrix{T}(M)
 Array(M::Tridiagonal) = Matrix(M)
 
-# For M<:Tridiagonal, similar(M[, neweltype]) should yield a Tridiagonal matrix.
-# On the other hand, similar(M, [neweltype,] shape...) should yield a sparse matrix.
-# The first method below effects the former, and the second the latter.
 similar(M::Tridiagonal, ::Type{T}) where {T} = Tridiagonal(similar(M.dl, T), similar(M.d, T), similar(M.du, T))
-# The method below is moved to SparseArrays for now
-# similar(M::Tridiagonal, ::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} = zeros(T, dims...)
 
 # Operations on Tridiagonal matrices
 copyto!(dest::Tridiagonal, src::Tridiagonal) = (copyto!(dest.dl, src.dl); copyto!(dest.d, src.d); copyto!(dest.du, src.du); dest)

stdlib/LinearAlgebra/src/uniformscaling.jl

@@ -391,7 +391,7 @@ end
 # so that the same promotion code can be used for hvcat. We pass the type T
 # so that we can re-use this code for sparse-matrix hcat etcetera.
 promote_to_arrays_(n::Int, ::Type, a::Number) = a
-promote_to_arrays_(n::Int, ::Type{Matrix}, J::UniformScaling{T}) where {T} = copyto!(Matrix{T}(undef, n,n), J)
+promote_to_arrays_(n::Int, ::Type{Matrix}, J::UniformScaling{T}) where {T} = Matrix(J, n, n)
 promote_to_arrays_(n::Int, ::Type, A::AbstractVecOrMat) = A
 promote_to_arrays(n,k, ::Type) = ()
 promote_to_arrays(n,k, ::Type{T}, A) where {T} = (promote_to_arrays_(n[k], T, A),)
@@ -407,7 +407,7 @@ for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"
  @eval begin
  @inline $f(A::Union{AbstractVecOrMat,UniformScaling}...) = $_f(A...)
  @inline $f(A::Union{AbstractVecOrMat,UniformScaling,Number}...) = $_f(A...)
- function $_f(A::Union{AbstractVecOrMat,UniformScaling,Number}...)
+ function $_f(A::Union{AbstractVecOrMat,UniformScaling,Number}...; array_type = promote_to_array_type(A))
  n = -1
  for a in A
  if !isa(a, UniformScaling)
@@ -420,14 +420,14 @@ for (f, _f, dim, name) in ((:hcat, :_hcat, 1, "rows"), (:vcat, :_vcat, 2, "cols"
  end
  end
  n == -1 && throw(ArgumentError($("$f of only UniformScaling objects cannot determine the matrix size")))
- return cat(promote_to_arrays(fill(n, length(A)), 1, promote_to_array_type(A), A...)..., dims=Val(3-$dim))
+ return cat(promote_to_arrays(fill(n, length(A)), 1, array_type, A...)..., dims=Val(3-$dim))
  end
  end
 end
 
 hvcat(rows::Tuple{Vararg{Int}}, A::Union{AbstractVecOrMat,UniformScaling}...) = _hvcat(rows, A...)
 hvcat(rows::Tuple{Vararg{Int}}, A::Union{AbstractVecOrMat,UniformScaling,Number}...) = _hvcat(rows, A...)
-function _hvcat(rows::Tuple{Vararg{Int}}, A::Union{AbstractVecOrMat,UniformScaling,Number}...)
+function _hvcat(rows::Tuple{Vararg{Int}}, A::Union{AbstractVecOrMat,UniformScaling,Number}...; array_type = promote_to_array_type(A))
  require_one_based_indexing(A...)
  nr = length(rows)
  sum(rows) == length(A) || throw(ArgumentError("mismatch between row sizes and number of arguments"))
@@ -479,14 +479,13 @@ function _hvcat(rows::Tuple{Vararg{Int}}, A::Union{AbstractVecOrMat,UniformScali
  j += rows[i]
  end
  end
- Atyp = promote_to_array_type(A)
- Amat = promote_to_arrays(n, 1, Atyp, A...)
+ Amat = promote_to_arrays(n, 1, array_type, A...)
  # We have two methods for promote_to_array_type, one returning Matrix and
  # another one returning SparseMatrixCSC (in SparseArrays.jl). In the dense
  # case, we cannot call hvcat for the promoted UniformScalings because this
  # causes a stack overflow. In the sparse case, however, we cannot call
  # typed_hvcat because we need a sparse output.
- if Atyp == Matrix
+ if array_type == Matrix
  return typed_hvcat(promote_eltype(Amat...), rows, Amat...)
  else
  return hvcat(rows, Amat...)

stdlib/LinearAlgebra/test/addmul.jl

@@ -6,7 +6,6 @@ using Base: rtoldefault
 using Test
 using LinearAlgebra
 using LinearAlgebra: AbstractTriangular
-using SparseArrays
 using Random
 
 _rand(::Type{T}) where {T <: AbstractFloat} = T(randn())

stdlib/LinearAlgebra/test/adjtrans.jl

@@ -2,7 +2,7 @@
 
 module TestAdjointTranspose
 
-using Test, LinearAlgebra, SparseArrays
+using Test, LinearAlgebra
 
 const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
 
@@ -354,14 +354,6 @@ end
  @test broadcast(+, Transpose(vec), 1, Transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == tvec + tvec .+ 1
  @test broadcast(+, Adjoint(vec), 1im, Adjoint(vec))::Adjoint{Complex{Int},Vector{Complex{Int}}} == avec + avec .+ 1im
  @test broadcast(+, Transpose(vec), 1im, Transpose(vec))::Transpose{Complex{Int},Vector{Complex{Int}}} == tvec + tvec .+ 1im
- # 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}}
 end
 
 @testset "Adjoint/Transpose-wrapped vector multiplication" begin

stdlib/LinearAlgebra/test/ambiguous_exec.jl

@@ -1,6 +1,6 @@
 # This file is a part of Julia. License is MIT: https://julialang.org/license
 
-using Test, LinearAlgebra, SparseArrays
+using Test, LinearAlgebra
 let ambig = detect_ambiguities(LinearAlgebra; recursive=true)
  @test isempty(ambig)
  ambig = Set{Any}(((m1.sig, m2.sig) for (m1, m2) in ambig))

stdlib/LinearAlgebra/test/bidiag.jl

@@ -2,7 +2,7 @@
 
 module TestBidiagonal
 
-using Test, LinearAlgebra, SparseArrays, Random
+using Test, LinearAlgebra, Random
 using LinearAlgebra: BlasReal, BlasFloat
 
 const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
@@ -98,8 +98,8 @@ Random.seed!(1)
  @test similar(ubd).uplo == ubd.uplo
  @test isa(similar(ubd, Int), Bidiagonal{Int})
  @test similar(ubd, Int).uplo == ubd.uplo
- @test isa(similar(ubd, (3, 2)), SparseMatrixCSC)
- @test isa(similar(ubd, Int, (3, 2)), SparseMatrixCSC{Int})
+ @test isa(similar(ubd, (3, 2)), Matrix)
+ @test isa(similar(ubd, Int, (3, 2)), Matrix{Int})
 
  # setindex! when off diagonal is zero bug
  Bu = Bidiagonal(rand(elty, 10), zeros(elty, 9), 'U')
@@ -432,9 +432,7 @@ using LinearAlgebra: fillstored!, UnitLowerTriangular
  exotic_arrays = Any[Tridiagonal(randn(3), randn(4), randn(3)),
  Bidiagonal(randn(3), randn(2), rand([:U,:L])),
  SymTridiagonal(randn(3), randn(2)),
- sparse(randn(3,4)),
  Diagonal(randn(5)),
- sparse(rand(3)),
  # LowerTriangular(randn(3,3)), # AbstractTriangular fill! deprecated, see below
  # UpperTriangular(randn(3,3)) # AbstractTriangular fill! deprecated, see below
  ]