Arrays and performance

Allocate or reuse a buffer

using PureRNGs, Random

rng = Philox4x32(7)
matrix, rng = rand_next(rng, Float32, 3, 4)

buffer = similar(matrix)
result, rng = rand_next!(rng, buffer)
@assert result === buffer

Dimensions are positional integers or one Dims tuple, as in Random. For uniform draws, put a result type before tuple dimensions: rand_next(rng, Float64, (3, 4)). A lone tuple is a population to pick from. Arrays use Julia's native column-major order.

A bang modifies the destination, not the immutable generator. rand_next! returns the destination and the advanced generator.

Use randn_next! and randexp_next! for normal and exponential fills. The ordinary Random fills require StatefulRNG.

Destinations that take the packed CPU path

Any destination is filled with the same stream, but the speed depends on how it indexes. A CPU destination whose IndexStyle is IndexLinear takes the packed fill, which decodes whole generator blocks at once: an Array, a contiguous or strided view of one, and a reshape of either. A destination whose style is IndexCartesian, such as a transpose or a view with a strided second index, falls back to one draw per element and runs several times slower. Copy such a destination into an Array, fill the Array, and copy back when the fill is on a hot path.

Process a long stream in chunks

function sum_batches(rng, buffer, batches)
    total = 0.0
    for _ in 1:batches
        _, rng = rand_next!(rng, buffer)
        total += sum(buffer)
    end
    return rng, total
end

buffer = Vector{Float64}(undef, 1024)
rng, total = sum_batches(Philox4x32(7), buffer, 4)
total

2067.5739637771267

Buffer reuse avoids repeated allocation. Continuing the returned state preserves the stream across chunk boundaries.

Control CPU threading

CPU fills run serially on the calling task by default, as Random fills do. Pass threaded=true to split a large fill into chunks across Julia's threads:

_, rng = rand_next!(rng, buffer; threaded=true)

([0.5731661621609925, 0.5856828351875951, 0.9045124819914252, 0.9967330155241246, 0.9644495975340932, 0.5222628162938953, 0.011823924712811684, 0.005738844092029294, 0.963658682913813, 0.08707696111753016  …  0.10399458004029005, 0.33213879957555437, 0.6262134002731766, 0.19805431075574376, 0.18541126914573924, 0.5645472322276337, 0.812014487245582, 0.9185997906193446, 0.9333417637691038, 0.2587995672333937], Philox4x32((0x00000007, 0x00000000), 271360))

Destination fills, allocating draws, sampling, and splitrng(rng, n) all take the keyword. It never changes the values: a threaded fill writes the same stream as the serial one. It affects CPU execution only. Fills too small to split stay on the calling task. Leave it off inside your own threaded loop, so the outer loop owns the parallelism.

Each internally threaded call creates tasks and synchronization objects. Reusing the destination does not remove those allocations. Large fills can amortize this cost. Serial fills avoid the task allocation.

See Parallel jobs for a complete outer-threaded workflow.

Packed booleans

Boolean draws use individual random bits. An ordinary Array{Bool} still stores one byte per element.

packed = falses(128)
_, rng = rand_next!(rng, packed)
packed

128-element BitVector:
 1
 1
 0
 1
 1
 0
 0
 1
 0
 1
 ⋮
 1
 0
 1
 0
 0
 1
 1
 1
 1

BitArray destinations are CPU-only. GPU Boolean arrays use their backend's array storage.

Measure the workload

Use the repository's benchmark/throughput.jl from a Julia session. Select the generator, result type, and device there.

Warm compilation before timing. Compare allocating draws separately from buffer fills. For GPU measurements, include synchronization and distinguish device generation from host transfers.

Scalar draws inside a kernel and bulk array fills are different workloads. Neither timing predicts the other.

See GPU kernels for asynchronous fills and Particle filter for a complete workload.