julia - 在 Julia GPU 代码中生成泊松分布随机数的正确方法？

Question

对于将在 GPU 上运行的随机求解器，我目前正在尝试绘制服从泊松分布的随机数。我需要为大型数组的每个条目分配一个数字。该数组存在于设备内存中，之后也会进行确定性更新。我面临的问题是分布的均值取决于条目的旧值。因此，我将不得不天真地做类似的事情:

CUDA.rand_poisson!(lambda=array*constant)

或:

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

这两个都不行，这并没有让我感到惊讶，但也许我只是需要更好地了解广播？然后我尝试编写一个如下所示的内核:

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

以及上述的许多细微变化。我收到大量关于动态函数调用的错误，这与我正在从我的内核中调用用于数组的函数这一事实有关。 2. 使用 rand() 的变体仅在没有 Poisson 参数的情况下有效(我猜它使用 Distributions 包？)执行此操作的正确方法是什么？

Answer 1

您可能需要 CURAND.jl，它提供 curand_poisson。

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