GPU（图形处理单元）最后是为计算机图形开发的，然而当初它们简直在所有须要高计算吞吐量的畛域无处不在。这一倒退是由 GPGPU(通用 GPU)接口的开发实现的，它容许咱们应用 GPU 进行通用计算编程。这些接口中最常见的是 CUDA，其次是 OpenCL 和最近刚呈现的 HIP。
Python 中应用 CUDA
CUDA 最后被设计为与 C 兼容起初的版本将其扩大到 c ++ 和 Fortran。在 Python 中应用 CUDA 的一种办法是通过 Numba，这是一种针对 Python 的即时 (JIT) 编译器，能够针对 gpu(它也针对 cpu，但这不在咱们探讨的范畴内)。Numba 为咱们提供了一个能够间接应用 Python 子集，Numba 将动静编译 Python 代码并运行它。尽管它没有实现残缺的 CUDA API，但与 cpu 相比它反对的个性曾经能够帮忙咱们进行并行计算的减速。
Numba 并不是惟一的抉择。CuPy 提供了通过基于 CUDA 的并且兼容 Numpy 的高级函数，PyCUDA 提供了对 CUDA API 更细粒度的管制，英伟达也公布了官网 CUDA Python。
本文不是 CUDA 或 Numba 的综合指南，本文的指标是通过用 Numba 和 CUDA 编写一些简略的示例，这样能够让你理解更多 GPU 相干的常识，无论是是不是应用 Python，甚至 C 编写代码，它都是一个很好的入门资源。
GPU 的并行编程简介
GPU 绝对于 CPU 的最大劣势是它们可能并行执行雷同的指令。单个 CPU 内核将一个接一个地串行运行指令。在 CPU 上进行并行化须要同时应用其多个内核（物理或虚构）。例如个别的计算机有 4-8 个内核，而 GPU 领有数千个计算外围。无关这两者的比拟，请参见上面的图 1。GPU 内核通常速度较慢，且只能执行简略的指令，但它们的数量通常能够补救这些毛病。
GPU 编程有四个次要方面问题：
1、了解如何思考和设计并行的算法。因为一些算法是串行设计的，把这些算法并行化可能是很艰难的。
2、学习如何将 CPU 上的构造（例如向量和图像）映射到 GPU 上例如线程和块。循环模式和辅助函数能够帮忙咱们解决这个问题。
3、了解驱动 GPU 编程的异步执行模型。不仅 GPU 和 CPU 互相独立地执行指令，GPU 的流还容许多个解决流在同一个 GPU 上运行，这种异步性在设计最佳解决流时十分重要。
4、抽象概念和具体代码之间的关系：这是通过学习 API 及其细微差别来实现的。
上图为简化的 CPU 架构 (左) 和 GPU 架构 (右)。计算产生在 ALU(算术逻辑单元) 中，DRAM 保留数据，缓存保留的数据能够更快地拜访，但通常容量更小。
开始编写代码
这里的环境要求是:Numba 版本 > 0.55 和一个 GPU。
import numpy as np
import numba
from numba import cuda
print(np.__version__)
print(numba.__version__)
cuda.detect()
# 1.21.6
# 0.55.2
#
# Found 1 CUDA devices
# id 0 b'Tesla T4' [SUPPORTED]
# Compute Capability: 7.5
# PCI Device ID: 4
# PCI Bus ID: 0
# UUID: GPU-e0b8547a-62e9-2ea2-44f6-9cd43bf7472d
# Watchdog: Disabled
# FP32/FP64 Performance Ratio: 32
# Summary:
# 1/1 devices are supportedNumba CUDA 的次要操作时是 CUDA.jit 的装璜器，它定义函数将在 GPU 中运行。
咱们首先写一个简略的函数，它承受两个数字相加而后将它们存储在第三个参数的第一个元素上。
# Example 1.1: Add scalars
@cuda.jit
def add_scalars(a, b, c):
c[0] = a + b
dev_c = cuda.device_array((1,), np.float32)
add_scalars[1, 1](2.0, 7.0, dev_c)
c = dev_c.copy_to_host()
print(f"2.0 + 7.0 = {c[0]}")
# 2.0 + 7.0 = 9.0因为 GPU 只能解决简略的工作，所以咱们的内核只运行一个函数，前面会将这个函数称之为内核。
第一个须要留神的是内核 (启动线程的 GPU 函数) 不能返回值。所以须要通过传递输出和输入来解决这个问题。这是 C 中常见的模式，但在 Python 中并不常见。
在调用内核之前，须要首先在设施上创立一个数组。如果想要显示返回值则须要将它复制回 CPU。这里就有一个隐形的问题：为什么抉择 float32(单精度浮点数)？这是因为尽管大多数 GPU 都反对双精度运算，但双精度运算的工夫可能是单精度运算的 4 倍甚至更长。所以最好习惯应用 np.float32 和 np.complex64 而不是 float / np.float64 和 complex / np.complex128
咱们的函数定义与一般的函数定定义雷同，但调用却略有不同。它在参数之前有方括号:add_scalars1, 1
这些方括号别离示意网格中的块数和块中的线程数，上面应用 CUDA 进行并行化时，会进一步探讨。
应用 CUDA 进行并行化编程
CUDA 网格
当内核启动时它会失去一个与之关联的网格，网格由块组成; 块由线程组成。下图 2 显示了一维 CUDA 网格。图中的网格有 4 个块。网格中的块数保留在一个非凡的变量中，该变量能够在内核中通过 gridDim.x 间接拜访，这里 x 是指网格的第一维度(在本例中是惟一的维度)。二维网格也有通过 y 还有三维网格 z 变量来拜访。到目前 2022 年，还没有四维网格或更高的网格。在内核外部能够通过应用 blockIdx.x 找出正在执行的块，例如咱们这个例子它将从 0 运行到 3。
每个块都有肯定数量的线程，保留在变量 blockDim.x 中。线程索引保留在变量 threadIdx.x 中，在这个示例中变量将从 0 运行到 7。
不同块中的线程被安顿以不同的形式运行，拜访不同的内存区域并在其余一些方面有所不同，本文次要介绍简略的入门所以咱们将跳过这些细节。
当咱们在第一个示例中应用参数 [1,1] 启动内核时，咱们通知 CUDA 用一个线程运行一个块。通过批改这两个值能够应用多个块和多现线程屡次运行内核。threadIdx.x 和 blockIdx.x 每个线程的惟一标识。
上面咱们对两个数组求和，这比对两个数字求和简单：假如每个数组都有 20 个元素。如上图所示，咱们能够用每个块 8 个线程启动内核。如果咱们心愿每个线程只解决一个数组元素，那么咱们至多须要 4 个块。启动 4 个块，每个块 8 个线程，咱们的网格将启动 32 个线程。
对于多线程解决，最须要弄清楚是如何将线程下标映射到数组下标（因为每个线程要独立解决局部数据）。threadIdx.x 从 0 运行到 7，因而它们本人无奈索引咱们的数组，不同的块也具备雷同的 threadIdx.x。然而他们有不同的 blockIdx.x。为了取得每个线程的惟一索引，咱们须要组合这些变量：
i = threadIdx.x + blockDim.x * blockIdx.x对于第一个块，blockIdx.x = 0，i 将从 0 运行到 7。对于第二个块，blockIdx.x = 1。因为 blockDim.x = 8，将失去 8 运行到 15。相似地，对于 blockIdx. x = 2，将从 16 跑到 23。在第四个也是最初一个区块中，将从 24 跑到 31。见下表 1。
这样尽管将每个线程映射到数组中的每个元素……然而当初咱们遇到了一些线程会溢出数组的问题，因为数组有 20 个元素，而 i 的最大值是 32-1。解决方案很简略：对于那些溢出线程，不要做任何事件！
# Example 1.2: Add arrays
@cuda.jit
def add_array(a, b, c):
i = cuda.threadIdx.x + cuda.blockDim.x * cuda.blockIdx.x
if i < a.size:
c[i] = a[i] + b[i]
N = 20
a = np.arange(N, dtype=np.float32)
b = np.arange(N, dtype=np.float32)
dev_c = cuda.device_array_like(a)
add_array[4, 8](a, b, dev_c)
c = dev_c.copy_to_host()
print(c)
# [0. 2. 4. 6. 8. 10. 12. 14. 16. 18. 20. 22. 24. 26. 28. 30. 32. 34. 36. 38.]在较新版本的 Numba 中可能会会收到一条正告，指出咱们应用内核应用了非设施上的数据。这条正告的产生的起因是将数据从主机挪动到设施十分慢，咱们应该在所有参数中应用设施数组调用内核。所以咱们须要事后将数组从主机挪动到设施：
dev_a = cuda.to_device(a)
dev_b = cuda.to_device(b)每个线程惟一索引的计算可能很快就会过期，Numba 提供了非常简单的包装器 cuda.grid，它以网格维度作为惟一参数调用。所以咱们新的内核将如下：
# Example 1.3: Add arrays with cuda.grid
@cuda.jit
def add_array(a, b, c):
i = cuda.grid(1)
if i < a.size:
c[i] = a[i] + b[i]
add_array[4, 8](dev_a, dev_b, dev_c)
c = dev_c.copy_to_host()
print(c)
# [0. 2. 4. 6. 8. 10. 12. 14. 16. 18. 20. 22. 24. 26. 28. 30. 32. 34. 36. 38.]如果咱们扭转数组的大小时会产生什么？咱们这里不扭转函数而更改网格参数（块数和每个块的线程数），这样就相当于启动至多与数组中的元素一样多的线程。
设置这些参数有一些”迷信“和一些”艺术“。对于“迷信”：（a）它们应该是 2 的倍数，通常在 32 到 1024 之间，并且（b）它们应该都被应用以最大化占用率（有多少线程同时处于活动状态）。所以 Nvidia 提供了一个能够帮忙计算这些的电子表格。（https://docs.nvidia.com/cuda/…）对于“艺术”而言，没有什么能够预测内核的行为，因而如果真的想优化这些参数，须要依据不同的输出来剖析代码。然而依据教训个别 GPU 的“正当”线程数是 256。
N = 1_000_000
a = np.arange(N, dtype=np.float32)
b = np.arange(N, dtype=np.float32)
dev_a = cuda.to_device(a)
dev_b = cuda.to_device(b)
dev_c = cuda.device_array_like(a)
threads_per_block = 256
blocks_per_grid = (N + (threads_per_block - 1)) // threads_per_block
# Note that
# blocks_per_grid == ceil(N / threads_per_block)
# ensures that blocks_per_grid * threads_per_block >= N
add_array[blocks_per_grid, threads_per_block](dev_a, dev_b, dev_c)
c = dev_c.copy_to_host()
np.allclose(a + b, c)
# True代码方面是没有问题了，然而因为硬件限度 GPU 不能运行任意数量的线程和块。通常每个块不能超过 1024 个线程，一个网格不能超过 2¹⁶ – 1 = 65535 个块。然而这并不是说能够启动 1024 × 65535 个线程……因为还须要依据寄存器占用的内存量以及其余的因素思考。还有一点就是解决不适宜 GPU RAM 的大型数组（也就是 OOM）。在这些状况下，能够应用多个 GPU 分段解决数组（单机多卡）。
在 Python 中，硬件限度能够通过 Nvidia 的 cuda-python 库的函数 cuDeviceGetAttribute 取得，具体请查看该函数阐明。
Grid-stride 循环
在每个网格的块数超过硬件限度但显存中能够包容残缺数组的状况下，能够应用一个线程来解决数组中的多个元素，这种办法被称为 Grid-stride。除了克服硬件限度之外，Grid-stride 还受害于重用线程，这样能够将线程创立 / 销毁开销降到最低。Mark Harris 的文章 CUDA Pro Tip: Write Flexible Kernels with Grid-Stride Loops 具体介绍这个办法，这里就不再介绍了。
在 CUDA 内核中增加一个循环来解决多个输出元素，这个循环的步幅等于网格中的线程数。这样如果网格中的线程总数 (threads_per_grid = blockDim.x * gridDim.x) 小于数组的元素数，则内核解决完索引 cuda.grid(1)它将解决索引 cuda.grid(1) + threads_per_grid，直到解决完所有数组元素，咱们来看代码。
# Example 1.4: Add arrays with grid striding
@cuda.jit
def add_array_gs(a, b, c):
i_start = cuda.grid(1)
threads_per_grid = cuda.blockDim.x * cuda.gridDim.x
for i in range(i_start, a.size, threads_per_grid):
c[i] = a[i] + b[i]
threads_per_block = 256
blocks_per_grid_gs = 32 * 80 # Use 32 * multiple of streaming multiprocessors
# 32 * 80 * 256 < 1_000_000 so one thread will process more than one array element
add_array_gs[blocks_per_grid_gs, threads_per_block](dev_a, dev_b, dev_c)
c = dev_c.copy_to_host()
np.allclose(a + b, c)
# True这段代码与下面的代码十分类似，不同之处在于从 cuda.grid(1)开始，然而执行更多的元素，每个 threads_per_grid 执行一个元素，直到达到数组的开端。
那么那种形式更快呢?
CUDA 内核的计算工夫
GPU 编程的指标就是进步速度。因而精确测量代码执行工夫十分重要。
CUDA 内核是由主机 (CPU) 启动的设施函数但它们是在 GPU 上执行的，GPU 和 CPU 不通信（除非咱们让它们通信）。因而当 GPU 内核被启动时，CPU 将简略地持续运行后续指令，不论它们是启动更多的内核还是执行其余 CPU 函数。所以如果在内核启动前后别离调用 time.time()，则只取得了内核启动所需的工夫，而不是计算运行所需的工夫。
所以这里就须要进行同步，也就是调用 cuda.synchronize()函数，这个函数将进行主机执行任何其余代码，直到 GPU 实现已在其中启动的每个内核的执行。
为了计算内核执行工夫，能够简略地计算内核运行而后同步所需的工夫，然而须要应用 time.perf_counter() 或 time.perf_counter_ns() 而不是 time.time()。
在应用 Numba 时，咱们还有一个细节须要留神：Numba 是一个 Just-In-Time 编译器，这意味着函数只有在被调用时才会被编译。因而计时函数的第一次调用也会计时编译步骤，这通常要慢得多。所以必须首先通过启动内核而后对其进行同步来编译代码，确保了下一个内核无需编译即可立刻运行，这样失去的工夫才是精确的。
from time import perf_counter_ns
# Compile and then clear GPU from tasks
add_array[blocks_per_grid, threads_per_block](dev_a, dev_b, dev_c)
cuda.synchronize()
timing = np.empty(101)
for i in range(timing.size):
tic = perf_counter_ns()
add_array[blocks_per_grid, threads_per_block](dev_a, dev_b, dev_c)
cuda.synchronize()
toc = perf_counter_ns()
timing[i] = toc - tic
timing *= 1e-3 # convert to μs
print(f"Elapsed time: {timing.mean():.0f} ± {timing.std():.0f} μs")
# Elapsed time: 201 ± 109 μs
# Compile and then clear GPU from tasks
add_array_gs[blocks_per_grid_gs, threads_per_block](dev_a, dev_b, dev_c)
cuda.synchronize()
timing_gs = np.empty(101)
for i in range(timing_gs.size):
tic = perf_counter_ns()
add_array_gs[blocks_per_grid_gs, threads_per_block](dev_a, dev_b, dev_c)
cuda.synchronize()
toc = perf_counter_ns()
timing_gs[i] = toc - tic
timing_gs *= 1e-3 # convert to μs
print(f"Elapsed time: {timing_gs.mean():.0f} ± {timing_gs.std():.0f} μs")
# Elapsed time: 178 ± 141 μs对于简略的内核，还能够测量算法的整个过程取得每秒浮点运算的数量。通常以 GFLOP/s (giga-FLOP /s)为度量单位。加法操作只蕴含一个触发器: 加法。因而，吞吐量由:
# G * FLOP / timing in s
gflops = 1e-9 * dev_a.size * 1e6 / timing.mean()
gflops_gs = 1e-9 * dev_a.size * 1e6 / timing_gs.mean()
print(f"GFLOP/s (algo 1): {gflops:.2f}")
print(f"GFLOP/s (algo 2): {gflops_gs:.2f}")
# GFLOP/s (algo 1): 4.99
# GFLOP/s (algo 2): 5.63一个二维的例子
上面应用一个更简单的例子，咱们创立一个 2D 内核来对图像利用对数校对。
给定一个值在 0 到 1 之间的图像 I(x, y)，对数校对图像由
Iᵪ(x, y) = γ log₂(1 + I(x, y))
首先让咱们获取一些数据!
import matplotlib.pyplot as plt
from skimage import data
moon = data.moon().astype(np.float32) / 255.
fig, ax = plt.subplots()
im = ax.imshow(moon, cmap="gist_earth")
ax.set_xticks([])
ax.set_yticks([])
ax.set_xticklabels([])
ax.set_yticklabels([])
fig.colorbar(im)
fig.tight_layout()数据在较低端饱和了。简直没有高于 0.6 的值。
当初编写核函数。
import math
# Example 1.5: 2D kernel
@cuda.jit
def adjust_log(inp, gain, out):
ix, iy = cuda.grid(2) # The first index is the fastest dimension
threads_per_grid_x, threads_per_grid_y = cuda.gridsize(2) # threads per grid dimension
n0, n1 = inp.shape # The last index is the fastest dimension
# Stride each dimension independently
for i0 in range(iy, n0, threads_per_grid_y):
for i1 in range(ix, n1, threads_per_grid_x):
out[i0, i1] = gain * math.log2(1 + inp[i0, i1])
threads_per_block_2d = (16, 16) # 256 threads total
blocks_per_grid_2d = (64, 64)
moon_gpu = cuda.to_device(moon)
moon_corr_gpu = cuda.device_array_like(moon_gpu)
adjust_log[blocks_per_grid_2d, threads_per_block_2d](moon_gpu, 1.0, moon_corr_gpu)
moon_corr = moon_corr_gpu.copy_to_host()
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.imshow(moon, cmap="gist_earth")
ax2.imshow(moon_corr, cmap="gist_earth")
ax1.set(title="Original image")
ax2.set(title="Log-corrected image")
for ax in (ax1, ax2):
ax.set_xticks([])
ax.set_yticks([])
ax.set_xticklabels([])
ax.set_yticklabels([])
fig.tight_layout()第一个 for 循环从 iy 开始，第二个最外面的循环从 ix 开始。为什么抉择这个程序呢? 这种抉择的内存拜访模式更无效。因为第一个网格索引是最快的，所以咱们想让它匹配最快的维度: 最初一个维度。
后果如下:
总结
本文中介绍了 Numba 和 CUDA 的基础知识，咱们能够创立简略的 CUDA 内核，并将其从内存挪动到 GPU 的显存来应用它们。还介绍了如何应用 Grid-stride 技术在 1D 和 2D 数组上迭代。
附录: 应用 Nvidia 的 cuda-python 来查看设施属性
要对 GPU 的确切属性进行细粒度管制，能够应用 Nvidia 提供的较低级别的官网 CUDA Python 包，代码如下：
# Need to: pip install --upgrade cuda-python
from cuda.cuda import CUdevice_attribute, cuDeviceGetAttribute, cuDeviceGetName, cuInit
# Initialize CUDA Driver API
(err,) = cuInit(0)
# Get attributes
err, DEVICE_NAME = cuDeviceGetName(128, 0)
DEVICE_NAME = DEVICE_NAME.decode("ascii").replace("\x00", "")
err, MAX_THREADS_PER_BLOCK = cuDeviceGetAttribute(CUdevice_attribute.CU_DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK, 0)
err, MAX_BLOCK_DIM_X = cuDeviceGetAttribute(CUdevice_attribute.CU_DEVICE_ATTRIBUTE_MAX_BLOCK_DIM_X, 0)
err, MAX_GRID_DIM_X = cuDeviceGetAttribute(CUdevice_attribute.CU_DEVICE_ATTRIBUTE_MAX_GRID_DIM_X, 0)
err, SMs = cuDeviceGetAttribute(CUdevice_attribute.CU_DEVICE_ATTRIBUTE_MULTIPROCESSOR_COUNT, 0)
print(f"Device Name: {DEVICE_NAME}")
print(f"Maximum number of multiprocessors: {SMs}")
print(f"Maximum number of threads per block: {MAX_THREADS_PER_BLOCK:10}")
print(f"Maximum number of blocks per grid: {MAX_BLOCK_DIM_X:10}")
print(f"Maximum number of threads per grid: {MAX_GRID_DIM_X:10}")
# Device Name: Tesla T4
# Maximum number of multiprocessors: 40
# Maximum number of threads per block: 1024
# Maximum number of blocks per grid: 1024
# Maximum number of threads per grid: 2147483647https://avoid.overfit.cn/post/9c18d13573384744934f52ad9361f55c
作者：Carlos Costa