numba cuda 不会通过 += 产生正确的结果(需要减少 gpu?)

numba cuda does not produce correct result with += (gpu reduction needed?)

我正在使用 numba cuda 来计算一个函数。

代码只是简单地将所有值加起来成为一个结果,但是 numba cuda 给了我一个与 numpy 不同的结果。

numba 代码

 import math
 def numba_example(number_of_maximum_loop,gs,ts,bs):
        from numba import cuda
        result = cuda.device_array([3,])

        @cuda.jit(device=True)

    def BesselJ0(x):
        return math.sqrt(2/math.pi/x)

    @cuda.jit
    def cuda_kernel(number_of_maximum_loop,result,gs,ts,bs):
        i = cuda.grid(1)
        if i < number_of_maximum_loop:
            result[0] += BesselJ0(i/100+gs)
            result[1] += BesselJ0(i/100+ts)
            result[2] += BesselJ0(i/100+bs)

    # Configure the blocks
    threadsperblock = 128
    blockspergrid = (number_of_maximum_loop + (threadsperblock - 1)) // threadsperblock

    # Start the kernel 
    cuda_kernel[blockspergrid, threadsperblock](number_of_maximum_loop,result,gs,ts,bs) 

    return result.copy_to_host()

numba_example(1000,20,20,20) 

输出:

array([ 0.17770302,  0.34166728,  0.35132036])

numpy 代码

import math
def numpy_example(number_of_maximum_loop,gs,ts,bs):
    import numpy as np
    result = np.zeros([3,])

    def BesselJ0(x):
        return math.sqrt(2/math.pi/x)

    for i in range(number_of_maximum_loop):
        result[0] += BesselJ0(i/100+gs)
        result[1] += BesselJ0(i/100+ts)
        result[2] += BesselJ0(i/100+bs)

    return result

numpy_example(1000,20,20,20) 

输出:

array([ 160.40546935,  160.40546935,  160.40546935])

我不知道我哪里错了。我想我可能会使用减少。但是用一个cuda内核似乎不可能完成。

是的,需要适当的并行缩减来将来自多个 GPU 线程的数据汇总到单个变量。

这是一个简单的例子,说明如何从单个内核完成它:

$ cat t23.py
import math
def numba_example(number_of_maximum_loop,gs,ts,bs):
    from numba import cuda
    result = cuda.device_array([3,])

    @cuda.jit(device=True)
    def BesselJ0(x):
        return math.sqrt(2/math.pi/x)

    @cuda.jit
    def cuda_kernel(number_of_maximum_loop,result,gs,ts,bs):
        i = cuda.grid(1)
        if i < number_of_maximum_loop:
            cuda.atomic.add(result, 0, BesselJ0(i/100+gs))
            cuda.atomic.add(result, 1, BesselJ0(i/100+ts))
            cuda.atomic.add(result, 2, BesselJ0(i/100+bs))

# Configure the blocks
    threadsperblock = 128
    blockspergrid = (number_of_maximum_loop + (threadsperblock - 1)) // threadsperblock

 # Start the kernel
    init = [0.0,0.0,0.0]
    result = cuda.to_device(init)
    cuda_kernel[blockspergrid, threadsperblock](number_of_maximum_loop,result,gs,ts,bs)

    return result.copy_to_host()

print(numba_example(1000,20,20,20))
$ python t23.py
[ 162.04299487  162.04299487  162.04299487]
$

您也可以直接使用 reduce 装饰器对 numba 进行适当的缩减,如 here 所述,尽管我不确定这样可以在单个内核中完成 3 次缩减。

最后,您可以使用 numba cuda 编写一个普通的 cuda 并行缩减,如 所示。我认为将其扩展到在单个内核中执行 3 个缩减应该不难。

当然,这 3 种不同的方法可能会有性能差异。

顺便说一句,如果您想知道我上面的代码与问题中的 python 代码之间的结果差异,我无法解释。当我 运行 你的 python 代码时,我得到的结果与我的答案中的 numba cuda 代码匹配:

$ cat t24.py
import math
def numpy_example(number_of_maximum_loop,gs,ts,bs):
    import numpy as np
    result = np.zeros([3,])

    def BesselJ0(x):
        return math.sqrt(2/math.pi/x)

    for i in range(number_of_maximum_loop):
        result[0] += BesselJ0(i/100+gs)
        result[1] += BesselJ0(i/100+ts)
        result[2] += BesselJ0(i/100+bs)

    return result

print(numpy_example(1000,20,20,20))
$ python t24.py
[ 162.04299487  162.04299487  162.04299487]
$

每个线程执行多个(甚至 1 个)原子添加很可能会产生不良性能。您可能会考虑采用平铺方法(使用共享数组对来自块的贡献或使用 warp 随机播放的 warp 求和),然后对每个 warp 或块执行一次原子操作。

顺便说一下,您发布的数值与 Robert Crovella 回复中的数值之间的差异似乎是 python 版本问题(特别是在 [=15= 中引入 // 用于整数除法) ] 3).在 Crovella 版本中将“i/100”更改为“i//100”并在 python 3 中执行会重现您的值。