Fix allocations in Base.:+ and add Base.sum method

Project.toml

@@ -1,6 +1,6 @@
 name = "ExtendableSparse"
 uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
-version = "2.0.1"
+version = "2.1"
 authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>", "Daniel Runge"]
 
 [deps]
@@ -35,20 +35,20 @@ ILUZero = "0.2"
 IncompleteLU = "^0.2.1"
 InteractiveUtils = "1.11.0"
 IterativeSolvers = "0.9"
-LinearAlgebra = "1.9"
+LinearAlgebra = "1.10"
 LinearSolve = "2.36.0, 3.7.1"
 Metis = "1"
 MultiFloats = "1, 2"
 OhMyThreads = "0.6, 0.7, 0.8"
-Printf = "1.9"
-Random = "1.9"
+Printf = "1.10"
+Random = "1.10"
 RecursiveFactorization = "0.2"
 SciMLPublic = "1.0.1"
-SparseArrays = "1.9"
+SparseArrays = "1.10"
 Sparspak = "0.3.6"
 StaticArrays = "1"
-Test = "1.9"
-julia = "1.9"
+Test = "1.10"
+julia = "1.10"
 
 [extras]
 AMGCLWrap = "4f76b812-4ba5-496d-b042-d70715554288"

src/ExtendableSparse.jl

@@ -8,7 +8,7 @@ module ExtendableSparse
 using DocStringExtensions: DocStringExtensions, SIGNATURES, TYPEDEF, TYPEDFIELDS, TYPEDSIGNATURES
 using ILUZero: ILUZero
 using LinearAlgebra: LinearAlgebra, Diagonal, Hermitian, Symmetric, Tridiagonal, mul!, ldiv!
-using SparseArrays: SparseArrays, AbstractSparseMatrix, AbstractSparseMatrixCSC, SparseMatrixCSC
+using SparseArrays: SparseArrays, AbstractSparseMatrix, AbstractSparseMatrixCSC, SparseMatrixCSC, indtype
 using SparseArrays: dropzeros!, findnz, nzrange, sparse, spzeros, rowvals, getcolptr, nonzeros, nnz
 using Sparspak: sparspaklu
 using SciMLPublic: @public

src/abstractextendablesparsematrixcsc.jl

@@ -106,6 +106,12 @@ Return element type.
 """
 Base.eltype(::AbstractExtendableSparseMatrixCSC{Tv, Ti}) where {Tv, Ti} = Tv
 
+"""
+    SparseArrays.indtype(A::AbstractExtendableSparseMatrixCSC{Tv, Ti})
+
+Return element type.
+"""
+SparseArrays.indtype(A::AbstractExtendableSparseMatrixCSC{Tv, Ti}) where {Tv, Ti} = Ti
 
 """
     SparseArrays.SparseMatrixCSC(A::AbstractExtendableSparseMatrixCSC)
@@ -360,6 +366,37 @@ function Base.:-(csc::SparseMatrixCSC, ext::AbstractExtendableSparseMatrixCSC)
     return csc - sparse(ext)
 end
 
+"""
+    Base.sum(M::AbstractVector{TM}) where TM<:AbstractExtendableSparseMatrixCSC
+
+Efficient sum of ExtendableSparse matrices.
+"""
+function Base.sum(M::AbstractVector{TM}) where {TM <: AbstractExtendableSparseMatrixCSC}
+    Ti = promote_type(indtype.(M)...)
+    Tv = promote_type(eltype.(M)...)
+    return _sum(M, Tv, Ti)
+end
+
+function _sum(M, ::Type{Tv}, ::Type{Ti}) where {Tv, Ti}
+    l = sum(nnz, M)
+    I = Vector{Ti}(undef, l)
+    J = Vector{Ti}(undef, l)
+    V = Vector{Tv}(undef, l)
+
+    i = 1
+    @time for m in M
+        (; colptr, nzval, rowval) = sparse(m)
+        for icsc in 1:(length(colptr) - 1)
+            for j in colptr[icsc]:(colptr[icsc + 1] - 1)
+                I[i] = rowval[j]
+                J[i] = icsc
+                V[i] = nzval[j]
+                i = i + 1
+            end
+        end
+    end
+    return SparseArrays.sparse!(I, J, V)
+end
 
 """
 $(TYPEDSIGNATURES)

src/genericextendablesparsematrixcsc.jl

@@ -137,12 +137,12 @@ end
 """
 function rawupdateindex!(
         ext::GenericExtendableSparseMatrixCSC,
-        op,
+        op::Op,
         v,
         i,
         j,
         part = 1
-    )
+    ) where {Op}
     k = findindex(ext.cscmatrix, i, j)
     return if k > 0
         ext.cscmatrix.nzval[k] = op(ext.cscmatrix.nzval[k], v)
