ExtendableSparse.jl

ExtendableSparse.jl copied to clipboard

extendable.jl provides several mul! methods and a Base.:+ override returning a sparse result when adding e.g. a sparse matrix to an ExtendableSparseMatrix. However, if an ExtendableSparseMatrix is multiplied by a scalar or added to another ExtendableSparseMatrix, i.e. running

using SparseArrays, LinearAlgebra, ExtendableSparse
A = ExtendableSparseMatrix(sprand(100,100,.5));
K = ExtendableSparseMatrix(sprand(100,100,.5));
A + 0.4 * K

the default Broadcasting machinery reverts to the Broadcast.DefaultArrayStyle type. This means that a dense destination container of type Array{Float64,2} is allocated and the existing SparseArrays machinery for sparse array maths is skipped completely leading to a significant performance and memory hit.

If there are no plans to implement custom broadcasting behaviour for the ExtendableSparseMatrix type, couldn't we provide new methods such as

function Base.:*(s::Tv, ext::ExtendableSparseMatrix{Tv,Ti}) where {Tv<:Number,Ti<:Integer}
       @inbounds flush!(ext)
       return s * ext.cscmatrix
end

function Base.:+(ext1::ExtendableSparseMatrix{Tv,Ti}, ext2::ExtendableSparseMatrix{Tv,Ti}) where {Tv<:Number,Ti<:Integer}
       @inbounds flush!(ext1)
       @inbounds flush!(ext2)
       return ext1.cscmatrix + ext2.cscmatrix
end

to guard against such unpleasant surprises when handling ExtendableSparseMatrix?