Setindex! for sparse arrays

Specific Domains
1

The question that follows was raised earlier at e.g. performance-enhancement-for-setindex-for-sparse-array.

I remain confused regarding the use of sparse(Is,Js,As,[n,n]). profview shows that computation time is spend in setindex!.

The code that follows

"""
    genLocStiffMat(element::Element)

Generates local stiffness matrix for single linear triangular element. 
"""
function genLocStiffMat(element::Element)
    v    = SVector(element.e1, element.e2, element.e3)   
    Emat = copy(element.Emat) 
    area = element.area 
    reluctivity = element.physics.reluctivity
    Iloc = SVector{9}(v[i] for j=1:3, i=1:3)
    Jloc = SVector{9}(v[i] for i=1:3, j=1:3) 
    Emat[3,:] .= 0.;  
    Amat = SMatrix{3,3}(area*reluctivity*(transpose(Emat)*Emat))
    Aloc = vec(Amat)
    return Iloc, Jloc, Aloc
end

"""
    genStiffMat(mesh::Mesh)

Generates global stiffness matrix on mesh of linear triangular elements. 
"""
function genStiffMat(mesh::Mesh)
 
    #..recover number of elements  
    nnodes      = mesh.nnodes 
    nelems      = mesh.nelems 
    dofPerElem  = mesh.dofPerElem
    dofPerElem2 = dofPerElem^2

    #..set range vector 
    irng = SVector{9}(1:9)

    #..preallocate the memory for local matrix contributions 
    Avalues = zeros(Float64,dofPerElem2*nelems)
    I = zeros(Int64,length(Avalues))
    J = zeros(Int64,length(Avalues))

    #..loop over number of elements..
    for element in mesh.Elements 
        Iloc, Jloc, Aloc = genLocStiffMat(element) 
        I[irng] .= Iloc 
        J[irng] .= Jloc 
        Avalues[irng] .= Aloc   
        irng = irng.+dofPerElem2
    end
 
    #..form sparse matrix        
    A = sparse(I,J,Avalues,nnodes,nnodes)
   
    return A; 
end

indeed has a Profvview as shown below.

genStiffMat_profile_view
genStiffMat_profile_view1546×1192 49.3 KB

Thanks.

2

This is only slowing down your code, just use irng = 1:9.

As for your question, constructing a sparse matrix is an expensive operation. You can accelerate it a bit by using instead A = SparseArrays.sparse!(I,J,Avalues,nnodes,nnodes), but it’s not going to make much of a difference.

What would really speed it up is to manually construct colptr, rowval, and nzval, and use the constructor directly: A = SparseMatrixCSC{Float64,Int}(nnodes, nnodes, colptr, rowval, nzval). That’s a huge pain in the ass, though.

1 Like
3

Many thanks for your input. Much appreciated.

I reduced the problem to the following MWE for sparse()

# A MWE for assembly of 1D linear Lagrange shape functions 
nelems = 10; h = 1/nelems; nnodes = nelems+1 
II = zeros(Int64,4*nelems)
JJ = zeros(Int64,4*nelems)
Avalues = zeros(4*nelems)

for k = 1:nelems 
    II[4*(k-1)+1:4*k]      .= [k, k, k+1, k+1]
    JJ[4*(k-1)+1:4*k]      .= [k, k+1, k, k+1]
    Avalues[4*(k-1)+1:4*k] .= [h, -h, -h, h]
end 

@code_warntype sparse(II,JJ,Avalues,nnodes,nnodes,+)

yielding the output with red font. Is this red font output entirely harmless?


MethodInstance for SparseArrays.sparse(::Vector{Int64}, ::Vector{Int64}, ::Vector{Float64}, ::Int64, ::Int64, ::typeof(+))
  from sparse(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector{Tv}, m::Integer, n::Integer, combine) where {Tv, Ti<:Integer} @ SparseArrays ~/.julia/juliaup/julia-1.10.9+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/SparseArrays/src/sparsematrix.jl:1050
