rand is slow

[jw3126]

Random number generation is currently slow in some cases. An extreme example is the following:

using FixedSizeArrays

function f(N)
    for _ in 1:N
        rand(Mat{6,4, Int64}, 1:10)
    end
end

function g(N)
    for _ in 1:N
        rand(1:10, 6, 4)
    end
end
f(1)
g(1)
N = 10^6
@time f(N) # 1.481381 seconds (4 allocations: 160 bytes)
@time g(N) # 0.443486 seconds (2.00 M allocations: 274.658 MB, 2.93% gc time)

Despite having huge memory allocation, the Array variant is still faster.

[c42f]

Hmm odd, this is true (though less severe) even when using Vec{3, Int64}.

[SimonDanisch]

After some detective work: turns out that the rand function that does all the work doesn't get inlined and also, the current code seems to copy the range for every call to rand (whaat?!).

julia>using FixedSizeArrays, BenchmarkTools
julia> function tm()
       r = 1:10
       Mat{6,4,Int}(
       (rand(r), rand(r), rand(r) ,rand(r), rand(r), rand(r)),
       (rand(r), rand(r), rand(r) ,rand(r), rand(r), rand(r)),
       (rand(r), rand(r), rand(r) ,rand(r), rand(r), rand(r)),
       (rand(r), rand(r), rand(r) ,rand(r), rand(r), rand(r)),
       )
       end
tm (generic function with 1 method)

julia> function tm2()
       r = Base.Random.RangeGenerator(1:10)
       g = Base.GLOBAL_RNG
       Mat{6,4,Int}(
       (rand(g,r), rand(g,r), rand(g,r) ,rand(g,r), rand(g,r), rand(g,r)),
       (rand(g,r), rand(g,r), rand(g,r) ,rand(g,r), rand(g,r), rand(g,r)),
       (rand(g,r), rand(g,r), rand(g,r) ,rand(g,r), rand(g,r), rand(g,r)),
       (rand(g,r), rand(g,r), rand(g,r) ,rand(g,r), rand(g,r), rand(g,r)),
       )
       end

tm2 (generic function with 1 method)

julia> function test()
              x = 1:10
              b = @benchmark rand($x)
              r1 = median(b.times)
              rg = Base.Random.RangeGenerator(x)
              b = @benchmark rand($(Base.GLOBAL_RNG), $rg)
             r2 = median(b.times)
              b = @benchmark tm()
              fsa1 = median(b.times)
              b = @benchmark tm2()
              fsa2 = median(b.times)
             b = @benchmark rand($(Mat{6,4, Int}), $x)
             fsa3 = median(b.times)
              b = @benchmark rand($x, 6,4)
             a1 = median(b.times)
              println("pure rand ", r1)
              println("special rand ", r2)
              println("theoretical max with pure rand for current impl: ", fsa3/r1/(6*4), " handrolled ", fsa1/r1/(6*4))
              println("theoretical max with special rand for special handrolled ", fsa2/r2/(6*4))
              end
julia>test()
pure rand 41.0
special rand 22.0
theoretical max with pure rand for current impl: 2.0863821138211383 handrolled 0.9664634146341463
theoretical max with special rand for special handrolled 0.9526515151515151

[jw3126]

Thanks for investigating this! I am confused. If I understand your benchmark correctly, a max speed up of 2x (or 4x for special) could theoretically be gained?
But the stupid rand(1:10, 6, 4) version was already nearly 3x as fast?