Per thread RNG state in multi-threaded code

[ViralBShah]

On the threads branch, I want to be able to have a differently seeded RNG in each thread. I am using the threads with the @threads macro to parallelize for loops. There doesn't seem to be a concept of a long running thread with global state.

I am guessing I should have a global array of RNG states, and every thread then just does something like rand(MT(threadid())). Is there a way I can do something like this for now? The @threads only applies to the for loop and I would like to do some thread local setup before that?

This is not going to work well with a Cilk-style threading model in general - but I just want something that can work for now.

[ViralBShah]

Cc: @kpamnany @rfourquet

[ViralBShah]

I'd just like to have a bit of setup for the RNGs done, just after the threads are spawned, but before the loop starts executing.

[kpamnany]

function seed_rng()
  rng_state[threadid()] = <whatever>
end

@threads all begin
  seed_rng()
  <etc.>
end

<etc.>

@threads all for <etc.>

A library that has to store state requires thread-awareness to be thread-safe. A thread-aware RNG can be implemented as above. Won't this work for you?

I can push a vectorized and thread-aware Xorshuffle implementation. Or we can create wrappers for the Intel MKL RNGs.

I don't quite see how Cilk-style threading changes anything about this.

[kpamnany]

Oh, I think I understand the issue now. JULIA_NUM_THREADS threads are created when Julia starts--they are long-running threads and can store global and/or per-thread state. The @threads macro can be seen as a way to send the threads a do-this-now command, with an implicit barrier at the end after which thread 1 alone continues on. Hope this makes sense...?

[ViralBShah]

Thanks, this is good for now. I thought that @threads will only work with for loops, which is not the case.

[timholy]

In a work-stealing model, it presumably won't be possible to run computations with a reproducible sequence of random numbers, because what gets used where will depend on when individual threads finish their chunks. But it seems this is a small price to pay, since in situations where this matters people can just use a single thread.

[ivarne]

I'm somewhat confused.

Do you want thread local RNGs when you are passing AbstractRNG objects manually, or do you think this is a feature that the global RNG should implement?

[ViralBShah]

For now I want thread local RNGs - which as discussed above, is easy.

[dluciv]

I vote for this case too.

Simple example below fails:

a = zeros(Threads.nthreads())
Threads.@threads for i = 1:Threads.nthreads()
  a[i] = rand()
end

After executing this, sometimes (very often in fact) all or at least some of a array values are similar.

[stevengj]

It would be nice to make GLOBAL_RNG thread-local by default. The most obvious way to do this is to allocate an array of nthreads() RNGs and use threadid() to grab the right one on each thread.

However, I'm concerned about the overhead of this approach: currently, this seems to give a 10-25% slowdown on my machine in simple benchmarks.

using Random
const RNG = [Random.GLOBAL_RNG]
trand() = rand(@inbounds(RNG[Threads.threadid()]))

If this function is called a bunch of times in a loop or something, it would be nice if RNG[Threads.threadid()] could be fetched only once, but I'm guessing the compiler is not smart enough to do this?

It would be nicer if Julia had better low-level support for thread-local variables, so that you could just stick @threadlocal in front of a global constant and it would make one copy per thread that is accessed transparently.