Advice for building sparsematrix

General Usage Performance
1

Hi,

I feel that my problem must be shared by some of you.

I need to build the sparse matrix associated to a 2d Finite Element method. Following the advice of Viral, I built 3 lists I,J,values containing the indices of nonzeros elements. So each element of I or J contains the linear index associated to a tuple (i,j) because my problem is 2d. However, this is very inconvenient.

Basically, I need to add values to a sparse matrix and this is slow.

Hence, I am looking for a convenient way to build the indices I,J,V so that I can call sparse(I,J,V) in the end.

By convenient way, I mean that I want to operations +,-,*,... to be defined somehow.

For example, would it be possible to have a sparse matrix defined as

struct spmat
I::Int64[]
J::Int64[]
val::Float64[]
end

and the plus operation would just push to I,J,val. Am I re-inventing the wheel?

Thank you for your suggestions.

2

In JuAFEM we have an assembler type (JuAFEM.jl/assembler.jl at 2d08f55b4410fc3415b6be2c53fc1be3251f3714 · KristofferC/JuAFEM.jl · GitHub). It does not support + per se since we call assemble! but might be useful anyway?. If you really want the + interface, here is a modifed version:

using SparseArrays

struct Assembler{T}
    I::Vector{Int}
    J::Vector{Int}
    V::Vector{T}
    Assembler(N = 0) = Assembler{Float64}(N)
    function Assembler{T}(N = 0) where T
        I, J, V = Int[], Int[], T[]
        sizehint!(I, N)
        sizehint!(J, N)
        sizehint!(V, N)
        return new{T}(I, J, V)
    end
end
function Base.setindex!(a::Assembler, v, i::Integer, j::Integer)
    push!(a.I, i)
    push!(a.J, j)
    push!(a.V, v)
    return v
end
# dummy method to make + work
function Base.getindex(a::Assembler{T}, i::Integer, j::Integer) where T
    return zero(T)
end
finalize(a::Assembler) = sparse(a.I, a.J, a.V)

with example usage:

julia> a = Assembler();

julia> a[1, 2] += 12;

julia> a[1, 2] += 12;

julia> a[2, 2] += 5;

julia> finalize(a)
2Ă—2 SparseMatrixCSC{Float64,Int64} with 2 stored entries:
  [1, 2]  =  24.0
  [2, 2]  =  5.0
2 Likes
3

If you “FEM” guys have to do something like this, it means a more abstract package would be useful…

Thank you!!

4

Also sparse sums (see combine keyword) for you:

julia> i = [1,2,2]; j = [1,2,2]; v=[1,1,1]; Array(sparse(i,j,v))
2Ă—2 Array{Int64,2}:
 1  0
 0  2

This might be enough for your purposes?

5

Yes but that’s not the question. Let us say that we have the following 1D linear operator


function rhs(out,x,V)
   for ii in eachindex(x)
       fp = 0.
       fm = 0.
       if ii>1 and ii<length(x)
           fl = V[ii-1] * x[ii-1] - V[ii]*x[ii]
           fp = V[ii+1] * x[ii+1] - V[ii]*x[ii]
       end
       out[ii] = 0.5*fp + 0.7*fl
   end
end

Can we build its sparse matrix? A very slow way is the following

function rhs_sp(out,x,V)
    A = spzeros(length(x),length(x))
   for ii in eachindex(x)
       fp = 0.
       fm = 0.
       if ii>1 and ii<length(x)
           A[ii,ii]   -= V[ii] *0.7
           A[ii,ii-1] +=  V[ii-1] *0.7

           A[ii,ii+1] += V[ii+1] *0.5
           A[ii,ii] -=  V[ii] *0.5
       end
       
   end
end

@fredrikekre suggestion is much better but still need a polish for easy use

1 Like
6

That sounded like the question to me too, and if I’m not mistaken is pretty much what Fredrik’s solution does (adds a new entry for each addition). If you need to obtain previously stored values, Fredrik’s solution won’t work (without first finalizing the matrix), and I think your problem becomes a bit more complicated.

Do you know beforehand what the structure of the non-zero entries will be? For example, in your last post, you could easily exploit the diagonal structure.

What’s the dimension of the matrix, and number of non-zeros?

One option would be a dictionary from an index tuple (probably combined into an Int for better performance) to its value, but that might not be performant enough.

7

Sorry, I was a bit quick.

The dimensions of the matrix is known beforehand but not the sparsity pattern.

8

To build a general sparse matrix, all the suggestions boil down to just pushing to I, J, V and then call sparse in the end. If you want to write some fancy interface on top of that or just do the pushing inside the function is just a matter of taste.

2 Likes
9

We all agree here. I guess we also all agree when I am saying that it is relatively easy to code the above function rhs in the case of FEM or FVM. What I am looking for is a generic way to transform rhs so that it builds I,J,V.

10

Off the top of my head, I see a few options.

If you only need to support addition and subtraction, and don’t mind the repeated indices (i.e. won’t be updating existing indices too much), just create I, J, V naively like in @mauro3 and @fredrikekre’s examples, and let sparse add them all up for you.

If you need to support other operations, like multiplication, or want to avoid the repeated indices, the easiest solution that comes to mind is to use a dictionary to map (row,col) indices to the corresponding index in I, J, V. That seems pretty trivial to code, but performance might suffer a bit since there will be many dictionary lookups, and there will be additional (temporary) storage for the dictionary.

Another option would be to call sparse with a custom combine function instead of the default +. That sounds a bit more complicated since you’d probably need to encode the operation within the value itself then, in which case you wouldn’t get the output in the format you want. But there’s probably a way around that. Here’s an example to give you an idea of what I mean; we start with the value 3, add 5, then multiply with 7:

julia> sparse([1,1,1], [1,1,1], [(3,~),(5,+),(7,*)], 1, 1, ((val1,fun1),(val2,fun2)) -> (fun2(val1, val2), fun1))
1Ă—1 SparseMatrixCSC{Tuple{Int64,Function},Int64} with 1 stored entry:
[1, 1] = (56, ~)
11

@fredrikekre Perhaps something like Assembler can be added in SparseArrays? Or a package of its own? It’s easier to build the sparse pattern item by item using a[i,j] = v, instead of making three push! calls for each new index that I encounter.