dSFMT: use fill_array_* API instead of genrand_* API

[rfourquet]

This is the same as PR #8808 (but on a branch in the main repo) updated to work also with the global RNG, which becomes an instance of MersenneTwister.

Consequently, rand[!] functions involving the global RNG become faster when filling arrays of Float64, but slower when producing Float64 one by one (because of the use of a global Julia object), and, worse, when filling arrays of other types: the problem is that there is a call to rand(T) at each iteration. The solution is to add a rand!{T}(mt::MersenneTwister, A::Array{T}) (which implies having rand(mt, T) methods) and have rand!{T}(A::Array{T}) = rand!(GLOBAL_RNG, A). So should I do this here, or leave that to (a re-write of) PR #8380 (hence keeping the current PR minimal)?

WRT points raised by @ViralBShah in #8808:

• there is already a one line test concerning issue Split Mersenne-Twister states between rand() and rand(N>382) #6573, I'm not sure what other test I could add. Also should the wrapping of dsfmt_genrand_* functions in dSFMT.jl be removed, or left for users who know what they do (a warning could be added in this file)?
• "make sure that Rmath-julia is also consistent with this" : I don't know what to look for.

In order to preserve previous performances when generating unique values or filling small arrays, I had to add a bunch of @inline annotations (I'm not sure that all of them are strictly necessary though), which I find noisy, what do you think?
What about exporting a function globalRNG which returns the GLOBAL_RNG object?

[StefanKarpinski]

I would say do it in two parts, so yes, keep this change minimal.

[rfourquet]

Rebased to make the global RNG const, making the scalar version of rand() less than twice as slow than it was (cf. 0a88117#commitcomment-8344766). I also did minor refactoring (renamed rand_abstract_array_Float64! and moved around few functions ("Low level API for MersenneTwister")).

[ViralBShah]

I think this can be backported without breaking any APIs.

Cc: @JuliaBackports @ivarne

[ViralBShah]

I just realized - we can remove all the _gv_ methods in dSFMT.jl. @rfourquet Would be great if you could do that since you are in the thick of this code.

[rfourquet]

Yes I'll do. How about dsfmt_genrand_* methods? (they could still be used by users who are very careful).

[ViralBShah]

Yes, we should remove those too.

Does the code become any faster if we specialize with methods like these, so that for the versions using the GLOBAL_RNG, one doesn't have to pass the MersenneTwitser object around across multiple function calls?

rand() = rand_close_open(GLOBAL_RNG)
rand(r::MersenneTwister) = rand_close_open(r)

[rfourquet]

I will test, but I had hoped that rand(r::MersenneTwister=GLOBAL_RNG) = rand_close_open(r) would be strictly equivalent (with inlining).
Also, in rand() = rand(GLOBAL_RNG), I thought the slowness due to the global object is not propagated to the subsequent calls to rand(r::MersenneTwister) and rand_close_open(r).

[Jutho]

Comparing the state before and after this merge, I find:
Before:

julia> versioninfo()
Julia Version 0.4.0-dev+1349
Commit da6087e* (2014-10-30 05:12 UTC)
Platform Info:
  System: Darwin (x86_64-apple-darwin14.0.0)
  CPU: Intel(R) Core(TM) i7-3720QM CPU @ 2.60GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

julia> @time A=rand(10,10,10,10,10,10,10,10);
elapsed time: 0.717336153 seconds (800390736 bytes allocated, 2.17% gc time)

julia> @time A=rand(10,10,10,10,10,10,10,10);
elapsed time: 0.712041278 seconds (800000336 bytes allocated)

julia> @time A=randn(10,10,10,10,10,10,10,10);
elapsed time: 1.879174375 seconds (800884100 bytes allocated, 3.66% gc time)

julia> @time A=randn(10,10,10,10,10,10,10,10);
elapsed time: 1.840849301 seconds (800000360 bytes allocated)

After

julia> versioninfo()
Julia Version 0.4.0-dev+1352
Commit 06fb0a3* (2014-10-30 08:19 UTC)
Platform Info:
  System: Darwin (x86_64-apple-darwin14.0.0)
  CPU: Intel(R) Core(TM) i7-3720QM CPU @ 2.60GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

julia> @time A=rand(10,10,10,10,10,10,10,10);
elapsed time: 0.457573658 seconds (801330504 bytes allocated, 3.45% gc time)

julia> @time A=rand(10,10,10,10,10,10,10,10);
elapsed time: 0.410389321 seconds (800000336 bytes allocated)

julia> @time A=randn(10,10,10,10,10,10,10,10);
elapsed time: 2.411182032 seconds (800948036 bytes allocated, 2.82% gc time)

julia> @time A=randn(10,10,10,10,10,10,10,10);
elapsed time: 2.37771637 seconds (800000360 bytes allocated)

So rand got faster, but randn got slower?

[rfourquet]

Yes, as a consequence of using a global Julia object for the global RNG. I will try to restore previous performances as much as I can (#8854 is a first step, and the same kind of optimization as in rand_AbstractArray_Float64! could be applied to other methods).

[ViralBShah]

Perhaps we should add a type annotation on GLOBAL_RNG too.

[rfourquet]

Adding a type annotation on GLOBAL_RNG doesn't seem to improve performances, nor what you suggested above (#8832 (comment)).

[StefanKarpinski]

There isn't yet any mechanism for specializing methods on argument values – even default argument values, although that is quite obviously one of the first things you'd want to do with such an optimization.

[rfourquet]

The global versions of randn can now not be faster than the ones working with a non-global RNG (I think it was the case before). So the use of a global variable for the Julia global RNG doesn't explain all the speed difference in @Jutho 's comment (#8832 (comment)): would you mind testing #8875 to see if this can improve things a bit?

[ivarne]

@rfourquet @ViralBShah
This is on my Backports list for 0.3.3, but I have lost overview over the RNG code. Lots of non-backported changes, makes review much harder.

[ViralBShah]

@rfourquet Would it be possible just to get the array fill improvements without any of the other API changes backported?

[rfourquet]

I could do this, but issue #9037 should be resolved first, so this seems difficult for 0.3.3.