Thread-safety of SparseArray

The setindex! op of SparseArray is not thread-safe. There is a MWE:

function foo(x::AbstractArray)
    y = similar(x)
    Threads.@threads for i = 1:length(x)
        y[i] = cos(x[i])
    end
end
using SparseArrays
x = sprandn(100, 100, 0.2)
foo(x)

This function would crash dua to thread-safety and @inbounds inner sparsematrix.

What’s the best pratice to write an generic function that support SparseArrays?

Hi!

First of all, in this exact MWE the output is better as a dense array, since cos(0) == 1.0. I’ll use sin for illustrative purposes below.

Speaking of writing generic functions, since sparse arrays behave a lot differently than dense ones, you have only two options — either use existing generic interface like map or write several specializations of your function for different array types.

In your case I would recommend using tmap from OhMyThreads.jl and specialize it something like that (assuming the output is not a dense array and keeps the same sparsity pattern):

foo(x::SparseMatrixCSC) = SparseMatrixCSC(x.m, x.n, x.colptr, x.rowval, tmap(sin, x.nzval))

Calling setindex! on a SparseArray is also extremely slow if you are changing the sparsity pattern (inserting new nonzeros). You may need a different data structure. What are you actually trying to do?

In your example code, you start with a sparse array but then assign every element to a nonzero value, which is obviously not realistic. We need a more realistic context to suggest a better approach.

I have some generic function that use multithreads with AbstractArray, but crash happends if user call them with their sparse inputs. I’m worried that this pattern of code is not so safe, maybe I should just normalize the inputs to DenseArray, or not calling similar(x) to creat output, or just check if the inputs satisfy thread-safety. I’m not sure how to write the safe and generic code in this context.

Strictly speaking, any time you use an abstract type, like DenseArray, it’s possible that some package will define a non-threadsafe subtype. (However, the contract of DenseArray, which specifies a particular type of memory layout, seems like it could make this unlikely. You could also use StridedArray, which includes strided subarrays of DenseArray.)

However, in your case it sounds like you may be better off just creating your output array with Array, since you require a performant, thread-safe setindex! and it sounds like you may have dense outputs.

Thread-safety of inserts into irregular sparse arrays in general is not really feasible and will always exhibit poor performance. Each insert to a sparse array will involve an allocation and data reshuffling. You probably want to detect sparsity (either with the type tree or using something like ArrayInterface.jl) and then create task local accumulators before combining them into a “static” SparseMatrixCSC.

Even in sparse libraries that are highly parallel you typically do this by creating a task local SparseMatrixCSC and then combine them.