atomic_hist_fill.jl

# Run with
#
#    JULIA_NUM_THREADS=64 numactl -C 0-63 julia
#
# or something like
#
#    for n in 1 2 4 8 16 32 64 128; do JULIA_NUM_THREADS="${n}" numactl -C 0-"$((${n} - 1))" julia atomic_hist_fill.jl ; done


using Base.Threads

using Random, StatsBase


Base.@propagate_inbounds function atomic_addindex!(A::Array{Int64}, v::U, I...) where {U}
    T = Int64
    i = LinearIndices(size(A))[I...]
    v_conv = convert(T, v)
    ptr = pointer(A, i)
    Base.Threads.llvmcall("%ptr = inttoptr i64 %0 to i64*\n%rv = atomicrmw add i64* %ptr, i64 %1 acq_rel\nret i64 %rv\n", T, Tuple{Ptr{T},T}, ptr, v_conv)::T
    A
end


function atomic_push!(h::Histogram{T,N}, xs::NTuple{N,Real}, w::Integer = 1) where {T<:Integer,N}
    # h.isdensity && error("Density histogram must have float-type weights")
    idx = StatsBase.binindex(h, xs)
    if checkbounds(Bool, h.weights, idx...)
        @inbounds atomic_addindex!(h.weights, w, idx...)
    end
    h
end


function locked_push!(l::Base.AbstractLock, h::Histogram{T,N}, xs::NTuple{N,Real}, w::Integer = 1) where {T<:Integer,N}
    # h.isdensity && error("Density histogram must have float-type weights")
    idx = StatsBase.binindex(h, xs)
    if checkbounds(Bool, h.weights, idx...)
        lock(l)
        try
            @inbounds h.weights[idx...] += w
        finally
            unlock(l)
        end
    end
    h
end



function fill_randn_st!(rngs::AbstractVector{<:AbstractRNG}, h::Histogram, n::Integer)
    for i in 1:n
        rng = rngs[threadid()]
        push!(h, (randn(rng), randn(rng)))
    end
    h
end

function fill_randn_mt_atomic!(rngs::AbstractVector{<:AbstractRNG}, h::Histogram, n::Integer)
    @threads for i in 1:n
        rng = rngs[threadid()]
        atomic_push!(h, (randn(rng), randn(rng)))
    end
    h
end

function fill_randn_mt_lock!(rngs::AbstractVector{<:AbstractRNG}, h::Histogram, n::Integer)
    l = SpinLock()
    @threads for i in 1:n
        rng = rngs[threadid()]
        locked_push!(l, h, (randn(rng), randn(rng)))
    end
    h
end


function gen_randn_st(rngs::AbstractVector{<:AbstractRNG}, n::Integer)
    for i in 1:n
        rng = rngs[threadid()]
        (randn(rng), randn(rng))
    end
end

function gen_randn_mt(rngs::AbstractVector{<:AbstractRNG}, n::Integer)
    @threads for i in 1:n
        rng = rngs[threadid()]
        (randn(rng), randn(rng))
    end
end


h = Histogram((-2:0.002:2, -2:0.002:2))
n = 10^6

using Random123
rngs = [Philox4x() for i in 1:nthreads()]

gen_randn_st(rngs, n)
gen_randn_mt(rngs, n)

fill_randn_st!(rngs, h, n)
fill_randn_mt_atomic!(rngs, h, n)
fill_randn_mt_lock!(rngs, h, n)


using BenchmarkTools

#display(@benchmark gen_randn_st($rngs, $n))
#display(@benchmark gen_randn_mt($rngs, $n))

#display(@benchmark fill_randn_st!($rngs, $h, $n))
#display(@benchmark fill_randn_mt_atomic!($rngs, $h, $n))
#display(@benchmark fill_randn_mt_lock!($rngs, $h, $n))

println("Using $(nthreads()) threads (mean time):")

for func in (gen_randn_st, gen_randn_mt)
    trial = @benchmark $func($rngs, $n)
    println("* $func: $(BenchmarkTools.prettytime(time(mean(trial))))")
end

for func in (fill_randn_st!, fill_randn_mt_atomic!, fill_randn_mt_lock!)
    trial = @benchmark $func($rngs, $h, $n)
    println("* $func: $(BenchmarkTools.prettytime(time(mean(trial))))")
end