numba gpu：如何计算两个数组的最大相对误差？

Question

我想计算两个数组的相对误差。纯numpy代码是：

# a1, a2 are the two array
np.abs( 1-a2/a1 ).max()

如何使用numba.cuda 加速上述代码？

在我看来：

@cuda.jit
def calculate(a1, a2):
    start = cuda.blockDim.x*cuda.blockIdx.x + cuda.threadIdx.x

    grid = cuda.gridDim.x*cuda.blockDim.x
    for id in range(start, a1.size, grid):
        r = abs(1-a2[id]/a1[id])

ca1 = cuda.to_device(a1)
ca2 = cuda.to_device(a2)

但是，我如何比较不同线程之间的r？

Answer 1

一种可能的方法是编写自己的shared memory parallel reduction。

如 cmets 所示，另一种可能的方法是使用 numba 的内置 reduce decorator。

这是一个展示两者的示例：

$ cat t79.py
from numba import cuda, float32, vectorize
import numpy as np
from numpy import random

#values of 0..10 are legal here
TPBP2 = 9
TPB   = 2**TPBP2
TPBH  = TPB//2
ds = 4096
#method 1: standard cuda parallel max-finding reduction
@cuda.jit
def max_error(a1, a2, err):
    s = cuda.shared.array(shape=(TPB), dtype=float32)
    x = cuda.grid(1)
    st = cuda.gridsize(1)
    tx = cuda.threadIdx.x
    s[tx] = 0
    cuda.syncthreads()
    for i in range(x, a1.size, st):
        s[tx] = max(s[tx], abs(1-a2[i]/a1[i]))
    mid = TPBH
    for i in range(TPBP2):
        cuda.syncthreads()
        if tx < mid:
            s[tx] = max(s[tx], s[tx+mid])
            mid >>= 1
    if tx == 0:
        err[cuda.blockIdx.x] = s[0]

# data
# for best performance we should choose blocks based on GPU occupancy
# but for demonstration since we don't know the GPU:
blocks = (ds+TPB-1)//TPB
a1= np.random.rand(ds).astype(np.float32)
a1 += 1
a2= np.random.rand(ds).astype(np.float32)
err = np.zeros(blocks).astype(np.float32)

# Start the kernel
max_error[blocks, TPB](a1,a2, err)
# we could perform another stage of GPU reduction here, but for simplicity:
my_err = np.max(err)
print(my_err)

#method 2: using numba features
@vectorize(['float32(float32,float32)'], target = 'cuda')
def my_error(a1,a2):
    return abs(1-a2/a1)
@cuda.reduce
def max_reduce(a,b):
    return max(a,b)
r = my_error(a1,a2)
my_err = max_reduce(r)
print(my_err)

$ python t79.py
0.9999707
0.9999707
$