Misunderstanding Allocs for Monte Carlo summation

I have a simple test function here:

function tester(nreps::Int, n_x::Int)
    res = zeros(Float64, n_x)
    for _ = 1:nreps
        res += rand(Float64, n_x)
    end
    return res
end

I do not understand the behaviour /meaning of Allocs as reported by BenchmarkTools
@benchmark tester(10, 1000) -> allocs estimate 21
@benchmark tester(100, 1000) -> allocs estimate 201

Replacing rand by rand! helps a little bit:

function tester(nreps::Int, n_x::Int)
    res = zeros(Float64, n_x)
    tmp = zeros(Float64, n_x)
    for _ = 1:nreps
        rand!(tmp)
        res += tmp
    end
    return res
end

but again I find that increasing nreps by a factor of 10 also increases allocations by a factor of 10.

Equivalent code in e.g. C should require the same memory for all values of nreps, since I am simply mutating an array in-place. What is Julia doing here, and how do I avoid the behaviour?

rand(Float64, n_x) is allocating an array of n_x elements. You probably want res .+= rand.(Float64)

2 Likes

So in general, the fast way to write this (with the understanding that rand! will be replaced by some other nontrivial function) is

function tester(nreps::Int, n_x::Int)
    res = zeros(Float64, n_x)
    tmp = zeros(Float64, n_x)
    for _ = 1:nreps
        rand!(tmp)
        res .+= tmp
    end
    return res
end

exactly.