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