Non allocating rand!(A) (get rid of pointer_to_array)

[rfourquet]

I discovered that rand!(A::Array{Int}) was allocating memory in a debugging session, thinking that the allocations came from a type instability. I think this shouldn't, but this patch is quite not pretty. It would be less so if getindex was defined for pointers (I have no hope that this happens), and unneeded if pointer_to_array was modified to not allocate (no idea if this is possible).

Here are timings I get before/after:

julia> function arr{T}(::Type{T}, n, r=1)
    a = zeros(T, n)
    m = MersenneTwister()
    tic()
    a = @allocated for i in 1:r
        rand!(m, a)
    end
    toq(), a/r
end

julia> [arr(Int, 10^i, 10^(9-i)) for i in 1:8] # master:
8-element Array{(Float64,Any),1}:
 (24.119626243,144.0)
 (5.701144431,144.0)
 (1.941585068,160.0)
 (1.328115407,176.0)
 (1.345973947,176.0)
 (1.344547687,176.0)
 (2.675291505,176.0)
 (2.861680458,176.0)

julia> [arr(Int, 10^i, 10^(9-i)) for i in 1:8] # this PR
8-element Array{(Float64,Any),1}:
(6.939816844,0.0)
 (3.781428907,0.0)
 (1.62118118,0.0)
 (1.174099115,0.0)
 (1.241850295,0.0)
 (1.351900588,0.0)
 (2.648783596,0.0)
 (2.839158974,0.0)

[ViralBShah]

What makes 10x10^8 so much slower than 10^8x10?

[rfourquet]

I should have given more explanations. On master, 10x10^8 allocates 144*10^8 bytes (i.e. 13.4 GB, and 38.60 gc time), whereas 10^8x10 allocates 176*10 bytes, so the allocation on small arrays is more visible.
Once the allocation is fixed (this PR), besides the fact that "a small loop within a big one is slower than a big one within a small one" in this case, the dSFMT fill_array! function is used only for lengths greater that 382 (so on small arrays, data has to be copied).

[ViralBShah]

But, why do we have 2 loops? Can't we just iterate over the length of the array?

Can you create a small self-contained example about the memory allocation, which may be easier to work with? Cc: @JeffBezanson

[rfourquet]

By the "outer loop", I was referring to the loop i in 1:r in the function arr I gave as an example. The smallest example about memory allocation is A=[1]; @allocated pointer_to_array(pointer(A),1). This is a small allocation, which occurs within the rand! function.

[JeffBezanson]

Yes this is really ugly. Surely there is another way to write this. I really do not want us to keep chasing moderate speedups at any cost.

[ViralBShah]

I am guessing pointer_to_array allocates memory for the array metadata.

[rfourquet]

Personally, I'm not sure I would want this change if it was only for getting the modest speedup. But the allocations, which grow linearly with the number of times rand! is called, can/will distract users trying to find efficiency bugs (which often come from unexpected allocations) in their programs.

[rfourquet]

What happens here is really pointer arithmetic/dereferencing; the master version wraps this into arrays for a nicer syntax, but uses @inbounds, The PR version only looses the getindex/setindex! syntax. Maybe one better way to write this is to create a lightweight version of pointer_to_array returning an immutable holding only a pointer and a length and which has getindex (maybe such a type already exists?).

[ViralBShah]

Yes, that is essentially the issue here. I don't know if getindex and setindex for pointer has been discussed before. I suspect it is intentionally not provided to discourage such code.

I think this issue will apply for all small array operations, and is not specific to rand. Perhaps the point should be noted in other issues on small array performance.

[ViralBShah]

I wonder how subarray fares here, with the recent improvements.

[ViralBShah]

I am going to close this. Need a better solution, and the issue here is not rand specific. Best to have a new PR if we have a different solution that solves this problem.

[rfourquet]

With the following type, it is enough to change pointer_to_array to UnsafePtrView in random.jl to get the same benefits as this PR:

immutable UnsafePtrView{T} <: DenseArray{T}
    ptr::Ptr{T}
    len::Int
end

UnsafePtrView{T}(A::Array{T}) = UnsafePtrView(pointer(A), length(A))
length(a::UnsafePtrView) = a.len
getindex(a::UnsafePtrView, i::Int) = unsafe_load(a.ptr, i)
setindex!(a::UnsafePtrView, x, i::Int) = unsafe_store!(a.ptr, x, i)
pointer(a::UnsafePtrView) = a.ptr

Would it be an acceptable solution? (until a better general solution is found)

[ViralBShah]

I would be ok to go with such a solution, as it can be easily jettisoned in the future without affecting the rest of the codebase. The gains for small random arrays seem worthwhile. Would love to hear what @JeffBezanson has to say.

[timholy]

See #3224.

[rfourquet]

As noted by @JeffBezanson, rand(::BigInt) is a good example of client of "small" array filling. With this PR, filling arrays of BigInt (with one or two limbs) can be 3x to 5x times faster than on master.

[ViralBShah]

So do we yet have a better way of doing this? I'd rather get the ugly fix in and have an open issue to fix it.