From 7073d71458e8d88c01908c951a17361c395e0985 Mon Sep 17 00:00:00 2001 From: Klaus Crusius Date: Mon, 19 Nov 2018 14:06:37 +0100 Subject: [PATCH 1/5] added sprandn methods with Type --- stdlib/SparseArrays/src/sparsematrix.jl | 4 +++- stdlib/SparseArrays/src/sparsevector.jl | 2 ++ stdlib/SparseArrays/test/sparse.jl | 7 +++++++ stdlib/SparseArrays/test/sparsevector.jl | 5 +++++ 4 files changed, 17 insertions(+), 1 deletion(-) diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl index 02c975168c6ec..7fac079ac97d6 100644 --- a/stdlib/SparseArrays/src/sparsematrix.jl +++ b/stdlib/SparseArrays/src/sparsematrix.jl @@ -1471,7 +1471,7 @@ sprand(r::AbstractRNG, ::Type{Bool}, m::Integer, n::Integer, density::AbstractFl sprand(::Type{T}, m::Integer, n::Integer, density::AbstractFloat) where {T} = sprand(GLOBAL_RNG, T, m, n, density) """ - sprandn([rng], m[,n],p::AbstractFloat) + sprandn([rng],[,Type],m[,n],p::AbstractFloat) Create a random sparse vector of length `m` or sparse matrix of size `m` by `n` with the specified (independent) probability `p` of any entry being nonzero, @@ -1488,6 +1488,8 @@ julia> sprandn(2, 2, 0.75) """ sprandn(r::AbstractRNG, m::Integer, n::Integer, density::AbstractFloat) = sprand(r,m,n,density,randn,Float64) sprandn(m::Integer, n::Integer, density::AbstractFloat) = sprandn(GLOBAL_RNG,m,n,density) +sprandn(r::AbstractRNG, ::Type{T}, m::Integer, n::Integer, density::AbstractFloat) where T = sprand(r,m,n,density,(r,i) -> randn(r,T,i), T) +sprandn(::Type{T}, m::Integer, n::Integer, density::AbstractFloat) where T = sprandn(GLOBAL_RNG,T,m,n,density) LinearAlgebra.fillstored!(S::SparseMatrixCSC, x) = (fill!(nzvalview(S), x); S) diff --git a/stdlib/SparseArrays/src/sparsevector.jl b/stdlib/SparseArrays/src/sparsevector.jl index e7f0e1a2db26a..ce9fbbc35e947 100644 --- a/stdlib/SparseArrays/src/sparsevector.jl +++ b/stdlib/SparseArrays/src/sparsevector.jl @@ -506,6 +506,8 @@ sprand(::Type{T}, n::Integer, p::AbstractFloat) where {T} = sprand(GLOBAL_RNG, T sprandn(n::Integer, p::AbstractFloat) = sprand(GLOBAL_RNG, n, p, randn) sprandn(r::AbstractRNG, n::Integer, p::AbstractFloat) = sprand(r, n, p, randn) +sprandn(::Type{T}, n::Integer, p::AbstractFloat) where T = sprand(GLOBAL_RNG, n, p, (r, i) -> randn(r, T, i)) +sprandn(r::AbstractRNG, ::Type{T}, n::Integer, p::AbstractFloat) where T = sprand(r, n, p, (r, i) -> randn(r, T, i)) ## Indexing into Matrices can return SparseVectors diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index c8e5ebd8b5acb..9eb2069d3984c 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/stdlib/SparseArrays/test/sparse.jl @@ -2335,4 +2335,11 @@ end @test m2.module == SparseArrays end +@testset "sprandn with type $T" for T in (Float64, Float32, Float16, ComplexF64, ComplexF32, ComplexF16) + @test sprandn(T, 5, 5, 0.5) isa AbstractSparseMatrix{T} +end +@testset "sprandn with invalid type $T" for T in (AbstractFloat, BigFloat, Complex) + @test_throws MethodError sprandn(T, 5, 5, 0.5) +end + end # module diff --git a/stdlib/SparseArrays/test/sparsevector.jl b/stdlib/SparseArrays/test/sparsevector.jl index 6061b21a6721c..d8ada3963e2c5 100644 --- a/stdlib/SparseArrays/test/sparsevector.jl +++ b/stdlib/SparseArrays/test/sparsevector.jl @@ -154,6 +154,11 @@ end end end + let xr = sprandn(ComplexF64, 1000, 0.9) + @test isa(xr, SparseVector{ComplexF64,Int}) + @test length(xr) == 1000 + end + let xr = sprand(Bool, 1000, 0.9) @test isa(xr, SparseVector{Bool,Int}) @test length(xr) == 1000 From 18a88c5694da91f4dd1f515653ab9a9ebcfb2b45 Mon Sep 17 00:00:00 2001 From: Klaus Crusius Date: Mon, 31 Dec 2018 20:03:46 +0100 Subject: [PATCH 2/5] tests for wrapped sparse --- stdlib/SparseArrays/src/sparsematrix.jl | 315 ++++++++++++++++++++++++ stdlib/SparseArrays/test/sparse.jl | 8 + 2 files changed, 323 insertions(+) diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl index 7324b368b05a2..52a19ec1b0f1a 100644 --- a/stdlib/SparseArrays/src/sparsematrix.jl +++ b/stdlib/SparseArrays/src/sparsematrix.jl @@ -412,6 +412,321 @@ SparseMatrixCSC{Tv}(M::Transpose{<:Any,SparseMatrixCSC}) where {Tv} = SparseMatr SparseMatrixCSC{Tv,Ti}(M::Adjoint{<:Any,SparseMatrixCSC}) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(copy(M)) SparseMatrixCSC{Tv,Ti}(M::Transpose{<:Any,SparseMatrixCSC}) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(copy(M)) +import LinearAlgebra: AbstractTriangular + +""" + SparseMatrixCSCSymmHerm + +`Symmetric` or `Hermitian` of a `SparseMatrixCSC` or `SparseMatrixCSCView`. +""" +const SparseMatrixCSCSymmHerm{Tv,Ti} = Union{Symmetric{Tv,<:SparseMatrixCSCUnion{Tv,Ti}}, + Hermitian{Tv,<:SparseMatrixCSCUnion{Tv,Ti}}} + +const AbstractTriangularSparse{Tv,Ti} = AbstractTriangular{Tv,<:SparseMatrixCSCUnion{Tv,Ti}} + +# converting Symmetric/Hermitian/AbstractTriangular/SubArray of SparseMatrixCSC +# and Transpose/Adjoint of AbstractTriangular of SparseMatrixCSC to SparseMatrixCSC +for wr in (Symmetric, Hermitian, Transpose, Adjoint, + UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular, + SubArray) + + @eval SparseMatrixCSC(A::$wr) = _sparsem(A) + @eval SparseMatrixCSC{Tv}(A::$wr{Tv}) where Tv = _sparsem(A) + @eval SparseMatrixCSC{Tv}(A::$wr) where Tv = SparseMatrixCSC{Tv}(_sparsem(A)) + @eval SparseMatrixCSC{Tv,Ti}(A::$wr) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(_sparsem(A)) +end + +""" + iswrsparse(::S) + iswrsparse(::Type{S}) + +Returns `true` if type `S` is backed by a sparse array, and `false` otherwise. +""" +iswrsparse(::T) where T<:AbstractArray = iswrsparse(T) +iswrsparse(::Type) = false +iswrsparse(::Type{T}) where T<:AbstractSparseArray = true + +""" + depth(::Type{S}) + +Returns 0 for unwrapped S, and nesting depth for wrapped (nested) abstract arrays. +""" +depth(::T) where T = depth(T) +depth(::Type{T}) where T<:AbstractArray = 0 + +for wr in (Symmetric, Hermitian, + LowerTriangular, UnitLowerTriangular, UpperTriangular, UnitUpperTriangular, + Transpose, Adjoint, SubArray, + Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal) + + pl = wr === SubArray ? :($wr{<:Any,<:Any,T}) : :($wr{<:Any,T}) + @eval iswrsparse(::Type{<:$pl}) where T = iswrsparse(T) + @eval depth(::Type{<:$pl}) where T = depth(T) + 1 +end + +# convert parent and re-wrap in same wrapper +_sparsewrap(A::Symmetric) = Symmetric(_sparsem(parent(A)), A.uplo == 'U' ? :U : :L) +_sparsewrap(A::Hermitian) = Hermitian(_sparsem(parent(A)), A.uplo == 'U' ? :U : :L) +_sparsewrap(A::SubArray) = SubArray(_sparsem(parent(A)), A.indices) +for ty in ( LowerTriangular, UnitLowerTriangular, + UpperTriangular, UnitUpperTriangular, + Transpose, Adjoint) + + @eval _sparsewrap(A::$ty) = $ty(_sparsem(parent(A))) +end +function _sparsewrap(A::Union{Diagonal,Bidiagonal,Tridiagonal,SymTridiagonal}) + dropzeros!(sparse(A)) +end + +""" + unwrap(A::AbstractMatrix) + +In case A is a wrapper type (`SubArray, Symmetric, Adjoint, SubArray, Triangular, Tridiagonal`, etc.) +convert to `Matrix` or `SparseMatrixCSC`, depending on final storage type of A. +For other types return A itself. +""" +unwrap(A::Any) = A +unwrap(A::AbstractMatrix) = iswrsparse(A) ? convert(SparseMatrixCSC, A) : convert(Array, A) + +import Base.copy +copy(A::SubArray) = getindex(unwrap(parent(A)), A.indices...) + +# For pure sparse matrices and vectors return A. +# For wrapped sparse matrices or vectors convert to SparseMatrixCSC. +# Handle nested wrappers properly. +# Use abstract matrix fallback if A is not sparse. +function _sparsem(@nospecialize A::AbstractArray{Tv}) where Tv + if iswrsparse(A) + if depth(A) >= 1 + _sparsem(_sparsewrap(A)) + else + A + end + else + # explicitly call abstract matrix fallback using getindex(A,...) + invoke(SparseMatrixCSC{Tv,Int}, Tuple{AbstractMatrix}, A) + end +end + +_sparsem(A::AbstractSparseMatrix) = A +_sparsem(A::AbstractSparseVector) = A + +# Transpose/Adjoint of sparse vector (returning sparse matrix) +function _sparsem(A::Union{Transpose{<:Any,<:AbstractSparseVector},Adjoint{<:Any,<:AbstractSparseVector}}) + B = parent(A) + n = length(B) + Ti = eltype(B.nzind) + fadj = A isa Transpose ? transpose : adjoint + colptr = fill!(Vector{Ti}(undef, n + 1), 0) + colptr[1] = 1 + colptr[B.nzind .+ 1] .= 1 + cumsum!(colptr, colptr) + rowval = fill!(similar(B.nzind), 1) + nzval = fadj.(nonzeros(B)) + SparseMatrixCSC(1, n, colptr, rowval, nzval) +end + +function _sparsem(A::Union{Transpose{<:Any,<:SparseMatrixCSC},Adjoint{<:Any,<:SparseMatrixCSC}}) + ftranspose(parent(A), A isa Transpose ? transpose : adjoint) +end + +# Symmetric/Hermitian of sparse matrix +_sparsem(A::SparseMatrixCSCSymmHerm) = _sparsem(A.uplo == 'U' ? nzrangeup : nzrangelo, A) +# Triangular of sparse matrix +_sparsem(A::UpperTriangular{T,<:AbstractSparseMatrix}) where T = triu(A.data) +_sparsem(A::LowerTriangular{T,<:AbstractSparseMatrix}) where T = tril(A.data) +# view of sparse matrix +_sparsem(S::SubArray{<:Any,2,<:SparseMatrixCSC}) = getindex(S.parent,S.indices...) + +# 4 cases: (Symmetric|Hermitian) variants (:U|:L) +function _sparsem(fnzrange::Function, sA::SparseMatrixCSCSymmHerm{Tv}) where {Tv} + A = sA.data + rowval = rowvals(A) + nzval = nonzeros(A) + m, n = size(A) + Ti = eltype(rowval) + fadj = sA isa Symmetric ? transpose : adjoint + newcolptr = Vector{Ti}(undef, n+1) + diagmap = fadj == transpose ? identity : real + + newcolptr[1] = 1 + colrange = fnzrange === nzrangeup ? (1:n) : (n:-1:1) + @inbounds for j = colrange + r = fnzrange(A, j); r1 = r.start; r2 = r.stop + newcolptr[j+1] = r2 - r1 + 1 + for k = r1:r2 + row = rowval[k] + if row != j + newcolptr[row+1] += 1 + end + end + end + cumsum!(newcolptr, newcolptr) + nz = newcolptr[n+1] - 1 + newrowval = Vector{Ti}(undef, nz) + newnzval = Vector{Tv}(undef, nz) + @inbounds for j = 1:n + newk = newcolptr[j] + for k = fnzrange(A, j) + i = rowval[k] + nzv = nzval[k] + if i != j + newrowval[newk] = i + newnzval[newk] = nzv + newk += 1 + ni = newcolptr[i] + newrowval[ni] = j + newnzval[ni] = fadj(nzv) + newcolptr[i] = ni + 1 + else + newrowval[newk] = i + newnzval[newk] = diagmap(nzv) + newk += 1 + end + end + newcolptr[j] = newk + end + _sparse_gen(m, n, newcolptr, newrowval, newnzval) +end + +# 2 cases: Unit(Upper|Lower)Triangular{Tv,SparseMatrixCSC} +function _sparsem(A::AbstractTriangularSparse{Tv}) where Tv + S = A.data + rowval = rowvals(S) + nzval = nonzeros(S) + m, n = size(S) + Ti = eltype(rowval) + fnzrange = A isa Union{UpperTriangular,UnitUpperTriangular} ? nzrangeup : nzrangelo + unit = A isa Union{UnitUpperTriangular,UnitLowerTriangular} + newcolptr = Vector{Ti}(undef, n+1) + newrowval = Vector{Ti}(undef, nnz(S)) + newnzval = Vector{Tv}(undef, nnz(S)) + newcolptr[1] = 1 + uplo = fnzrange == nzrangeup + newk = 1 + @inbounds for j = 1:n + newkk = newk + if unit + newk += !uplo + end + r = fnzrange(S, j); r1 = r.start; r2 = r.stop + for k = r1:r2 + i = rowval[k] + if i != j || i == j && !unit + newrowval[newk] = i + newnzval[newk] = nzval[k] + newk += 1 + end + end + if unit + uplo && (newkk = newk) + newrowval[newkk] = j + newnzval[newkk] = one(Tv) + newk += uplo + end + newcolptr[j+1] = newk + end + nz = newcolptr[n+1] - 1 + resize!(newrowval, nz) + resize!(newnzval, nz) + SparseMatrixCSC(m, n, newcolptr, newrowval, newnzval) +end + +# 8 cases: (Transpose|Adjoint){Tv,[Unit](Upper|Lower)Triangular} +function _sparsem(taA::Union{Transpose{Tv,<:AbstractTriangularSparse}, + Adjoint{Tv,<:AbstractTriangularSparse}}) where {Tv} + + sA = taA.parent + A = sA.data + rowval = rowvals(A) + nzval = nonzeros(A) + m, n = size(A) + Ti = eltype(rowval) + fnzrange = sA isa Union{UpperTriangular,UnitUpperTriangular} ? nzrangeup : nzrangelo + fadj = taA isa Transpose ? transpose : adjoint + unit = sA isa Union{UnitUpperTriangular,UnitLowerTriangular} + uplo = A isa Union{UpperTriangular,UnitUpperTriangular} + + newcolptr = Vector{Ti}(undef, n+1) + fill!(newcolptr, 1unit) + newcolptr[1] = 1 + @inbounds for j = 1:n + for k = fnzrange(A, j) + i = rowval[k] + if i != j || i == j && !unit + newcolptr[i+1] += 1 + end + end + end + cumsum!(newcolptr, newcolptr) + nz = newcolptr[n+1] - 1 + newrowval = Vector{Ti}(undef, nz) + newnzval = Vector{Tv}(undef, nz) + + @inbounds for j = 1:n + if !uplo && unit + ni = newcolptr[j] + newrowval[ni] = j + newnzval[ni] = fadj(one(Tv)) + newcolptr[j] = ni + 1 + end + for k = fnzrange(A, j) + i = rowval[k] + nzv = nzval[k] + if i != j || i == j && !unit + ni = newcolptr[i] + newrowval[ni] = j + newnzval[ni] = fadj(nzv) + newcolptr[i] = ni + 1 + end + end + if uplo && unit + ni = newcolptr[j] + newrowval[ni] = j + newnzval[ni] = fadj(one(Tv)) + newcolptr[j] = ni + 1 + end + end + _sparse_gen(n, m, newcolptr, newrowval, newnzval) +end + +function _sparse_gen(m, n, newcolptr, newrowval, newnzval) + @inbounds for j = n:-1:1 + newcolptr[j+1] = newcolptr[j] + end + newcolptr[1] = 1 + SparseMatrixCSC(m, n, newcolptr, newrowval, newnzval) +end + +""" + nzrangeup(A::AbstractMatrix, i::Integer) +restrict `nzrange(A, i)` at diagonal as end (if only upper part matters) +""" +nzrangeup(A, i) = nzrangeup(A, nzrange(A, i), i) +""" + nzrangelo(A::AbstractMatrix, i::Integer) +restrict `nzrange(A, i)` at diagonal as start (if only lower part matters) +""" +nzrangelo(A, i) = nzrangelo(A, nzrange(A, i), i) +function nzrangeup(A, r::AbstractUnitRange, i) + r1 = r.start; r2 = r.stop + rv = rowvals(A) + r1:searchsortedlast(rowvals(A), i, r1, r2, Forward) + # @inbounds r2 < r1 || rv[r2] <= i ? r : r1:searchsortedlast(rv, i, r1, r2, Forward) +end +function nzrangeup(A, r::AbstractVector{<:Integer}, i) + view(r, 1:searchsortedlast(view(rowvals(A), r), i)) +end +function nzrangelo(A, r::AbstractUnitRange, i) + r1 = r.start; r2 = r.stop + rv = rowvals(A) + # searchsortedfirst(rv, i, r1, r2, Forward):r2 + @inbounds r2 < r1 || rv[r1] >= i ? r : searchsortedfirst(rv, i, r1, r2, Forward):r2 +end +function nzrangelo(A, r::AbstractVector{<:Integer}, i) + view(r, searchsortedfirst(view(rowvals(A), r), i):length(r)) +end + # converting from SparseMatrixCSC to other matrix types function Matrix(S::SparseMatrixCSC{Tv}) where Tv # Handle cases where zero(Tv) is not defined but the array is dense. diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index 7784c67baf078..ee6f321d04025 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/stdlib/SparseArrays/test/sparse.jl @@ -2446,4 +2446,12 @@ end end end +@testset "wrappers of sparse" begin + m, n = 10, 10 + A = sprand(m, n, 0.5) + B = Matrix(A) + @test sparse(UpperTriangular(A')) == UpperTriangular(B') + @test sparse(Adjoint(UpperTriangular(A'))) == Adjoint(UpperTriangular(B')) +end + end # module From 5a1fd7142d512e2b849b8096870cb2a888622f78 Mon Sep 17 00:00:00 2001 From: Klaus Crusius Date: Tue, 1 Jan 2019 19:12:54 +0100 Subject: [PATCH 3/5] Fast conversion of wrapped types to SparseMatrixCSC --- stdlib/SparseArrays/test/sparse.jl | 32 ++++++++++++++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index ee6f321d04025..5c979996c3772 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/stdlib/SparseArrays/test/sparse.jl @@ -2447,9 +2447,37 @@ end end @testset "wrappers of sparse" begin - m, n = 10, 10 - A = sprand(m, n, 0.5) + m = n = 10 + A = spzeros(ComplexF64, m, n) + A[:,1] = 1:m + A[:,2] = [1 3 0 0 0 0 0 0 0 0]' + A[:,3] = [2 4 0 0 0 0 0 0 0 0]' + A[:,4] = [0 0 0 0 5 3 0 0 0 0]' + A[:,5] = [0 0 0 0 6 2 0 0 0 0]' + A[:,6] = [0 0 0 0 7 4 0 0 0 0]' + A[:,7:n] = rand(ComplexF64, m, n-6) B = Matrix(A) + dowrap(wr, A) = wr(A) + dowrap(wr::Tuple, A) = (wr[1])(A, wr[2:end]...) + + @testset "sparse($wr(A))" for wr in ( + Symmetric, (Symmetric, :L), Hermitian, (Hermitian, :L), + Transpose, Adjoint, + UpperTriangular, LowerTriangular, + UnitUpperTriangular, UnitLowerTriangular, + (view, 3:6, 2:5)) + + @test SparseMatrixCSC(dowrap(wr, A)) == Matrix(dowrap(wr, B)) + end + + @testset "sparse($at($wr))" for at = (Transpose, Adjoint), wr = + (UpperTriangular, LowerTriangular, + UnitUpperTriangular, UnitLowerTriangular) + + @test SparseMatrixCSC(at(wr(A))) == Matrix(at(wr(B))) + end + + @test sparse([1,2,3,4,5]') == SparseMatrixCSC([1 2 3 4 5]) @test sparse(UpperTriangular(A')) == UpperTriangular(B') @test sparse(Adjoint(UpperTriangular(A'))) == Adjoint(UpperTriangular(B')) end From 1244d885da358ca87f2d538135c816cfbb6e12f0 Mon Sep 17 00:00:00 2001 From: Klaus Crusius Date: Thu, 3 Jan 2019 11:18:12 +0100 Subject: [PATCH 4/5] uncommented small optimization --- stdlib/SparseArrays/src/sparsematrix.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl index 52a19ec1b0f1a..bad939dc0ce1f 100644 --- a/stdlib/SparseArrays/src/sparsematrix.jl +++ b/stdlib/SparseArrays/src/sparsematrix.jl @@ -711,8 +711,7 @@ nzrangelo(A, i) = nzrangelo(A, nzrange(A, i), i) function nzrangeup(A, r::AbstractUnitRange, i) r1 = r.start; r2 = r.stop rv = rowvals(A) - r1:searchsortedlast(rowvals(A), i, r1, r2, Forward) - # @inbounds r2 < r1 || rv[r2] <= i ? r : r1:searchsortedlast(rv, i, r1, r2, Forward) + @inbounds r2 < r1 || rv[r2] <= i ? r : r1:searchsortedlast(rv, i, r1, r2, Forward) end function nzrangeup(A, r::AbstractVector{<:Integer}, i) view(r, 1:searchsortedlast(view(rowvals(A), r), i)) @@ -720,7 +719,6 @@ end function nzrangelo(A, r::AbstractUnitRange, i) r1 = r.start; r2 = r.stop rv = rowvals(A) - # searchsortedfirst(rv, i, r1, r2, Forward):r2 @inbounds r2 < r1 || rv[r1] >= i ? r : searchsortedfirst(rv, i, r1, r2, Forward):r2 end function nzrangelo(A, r::AbstractVector{<:Integer}, i) From 86ae00e9f6eafbd7c3496d7a565fa1ca479849f5 Mon Sep 17 00:00:00 2001 From: Klaus Crusius Date: Wed, 27 Feb 2019 11:02:56 +0100 Subject: [PATCH 5/5] moved sparse conversions to sparseconvert.jl --- stdlib/SparseArrays/src/SparseArrays.jl | 3 + stdlib/SparseArrays/src/sparseconvert.jl | 288 +++++++++++++++++++++ stdlib/SparseArrays/src/sparsematrix.jl | 313 ----------------------- 3 files changed, 291 insertions(+), 313 deletions(-) create mode 100644 stdlib/SparseArrays/src/sparseconvert.jl diff --git a/stdlib/SparseArrays/src/SparseArrays.jl b/stdlib/SparseArrays/src/SparseArrays.jl index f43ff255e003b..55f112048ee05 100644 --- a/stdlib/SparseArrays/src/SparseArrays.jl +++ b/stdlib/SparseArrays/src/SparseArrays.jl @@ -36,8 +36,11 @@ export AbstractSparseArray, AbstractSparseMatrix, AbstractSparseVector, issparse, nonzeros, nzrange, rowvals, sparse, sparsevec, spdiagm, sprand, sprandn, spzeros, nnz, permute, findnz +export unwrap, iswrsparse + include("abstractsparse.jl") include("sparsematrix.jl") +include("sparseconvert.jl") include("sparsevector.jl") include("higherorderfns.jl") include("linalg.jl") diff --git a/stdlib/SparseArrays/src/sparseconvert.jl b/stdlib/SparseArrays/src/sparseconvert.jl new file mode 100644 index 0000000000000..bdadc3a4729ff --- /dev/null +++ b/stdlib/SparseArrays/src/sparseconvert.jl @@ -0,0 +1,288 @@ +# This file is a part of Julia. License is MIT: https://julialang.org/license + +import LinearAlgebra: AbstractTriangular + +""" + SparseMatrixCSCSymmHerm + +`Symmetric` or `Hermitian` of a `SparseMatrixCSC` or `SparseMatrixCSCView`. +""" +const SparseMatrixCSCSymmHerm{Tv,Ti} = Union{Symmetric{Tv,<:SparseMatrixCSCUnion{Tv,Ti}}, + Hermitian{Tv,<:SparseMatrixCSCUnion{Tv,Ti}}} + +const AbstractTriangularSparse{Tv,Ti} = AbstractTriangular{Tv,<:SparseMatrixCSCUnion{Tv,Ti}} + +# converting Symmetric/Hermitian/AbstractTriangular/SubArray of SparseMatrixCSC +# and Transpose/Adjoint of AbstractTriangular of SparseMatrixCSC to SparseMatrixCSC +for wr in (Symmetric, Hermitian, Transpose, Adjoint, + UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular, + SubArray) + + @eval SparseMatrixCSC(A::$wr) = _sparsem(A) + @eval SparseMatrixCSC{Tv}(A::$wr{Tv}) where Tv = _sparsem(A) + @eval SparseMatrixCSC{Tv}(A::$wr) where Tv = SparseMatrixCSC{Tv}(_sparsem(A)) + @eval SparseMatrixCSC{Tv,Ti}(A::$wr) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(_sparsem(A)) +end + +""" + iswrsparse(::S) + iswrsparse(::Type{S}) + +Returns `true` if type `S` is backed by a sparse array, and `false` otherwise. +""" +iswrsparse(::T) where T<:AbstractArray = iswrsparse(T) +iswrsparse(::Type) = false +iswrsparse(::Type{T}) where T<:AbstractSparseArray = true + +""" + depth(::Type{S}) + +Returns 0 for unwrapped S, and nesting depth for wrapped (nested) abstract arrays. +""" +depth(::T) where T = depth(T) +depth(::Type{T}) where T<:AbstractArray = 0 + +for wr in (Symmetric, Hermitian, + LowerTriangular, UnitLowerTriangular, UpperTriangular, UnitUpperTriangular, + Transpose, Adjoint, SubArray, + Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal) + + pl = wr === SubArray ? :($wr{<:Any,<:Any,T}) : :($wr{<:Any,T}) + @eval iswrsparse(::Type{<:$pl}) where T = iswrsparse(T) + @eval depth(::Type{<:$pl}) where T = depth(T) + 1 +end + +# convert parent and re-wrap in same wrapper +_sparsewrap(A::Symmetric) = Symmetric(_sparsem(parent(A)), A.uplo == 'U' ? :U : :L) +_sparsewrap(A::Hermitian) = Hermitian(_sparsem(parent(A)), A.uplo == 'U' ? :U : :L) +_sparsewrap(A::SubArray) = SubArray(_sparsem(parent(A)), A.indices) +for ty in ( LowerTriangular, UnitLowerTriangular, + UpperTriangular, UnitUpperTriangular, + Transpose, Adjoint) + + @eval _sparsewrap(A::$ty) = $ty(_sparsem(parent(A))) +end +function _sparsewrap(A::Union{Diagonal,Bidiagonal,Tridiagonal,SymTridiagonal}) + dropzeros!(sparse(A)) +end + +""" + unwrap(A::AbstractMatrix) + +In case A is a wrapper type (`SubArray, Symmetric, Adjoint, SubArray, Triangular, Tridiagonal`, etc.) +convert to `Matrix` or `SparseMatrixCSC`, depending on final storage type of A. +For other types return A itself. +""" +unwrap(A::Any) = A +unwrap(A::AbstractMatrix) = iswrsparse(A) ? convert(SparseMatrixCSC, A) : convert(Array, A) + +import Base.copy +copy(A::SubArray{T,2}) where T = getindex(unwrap(parent(A)), A.indices...) + +# For pure sparse matrices and vectors return A. +# For wrapped sparse matrices or vectors convert to SparseMatrixCSC. +# Handle nested wrappers properly. +# Use abstract matrix fallback if A is not sparse. +function _sparsem(@nospecialize A::AbstractArray{Tv}) where Tv + if iswrsparse(A) + if depth(A) >= 1 + _sparsem(_sparsewrap(A)) + else + A + end + else + # explicitly call abstract matrix fallback using getindex(A,...) + invoke(SparseMatrixCSC{Tv,Int}, Tuple{AbstractMatrix}, A) + end +end + +_sparsem(A::AbstractSparseMatrix) = A +_sparsem(A::AbstractSparseVector) = A + +# Transpose/Adjoint of sparse vector (returning sparse matrix) +function _sparsem(A::Union{Transpose{<:Any,<:AbstractSparseVector},Adjoint{<:Any,<:AbstractSparseVector}}) + B = parent(A) + n = length(B) + Ti = eltype(B.nzind) + fadj = A isa Transpose ? transpose : adjoint + colptr = fill!(Vector{Ti}(undef, n + 1), 0) + colptr[1] = 1 + colptr[B.nzind .+ 1] .= 1 + cumsum!(colptr, colptr) + rowval = fill!(similar(B.nzind), 1) + nzval = fadj.(nonzeros(B)) + SparseMatrixCSC(1, n, colptr, rowval, nzval) +end + +function _sparsem(A::Union{Transpose{<:Any,<:SparseMatrixCSC},Adjoint{<:Any,<:SparseMatrixCSC}}) + ftranspose(parent(A), A isa Transpose ? transpose : adjoint) +end + +# Symmetric/Hermitian of sparse matrix +_sparsem(A::SparseMatrixCSCSymmHerm) = _sparsem(A.uplo == 'U' ? nzrangeup : nzrangelo, A) +# Triangular of sparse matrix +_sparsem(A::UpperTriangular{T,<:AbstractSparseMatrix}) where T = triu(A.data) +_sparsem(A::LowerTriangular{T,<:AbstractSparseMatrix}) where T = tril(A.data) +# view of sparse matrix +_sparsem(S::SubArray{<:Any,2,<:SparseMatrixCSC}) = getindex(S.parent,S.indices...) + +# 4 cases: (Symmetric|Hermitian) variants (:U|:L) +function _sparsem(fnzrange::Function, sA::SparseMatrixCSCSymmHerm{Tv}) where {Tv} + A = sA.data + rowval = rowvals(A) + nzval = nonzeros(A) + m, n = size(A) + Ti = eltype(rowval) + fadj = sA isa Symmetric ? transpose : adjoint + newcolptr = Vector{Ti}(undef, n+1) + diagmap = fadj == transpose ? identity : real + + newcolptr[1] = 1 + colrange = fnzrange === nzrangeup ? (1:n) : (n:-1:1) + @inbounds for j = colrange + r = fnzrange(A, j); r1 = r.start; r2 = r.stop + newcolptr[j+1] = r2 - r1 + 1 + for k = r1:r2 + row = rowval[k] + if row != j + newcolptr[row+1] += 1 + end + end + end + cumsum!(newcolptr, newcolptr) + nz = newcolptr[n+1] - 1 + newrowval = Vector{Ti}(undef, nz) + newnzval = Vector{Tv}(undef, nz) + @inbounds for j = 1:n + newk = newcolptr[j] + for k = fnzrange(A, j) + i = rowval[k] + nzv = nzval[k] + if i != j + newrowval[newk] = i + newnzval[newk] = nzv + newk += 1 + ni = newcolptr[i] + newrowval[ni] = j + newnzval[ni] = fadj(nzv) + newcolptr[i] = ni + 1 + else + newrowval[newk] = i + newnzval[newk] = diagmap(nzv) + newk += 1 + end + end + newcolptr[j] = newk + end + _sparse_gen(m, n, newcolptr, newrowval, newnzval) +end + +# 2 cases: Unit(Upper|Lower)Triangular{Tv,SparseMatrixCSC} +function _sparsem(A::AbstractTriangularSparse{Tv}) where Tv + S = A.data + rowval = rowvals(S) + nzval = nonzeros(S) + m, n = size(S) + Ti = eltype(rowval) + fnzrange = A isa Union{UpperTriangular,UnitUpperTriangular} ? nzrangeup : nzrangelo + unit = A isa Union{UnitUpperTriangular,UnitLowerTriangular} + newcolptr = Vector{Ti}(undef, n+1) + newrowval = Vector{Ti}(undef, nnz(S)) + newnzval = Vector{Tv}(undef, nnz(S)) + newcolptr[1] = 1 + uplo = fnzrange == nzrangeup + newk = 1 + @inbounds for j = 1:n + newkk = newk + if unit + newk += !uplo + end + r = fnzrange(S, j); r1 = r.start; r2 = r.stop + for k = r1:r2 + i = rowval[k] + if i != j || i == j && !unit + newrowval[newk] = i + newnzval[newk] = nzval[k] + newk += 1 + end + end + if unit + uplo && (newkk = newk) + newrowval[newkk] = j + newnzval[newkk] = one(Tv) + newk += uplo + end + newcolptr[j+1] = newk + end + nz = newcolptr[n+1] - 1 + resize!(newrowval, nz) + resize!(newnzval, nz) + SparseMatrixCSC(m, n, newcolptr, newrowval, newnzval) +end + +# 8 cases: (Transpose|Adjoint){Tv,[Unit](Upper|Lower)Triangular} +function _sparsem(taA::Union{Transpose{Tv,<:AbstractTriangularSparse}, + Adjoint{Tv,<:AbstractTriangularSparse}}) where {Tv} + + sA = taA.parent + A = sA.data + rowval = rowvals(A) + nzval = nonzeros(A) + m, n = size(A) + Ti = eltype(rowval) + fnzrange = sA isa Union{UpperTriangular,UnitUpperTriangular} ? nzrangeup : nzrangelo + fadj = taA isa Transpose ? transpose : adjoint + unit = sA isa Union{UnitUpperTriangular,UnitLowerTriangular} + uplo = A isa Union{UpperTriangular,UnitUpperTriangular} + + newcolptr = Vector{Ti}(undef, n+1) + fill!(newcolptr, 1unit) + newcolptr[1] = 1 + @inbounds for j = 1:n + for k = fnzrange(A, j) + i = rowval[k] + if i != j || i == j && !unit + newcolptr[i+1] += 1 + end + end + end + cumsum!(newcolptr, newcolptr) + nz = newcolptr[n+1] - 1 + newrowval = Vector{Ti}(undef, nz) + newnzval = Vector{Tv}(undef, nz) + + @inbounds for j = 1:n + if !uplo && unit + ni = newcolptr[j] + newrowval[ni] = j + newnzval[ni] = fadj(one(Tv)) + newcolptr[j] = ni + 1 + end + for k = fnzrange(A, j) + i = rowval[k] + nzv = nzval[k] + if i != j || i == j && !unit + ni = newcolptr[i] + newrowval[ni] = j + newnzval[ni] = fadj(nzv) + newcolptr[i] = ni + 1 + end + end + if uplo && unit + ni = newcolptr[j] + newrowval[ni] = j + newnzval[ni] = fadj(one(Tv)) + newcolptr[j] = ni + 1 + end + end + _sparse_gen(n, m, newcolptr, newrowval, newnzval) +end + +function _sparse_gen(m, n, newcolptr, newrowval, newnzval) + @inbounds for j = n:-1:1 + newcolptr[j+1] = newcolptr[j] + end + newcolptr[1] = 1 + SparseMatrixCSC(m, n, newcolptr, newrowval, newnzval) +end + diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl index 05a7e68ad6a00..b4f02b2d0f163 100644 --- a/stdlib/SparseArrays/src/sparsematrix.jl +++ b/stdlib/SparseArrays/src/sparsematrix.jl @@ -431,319 +431,6 @@ SparseMatrixCSC{Tv}(M::Transpose{<:Any,SparseMatrixCSC}) where {Tv} = SparseMatr SparseMatrixCSC{Tv,Ti}(M::Adjoint{<:Any,SparseMatrixCSC}) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(copy(M)) SparseMatrixCSC{Tv,Ti}(M::Transpose{<:Any,SparseMatrixCSC}) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(copy(M)) -import LinearAlgebra: AbstractTriangular - -""" - SparseMatrixCSCSymmHerm - -`Symmetric` or `Hermitian` of a `SparseMatrixCSC` or `SparseMatrixCSCView`. -""" -const SparseMatrixCSCSymmHerm{Tv,Ti} = Union{Symmetric{Tv,<:SparseMatrixCSCUnion{Tv,Ti}}, - Hermitian{Tv,<:SparseMatrixCSCUnion{Tv,Ti}}} - -const AbstractTriangularSparse{Tv,Ti} = AbstractTriangular{Tv,<:SparseMatrixCSCUnion{Tv,Ti}} - -# converting Symmetric/Hermitian/AbstractTriangular/SubArray of SparseMatrixCSC -# and Transpose/Adjoint of AbstractTriangular of SparseMatrixCSC to SparseMatrixCSC -for wr in (Symmetric, Hermitian, Transpose, Adjoint, - UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular, - SubArray) - - @eval SparseMatrixCSC(A::$wr) = _sparsem(A) - @eval SparseMatrixCSC{Tv}(A::$wr{Tv}) where Tv = _sparsem(A) - @eval SparseMatrixCSC{Tv}(A::$wr) where Tv = SparseMatrixCSC{Tv}(_sparsem(A)) - @eval SparseMatrixCSC{Tv,Ti}(A::$wr) where {Tv,Ti} = SparseMatrixCSC{Tv,Ti}(_sparsem(A)) -end - -""" - iswrsparse(::S) - iswrsparse(::Type{S}) - -Returns `true` if type `S` is backed by a sparse array, and `false` otherwise. -""" -iswrsparse(::T) where T<:AbstractArray = iswrsparse(T) -iswrsparse(::Type) = false -iswrsparse(::Type{T}) where T<:AbstractSparseArray = true - -""" - depth(::Type{S}) - -Returns 0 for unwrapped S, and nesting depth for wrapped (nested) abstract arrays. -""" -depth(::T) where T = depth(T) -depth(::Type{T}) where T<:AbstractArray = 0 - -for wr in (Symmetric, Hermitian, - LowerTriangular, UnitLowerTriangular, UpperTriangular, UnitUpperTriangular, - Transpose, Adjoint, SubArray, - Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal) - - pl = wr === SubArray ? :($wr{<:Any,<:Any,T}) : :($wr{<:Any,T}) - @eval iswrsparse(::Type{<:$pl}) where T = iswrsparse(T) - @eval depth(::Type{<:$pl}) where T = depth(T) + 1 -end - -# convert parent and re-wrap in same wrapper -_sparsewrap(A::Symmetric) = Symmetric(_sparsem(parent(A)), A.uplo == 'U' ? :U : :L) -_sparsewrap(A::Hermitian) = Hermitian(_sparsem(parent(A)), A.uplo == 'U' ? :U : :L) -_sparsewrap(A::SubArray) = SubArray(_sparsem(parent(A)), A.indices) -for ty in ( LowerTriangular, UnitLowerTriangular, - UpperTriangular, UnitUpperTriangular, - Transpose, Adjoint) - - @eval _sparsewrap(A::$ty) = $ty(_sparsem(parent(A))) -end -function _sparsewrap(A::Union{Diagonal,Bidiagonal,Tridiagonal,SymTridiagonal}) - dropzeros!(sparse(A)) -end - -""" - unwrap(A::AbstractMatrix) - -In case A is a wrapper type (`SubArray, Symmetric, Adjoint, SubArray, Triangular, Tridiagonal`, etc.) -convert to `Matrix` or `SparseMatrixCSC`, depending on final storage type of A. -For other types return A itself. -""" -unwrap(A::Any) = A -unwrap(A::AbstractMatrix) = iswrsparse(A) ? convert(SparseMatrixCSC, A) : convert(Array, A) - -import Base.copy -copy(A::SubArray) = getindex(unwrap(parent(A)), A.indices...) - -# For pure sparse matrices and vectors return A. -# For wrapped sparse matrices or vectors convert to SparseMatrixCSC. -# Handle nested wrappers properly. -# Use abstract matrix fallback if A is not sparse. -function _sparsem(@nospecialize A::AbstractArray{Tv}) where Tv - if iswrsparse(A) - if depth(A) >= 1 - _sparsem(_sparsewrap(A)) - else - A - end - else - # explicitly call abstract matrix fallback using getindex(A,...) - invoke(SparseMatrixCSC{Tv,Int}, Tuple{AbstractMatrix}, A) - end -end - -_sparsem(A::AbstractSparseMatrix) = A -_sparsem(A::AbstractSparseVector) = A - -# Transpose/Adjoint of sparse vector (returning sparse matrix) -function _sparsem(A::Union{Transpose{<:Any,<:AbstractSparseVector},Adjoint{<:Any,<:AbstractSparseVector}}) - B = parent(A) - n = length(B) - Ti = eltype(B.nzind) - fadj = A isa Transpose ? transpose : adjoint - colptr = fill!(Vector{Ti}(undef, n + 1), 0) - colptr[1] = 1 - colptr[B.nzind .+ 1] .= 1 - cumsum!(colptr, colptr) - rowval = fill!(similar(B.nzind), 1) - nzval = fadj.(nonzeros(B)) - SparseMatrixCSC(1, n, colptr, rowval, nzval) -end - -function _sparsem(A::Union{Transpose{<:Any,<:SparseMatrixCSC},Adjoint{<:Any,<:SparseMatrixCSC}}) - ftranspose(parent(A), A isa Transpose ? transpose : adjoint) -end - -# Symmetric/Hermitian of sparse matrix -_sparsem(A::SparseMatrixCSCSymmHerm) = _sparsem(A.uplo == 'U' ? nzrangeup : nzrangelo, A) -# Triangular of sparse matrix -_sparsem(A::UpperTriangular{T,<:AbstractSparseMatrix}) where T = triu(A.data) -_sparsem(A::LowerTriangular{T,<:AbstractSparseMatrix}) where T = tril(A.data) -# view of sparse matrix -_sparsem(S::SubArray{<:Any,2,<:SparseMatrixCSC}) = getindex(S.parent,S.indices...) - -# 4 cases: (Symmetric|Hermitian) variants (:U|:L) -function _sparsem(fnzrange::Function, sA::SparseMatrixCSCSymmHerm{Tv}) where {Tv} - A = sA.data - rowval = rowvals(A) - nzval = nonzeros(A) - m, n = size(A) - Ti = eltype(rowval) - fadj = sA isa Symmetric ? transpose : adjoint - newcolptr = Vector{Ti}(undef, n+1) - diagmap = fadj == transpose ? identity : real - - newcolptr[1] = 1 - colrange = fnzrange === nzrangeup ? (1:n) : (n:-1:1) - @inbounds for j = colrange - r = fnzrange(A, j); r1 = r.start; r2 = r.stop - newcolptr[j+1] = r2 - r1 + 1 - for k = r1:r2 - row = rowval[k] - if row != j - newcolptr[row+1] += 1 - end - end - end - cumsum!(newcolptr, newcolptr) - nz = newcolptr[n+1] - 1 - newrowval = Vector{Ti}(undef, nz) - newnzval = Vector{Tv}(undef, nz) - @inbounds for j = 1:n - newk = newcolptr[j] - for k = fnzrange(A, j) - i = rowval[k] - nzv = nzval[k] - if i != j - newrowval[newk] = i - newnzval[newk] = nzv - newk += 1 - ni = newcolptr[i] - newrowval[ni] = j - newnzval[ni] = fadj(nzv) - newcolptr[i] = ni + 1 - else - newrowval[newk] = i - newnzval[newk] = diagmap(nzv) - newk += 1 - end - end - newcolptr[j] = newk - end - _sparse_gen(m, n, newcolptr, newrowval, newnzval) -end - -# 2 cases: Unit(Upper|Lower)Triangular{Tv,SparseMatrixCSC} -function _sparsem(A::AbstractTriangularSparse{Tv}) where Tv - S = A.data - rowval = rowvals(S) - nzval = nonzeros(S) - m, n = size(S) - Ti = eltype(rowval) - fnzrange = A isa Union{UpperTriangular,UnitUpperTriangular} ? nzrangeup : nzrangelo - unit = A isa Union{UnitUpperTriangular,UnitLowerTriangular} - newcolptr = Vector{Ti}(undef, n+1) - newrowval = Vector{Ti}(undef, nnz(S)) - newnzval = Vector{Tv}(undef, nnz(S)) - newcolptr[1] = 1 - uplo = fnzrange == nzrangeup - newk = 1 - @inbounds for j = 1:n - newkk = newk - if unit - newk += !uplo - end - r = fnzrange(S, j); r1 = r.start; r2 = r.stop - for k = r1:r2 - i = rowval[k] - if i != j || i == j && !unit - newrowval[newk] = i - newnzval[newk] = nzval[k] - newk += 1 - end - end - if unit - uplo && (newkk = newk) - newrowval[newkk] = j - newnzval[newkk] = one(Tv) - newk += uplo - end - newcolptr[j+1] = newk - end - nz = newcolptr[n+1] - 1 - resize!(newrowval, nz) - resize!(newnzval, nz) - SparseMatrixCSC(m, n, newcolptr, newrowval, newnzval) -end - -# 8 cases: (Transpose|Adjoint){Tv,[Unit](Upper|Lower)Triangular} -function _sparsem(taA::Union{Transpose{Tv,<:AbstractTriangularSparse}, - Adjoint{Tv,<:AbstractTriangularSparse}}) where {Tv} - - sA = taA.parent - A = sA.data - rowval = rowvals(A) - nzval = nonzeros(A) - m, n = size(A) - Ti = eltype(rowval) - fnzrange = sA isa Union{UpperTriangular,UnitUpperTriangular} ? nzrangeup : nzrangelo - fadj = taA isa Transpose ? transpose : adjoint - unit = sA isa Union{UnitUpperTriangular,UnitLowerTriangular} - uplo = A isa Union{UpperTriangular,UnitUpperTriangular} - - newcolptr = Vector{Ti}(undef, n+1) - fill!(newcolptr, 1unit) - newcolptr[1] = 1 - @inbounds for j = 1:n - for k = fnzrange(A, j) - i = rowval[k] - if i != j || i == j && !unit - newcolptr[i+1] += 1 - end - end - end - cumsum!(newcolptr, newcolptr) - nz = newcolptr[n+1] - 1 - newrowval = Vector{Ti}(undef, nz) - newnzval = Vector{Tv}(undef, nz) - - @inbounds for j = 1:n - if !uplo && unit - ni = newcolptr[j] - newrowval[ni] = j - newnzval[ni] = fadj(one(Tv)) - newcolptr[j] = ni + 1 - end - for k = fnzrange(A, j) - i = rowval[k] - nzv = nzval[k] - if i != j || i == j && !unit - ni = newcolptr[i] - newrowval[ni] = j - newnzval[ni] = fadj(nzv) - newcolptr[i] = ni + 1 - end - end - if uplo && unit - ni = newcolptr[j] - newrowval[ni] = j - newnzval[ni] = fadj(one(Tv)) - newcolptr[j] = ni + 1 - end - end - _sparse_gen(n, m, newcolptr, newrowval, newnzval) -end - -function _sparse_gen(m, n, newcolptr, newrowval, newnzval) - @inbounds for j = n:-1:1 - newcolptr[j+1] = newcolptr[j] - end - newcolptr[1] = 1 - SparseMatrixCSC(m, n, newcolptr, newrowval, newnzval) -end - -""" - nzrangeup(A::AbstractMatrix, i::Integer) -restrict `nzrange(A, i)` at diagonal as end (if only upper part matters) -""" -nzrangeup(A, i) = nzrangeup(A, nzrange(A, i), i) -""" - nzrangelo(A::AbstractMatrix, i::Integer) -restrict `nzrange(A, i)` at diagonal as start (if only lower part matters) -""" -nzrangelo(A, i) = nzrangelo(A, nzrange(A, i), i) -function nzrangeup(A, r::AbstractUnitRange, i) - r1 = r.start; r2 = r.stop - rv = rowvals(A) - @inbounds r2 < r1 || rv[r2] <= i ? r : r1:searchsortedlast(rv, i, r1, r2, Forward) -end -function nzrangeup(A, r::AbstractVector{<:Integer}, i) - view(r, 1:searchsortedlast(view(rowvals(A), r), i)) -end -function nzrangelo(A, r::AbstractUnitRange, i) - r1 = r.start; r2 = r.stop - rv = rowvals(A) - @inbounds r2 < r1 || rv[r1] >= i ? r : searchsortedfirst(rv, i, r1, r2, Forward):r2 -end -function nzrangelo(A, r::AbstractVector{<:Integer}, i) - view(r, searchsortedfirst(view(rowvals(A), r), i):length(r)) -end - # converting from SparseMatrixCSC to other matrix types function Matrix(S::SparseMatrixCSC{Tv}) where Tv # Handle cases where zero(Tv) is not defined but the array is dense.