diff --git a/Project.toml b/Project.toml index 6edec0d..8b106a4 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ExtendableSparse" uuid = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3" -version = "2.0.1" +version = "2.1" authors = ["Juergen Fuhrmann ", "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" diff --git a/src/ExtendableSparse.jl b/src/ExtendableSparse.jl index 52b173c..6857d96 100644 --- a/src/ExtendableSparse.jl +++ b/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 diff --git a/src/abstractextendablesparsematrixcsc.jl b/src/abstractextendablesparsematrixcsc.jl index 2bcf273..41ca7ee 100644 --- a/src/abstractextendablesparsematrixcsc.jl +++ b/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) diff --git a/src/genericextendablesparsematrixcsc.jl b/src/genericextendablesparsematrixcsc.jl index 5d4ce96..52718d7 100644 --- a/src/genericextendablesparsematrixcsc.jl +++ b/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) diff --git a/src/genericmtextendablesparsematrixcsc.jl b/src/genericmtextendablesparsematrixcsc.jl index f376783..a97658d 100644 --- a/src/genericmtextendablesparsematrixcsc.jl +++ b/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) diff --git a/src/sparsematrixcsc.jl b/src/sparsematrixcsc.jl index 8065610..16f161b 100644 --- a/src/sparsematrixcsc.jl +++ b/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) diff --git a/src/sparsematrixdict.jl b/src/sparsematrixdict.jl index 6d8c5d6..00f1c68 100644 --- a/src/sparsematrixdict.jl +++ b/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 diff --git a/src/sparsematrixdilnkc.jl b/src/sparsematrixdilnkc.jl index 081bd29..2c7ae15 100644 --- a/src/sparsematrixdilnkc.jl +++ b/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 diff --git a/src/sparsematrixlnk.jl b/src/sparsematrixlnk.jl index 83d897f..ff5cc39 100644 --- a/src/sparsematrixlnk.jl +++ b/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) diff --git a/test/alltests.jl b/test/alltests.jl index 3a784c8..0204826 100644 --- a/test/alltests.jl +++ b/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