@@ -151,16 +151,17 @@ function rawupdateindex!(
     end
 end
 
+
 """
     $(TYPEDSIGNATURES)
 """
 function updateindex!(
         ext::GenericExtendableSparseMatrixCSC,
-        op,
+        op::Op,
         v,
         i,
         j
-    )
+    ) where {Op}
     k = findindex(ext.cscmatrix, i, j)
     return if k > 0
         ext.cscmatrix.nzval[k] = op(ext.cscmatrix.nzval[k], v)

src/genericmtextendablesparsematrixcsc.jl

@@ -148,12 +148,12 @@ nnznew(ext::GenericMTExtendableSparseMatrixCSC) = sum(nnz, ext.xmatrices)
 """
 function rawupdateindex!(
         ext::GenericMTExtendableSparseMatrixCSC,
-        op,
+        op::Op,
         v,
         i,
         j,
         tid = 1
-    )
+    ) where {Op}
     k = findindex(ext.cscmatrix, i, j)
     return if k > 0
         ext.cscmatrix.nzval[k] = op(ext.cscmatrix.nzval[k], v)
@@ -168,12 +168,12 @@ end
 """
 function updateindex!(
         ext::GenericMTExtendableSparseMatrixCSC,
-        op,
+        op::Op,
         v,
         i,
         j,
         tid = 1
-    )
+    ) where {Op}
     k = findindex(ext.cscmatrix, i, j)
     return if k > 0
         ext.cscmatrix.nzval[k] = op(ext.cscmatrix.nzval[k], v)

src/sparsematrixcsc.jl

