Add docs on task-specific buffering using multithreading

[IanButterworth]

It's common to see people using threadid()-based buffers, for instance in a toy sum example

function sum_multi(a)
     buffers = zeros(eltype(a), Threads.nthreads())
     Threads.@threads for i in eachindex(a)
         buffers[Threads.threadid()] += a[i]
     end
     s = 0
     for b in buffers
         s += b
     end
     return s
 end

but after task migration #40715 (which is undocumented AFAICT) this practice is not race-safe.

This is an attempt to document some best practice.

I expected the task_local_storage api to be the right answer here, but I couldn't figure out how to reduce it after the @threads loop returned.

Help appreciated in making this correct

[vtjnash]

Thanks for tackling this

[vtjnash]

Also because nthreads is soft-deprecated (no warning, but may return incorrect answers now)

[IanButterworth]

I didn't add a note on this. I don't understand whether the docs should or shouldn't recommend nthreads()

[torfjelde]

I'm not an expert on multi-threaded parallelism, so I might just be wrong.

But this doesn't seem like it'll be particularly performant as the values of the WeakKeyDict won't be contiguous in memory, and so the reduction at the needs to gather these from likely different parts of the heap => slow reduction. In fact, when I tried running something like this on an example of mine it was incredibly slow; probably partially due to what I just mentioned, and partially due to the slowness of WeakKeyDict (I believe you can also use an IdDict which at least should be faster that WeakKeyDict).

An alternative is Atomic as mentioned below, or, even better, you can do a Vector{Atomic{T}} of some buffer_length. In the latter scenario you can then just randomly pick an index for each Task, and act on the corresponding Atomic{T} atomically, e.g. atomic_add!(buffer[rand(1:buffer_size)], i). If the work within each of the tasks is fairly uniform, then this random picking of index to add to should result in fairly even congestion for the different buffers.

I'm sure there are much, much better ways of doing this, but the above is a quick-and-easy way to implement a much faster version of this kind of map-reduce approach using a buffer.

[Moelf]

just to move some of the comments here to provide some context as to why this PR is very much needed.

Problem with task local storage

• it's not type stable out of the box
• if what you're buffering involves large RAM or slow I/O and the result is used across task boundary, tasks may over spawn and overlap with each other (see this for example)

Our ecosystem contains many silently incorrect code now without this

• https://github.com/mkarikom/AdditiveCellCom.jl/blob/v0.3.11/src/features.jl#L71
• https://github.com/apache/arrow-julia/blob/6fe4ec0058657a5fe9a18b226a54a03a56d5bf81/src/Arrow.jl#L103-L106
• https://github.com/JuliaData/Parsers.jl/blob/41c1bdd31323ec937680b25f723c8a91b190ba7e/src/floats.jl#L11-L21
• https://github.com/pagnani/ArDCA.jl/blob/v0.5.1/src/dca.jl#L20
• https://github.com/ivanslapnicar/Arrowhead.jl/blob/v1.0.0/src/arrowhead4.jl#L491
• https://github.com/oxinabox/AutoPreallocation.jl/blob/v0.3.1/src/preallocate.jl#L13 (?)
• https://github.com/ITensor/ITensorParallel.jl/blob/v0.1.0/src/references/ExtendedProjMPO.jl#L95
• https://github.com/SciML/OrdinaryDiffEq.jl/blob/b7bb73f46669e6f35c33baf2b4a1a1dc9892a091/src/perform_step/extrapolation_perform_step.jl#L52-L55

The easiest way for these to be fixed is probably to provide a replacement for what used to work, thus this PR is pretty important.

[vtjnash]

At least most of those had the decency to avoid having inbounds annotations, which can quickly (and often) turns mild wrong code into drastic wrong answers

[vtjnash]

Those seem more impetus to continue looking into #48543 perhaps?

[KristofferC]

I think many of those mimiced the implementation of the old default rng.

[IanButterworth]

Perhaps this is ready for final review? @MasonProtter perhaps?

[MasonProtter]

Sorry for suggested that code and now I'm suggesting changes, but I guess it would be nice and concise here to just re-use the single threaded version of the function and run it on each of the spawned chunks, and also use it for reducing the returned data?

[IanButterworth]

@MasonProtter I absorbed your suggestions. Can you give this a last review please. Thanks

[MasonProtter]

I like it!