diff --git a/NEWS.md b/NEWS.md index 8321cb599411a..d048c6f2cccd2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -113,6 +113,7 @@ Standard library changes #### SparseArrays +* new `sizehint!(::SparseMatrixCSC, ::Integer)` method ([#30676]). #### Dates diff --git a/stdlib/SparseArrays/src/higherorderfns.jl b/stdlib/SparseArrays/src/higherorderfns.jl index 383211267ee3b..a5941da764883 100644 --- a/stdlib/SparseArrays/src/higherorderfns.jl +++ b/stdlib/SparseArrays/src/higherorderfns.jl @@ -8,7 +8,7 @@ import Base: map, map!, broadcast, copy, copyto!, _extrema_dims, _extrema_itr using Base: front, tail, to_shape using ..SparseArrays: SparseVector, SparseMatrixCSC, AbstractSparseVector, AbstractSparseMatrixCSC, - AbstractSparseMatrix, AbstractSparseArray, indtype, nnz, nzrange, + AbstractSparseMatrix, AbstractSparseArray, indtype, nnz, nzrange, spzeros, SparseVectorUnion, AdjOrTransSparseVectorUnion, nonzeroinds, nonzeros, rowvals, getcolptr, widelength using Base.Broadcast: BroadcastStyle, Broadcasted, flatten using LinearAlgebra @@ -132,12 +132,17 @@ function trimstorage!(A::SparseVecOrMat, maxstored) resize!(storedvals(A), maxstored) return maxstored end + function expandstorage!(A::SparseVecOrMat, maxstored) - length(storedinds(A)) < maxstored && resize!(storedinds(A), maxstored) - length(storedvals(A)) < maxstored && resize!(storedvals(A), maxstored) + if length(storedinds(A)) < maxstored + resize!(storedinds(A), maxstored) + resize!(storedvals(A), maxstored) + end return maxstored end +_checkbuffers(S::SparseMatrixCSC) = (@assert length(getcolptr(S)) == size(S, 2) + 1 && getcolptr(S)[end] - 1 == length(rowvals(S)) == length(nonzeros(S)); S) +_checkbuffers(S::SparseVector) = (@assert length(storedvals(S)) == length(storedinds(S)); S) # (2) map[!] entry points map(f::Tf, A::SparseVector) where {Tf} = _noshapecheck_map(f, A) @@ -181,7 +186,7 @@ copy(bc::SpBroadcasted1) = _noshapecheck_map(bc.f, bc.args[1]) storedvals(C)[1] = fofnoargs broadcast!(f, view(storedvals(C), 2:length(storedvals(C)))) end - return C + return _checkbuffers(C) end function _diffshape_broadcast(f::Tf, A::SparseVecOrMat, Bs::Vararg{SparseVecOrMat,N}) where {Tf,N} @@ -224,22 +229,17 @@ _maxnnzfrom(shape::NTuple{2}, A::AbstractSparseMatrixCSC) = nnz(A) * div(shape[1 @inline _unchecked_maxnnzbcres(shape, As...) = _unchecked_maxnnzbcres(shape, As) @inline _checked_maxnnzbcres(shape::NTuple{1}, As...) = shape[1] != 0 ? _unchecked_maxnnzbcres(shape, As) : 0 @inline _checked_maxnnzbcres(shape::NTuple{2}, As...) = shape[1] != 0 && shape[2] != 0 ? _unchecked_maxnnzbcres(shape, As) : 0 -@inline function _allocres(shape::NTuple{1}, indextype, entrytype, maxnnz) - storedinds = Vector{indextype}(undef, maxnnz) - storedvals = Vector{entrytype}(undef, maxnnz) - return SparseVector(shape..., storedinds, storedvals) -end -@inline function _allocres(shape::NTuple{2}, indextype, entrytype, maxnnz) - pointers = ones(indextype, shape[2] + 1) - storedinds = Vector{indextype}(undef, maxnnz) - storedvals = Vector{entrytype}(undef, maxnnz) - return SparseMatrixCSC(shape..., pointers, storedinds, storedvals) +@inline function _allocres(shape::Union{NTuple{1},NTuple{2}}, indextype, entrytype, maxnnz) + X = spzeros(entrytype, indextype, shape) + resize!(storedinds(X), maxnnz) + resize!(storedvals(X), maxnnz) + return X end # (4) _map_zeropres!/_map_notzeropres! specialized for a single sparse vector/matrix "Stores only the nonzero entries of `map(f, Array(A))` in `C`." function _map_zeropres!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat) where Tf - spaceC::Int = min(length(storedinds(C)), length(storedvals(C))) + spaceC::Int = length(nonzeros(C)) Ck = 1 @inbounds for j in columns(C) setcolptr!(C, j, Ck) @@ -255,7 +255,7 @@ function _map_zeropres!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat) where Tf end @inbounds setcolptr!(C, numcols(C) + 1, Ck) trimstorage!(C, Ck - 1) - return C + return _checkbuffers(C) end """ Densifies `C`, storing `fillvalue` in place of each unstored entry in `A` and @@ -274,7 +274,7 @@ function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseVecOrMa end # NOTE: Combining the fill! above into the loop above to avoid multiple sweeps over / # nonsequential access of storedvals(C) does not appear to improve performance. - return C + return _checkbuffers(C) end # helper functions for these methods and some of those below @inline _densecoloffsets(A::SparseVector) = 0 @@ -297,7 +297,7 @@ end # (5) _map_zeropres!/_map_notzeropres! specialized for a pair of sparse vectors/matrices function _map_zeropres!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat, B::SparseVecOrMat) where Tf - spaceC::Int = min(length(storedinds(C)), length(storedvals(C))) + spaceC::Int = length(nonzeros(C)) rowsentinelA = convert(indtype(A), numrows(C) + 1) rowsentinelB = convert(indtype(B), numrows(C) + 1) Ck = 1 @@ -336,7 +336,7 @@ function _map_zeropres!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat, B::SparseVe end @inbounds setcolptr!(C, numcols(C) + 1, Ck) trimstorage!(C, Ck - 1) - return C + return _checkbuffers(C) end function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseVecOrMat, B::SparseVecOrMat) where Tf # Build dense matrix structure in C, expanding storage if necessary @@ -368,13 +368,13 @@ function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseVecOrMa Cx != fillvalue && (storedvals(C)[jo + Ci] = Cx) end end - return C + return _checkbuffers(C) end # (6) _map_zeropres!/_map_notzeropres! for more than two sparse matrices / vectors function _map_zeropres!(f::Tf, C::SparseVecOrMat, As::Vararg{SparseVecOrMat,N}) where {Tf,N} - spaceC::Int = min(length(storedinds(C)), length(storedvals(C))) + spaceC::Int = length(nonzeros(C)) rowsentinel = numrows(C) + 1 Ck = 1 stopks = _colstartind_all(1, As) @@ -398,7 +398,7 @@ function _map_zeropres!(f::Tf, C::SparseVecOrMat, As::Vararg{SparseVecOrMat,N}) end @inbounds setcolptr!(C, numcols(C) + 1, Ck) trimstorage!(C, Ck - 1) - return C + return _checkbuffers(C) end function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, As::Vararg{SparseVecOrMat,N}) where {Tf,N} # Build dense matrix structure in C, expanding storage if necessary @@ -421,7 +421,7 @@ function _map_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, As::Vararg{Spars activerow = min(rows...) end end - return C + return _checkbuffers(C) end # helper methods for map/map! methods just above @@ -462,7 +462,7 @@ end # (7) _broadcast_zeropres!/_broadcast_notzeropres! specialized for a single (input) sparse vector/matrix function _broadcast_zeropres!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat) where Tf isempty(C) && return _finishempty!(C) - spaceC::Int = min(length(storedinds(C)), length(storedvals(C))) + spaceC::Int = length(nonzeros(C)) # C and A cannot have the same shape, as we directed that case to map in broadcast's # entry point; here we need efficiently handle only heterogeneous C-A combinations where # one or both of C and A has at least one singleton dimension. @@ -509,7 +509,7 @@ function _broadcast_zeropres!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat) where end @inbounds setcolptr!(C, numcols(C) + 1, Ck) trimstorage!(C, Ck - 1) - return C + return _checkbuffers(C) end function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseVecOrMat) where Tf # For information on this code, see comments in similar code in _broadcast_zeropres! above @@ -540,14 +540,14 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV end end end - return C + return _checkbuffers(C) end # (8) _broadcast_zeropres!/_broadcast_notzeropres! specialized for a pair of (input) sparse vectors/matrices function _broadcast_zeropres!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat, B::SparseVecOrMat) where Tf isempty(C) && return _finishempty!(C) - spaceC::Int = min(length(storedinds(C)), length(storedvals(C))) + spaceC::Int = length(nonzeros(C)) rowsentinelA = convert(indtype(A), numrows(C) + 1) rowsentinelB = convert(indtype(B), numrows(C) + 1) # C, A, and B cannot all have the same shape, as we directed that case to map in broadcast's @@ -711,7 +711,7 @@ function _broadcast_zeropres!(f::Tf, C::SparseVecOrMat, A::SparseVecOrMat, B::Sp end @inbounds setcolptr!(C, numcols(C) + 1, Ck) trimstorage!(C, Ck - 1) - return C + return _checkbuffers(C) end function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseVecOrMat, B::SparseVecOrMat) where Tf # For information on this code, see comments in similar code in _broadcast_zeropres! above @@ -810,7 +810,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, A::SparseV end end end - return C + return _checkbuffers(C) end _finishempty!(C::SparseVector) = C _finishempty!(C::AbstractSparseMatrixCSC) = (fill!(getcolptr(C), 1); C) @@ -861,7 +861,7 @@ end # (9) _broadcast_zeropres!/_broadcast_notzeropres! for more than two (input) sparse vectors/matrices function _broadcast_zeropres!(f::Tf, C::SparseVecOrMat, As::Vararg{SparseVecOrMat,N}) where {Tf,N} isempty(C) && return _finishempty!(C) - spaceC::Int = min(length(storedinds(C)), length(storedvals(C))) + spaceC::Int = length(nonzeros(C)) expandsverts = _expandsvert_all(C, As) expandshorzs = _expandshorz_all(C, As) rowsentinel = numrows(C) + 1 @@ -909,7 +909,7 @@ function _broadcast_zeropres!(f::Tf, C::SparseVecOrMat, As::Vararg{SparseVecOrMa end @inbounds setcolptr!(C, numcols(C) + 1, Ck) trimstorage!(C, Ck - 1) - return C + return _checkbuffers(C) end function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, As::Vararg{SparseVecOrMat,N}) where {Tf,N} isempty(C) && return _finishempty!(C) @@ -950,7 +950,7 @@ function _broadcast_notzeropres!(f::Tf, fillvalue, C::SparseVecOrMat, As::Vararg end end end - return C + return _checkbuffers(C) end # helper method for broadcast/broadcast! methods just above diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl index 0ac2806b77a04..73fc0c7dee03a 100644 --- a/stdlib/SparseArrays/src/linalg.jl +++ b/stdlib/SparseArrays/src/linalg.jl @@ -935,16 +935,15 @@ function triu(S::AbstractSparseMatrixCSC{Tv,Ti}, k::Integer=0) where {Tv,Ti} end rowval = Vector{Ti}(undef, nnz) nzval = Vector{Tv}(undef, nnz) - A = SparseMatrixCSC(m, n, colptr, rowval, nzval) for col = max(k+1,1) : n c1 = getcolptr(S)[col] - for c2 in nzrange(A, col) - rowvals(A)[c2] = rowvals(S)[c1] - nonzeros(A)[c2] = nonzeros(S)[c1] + for c2 in colptr[col]:colptr[col+1]-1 + rowval[c2] = rowvals(S)[c1] + nzval[c2] = nonzeros(S)[c1] c1 += 1 end end - A + SparseMatrixCSC(m, n, colptr, rowval, nzval) end function tril(S::AbstractSparseMatrixCSC{Tv,Ti}, k::Integer=0) where {Tv,Ti} @@ -965,17 +964,16 @@ function tril(S::AbstractSparseMatrixCSC{Tv,Ti}, k::Integer=0) where {Tv,Ti} end rowval = Vector{Ti}(undef, nnz) nzval = Vector{Tv}(undef, nnz) - A = SparseMatrixCSC(m, n, colptr, rowval, nzval) for col = 1 : min(n, m+k) c1 = getcolptr(S)[col+1]-1 - l2 = getcolptr(A)[col+1]-1 - for c2 = 0 : l2 - getcolptr(A)[col] - rowvals(A)[l2 - c2] = rowvals(S)[c1] - nonzeros(A)[l2 - c2] = nonzeros(S)[c1] + l2 = colptr[col+1]-1 + for c2 = 0 : l2 - colptr[col] + rowval[l2 - c2] = rowvals(S)[c1] + nzval[l2 - c2] = nonzeros(S)[c1] c1 -= 1 end end - A + SparseMatrixCSC(m, n, colptr, rowval, nzval) end ## diff @@ -1340,7 +1338,6 @@ end ## kron @inline function kron!(C::SparseMatrixCSC, A::AbstractSparseMatrixCSC, B::AbstractSparseMatrixCSC) - nnzC = nnz(A)*nnz(B) mA, nA = size(A); mB, nB = size(B) mC, nC = mA*mB, nA*nB @@ -1348,11 +1345,9 @@ end nzvalC = nonzeros(C) colptrC = getcolptr(C) - @boundscheck begin - length(colptrC) == nC+1 || throw(DimensionMismatch("expect C to be preallocated with $(nC+1) colptrs ")) - length(rowvalC) == nnzC || throw(DimensionMismatch("expect C to be preallocated with $(nnzC) rowvals")) - length(nzvalC) == nnzC || throw(DimensionMismatch("expect C to be preallocated with $(nnzC) nzvals")) - end + nnzC = nnz(A)*nnz(B) + resize!(nzvalC, nnzC) + resize!(rowvalC, nnzC) col = 1 @inbounds for j = 1:nA @@ -1381,16 +1376,13 @@ end end @inline function kron!(z::SparseVector, x::SparseVector, y::SparseVector) - nnzx = nnz(x); nnzy = nnz(y); nnzz = nnz(z); + nnzx = nnz(x); nnzy = nnz(y); nzind = nonzeroinds(z) nzval = nonzeros(z) - @boundscheck begin - nnzval = length(nzval); nnzind = length(nzind) - nnzz = nnzx*nnzy - nnzval == nnzz || throw(DimensionMismatch("expect z to be preallocated with $nnzz nonzeros")) - nnzind == nnzz || throw(DimensionMismatch("expect z to be preallocated with $nnzz nonzeros")) - end + nnzz = nnzx*nnzy + resize!(nzind, nnzz) + resize!(nzval, nnzz) @inbounds for i = 1:nnzx, j = 1:nnzy this_ind = (i-1)*nnzy+j @@ -1402,17 +1394,12 @@ end # sparse matrix ⊗ sparse matrix function kron(A::AbstractSparseMatrixCSC{T1,S1}, B::AbstractSparseMatrixCSC{T2,S2}) where {T1,S1,T2,S2} - nnzC = nnz(A)*nnz(B) mA, nA = size(A); mB, nB = size(B) mC, nC = mA*mB, nA*nB Tv = typeof(one(T1)*one(T2)) Ti = promote_type(S1,S2) - colptrC = Vector{Ti}(undef, nC+1) - rowvalC = Vector{Ti}(undef, nnzC) - nzvalC = Vector{Tv}(undef, nnzC) - colptrC[1] = 1 - # skip sparse_check - C = SparseMatrixCSC{Tv, Ti}(mC, nC, colptrC, rowvalC, nzvalC) + C = spzeros(Tv, Ti, mC, nC) + sizehint!(C, nnz(A)*nnz(B)) return @inbounds kron!(C, A, B) end diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl index 888db226ddd14..58dfd88abd339 100644 --- a/stdlib/SparseArrays/src/sparsematrix.jl +++ b/stdlib/SparseArrays/src/sparsematrix.jl @@ -25,10 +25,9 @@ struct SparseMatrixCSC{Tv,Ti<:Integer} <: AbstractSparseMatrixCSC{Tv,Ti} function SparseMatrixCSC{Tv,Ti}(m::Integer, n::Integer, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}) where {Tv,Ti<:Integer} - @noinline throwsz(str, lbl, k) = - throw(ArgumentError("number of $str ($lbl) must be ≥ 0, got $k")) - m < 0 && throwsz("rows", 'm', m) - n < 0 && throwsz("columns", 'n', n) + sparse_check_Ti(m, n, Ti) + _goodbuffers(Int(m), Int(n), colptr, rowval, nzval) || + throw(ArgumentError("Illegal buffers for SparseMatrixCSC construction $n $colptr $rowval $nzval")) new(Int(m), Int(n), colptr, rowval, nzval) end end @@ -80,6 +79,15 @@ end size(S::SparseMatrixCSC) = (getfield(S, :m), getfield(S, :n)) +_goodbuffers(S::SparseMatrixCSC) = _goodbuffers(size(S)..., getcolptr(S), getrowval(S), nonzeros(S)) +_checkbuffers(S::SparseMatrixCSC) = (@assert _goodbuffers(S); S) + +function _goodbuffers(m, n, colptr, rowval, nzval) + (length(colptr) == n + 1 && colptr[end] - 1 == length(rowval) == length(nzval)) + # stronger check for debugging purposes + # && all(issorted(@view rowval[colptr[i]:colptr[i+1]-1]) for i=1:n) +end + # Define an alias for views of a SparseMatrixCSC which include all rows and a unit range of the columns. # Also define a union of SparseMatrixCSC and this view since many methods can be defined efficiently for # this union by extracting the fields via the get function: getcolptr, getrowval, and getnzval. The key @@ -220,6 +228,7 @@ Base.replace_in_print_matrix(A::AbstractSparseMatrix, i::Integer, j::Integer, s: Base.isstored(A, i, j) ? s : Base.replace_with_centered_mark(s) function Base.show(io::IO, ::MIME"text/plain", S::AbstractSparseMatrixCSC) + _checkbuffers(S) xnnz = nnz(S) m, n = size(S) print(io, m, "×", n, " ", typeof(S), " with ", xnnz, " stored ", @@ -295,6 +304,7 @@ function _show_with_braille_patterns(io::IOContext, S::AbstractSparseMatrixCSC) end function Base.show(io::IOContext, S::AbstractSparseMatrixCSC) + _checkbuffers(S) if max(size(S)...) < 16 && !(get(io, :compact, false)::Bool) ioc = IOContext(io, :compact => true) println(ioc) @@ -307,7 +317,7 @@ end ## Reshape -function sparse_compute_reshaped_colptr_and_rowval(colptrS::Vector{Ti}, rowvalS::Vector{Ti}, +function sparse_compute_reshaped_colptr_and_rowval!(colptrS::Vector{Ti}, rowvalS::Vector{Ti}, mS::Int, nS::Int, colptrA::Vector{Ti}, rowvalA::Vector{Ti}, mA::Int, nA::Int) where Ti lrowvalA = length(rowvalA) @@ -350,7 +360,7 @@ function copy(ra::ReshapedArray{<:Any,2,<:AbstractSparseMatrixCSC}) rowval = similar(rowvals(a)) nzval = copy(nonzeros(a)) - sparse_compute_reshaped_colptr_and_rowval(colptr, rowval, mS, nS, getcolptr(a), rowvals(a), mA, nA) + sparse_compute_reshaped_colptr_and_rowval!(colptr, rowval, mS, nS, getcolptr(a), rowvals(a), mA, nA) return SparseMatrixCSC(mS, nS, colptr, rowval, nzval) end @@ -377,7 +387,7 @@ function copyto!(A::AbstractSparseMatrixCSC, B::AbstractSparseMatrixCSC) copyto!(rowvals(A), rowvals(B)) else # This is like a "reshape B into A". - sparse_compute_reshaped_colptr_and_rowval(getcolptr(A), rowvals(A), size(A, 1), size(A, 2), getcolptr(B), rowvals(B), size(B, 1), size(B, 2)) + sparse_compute_reshaped_colptr_and_rowval!(getcolptr(A), rowvals(A), size(A, 1), size(A, 2), getcolptr(B), rowvals(B), size(B, 1), size(B, 2)) end else widelength(A) >= widelength(B) || throw(BoundsError()) @@ -407,10 +417,10 @@ function copyto!(A::AbstractSparseMatrixCSC, B::AbstractSparseMatrixCSC) @inbounds for i in 2:length(getcolptr(A)) getcolptr(A)[i] += nnzB - lastmodptrA end - sparse_compute_reshaped_colptr_and_rowval(getcolptr(A), rowvals(A), size(A, 1), lastmodcolA-1, getcolptr(B), rowvals(B), size(B, 1), size(B, 2)) + sparse_compute_reshaped_colptr_and_rowval!(getcolptr(A), rowvals(A), size(A, 1), lastmodcolA-1, getcolptr(B), rowvals(B), size(B, 1), size(B, 2)) end copyto!(nonzeros(A), nonzeros(B)) - return A + return _checkbuffers(A) end copyto!(A::AbstractMatrix, B::AbstractSparseMatrixCSC) = _sparse_copyto!(A, B) @@ -471,7 +481,7 @@ function _sparsesimilar(S::AbstractSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew} end # parent methods for similar that preserves only storage space (for when new dims are 2-d) _sparsesimilar(S::AbstractSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{2}) where {TvNew,TiNew} = - SparseMatrixCSC(dims..., fill(one(TiNew), last(dims)+1), similar(rowvals(S), TiNew, 0), similar(nonzeros(S), TvNew, 0)) + sizehint!(spzeros(TvNew, TiNew, dims...), length(nonzeros(S))) # parent method for similar that allocates an empty sparse vector (for when new dims are 1-d) _sparsesimilar(S::AbstractSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{1}) where {TvNew,TiNew} = SparseVector(dims..., similar(rowvals(S), TiNew, 0), similar(nonzeros(S), TvNew, 0)) @@ -496,6 +506,12 @@ similar(S::AbstractSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew}, m::Integer) wh similar(S::AbstractSparseMatrixCSC, ::Type{TvNew}, ::Type{TiNew}, m::Integer, n::Integer) where {TvNew,TiNew} = _sparsesimilar(S, TvNew, TiNew, (m, n)) +function Base.sizehint!(S::SparseMatrixCSC, n::Integer) + nhint = min(n, widelength(S)) + sizehint!(getrowval(S), nhint) + sizehint!(nonzeros(S), nhint) + return S +end # converting between SparseMatrixCSC types SparseMatrixCSC(S::AbstractSparseMatrixCSC) = copy(S) @@ -658,6 +674,7 @@ SparseMatrixCSC{Tv,Ti}(M::Transpose{<:Any,<:AbstractSparseMatrixCSC}) where {Tv, # converting from SparseMatrixCSC to other matrix types function Matrix(S::AbstractSparseMatrixCSC{Tv}) where Tv + _checkbuffers(S) A = Matrix{Tv}(undef, size(S, 1), size(S, 2)) copyto!(A, S) return A @@ -1016,6 +1033,8 @@ respectively. Simultaneously fixes the one-position-forward shift in `getcolptr( """ @noinline function _distributevals_halfperm!(X::AbstractSparseMatrixCSC{Tv,Ti}, A::AbstractSparseMatrixCSC{TvA,Ti}, q::AbstractVector{<:Integer}, f::Function) where {Tv,TvA,Ti} + resize!(nonzeros(X), nnz(A)) + resize!(rowvals(X), nnz(A)) @inbounds for Xi in 1:size(A, 2) Aj = q[Xi] for Ak in nzrange(A, Aj) @@ -1035,16 +1054,8 @@ function ftranspose!(X::AbstractSparseMatrixCSC{Tv,Ti}, A::AbstractSparseMatrixC throw(DimensionMismatch(string("destination argument `X`'s column count, ", "`size(X, 2) (= $(size(X, 2)))`, must match source argument `A`'s row count, `size(A, 1) (= $(size(A, 1)))`"))) elseif size(X, 1) != size(A, 2) - throw(DimensionMismatch(string("destination argument `X`'s row count, - `size(X, 1) (= $(size(X, 1)))`, must match source argument `A`'s column count, `size(A, 2) (= $(size(A, 2)))`"))) - elseif length(rowvals(X)) < nnz(A) - throw(ArgumentError(string("the length of destination argument `X`'s `rowval` ", - "array, `length(rowvals(X)) (= $(length(rowvals(X))))`, must be greater than or ", - "equal to source argument `A`'s allocated entry count, `nnz(A) (= $(nnz(A)))`"))) - elseif length(nonzeros(X)) < nnz(A) - throw(ArgumentError(string("the length of destination argument `X`'s `nzval` ", - "array, `length(nonzeros(X)) (= $(length(nonzeros(X))))`, must be greater than or ", - "equal to source argument `A`'s allocated entry count, `nnz(A) (= $(nnz(A)))`"))) + throw(DimensionMismatch(string("destination argument `X`'s row count, ", + "`size(X, 1) (= $(size(X, 1)))`, must match source argument `A`'s column count, `size(A, 2) (= $(size(A, 2)))`"))) end halfperm!(X, A, 1:size(A, 2), f) end @@ -1055,8 +1066,9 @@ adjoint!(X::AbstractSparseMatrixCSC{Tv,Ti}, A::AbstractSparseMatrixCSC{Tv,Ti}) w function ftranspose(A::AbstractSparseMatrixCSC{TvA,Ti}, f::Function, eltype::Type{Tv} = TvA) where {Tv,TvA,Ti} X = SparseMatrixCSC(size(A, 2), size(A, 1), ones(Ti, size(A, 1)+1), - Vector{Ti}(undef, nnz(A)), - Vector{Tv}(undef, nnz(A))) + Vector{Ti}(undef, 0), + Vector{Tv}(undef, 0)) + sizehint!(X, nnz(A)) halfperm!(X, A, 1:size(A, 2), f) end adjoint(A::AbstractSparseMatrixCSC) = Adjoint(A) @@ -1274,10 +1286,11 @@ function permute!(X::AbstractSparseMatrixCSC{Tv,Ti}, A::AbstractSparseMatrixCSC{ p::AbstractVector{<:Integer}, q::AbstractVector{<:Integer}) where {Tv,Ti} _checkargs_sourcecompatdest_permute!(A, X) _checkargs_sourcecompatperms_permute!(A, p, q) - C = SparseMatrixCSC(size(A, 2), size(A, 1), - ones(Ti, size(A, 1) + 1), - Vector{Ti}(undef, nnz(A)), - Vector{Tv}(undef, nnz(A))) + # bypass strict buffer checking + C = spzeros(Tv, Ti, size(A,2), size(A,1)) + resize!(getrowval(C), nnz(A)) + resize!(getnzval(C), nnz(A)) + _checkargs_permutationsvalid_permute!(p, getcolptr(C), q, getcolptr(X)) unchecked_noalias_permute!(X, A, p, q, C) end @@ -1293,10 +1306,9 @@ end function permute!(A::AbstractSparseMatrixCSC{Tv,Ti}, p::AbstractVector{<:Integer}, q::AbstractVector{<:Integer}) where {Tv,Ti} _checkargs_sourcecompatperms_permute!(A, p, q) - C = SparseMatrixCSC(size(A, 2), size(A, 1), - ones(Ti, size(A, 1) + 1), - Vector{Ti}(undef, nnz(A)), - Vector{Tv}(undef, nnz(A))) + C = spzeros(Tv, Ti, size(A,2), size(A,1)) + resize!(getrowval(C), nnz(A)) + resize!(getnzval(C), nnz(A)) workcolptr = Vector{Ti}(undef, size(A, 2) + 1) _checkargs_permutationsvalid_permute!(p, getcolptr(C), q, workcolptr) unchecked_aliasing_permute!(A, p, q, C, workcolptr) @@ -1355,14 +1367,14 @@ julia> permute(A, [1, 2, 3, 4], [4, 3, 2, 1]) function permute(A::AbstractSparseMatrixCSC{Tv,Ti}, p::AbstractVector{<:Integer}, q::AbstractVector{<:Integer}) where {Tv,Ti} _checkargs_sourcecompatperms_permute!(A, p, q) - X = SparseMatrixCSC(size(A, 1), size(A, 2), - ones(Ti, size(A, 2) + 1), - Vector{Ti}(undef, nnz(A)), - Vector{Tv}(undef, nnz(A))) - C = SparseMatrixCSC(size(A, 2), size(A, 1), - ones(Ti, size(A, 1) + 1), - Vector{Ti}(undef, nnz(A)), - Vector{Tv}(undef, nnz(A))) + # bypass strict buffer checking + X = spzeros(Tv, Ti, size(A,1), size(A,2)) + resize!(getrowval(X), nnz(A)) + resize!(getnzval(X), nnz(A)) + # bypass strict buffer checking + C = spzeros(Tv, Ti, size(A,2), size(A,1)) + resize!(getrowval(C), nnz(A)) + resize!(getnzval(C), nnz(A)) _checkargs_permutationsvalid_permute!(p, getcolptr(C), q, getcolptr(X)) unchecked_noalias_permute!(X, A, p, q, C) end @@ -1455,6 +1467,7 @@ For an out-of-place version, see [`dropzeros`](@ref). For algorithmic information, see `fkeep!`. """ dropzeros!(A::AbstractSparseMatrixCSC) = fkeep!(A, (i, j, x) -> !iszero(x)) + """ dropzeros(A::AbstractSparseMatrixCSC;) @@ -2480,7 +2493,7 @@ function permute_rows!(S::AbstractSparseMatrixCSC{Tv,Ti}, pI::Vector{Int}) where k += 1 end end - S + return _checkbuffers(S) end function getindex_general(A::AbstractSparseMatrixCSC, I::AbstractVector, J::AbstractVector) @@ -2634,6 +2647,7 @@ function Base.fill!(V::SubArray{Tv, <:Any, <:AbstractSparseMatrixCSC{Tv}, <:Tupl else _spsetnz_setindex!(A, convert(Tv, x), I, J) end + return _checkbuffers(A) end """ Helper method for immediately preceding fill! method. For all (i,j) such that i in I and @@ -2917,7 +2931,7 @@ function setindex!(A::AbstractSparseMatrixCSC{Tv,Ti}, V::AbstractVecOrMat, Ix::U deleteat!(rowvalA, colptrA[end]:length(rowvalA)) deleteat!(nzvalA, colptrA[end]:length(nzvalA)) - return A + return _checkbuffers(A) end # Logical setindex! @@ -3026,7 +3040,7 @@ function setindex!(A::AbstractSparseMatrixCSC, x::AbstractArray, I::AbstractMatr deleteat!(rowvalB, bidx:n) end end - A + return _checkbuffers(A) end function setindex!(A::AbstractSparseMatrixCSC, x::AbstractArray, Ix::AbstractVector{<:Integer}) @@ -3137,7 +3151,7 @@ function setindex!(A::AbstractSparseMatrixCSC, x::AbstractArray, Ix::AbstractVec deleteat!(rowvalB, bidx:n) end end - A + return _checkbuffers(A) end ## dropstored! methods @@ -3173,7 +3187,7 @@ function dropstored!(A::AbstractSparseMatrixCSC, i::Integer, j::Integer) @inbounds getcolptr(A)[m] -= 1 end end - return A + return _checkbuffers(A) end """ dropstored!(A::AbstractSparseMatrixCSC, I::AbstractVector{<:Integer}, J::AbstractVector{<:Integer}) @@ -3260,7 +3274,7 @@ function dropstored!(A::AbstractSparseMatrixCSC, deleteat!(rowvalA, rowidx:nnzA) deleteat!(nzvalA, rowidx:nnzA) end - return A + return _checkbuffers(A) end dropstored!(A::AbstractSparseMatrixCSC, i::Integer, J::AbstractVector{<:Integer}) = dropstored!(A, [i], J) dropstored!(A::AbstractSparseMatrixCSC, I::AbstractVector{<:Integer}, j::Integer) = dropstored!(A, I, [j]) @@ -3694,74 +3708,6 @@ function tr(A::AbstractSparseMatrixCSC{Tv}) where Tv return s end - -# Sort all the indices in each column of a CSC sparse matrix -# sortSparseMatrixCSC!(A, sortindices = :sortcols) # Sort each column with sort() -# sortSparseMatrixCSC!(A, sortindices = :doubletranspose) # Sort with a double transpose -function sortSparseMatrixCSC!(A::AbstractSparseMatrixCSC{Tv,Ti}; sortindices::Symbol = :sortcols) where {Tv,Ti} - if sortindices === :doubletranspose - nB, mB = size(A) - B = SparseMatrixCSC(mB, nB, Vector{Ti}(undef, nB+1), similar(rowvals(A)), similar(nonzeros(A))) - transpose!(B, A) - transpose!(A, B) - return A - end - - m, n = size(A) - colptr = getcolptr(A); rowval = rowvals(A); nzval = nonzeros(A) - - index = zeros(Ti, m) - row = zeros(Ti, m) - val = zeros(Tv, m) - - perm = Base.Perm(Base.ord(isless, identity, false, Base.Order.Forward), row) - - @inbounds for i = 1:n - nzr = nzrange(A, i) - numrows = length(nzr) - if numrows <= 1 - continue - elseif numrows == 2 - f = first(nzr) - s = f+1 - if rowval[f] > rowval[s] - rowval[f], rowval[s] = rowval[s], rowval[f] - nzval[f], nzval[s] = nzval[s], nzval[f] - end - continue - end - resize!(row, numrows) - resize!(index, numrows) - - jj = 1 - @simd for j = nzr - row[jj] = rowval[j] - val[jj] = nzval[j] - jj += 1 - end - - if numrows <= 16 - alg = Base.Sort.InsertionSort - else - alg = Base.Sort.QuickSort - end - - # Reset permutation - index .= 1:numrows - - sort!(index, alg, perm) - - jj = 1 - @simd for j = nzr - rowval[j] = row[index[jj]] - nzval[j] = val[index[jj]] - jj += 1 - end - end - - return A -end - ## rotations function rot180(A::AbstractSparseMatrixCSC) @@ -3837,7 +3783,7 @@ function circshift!(O::AbstractSparseMatrixCSC, X::AbstractSparseMatrixCSC, (r,c @inbounds for i=1:size(O, 2) subvector_shifter!(rowvals(O), nonzeros(O), getcolptr(O)[i], getcolptr(O)[i+1]-1, size(O, 1), r) end - return O + return _checkbuffers(O) end circshift!(O::AbstractSparseMatrixCSC, X::AbstractSparseMatrixCSC, (r,)::Base.DimsInteger{1}) = circshift!(O, X, (r,0)) diff --git a/stdlib/SparseArrays/src/sparsevector.jl b/stdlib/SparseArrays/src/sparsevector.jl index 52142553e5019..f893fb04eadf3 100644 --- a/stdlib/SparseArrays/src/sparsevector.jl +++ b/stdlib/SparseArrays/src/sparsevector.jl @@ -84,6 +84,13 @@ rowvals(x::SparseVectorUnion) = nonzeroinds(x) indtype(x::SparseColumnView) = indtype(parent(x)) indtype(x::SparseVectorView) = indtype(parent(x)) + +function Base.sizehint!(v::SparseVector, newlen::Integer) + sizehint!(nonzeroinds(v), newlen) + sizehint!(nonzeros(v), newlen) + return v +end + ## similar # # parent method for similar that preserves stored-entry structure (for when new and old dims match) @@ -92,10 +99,11 @@ _sparsesimilar(S::SparseVector, ::Type{TvNew}, ::Type{TiNew}) where {TvNew,TiNew # parent method for similar that preserves nothing (for when new dims are 1-d) _sparsesimilar(S::SparseVector, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{1}) where {TvNew,TiNew} = SparseVector(dims..., similar(nonzeroinds(S), TiNew, 0), similar(nonzeros(S), TvNew, 0)) -# parent method for similar that preserves storage space (for when new dims are 2-d) -_sparsesimilar(S::SparseVector, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{2}) where {TvNew,TiNew} = - SparseMatrixCSC(dims..., fill(one(TiNew), last(dims)+1), similar(nonzeroinds(S), TiNew, 0), similar(nonzeros(S), TvNew, 0)) - +# parent method for similar that preserves storage space (for old and new dims differ, and new is 2d) +function _sparsesimilar(S::SparseVector, ::Type{TvNew}, ::Type{TiNew}, dims::Dims{2}) where {TvNew,TiNew} + S1 = SparseMatrixCSC(dims..., fill(one(TiNew), last(dims)+1), similar(nonzeroinds(S), TiNew, 0), similar(nonzeros(S), TvNew, 0)) + return sizehint!(S1, min(widelength(S1), length(nonzeroinds(S)))) +end # The following methods hook into the AbstractArray similar hierarchy. The first method # covers similar(A[, Tv]) calls, which preserve stored-entry structure, and the latter # methods cover similar(A[, Tv], shape...) calls, which preserve nothing if the dims diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index a622d2bc8bcc0..f2754ba5ade81 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/stdlib/SparseArrays/test/sparse.jl @@ -48,7 +48,7 @@ end S = sparse(I, 3, 3) fill!(S, 0) @test iszero(S) # test success with stored zeros via fill! - @test iszero(SparseMatrixCSC(2, 2, [1,2,3], [1,2], [0,0,1])) # test success with nonzeros beyond data range + @test_throws ArgumentError iszero(SparseMatrixCSC(2, 2, [1,2,3], [1,2], [0,0,1])) # test failure with nonzeros beyond data range end @testset "isone specialization for SparseMatrixCSC" begin @test isone(sparse(I, 3, 3)) # test success @@ -673,8 +673,6 @@ end @testset "common error checking of [c]transpose! methods (ftranspose!)" begin @test_throws DimensionMismatch transpose!(A[:, 1:(smalldim - 1)], A) @test_throws DimensionMismatch transpose!(A[1:(smalldim - 1), 1], A) - @test_throws ArgumentError transpose!((B = similar(A); resize!(rowvals(B), nnz(A) - 1); B), A) - @test_throws ArgumentError transpose!((B = similar(A); resize!(nonzeros(B), nnz(A) - 1); B), A) end @testset "common error checking of permute[!] methods / source-perm compat" begin @test_throws DimensionMismatch permute(A, p[1:(end - 1)], q) @@ -2416,19 +2414,6 @@ end @test String(take!(io)) == "⠛⠛" end -@testset "check buffers" for n in 1:3 - local A - rowval = [1,2,3] - nzval1 = Int[] - nzval2 = [1,1,1] - A = SparseMatrixCSC(n, n, [1:n+1;], rowval, nzval1) - @test nnz(A) == n - @test_throws BoundsError A[n,n] - A = SparseMatrixCSC(n, n, [1:n+1;], rowval, nzval2) - @test nnz(A) == n - @test A == Matrix(I, n, n) -end - @testset "reverse search direction if step < 0 #21986" begin local A, B A = guardseed(1234) do @@ -2530,8 +2515,6 @@ end # 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])) @test_throws MethodError count(SparseVector(2, Int[1], Any[true])) - # count should run only over nonzeros(S)[1:nnz(S)], not nonzeros(S) in full - @test count(SparseMatrixCSC(2, 2, Int[1, 2, 3], Int[1, 2], Bool[true, true, true])) == 2 end @testset "sparse findprev/findnext operations" begin @@ -2599,15 +2582,6 @@ end @test sum(s, dims=2) == reshape([1, 2, 3], 3, 1) end -@testset "mapreduce of sparse matrices with trailing elements in nzval #26534" begin - B = SparseMatrixCSC{Int,Int}(2, 3, - [1, 3, 4, 5], - [1, 2, 1, 2, 999, 999, 999, 999], - [1, 2, 3, 6, 999, 999, 999, 999] - ) - @test maximum(B) == 6 -end - _length_or_count_or_five(::Colon) = 5 _length_or_count_or_five(x::AbstractVector{Bool}) = count(x) _length_or_count_or_five(x) = length(x) @@ -2883,19 +2857,6 @@ end @test sparse(deepwrap(A)) == Matrix(deepwrap(B)) end -@testset "unary operations on matrices where length(nzval)>nnz" begin - # this should create a sparse matrix with length(nzval)>nnz - A = SparseMatrixCSC(Complex{BigInt}[1+im 2+2im]')'[1:1, 2:2] - # ...ensure it does! If necessary, the test needs to be updated to use - # another mechanism to create a suitable A. - resize!(nonzeros(A), 2) - @assert length(nonzeros(A)) > nnz(A) - @test -A == fill(-2-2im, 1, 1) - @test conj(A) == fill(2-2im, 1, 1) - conj!(A) - @test A == fill(2-2im, 1, 1) -end - @testset "issue #31453" for T in [UInt8, Int8, UInt16, Int16, UInt32, Int32] i = Int[1, 2] j = Int[2, 1] @@ -2945,14 +2906,8 @@ end @test_throws ArgumentError SparseMatrixCSC(10, 3, [1,2,1,2], Int[], Float64[]) # rowwal (and nzval) short @test_throws ArgumentError SparseMatrixCSC(10, 3, [1,2,2,4], [1,2], Float64[]) - # nzval short - @test SparseMatrixCSC(10, 3, [1,2,2,4], [1,2,3], Float64[]) !== nothing - # length(rowval) >= typemax - @test_throws ArgumentError SparseMatrixCSC(5, 1, Int8[1,2], fill(Int8(1),127), Int[1,2,3]) - @test SparseMatrixCSC{Int,Int8}(5, 1, Int8[1,2], fill(Int8(1),127), Int[1,2,3]) != 0 # length(nzval) >= typemax - @test_throws ArgumentError SparseMatrixCSC(5, 1, Int8[1,2], Int8[1], fill(7, 127)) - @test SparseMatrixCSC{Int,Int8}(5, 1, Int8[1,2], Int8[1], fill(7, 127)) != 0 + @test_throws ArgumentError SparseMatrixCSC(5, 1, Int8[1,2], fill(Int8(1), 127), fill(7, 127)) # length(I) >= typemax @test_throws ArgumentError sparse(UInt8.(1:255), fill(UInt8(1), 255), fill(1, 255)) diff --git a/stdlib/SparseArrays/test/sparsevector.jl b/stdlib/SparseArrays/test/sparsevector.jl index 6d29633bb04e0..03aed95a43a83 100644 --- a/stdlib/SparseArrays/test/sparsevector.jl +++ b/stdlib/SparseArrays/test/sparsevector.jl @@ -1419,14 +1419,14 @@ end # test entry points to similar with entry type, index type, and non-Dims shape specification @test similar(A, Float32, Int8, 6, 6) == similar(A, Float32, Int8, (6, 6)) @test similar(A, Float32, Int8, 6) == similar(A, Float32, Int8, (6,)) - # test similar with Dims{2} specification (preserves storage space only, not stored-entry structure) + # test similar with Dims{2} specification (preserves allocated storage space only, not stored-entry structure) simA = similar(A, (6,6)) @test typeof(simA) == SparseMatrixCSC{eltype(nonzeros(A)),eltype(nonzeroinds(A))} @test size(simA) == (6,6) @test getcolptr(simA) == fill(1, 6+1) @test length(rowvals(simA)) == 0 @test length(nonzeros(simA)) == 0 - # test similar with entry type and Dims{2} specification (preserves storage space only) + # test similar with entry type and Dims{2} specification (preserves allocated storage space only) simA = similar(A, Float32, (6,6)) @test typeof(simA) == SparseMatrixCSC{Float32,eltype(nonzeroinds(A))} @test size(simA) == (6,6) diff --git a/stdlib/SuiteSparse/src/cholmod.jl b/stdlib/SuiteSparse/src/cholmod.jl index af0dff10ed473..d8dafa853702d 100644 --- a/stdlib/SuiteSparse/src/cholmod.jl +++ b/stdlib/SuiteSparse/src/cholmod.jl @@ -1094,60 +1094,46 @@ function Vector{T}(D::Dense{T}) where T end Vector(D::Dense{T}) where {T} = Vector{T}(D) +function _extract_args(s) + return (s.nrow, s.ncol, increment(unsafe_wrap(Array, s.p, (s.ncol + 1,), own = false)), + increment(unsafe_wrap(Array, s.i, (s.nzmax,), own = false)), + copy(unsafe_wrap(Array, s.x, (s.nzmax,), own = false))) +end + +# Trim extra elements in rowval and nzval left around sometimes by CHOLMOD rutines +function _trim_nz_builder!(m, n, colptr, rowval, nzval) + l = colptr[end] - 1 + resize!(rowval, l) + resize!(nzval, l) + SparseMatrixCSC(m, n, colptr, rowval, nzval) +end + function SparseMatrixCSC{Tv,SuiteSparse_long}(A::Sparse{Tv}) where Tv s = unsafe_load(pointer(A)) if s.stype != 0 throw(ArgumentError("matrix has stype != 0. Convert to matrix " * "with stype == 0 before converting to SparseMatrixCSC")) end - - B = SparseMatrixCSC(s.nrow, s.ncol, - increment(unsafe_wrap(Array, s.p, (s.ncol + 1,), own = false)), - increment(unsafe_wrap(Array, s.i, (s.nzmax,), own = false)), - copy(unsafe_wrap(Array, s.x, (s.nzmax,), own = false))) - - if s.sorted == 0 - return SparseArrays.sortSparseMatrixCSC!(B) - else - return B - end + args = _extract_args(s) + s.sorted == 0 && _sort_buffers!(args...); + return _trim_nz_builder!(args...) end function Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}}(A::Sparse{Float64}) s = unsafe_load(pointer(A)) - if !issymmetric(A) - throw(ArgumentError("matrix is not symmetric")) - end - - B = Symmetric(SparseMatrixCSC(s.nrow, s.ncol, - increment(unsafe_wrap(Array, s.p, (s.ncol + 1,), own = false)), - increment(unsafe_wrap(Array, s.i, (s.nzmax,), own = false)), - copy(unsafe_wrap(Array, s.x, (s.nzmax,), own = false))), s.stype > 0 ? :U : :L) - - if s.sorted == 0 - return SparseArrays.sortSparseMatrixCSC!(B.data) - else - return B - end + issymmetric(A) || throw(ArgumentError("matrix is not symmetric")) + args = _extract_args(s) + s.sorted == 0 && _sort_buffers!(args...) + Symmetric(_trim_nz_builder!(args...), s.stype > 0 ? :U : :L) end convert(T::Type{Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}}}, A::Sparse{Float64}) = T(A) function Hermitian{Tv,SparseMatrixCSC{Tv,SuiteSparse_long}}(A::Sparse{Tv}) where Tv<:VTypes s = unsafe_load(pointer(A)) - if !ishermitian(A) - throw(ArgumentError("matrix is not Hermitian")) - end - - B = Hermitian(SparseMatrixCSC(s.nrow, s.ncol, - increment(unsafe_wrap(Array, s.p, (s.ncol + 1,), own = false)), - increment(unsafe_wrap(Array, s.i, (s.nzmax,), own = false)), - copy(unsafe_wrap(Array, s.x, (s.nzmax,), own = false))), s.stype > 0 ? :U : :L) - - if s.sorted == 0 - return SparseArrays.sortSparseMatrixCSC!(B.data) - else - return B - end + ishermitian(A) || throw(ArgumentError("matrix is not Hermitian")) + args = _extract_args(s) + s.sorted == 0 && _sort_buffers!(args...) + Hermitian(_trim_nz_builder!(args...), s.stype > 0 ? :U : :L) end convert(T::Type{Hermitian{Tv,SparseMatrixCSC{Tv,SuiteSparse_long}}}, A::Sparse{Tv}) where {Tv<:VTypes} = T(A) @@ -1156,14 +1142,14 @@ function sparse(A::Sparse{Float64}) # Notice! Cannot be type stable because of s if s.stype == 0 return SparseMatrixCSC{Float64,SuiteSparse_long}(A) end - return Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}}(A) + Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}}(A) end function sparse(A::Sparse{ComplexF64}) # Notice! Cannot be type stable because of stype s = unsafe_load(pointer(A)) if s.stype == 0 return SparseMatrixCSC{ComplexF64,SuiteSparse_long}(A) end - return Hermitian{ComplexF64,SparseMatrixCSC{ComplexF64,SuiteSparse_long}}(A) + Hermitian{ComplexF64,SparseMatrixCSC{ComplexF64,SuiteSparse_long}}(A) end function sparse(F::Factor) s = unsafe_load(pointer(F)) @@ -1175,7 +1161,8 @@ function sparse(F::Factor) L, d = getLd!(LD) A = (L * Diagonal(d)) * L' end - SparseArrays.sortSparseMatrixCSC!(A) + # no need to sort buffers here, as A isa SparseMatrixCSC + # and it is taken care in sparse p = get_perm(F) if p != [1:s.n;] pinv = Vector{Int}(undef, length(p)) @@ -1968,4 +1955,57 @@ end (*)(A::SparseVecOrMat{Float64,Ti}, B::Hermitian{Float64,SparseMatrixCSC{Float64,Ti}}) where {Ti} = sparse(Sparse(A)*Sparse(B)) +# Sort all the indices in each column for the construction of a CSC sparse matrix +function _sort_buffers!(m, n, colptr::Vector{Ti}, rowval::Vector{Ti}, nzval::Vector{Tv}) where {Ti <: Integer, Tv} + index = Base.zeros(Ti, m) + row = Base.zeros(Ti, m) + val = Base.zeros(Tv, m) + + perm = Base.Perm(Base.ord(isless, identity, false, Base.Order.Forward), row) + + @inbounds for i = 1:n + nzr = colptr[i]:colptr[i+1]-1 + numrows = length(nzr) + if numrows <= 1 + continue + elseif numrows == 2 + f = first(nzr) + s = f+1 + if rowval[f] > rowval[s] + rowval[f], rowval[s] = rowval[s], rowval[f] + nzval[f], nzval[s] = nzval[s], nzval[f] + end + continue + end + resize!(row, numrows) + resize!(index, numrows) + + jj = 1 + @simd for j = nzr + row[jj] = rowval[j] + val[jj] = nzval[j] + jj += 1 + end + + if numrows <= 16 + alg = Base.Sort.InsertionSort + else + alg = Base.Sort.QuickSort + end + + # Reset permutation + index .= 1:numrows + + Base.sort!(index, alg, perm) + + jj = 1 + @simd for j = nzr + rowval[j] = row[index[jj]] + nzval[j] = val[index[jj]] + jj += 1 + end + end +end + + end #module diff --git a/stdlib/SuiteSparse/test/cholmod.jl b/stdlib/SuiteSparse/test/cholmod.jl index 501be8860436e..edc31f98c5b96 100644 --- a/stdlib/SuiteSparse/test/cholmod.jl +++ b/stdlib/SuiteSparse/test/cholmod.jl @@ -795,12 +795,13 @@ end end end -@testset "Check inputs to Sparse. Related to #20024" for A_ in ( - SparseMatrixCSC(2, 2, [1, 2, 3], CHOLMOD.SuiteSparse_long[1,2], Float64[]), - SparseMatrixCSC(2, 2, [1, 2, 3], CHOLMOD.SuiteSparse_long[1,2], Float64[1.0])) - args = (size(A_)..., getcolptr(A_) .- 1, rowvals(A_) .- 1, nonzeros(A_)) - @test_throws ArgumentError CHOLMOD.Sparse(args...) - @test_throws ArgumentError CHOLMOD.Sparse(A_) +@testset "Check inputs to Sparse. Related to #20024" for t_ in ( + (2, 2, [1, 2], CHOLMOD.SuiteSparse_long[], Float64[]), + (2, 2, [1, 2, 3], CHOLMOD.SuiteSparse_long[1], Float64[]), + (2, 2, [1, 2, 3], CHOLMOD.SuiteSparse_long[], Float64[1.0]), + (2, 2, [1, 2, 3], CHOLMOD.SuiteSparse_long[1], Float64[1.0])) + @test_throws ArgumentError SparseMatrixCSC(t_...) + @test_throws ArgumentError CHOLMOD.Sparse(t_[1], t_[2], t_[3] .- 1, t_[4] .- 1, t_[5]) end @testset "sparse right multiplication of Symmetric and Hermitian matrices #21431" begin