@@ -48,7 +48,7 @@ A
 ```
 
 """
-function updateindex!(csc::SparseMatrixCSC{Tv, Ti}, op, v, i, j) where {Tv, Ti <: Integer}
+function updateindex!(csc::SparseMatrixCSC{Tv, Ti}, op::Op, v, i, j) where {Tv, Ti <: Integer, Op}
     k = findindex(csc, i, j)
     if k > 0 # update existing value
         csc.nzval[k] = op(csc.nzval[k], v)

src/sparsematrixdict.jl

@@ -85,7 +85,7 @@ end
 """
     $(TYPEDSIGNATURES)
 """
-function rawupdateindex!(m::SparseMatrixDict{Tv, Ti}, op, v, i, j) where {Tv, Ti}
+function rawupdateindex!(m::SparseMatrixDict{Tv, Ti}, op::Op, v, i, j) where {Tv, Ti, Op}
     p = Pair(i, j)
     return m.values[p] = op(get(m.values, p, zero(Tv)), v)
 end
@@ -94,7 +94,7 @@ end
 """
     $(TYPEDSIGNATURES)
 """
-function updateindex!(m::SparseMatrixDict{Tv, Ti}, op, v, i, j) where {Tv, Ti}
+function updateindex!(m::SparseMatrixDict{Tv, Ti}, op::Op, v, i, j) where {Tv, Ti, Op}
     p = Pair(i, j)
     v1 = op(get(m.values, p, zero(Tv)), v)
     if !iszero(v1)
@@ -135,11 +135,7 @@ function SparseArrays.sparse(mat::SparseMatrixDict{Tv, Ti}) where {Tv, Ti}
         V[i] = v
         i = i + 1
     end
-    @static if VERSION >= v"1.10"
-        return SparseArrays.sparse!(I, J, V, size(mat)..., +)
-    else
-        return SparseArrays.sparse(I, J, V, size(mat)..., +)
-    end
+    return SparseArrays.sparse!(I, J, V, size(mat)..., +)
 end
 
 """
@@ -172,11 +168,7 @@ function Base.:+(dictmatrix::SparseMatrixDict{Tv, Ti}, cscmatrix::SparseMatrixCS
         end
 
         @assert l == i - 1
-        @static if VERSION >= v"1.10"
-            return SparseArrays.sparse!(I, J, V, m, n, +)
-        else
-            return SparseArrays.sparse(I, J, V, m, n, +)
-        end
+        return SparseArrays.sparse!(I, J, V, m, n, +)
     end
     return cscmatrix
 end

src/sparsematrixdilnkc.jl

@@ -218,7 +218,7 @@ Update element of the matrix  with operation `op`.
 It assumes that `op(0,0)==0`. If `v` is zero, no new 
 entry is created.
 """
-function updateindex!(lnk::SparseMatrixDILNKC{Tv, Ti}, op, v, i, j) where {Tv, Ti}
+function updateindex!(lnk::SparseMatrixDILNKC{Tv, Ti}, op::Op, v, i, j) where {Tv, Ti, Op}
     k, k0 = findindex(lnk, i, j)
     if k > 0
         lnk.nzval[k] = op(lnk.nzval[k], v)
@@ -238,7 +238,7 @@ Update element of the matrix  with operation `op`.
 It assumes that `op(0,0)==0`. If `v` is zero a new entry
 is created nevertheless.
 """
-function rawupdateindex!(lnk::SparseMatrixDILNKC{Tv, Ti}, op, v, i, j) where {Tv, Ti}
+function rawupdateindex!(lnk::SparseMatrixDILNKC{Tv, Ti}, op::Op, v, i, j) where {Tv, Ti, Op}
     k, k0 = findindex(lnk, i, j)
     if k > 0
         lnk.nzval[k] = op(lnk.nzval[k], v)
@@ -443,11 +443,7 @@ function Base.sum(lnkdictmatrices::Vector{SparseMatrixDILNKC{Tv, Ti}}, cscmatrix
             end
         end
         @assert l == i - 1
-        @static if VERSION >= v"1.10"
-            return SparseArrays.sparse!(I, J, V, m, n, +)
-        else
-            return SparseArrays.sparse(I, J, V, m, n, +)
-        end
+        return SparseArrays.sparse!(I, J, V, m, n, +)
     end
     return cscmatrix
 end

src/sparsematrixlnk.jl

@@ -18,7 +18,7 @@ can be conveniently updated via `push!`.  No copying of existing data is necessa
 
 Via the type aliases [`STExtendableSparseMatrixCSC`](@ref), [`ExtendableSparseMatrixCSC`](@ref),
 and [`ExtendableSparseMatrix`](@ref) this extension is used as default for handling
-scalar assembly.
+sequential assembly.
 
 
 $(TYPEDFIELDS)
@@ -216,7 +216,7 @@ Update element of the matrix  with operation `op`.
 It assumes that `op(0,0)==0`. If `v` is zero, no new 
 entry is created.
 """
-function updateindex!(lnk::SparseMatrixLNK{Tv, Ti}, op, v, i, j) where {Tv, Ti}
+function updateindex!(lnk::SparseMatrixLNK{Tv, Ti}, op::Op, v, i, j) where {Tv, Ti, Op}
     # Set the first  column entry if it was not yet set.
     if lnk.rowval[j] == 0 && !iszero(v)
         lnk.rowval[j] = i
@@ -243,7 +243,7 @@ Update element of the matrix  with operation `op`.
 It assumes that `op(0,0)==0`. If `v` is zero a new entry
 is created nevertheless.
 """
-function rawupdateindex!(lnk::SparseMatrixLNK{Tv, Ti}, op, v, i, j) where {Tv, Ti}
+function rawupdateindex!(lnk::SparseMatrixLNK{Tv, Ti}, op::Op, v, i, j) where {Tv, Ti, Op}
     # Set the first  column entry if it was not yet set.
     if lnk.rowval[j] == 0
         lnk.rowval[j] = i
@@ -278,13 +278,13 @@ SparseArrays.nnz(lnk::SparseMatrixLNK) = lnk.nnz
 
 # Struct holding pair of value and row
 # number, for sorting
-mutable struct ColEntry{Tv, Ti <: Integer}
+struct ColEntry{Tv, Ti <: Integer}
     rowval::Ti
     nzval::Tv
 end
 
 # Comparison method for sorting
-Base.isless(x::ColEntry, y::ColEntry) = (x.rowval < y.rowval)
+Base.isless(x::ColEntry{Tv, Ti}, y::ColEntry{Tv, Ti}) where {Tv, Ti} = (x.rowval < y.rowval)
 
 """
 $(TYPEDSIGNATURES)

test/alltests.jl

@@ -21,13 +21,13 @@ end
     @test ExplicitImports.check_no_implicit_imports(ExtendableSparse) === nothing
     @test ExplicitImports.check_all_explicit_imports_via_owners(ExtendableSparse) === nothing
     @static if VERSION >= v"1.11.0"
-        @test ExplicitImports.check_all_explicit_imports_are_public(ExtendableSparse, ignore = (:AbstractSparseMatrixCSC, :getcolptr)) === nothing
+        @test ExplicitImports.check_all_explicit_imports_are_public(ExtendableSparse, ignore = (:AbstractSparseMatrixCSC, :getcolptr, :indtype)) === nothing
     end
     @test ExplicitImports.check_no_stale_explicit_imports(ExtendableSparse) === nothing
     @test ExplicitImports.check_all_qualified_accesses_via_owners(ExtendableSparse) === nothing
     @test ExplicitImports.check_all_qualified_accesses_are_public(
         ExtendableSparse,
-        ignore = (:AbstractSparseMatrixCSC, :AbstractTriangular, :getcolptr, :Forward, :USE_GPL_LIBS, :_checkbuffers, :print_array, :sparse!)
+        ignore = (:AbstractSparseMatrixCSC, :AbstractTriangular, :getcolptr, :Forward, :USE_GPL_LIBS, :_checkbuffers, :print_array, :sparse!, :indtype)
     ) === nothing
     @test ExplicitImports.check_no_self_qualified_accesses(ExtendableSparse) === nothing
 end