编写CUDA内核
介绍
与用于CPU编程的传统顺序模型不同，CUDA具有执行模型。在CUDA中，编写的代码将同时由多个线程（通常成百上千个）执行。解决方案将通过定义网格，块和线程层次结构进行建模。
Numba的CUDA支持提供了用于声明和管理此线程层次结构的工具。这些功能与NVidia的CUDA C语言开放的功能非常相似。
Numba还开放了三种GPU内存：全局设备内存（连接到GPU本身的大型，相对较慢的片外内存），片上 共享内存和本地内存。对于除最简单算法以外的所有算法，务必仔细考虑如何使用和访问内存，以最大程度地减少带宽需求和争用，这一点很重要。
内核声明
一个核心功能是指从CPU代码（*）称为GPU功能。它具有两个基本特征：
- 内核无法显式返回值；所有结果数据都必须写入传递给函数的数组中（如果计算标量，则可能传递一个单元素数组）；
- 内核在调用时显式声明其线程层次结构：即线程块数和每个块的线程数（请注意，虽然内核仅编译一次，但可以使用不同的块大小或网格大小多次调用）。
用Numba编写CUDA内核看起来非常像为CPU编写JIT函数：
@cuda.jit
def increment_by_one(an_array):
"""
Increment all array elements by one.
"""
# code elided here; read further for different implementations
（*）注意：较新的CUDA支持设备端内核启动；此功能称为动态并行性，但Numba当前不支持它）
内核调用
通常以以下方式启动内核：
threadsperblock = 32
blockspergrid = (an_array.size + (threadsperblock - 1)) // threadsperblock
increment_by_one[blockspergrid, threadsperblock](an_array)
注意到两个步骤：
- 通过指定多个块（或“每个网格的块”）和每个块的线程数来实例化内核。两者的乘积将给出启动的线程总数。内核实例化是通过采用已编译的内核函数（在此处increment_by_one）并用整数元组对其进行索引来完成的。
- 通过将输入数组（如果需要，以及任何单独的输出数组）传递给内核来运行内核。内核异步运行：启动将其在设备上的执行排队，然后立即返回。可以 cuda.synchronize()用来等待所有先前的内核启动完成执行。
注意
传递驻留在主机内存中的数组，将隐式地导致将副本复制回主机，这将是同步的。在这种情况下，直到将数据复制回内核启动才会返回，因此似乎是同步执行的。
选择块大小
在声明内核所需的线程数时，具有两级层次结构似乎很奇怪。块大小（即每个块的线程数）通常很关键：
- 在软件方面，块大小确定多少线程共享内存的给定区域。
- 在硬件方面，块的大小必须足够大以完全占用执行单元。建议可在 CUDA C编程指南中找到。
多维块和网格
为了帮助处理多维数组，CUDA允许指定多维块和网格。在上面的示例中，可以使blockspergridandthreadsperblock元组为一个，两个或三个整数。与等效大小的一维声明相比，这不会改变所生成代码的效率或行为，但可以帮助以更自然的方式编写算法。
Thread线程定位
运行内核时，内核函数的代码由每个线程执行一次。因此，它必须知道它在哪个线程中，以便知道它负责哪个数组元素（复杂算法可以定义更复杂的职责，但是基本原理是相同的）。
一种方法是让线程确定其在网格和块中的位置，然后手动计算相应的数组位置：
@cuda.jit
def increment_by_one(an_array):
# Thread id in a 1D block
tx = cuda.threadIdx.x
# Block id in a 1D grid
ty = cuda.blockIdx.x
# Block width, i.e. number of threads per block
bw = cuda.blockDim.x
# Compute flattened index inside the array
pos = tx + ty * bw
if pos < an_array.size: # Check array boundaries
an_array[pos] += 1
注意
除非确定块大小和网格大小是阵列大小的除数，否则必须如上所述检查边界。
threadIdx，blockIdx，blockDim和gridDim 是由CUDA后端为知道Thread线程层次结构的几何形状和当前线程的该几何形状内的位置，唯一目的提供特殊对象。
这些对象可以是1D，2D或3D，具体取决于调用内核的方式 。在每个维度访问该值，可使用x，y并z分别这些对象的属性。
numba.cuda.threadIdx
当前线程块中的线程索引。对于1D块，索引（由x属性赋予）是一个整数，范围从0（包括）到numba.cuda.blockDim排除（exclusive）。当使用多个维度时，每个维度都存在类似的规则。
numba.cuda.blockDim
实例化内核时声明的线程块的形状。对于给定内核中的所有线程，即使属于不同的块（即，每个块“已满”），该值也相同。
numba.cuda.blockIdx
线程网格中的块索引启动了内核。对于一维网格，索引（由x属性赋予）是一个整数，范围从0（含）到numba.cuda.gridDim不包含（exclusive）。当使用多个维度时，每个维度都存在类似的规则。
numba.cuda.gridDim
实例化内核时，声明的块网格形状，即此内核调用启动的块总数。
绝对位置
简单的算法将倾向于总是以与上例相同的方式使用线程索引。Numba提供了其它工具来自动执行此类计算：
numba.cuda.grid（ndim ）
返回当前线程在整个块网格中的绝对位置。 ndim应该与实例化内核时声明的维数相对应。如果ndim为1，则返回一个整数。如果ndim为2或3，则返回给定整数的元组。
numba.cuda.gridsize（ndim ）
返回整个块网格中Thread线程的绝对尺寸（或形状）。 ndim与grid()上述含义相同。
使用这些功能，递增示例可以变成：
@cuda.jit
def increment_by_one(an_array):
pos = cuda.grid(1)
if pos < an_array.size:
an_array[pos] += 1
二维数组和线程网格的相同示例为：
@cuda.jit
def increment_a_2D_array(an_array):
x, y = cuda.grid(2)
if x < an_array.shape[0] and y < an_array.shape[1]:
an_array[x, y] += 1
注意，实例化内核时，网格计算仍必须手动完成，例如：
threadsperblock = (16, 16)
blockspergrid_x = math.ceil(an_array.shape[0] / threadsperblock[0])
blockspergrid_y = math.ceil(an_array.shape[1] / threadsperblock[1])
blockspergrid = (blockspergrid_x, blockspergrid_y)
increment_a_2D_array[blockspergrid, threadsperblock](an_array)
进一步阅读
请参阅《CUDA C编程指南》，以详细了解CUDA编程。