3

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?

physicsjoke
  • 31
  • 3

1 Answers1

1

You may want CURAND.jl, which provides curand_poisson.

using CURAND
n = 10
lambda = .5
curand_poisson(n, lambda)
Bill
  • 5,600
  • 15
  • 27
  • Hi, I'm having a hard time installing that package... It also seems to be "phased out" and the dependencies (CUDAdrv) report incompatibility. Which version are you running? – physicsjoke Jul 26 '22 at 09:17
  • 1
    CURAND is part of CUDA.jl nowadays, you shouldn't use the CURAND.jl package (it's old, and as you noticed probably not even installable). The equivalent function in CUDA.CURAND is `rand_poisson(!)`. – maleadt Jul 26 '22 at 14:33
  • Old code in a folder, not used lately was the source, physicsjoke. So anyway in that case have you tried rand_poisson already? – Bill Jul 27 '22 at 03:11
  • Yes, with similar not-results :( – physicsjoke Jul 27 '22 at 11:49