Static Parameters
  Tv = Float64
  Ti = Int64
Arguments
  #self#::Core.Const(SparseArrays.sparse)
  I::Vector{Int64}
  J::Vector{Int64}
  V::Vector{Float64}
  m::Int64
  n::Int64
  combine::Core.Const(+)
Locals
  cscnzval::Vector{Float64}
  cscrowval::Vector{Int64}
  klasttouch::Vector{Int64}
  csccolptr::Vector{Int64}
  csrnzval::Vector{Float64}
  csrcolval::Vector{Int64}
  csrrowptr::Vector{Int64}
  coolen::Int64
Body::SparseMatrixCSC{Float64, Int64}
1 ──       Core.NewvarNode(:(cscnzval))
│          Core.NewvarNode(:(cscrowval))
│          Core.NewvarNode(:(klasttouch))
│          Core.NewvarNode(:(csccolptr))
│          Core.NewvarNode(:(csrnzval))
│          Core.NewvarNode(:(csrcolval))
│          Core.NewvarNode(:(csrrowptr))
│          SparseArrays.require_one_based_indexing(I, J, V)
│          (coolen = SparseArrays.length(I))
│    %10 = SparseArrays.length(J)::Int64
│    %11 = (%10 != coolen)::Bool
└───       goto #3 if not %11
2 ──       goto #4
3 ── %14 = SparseArrays.length(V)::Int64
│    %15 = (%14 != coolen)::Bool
└───       goto #5 if not %15
4 ┄─ %17 = SparseArrays.length(I)::Int64
│    %18 = SparseArrays.length(J)::Int64
│    %19 = Base.string("length(I) (=", %17, ") == length(J) (= ", %18, ") == length(V) (= ")::String
│    %20 = SparseArrays.length(V)::Any
│    %21 = Base.string(%20, ")")::Any
│    %22 = SparseArrays.string("the first three arguments' lengths must match, ", %19, %21)::Any
│    %23 = SparseArrays.ArgumentError(%22)::Any
└───       SparseArrays.throw(%23)
5 ┄─ %25 = Base.hastypemax::Core.Const(Base.hastypemax)
│    %26 = (%25)($(Expr(:static_parameter, 2)))::Core.Const(true)
└───       goto #8 if not %26
6 ── %28 = coolen::Int64
│    %29 = SparseArrays.typemax($(Expr(:static_parameter, 2)))::Core.Const(9223372036854775807)
│    %30 = (%28 >= %29)::Bool
└───       goto #8 if not %30
7 ── %32 = Base.string("the index type ", $(Expr(:static_parameter, 2)), " cannot hold ", coolen, " elements; use a larger index type")::Any
│    %33 = SparseArrays.ArgumentError(%32)::Any
└───       SparseArrays.throw(%33)
8 ┄─ %35 = (m == 0)::Bool
└───       goto #10 if not %35
9 ──       goto #13
10 ─ %38 = (n == 0)::Bool
└───       goto #12 if not %38
11 ─       goto #13
12 ─ %41 = (coolen == 0)::Bool
└───       goto #19 if not %41
13 ┄ %43 = (coolen != 0)::Bool
└───       goto #18 if not %43
14 ─ %45 = (n == 0)::Bool
└───       goto #16 if not %45
15 ─ %47 = SparseArrays.ArgumentError("column indices J[k] must satisfy 1 <= J[k] <= n")::Any
│          SparseArrays.throw(%47)
└───       Core.Const(:(goto %54))
16 ┄ %50 = (m == 0)::Bool
└───       goto #18 if not %50
17 ─ %52 = SparseArrays.ArgumentError("row indices I[k] must satisfy 1 <= I[k] <= m")::Any
└───       SparseArrays.throw(%52)
18 ┄ %54 = SparseArrays.one($(Expr(:static_parameter, 2)))::Core.Const(1)
│    %55 = (n + 1)::Int64
│    %56 = SparseArrays.fill(%54, %55)::Vector{Int64}
│    %57 = Core.apply_type(SparseArrays.Vector, $(Expr(:static_parameter, 2)))::Core.Const(Vector{Int64})
│    %58 = (%57)()::Vector{Int64}
│    %59 = Core.apply_type(SparseArrays.Vector, $(Expr(:static_parameter, 1)))::Core.Const(Vector{Float64})
│    %60 = (%59)()::Vector{Float64}
│    %61 = SparseArrays.SparseMatrixCSC(m, n, %56, %58, %60)::SparseMatrixCSC{Float64, Int64}
└───       return %61
19 ─ %63 = Core.apply_type(SparseArrays.Vector, $(Expr(:static_parameter, 2)))::Core.Const(Vector{Int64})
│    %64 = SparseArrays.undef::Core.Const(UndefInitializer())
│    %65 = (m + 1)::Int64
│          (csrrowptr = (%63)(%64, %65))
│    %67 = Core.apply_type(SparseArrays.Vector, $(Expr(:static_parameter, 2)))::Core.Const(Vector{Int64})
│    %68 = SparseArrays.undef::Core.Const(UndefInitializer())
│          (csrcolval = (%67)(%68, coolen))
│    %70 = Core.apply_type(SparseArrays.Vector, $(Expr(:static_parameter, 1)))::Core.Const(Vector{Float64})
│    %71 = SparseArrays.undef::Core.Const(UndefInitializer())
│          (csrnzval = (%70)(%71, coolen))
│    %73 = Core.apply_type(SparseArrays.Vector, $(Expr(:static_parameter, 2)))::Core.Const(Vector{Int64})
│    %74 = SparseArrays.undef::Core.Const(UndefInitializer())
│    %75 = (n + 1)::Int64
│          (csccolptr = (%73)(%74, %75))
│    %77 = Core.apply_type(SparseArrays.Vector, $(Expr(:static_parameter, 2)))::Core.Const(Vector{Int64})
│    %78 = SparseArrays.undef::Core.Const(UndefInitializer())
│          (klasttouch = (%77)(%78, n))
│    %80 = Core.apply_type(SparseArrays.Vector, $(Expr(:static_parameter, 2)))::Core.Const(Vector{Int64})
│          (cscrowval = (%80)())
│    %82 = Core.apply_type(SparseArrays.Vector, $(Expr(:static_parameter, 1)))::Core.Const(Vector{Float64})
│          (cscnzval = (%82)())
│    %84 = SparseArrays.sparse!(I, J, V, m, n, combine, klasttouch, csrrowptr, csrcolval, csrnzval, csccolptr, cscrowval, cscnzval)::SparseMatrixCSC{Float64, Int64}
└───       return %84
4

That’s strange: │ %20 = SparseArrays.length(V)::Any
It definitely shouldn’t happen. I could reproduce it here with Julia 1.11, but not with 1.12. What version of Julia are you using?

5

I have to say I don’t understand why in your case setindex! should have been called at all. You built the sparse matrix from COO triples!

6

The method’s length(V) calls are all in a validation step:

function sparse(I::AbstractVector{Ti}, J::AbstractVector{Ti}, V::AbstractVector{Tv}, m::Integer, n::Integer, combine) where {Tv,Ti<:Integer}
    require_one_based_indexing(I, J, V)
    coolen = length(I)
    if length(J) != coolen || length(V) != coolen
        throw(ArgumentError(string("the first three arguments' lengths must match, ",
              "length(I) (=$(length(I))) == length(J) (= $(length(J))) == length(V) (= ",
              "$(length(V)))")))
    end

The first one in the condition was actually inferred well, it’s the interpolated call that gets inferred as ::Any, which is weird because the other interpolated length calls are inferred just fine.

I wouldn’t expect that to be a problem for performance because that branch wouldn’t be executed in normal practice, though extra compiled code could affect CPU code caches.