Correct way to generate Poisson-distributed random numbers in Julia GPU code?

Question

For a stochastic solver that will run on a GPU, I'm currently trying to draw Poisson-distributed random numbers. I will need one number for each entry of a large array. The array lives in device memory and will also be deterministically updated afterwards. The problem I'm facing is that the mean of the distribution depends on the old value of the entry. Therefore, I would have to do naively do something like:

CUDA.rand_poisson!(lambda=array*constant)

or:

array = CUDA.rand_poisson(lambda=array*constant)

Both of which don't work, which does not really surprise me, but maybe I just need to get a better understanding of broadcasting? Then I tried writing a kernel which looks like this:

function cu_draw_rho!(rho::CuDeviceVector{FloatType}, λ::FloatType)
    idx = (blockIdx().x - 1i32) * blockDim().x + threadIdx().x
    stride = gridDim().x * blockDim().x
    @inbounds for i=idx:stride:length(rho)
        l = rho[i]*λ
        # 1. variant
        rho[i] > 0.f0 && (rho[i] = FloatType(CUDA.rand_poisson(UInt32,1;lambda=l)))
        # 2. variant
        rho[i] > 0.f0 && (rho[i] = FloatType(rand(Poisson(lambda=l))))
    end
    return
end

And many slight variations of the above. I get tons of errors about dynamic function calls, which I connect to the fact that I'm calling functions that are meant for arrays from my kernels. the 2. variant of using rand() works only without the Poisson argument (which uses the Distributions package, I guess?) What is the correct way to do this?

Answer 1

You may want CURAND.jl, which provides curand_poisson.

using CURAND
n = 10
lambda = .5
curand_poisson(n, lambda)