使用 CUDA 通过 GPU 编程实现性能最大化
CPU 并不是唯一能进行并行计算的处理器，还有图形处理器（GPU）。这些处理器旨在极其快速高效地处理矢量数据，用于图像渲染、三维引擎和多边形基元的处理。它们的发展主要得益于游戏市场的大力推动，多年来取得了非凡的成果。有趣的是，这项技术和 GPU 也已成为科学计算领域非常强大的工具。

我们将讨论以下主题：

GPU 架构
GPU 编程
用于 CUDA 的 Numba
OpenCL
GPU 架构
实际上，GPU 可被视为由数百万个轻型内核并行组成，每个内核都能独立进行计算。这种高度并行的结构能够在极短的时间内对超大矢量和矩阵进行计算，这当然是 CPU 无法做到的。因此，这些特性最近引起了科学界的关注，导致了该领域开发模型和工具以及高度专业化库的发展。

让我们详细了解一下GPU的内部结构。这些概念将有助于更好地理解 GPU 编程模型。GPU 的内部结构分为两级并行。第一个层次由组成 GPU 的多个处理单元组成，称为流式多处理器（SM Streaming Multiprocessor）。这些单元各自独立工作，同时执行任务。第二级并行发生在 SM 本身内部。每个 SM 又被细分为更多的计算单元，称为流处理器 (SP Stream Processor)。每个 SP 都有一个能够执行线程的内核，该线程独立工作，并与其他 SP 同时工作。

SM 和 SP 的细分是结构性的，在此基础上可以创建一个逻辑组织，我们将把它用于 GPU 编程模型。SP 是 GPU 的基本计算单元，可分为多个逻辑块。这些逻辑块将用于执行特定的执行模式。也就是说，其中的所有内核将同时执行相同的指令，但指向不同的内存区域，各自处理分配给它们的特定部分数据。它们将访问 GPU 本身内部的内存，以便尽可能快速地访问内存中的值。然后，CPU 将再次传输并重新加载这些数据，以便进一步处理。

GPU 具有高度并行的特性，拥有数千或数百万个处理单元（或线程），能够进行同样多的并行计算，因此它的功能非常强大，但又与 CPU 完全不同。因此，有必要开发不同于传统编程方法的特定编程方法。您必须利用不同的驱动程序和编译器，这些驱动程序和编译器能够为 GPU 处理的特定语言编译部分代码，而 GPU 又是由不同的制造商生产的，因此需要针对每种 GPU 采用特定的方法。

目前，有一些更有效的解决方案可以让我们利用 GPU 的计算能力，并将其与传统的、通常是顺序执行的 CPU 相结合。在不同的显卡制造商中，英伟达（NVIDIA）公司可能是在计算能力和在 PC 上的普及方面走得最远的一家。对于这种类型的显卡，有一种专门的开发环境，称为 “计算统一设备架构（CUDA Compute Unified Device Architecture）”，是专为英伟达显卡设计的。这种开发环境主要以 C 语言为基础，利用特定的扩展功能可以使用 GPU 的特定功能，但也可以使用其他编程语言，如 Fortran、Java，尤其是 Python。

OpenCL 框架是另一项允许我们将 GPU 与 CPU 一起编程以挖掘其并行计算潜力的技术。它也是基于 C 和 C ++ 语言，但随着时间的推移，也可以将其集成到 Python 等其他编程语言中。与专为英伟达™（NVIDIA®）显卡创建的 CUDA 相比，这项技术更具通用性。事实上，由于发布了集成 OpenCL 的 SDK 工具包和特定驱动程序，OpenCL 不仅可以使用英伟达™（NVIDIA®）显卡，还可以使用其他制造商（英特尔、AMD 等）的 GPU 显卡。至于英伟达™（NVIDIA®）的 GPU 显卡，OpenCL 仍将依赖于完全兼容并可与之集成的 CUDA 技术。

Python GPU 编程
由于一系列库（通常是相应的 C 库在 Python 中的封装）的存在，GPU 编程在 Python 中得以实现。这些库的发展如此迅速，以至于其中一些库随着时间的推移或多或少地取得了成功，与这些技术的不断发展交替进行。目前，有两个库特别适用于 GPU 并行编程：

Numba
pyOpenCL
当然，它们并不是这方面的唯一库。还有许多其他库可用，它们通常使用相同的基本技术，如 CUDA 和 OpenCL，但其中一些库的使用相当困难，通常与其他库不兼容，需要用 C 语言编译基础库（安装编译器规范，必须符合特定版本）。通常情况下，使用这些库建立开发环境需要花费数天时间，更新的文档也很少，而且并不总能取得成功。Numba 和 pyOpenCL 与 Python 软件包管理环境（如 Pip 和 Anaconda）完美集成，因此安装快捷方便，可以立即开始工作。

Numba

Numba 是一个开源库，包含一个 JIT（Just-In-Time）编译器，能够使用 LLVM（Low Level Virtual Machine）技术将 Python 代码和 NumPy 对象的一部分直接翻译成特定的机器代码。这项技术由 llvmilite 软件包提供，已集成到 Numba 中。得益于这种集成编译系统，Numba让我们能够以一种简单的方式创建使用CPU和GPU特性的代码，而无需对Python语言做太多改动。

因此，Numba允许我们用Python语言编写经典代码，并在其中添加特殊符号（适当的Python装饰符），这些符号将标识出专门针对不同GPU卡编译的代码部分。

用于 CUDA 的 Numba

作为本书的一部分，我们将使用Numba来专门处理英伟达™（NVIDIA®）GPU，因此将依赖于CUDA。要使用这项技术，首先需要了解其编程模型和基本概念。只有了解了这些概念，我们才能很好地理解要想最大限度地发挥 GPU 的潜能所需要执行的各个步骤。

首先，在这个编程模型中，您需要确定运行用 Python 编写的程序的 CPU 及其连接的 GPU。CPU 将被识别为主机，因为它将承载程序执行的主流程，通常是顺序执行。GPU 将被识别为设备，因为它们可能是异构的，需要不同的驱动程序和编译器才能运行。专为同时在 CPU 和 GPU 上运行而编写的 Python 代码将在主机上以传统方式编码和执行。当发现有部分代码要在 GPU 上运行时，将根据其所针对的技术进行专门编译。由于我们考虑的是英伟达™（NVIDIA®）公司的 GPU 显卡，Numba 将使用 CUDA 技术及其内部的编译器来创建特定的编译包（称为内核），并将其发送到这些显卡上运行。我们将在下文中看到，这些内核在代码中很容易识别，其形式是由特定装饰器标识的函数。

除了内核，还有其他一些针对 GPU 的函数，称为设备函数。这些函数都是专门为在不同设备上运行而编译的，但内核在 GPU 上运行时会调用这些函数（而不是主机）。

在内存方面，也有必要做一些说明。与所有传统 Python 程序一样，在 CPU 上执行的程序部分可以访问系统内存，我们将其定义为主机内存。至于内核，它们将在 GPU 上执行，因此必须访问这些卡上的特定内存，我们将其定义为设备内存。

在 GPU 编程模型中，你必须明确分配这种类型的内存（主机内存和设备内存），并管理它们之间的传输。由于 GPU 能够以最佳方式处理矢量，因此将在基于 CUDA 的 NumPy 数组模型中使用矢量，该模型专门针对 GPU 上的这些传输和分配操作进行了优化。

通过重新列出使用 CUDA 的 GPU 编程模型的主角，我们将得到以下列表：

主机Host（CPU）
设备Device（GPU）
内核Kernel（由主机启动并在设备上执行的 GPU 功能）
设备函数Device function（由设备启动并在设备上执行的 GPU 函数）
主机内存（系统的主内存）
设备内存（GPU 卡上的内存）
GPU 编程模型的逻辑层次结构

正如我们所见，CUDA 采用的执行模式与 CPU 编程采用的传统顺序模式不同。使用 CUDA 时，代码将由多个线程（通常是数百或数千个）同时执行。除了包括 SM 和 SP 单元在内的多级物理结构外，CUDA 编程模型还提供了通过线程层次结构（按网格和区块分组）组织的逻辑细分。逻辑层次结构的代表图如图所示：

Numba 提供了一整套函数、类和属性，让您能够管理这种层次结构和组成它的元素。在每个层次结构中，都对应着 GPU 中特定类型的内存：

全局内存(Global memory)
共享内存(Shared memory)
本地内存(Local memory)
在使用 Numba（以及其他使用 CUDA 的程序库）执行程序时，必须始终牢记这些不同类型的内存。要创建高效算法，就必须优化使用这些内存，尽量减少访问和数据交换，以及需要时间资源的操作。
安装 CUDA

首先要做的是从互联网上安装 CUDA 工具包。它可以在英伟达™（NVIDIA®）开发者网站的专门页面（https://developer.nvidia.com/cuda-downloads）上免费下载。

在我们的例子中，选择系统特性（Windows 11, x86_64）后，会自动下载以下版本：cuda_11.7.1_516.94_windows.exe

当然，就您的情况而言，您将获得最新的版本。下载可执行文件后，请将其安装到计算机上。按照安装过程中的提示进行操作。

安装 Numba for CUDA

在Anaconda平台上安装Numba非常简单。要安装Numba，只需打开为使用GPU而专门准备的虚拟环境控制台，然后输入以下命令：conda install numba cudatoolkit

如您所见，还必须安装 cudatoolkit 软件包。许多使用英伟达™（NVIDIA®）GPU 的库都已经包含了该软件包的依赖项，但 Numba 是一个用途广泛的库，并不总是必须使用该类型的 GPU。因此，有必要单独安装 cudatoolkit。

声明和调用内核

在启动程序之前，我们看到用于在 GPU 上执行的代码部分被定义为内核函数。这类函数有两个基本特征，区别于 Python 中的普通函数。首先，内核函数是不能显式返回值的函数，也就是说，它们不能使用返回值来返回函数之外的值。这些函数只对作为参数传递给函数的数组中包含的数据工作，因此这些数组中必须有一个包含计算结果。在内核声明中，还必须考虑到调用时将设置的线程层次结构。因此，必须指定参与计算的线程块数量以及每个线程块中的线程数量。

在 Python 代码中，通过使用特定的装饰器，可以轻松识别实现内核函数的代码部分： @cuda.jit：

@cuda.jit

def kernel_func(array):

“””

Code here

“””
一旦定义了内核，就必须在代码需要时调用它。不过，在此之前，正如我们所说，有必要表达执行内核时必须设置和使用的线程层次结构。然后，我们将设置每个区块的线程数，以及每个网格的区块数：

threadsperblock = 32
blockspergrid = (array.size // threadsperblock) + 1
很明显，通过这样定义这两个参数，我们可以确保每次都考虑到所使用数组的大小，并将相应数量的线程按照适合以最优方式计算的层次进行排列。

然后，我们将调用内核函数，将这些参数指定为函数的索引。这两个值的乘积将实际匹配计算中使用的线程：

kernel_funcblockspergrid, threadsperblock
一旦调用，内核函数将被编译并作为内核发送到相应的设备，由所有被调用的线程执行。所有内核都将以异步方式并行运行，每个内核都独立运行。可以通过添加 cuda.synchronize()来同步计算。在这种情况下，代码将等待所有内核执行完毕后再继续执行代码。

我们将在下面的示例中看到，选择正确的线程层次结构对于在 GPU 上实现最佳计算性能至关重要。区块大小决定了有多少线程将共享某个内存区域（共享内存），同时还必须足够大，以便能够在一个步骤中完成输入数据的计算。

此外，在内核函数声明中，我们必须对相关线程和区块的索引进行优化管理，以确保每个线程都能准确知道要处理输入数组中的哪个元素。在这方面，CUDA 提供了一些特殊对象。例如，threadIdx、threadIdy 和 threadIdz 指定了一个区块中使用的各线程在三个维度上的索引。而 blockIdx、blockIdy 和 blockIdz 则指定了网格内各个区块的索引。区块的大小可以通过 blockDim.x、blockDim.y 和 blocDim.z 来指定。

有鉴于此，举例来说，如果你想实现一个内核，只对数组中的所有元素做 1 的递增，我们将不得不使用之前的索引对象，如下所示：

@cuda.jit
def add_one(a):
tx = cuda.threadIdx.x
ty = cuda.blockIdx.x
dim = cuda.blockDim.x

pos = tx + ty * dim
if pos < a.size:  
    a[pos] += 1

由于我们考虑的数组是一维的，因此我们可以在内核中考虑只有 x 存在的一维几何图形。因此，计算中涉及的所有索引之间的迭代可以用表达式与数组元素的位置相关联：pos = tx + ty* dim

此外，在计算过程中还可以设置一个条件，以确保所获得的索引包含在数组的维数中。我们添加以下条件 if pos < array.size：

最终，所有以数组中每个元素的 pos 为索引的内核将以并行和异步的方式同时增加一个单位的值。

设备函数

至于设备函数，与内核函数相比，它们在声明方面略有不同。只需在装饰器中添加 device = True 即可：

@cuda.jit(device=True)
def device_func(a, b):
return a + b
与内核函数不同，设备函数可以像所有普通函数一样返回值。

使用 Numba 的编程示例

到目前为止，我们所积累的有关 Numba 的概念已经足够让我们将其付诸实践，并编译一小段示例代码。我们使用之前的内核，将数组元素的值增加一个单位。

让我们编写以下代码：

import numpy as np
from numba import cuda

@cuda.jit
def add_one(a):
tx = cuda.threadIdx.x
ty = cuda.blockIdx.x
dim = cuda.blockDim.x

pos = tx + ty * dim
if pos < a.size:  
    a[pos] += 1

n = 10
a_host = np.random.random(n)
print(“Vector a: %s” %a_host)
a_dev = cuda.to_device(a_host)

threadsperblock = 128
blockspergrid = (a_host.size // threadsperblock) + 1
add_onethreadsperblock, blockspergrid

a_host = a_dev.copy_to_host()
print(“New Vector a: %s” %a_host)

通过使用函数 cuda.to_device()，我们将这个数组转移到了设备内存中，该函数在设备内存中分配了一个与向量 a_host 相同的内存空间，然后复制了与数组 a_dev 标识的内容（因为它位于设备内存中）。

下一步是定义内核运行所涉及的线程层次。在这里，我们也复制了之前看到的定义, 我们插入了大量线程（n = 128），以避免在执行过程中输出警告信息

NumbaPerformanceWarning: Grid size 10 will likely result in GPU under-utilization due to low occupancy.
我们之所以看到这条错误信息，是因为程序建议使用的网格大小应与 GPU 中的 SP 数量相当，且不能太小，以免无法充分发挥 GPU 硬件的潜力。

在为线程定义了适当的层次结构后，内核函数会在 a_dev 数组上调用刚刚定义的参数：

add_onethreadsperblock,blockpergrid
调用内核后，就必须管理计算结果，即元素增加了一个单位的向量 a_dev。但该向量位于 GPU 上，因此需要将其传回主机内存。为此，请使用 copy_to_host()函数将 GPU（设备）上的内存区域复制到主机上：a_host = a_dev.copy_to_host()

现在变量 a_host 将包含更新后的值。我们可以在输出中打印它，并将结果与计算前的值进行比较：print(“New Vector a: %s” %a_host)

完整运行代码后，我们将得到类似下面的结果。

python numba_01.py

Vector a: [0.84178191 0.52956292 0.17327604 0.93608128 0.51674129 0.34798135
0.73725271 0.02645266 0.77222292 0.80086225]
New Vector a: [1.84178191 1.52956292 1.17327604 1.93608128 1.51674129 1.34798135
1.73725271 1.02645266 1.77222292 1.80086225]

可以看出，计算结果是正确的。

进一步修改

让我们继续修改前面的代码，同时说明其他一些有用的概念。让我们来看看线程几何的定义。我们已经看到，线程几何的管理可能很复杂，而且难以理解。对我们有利的是，还有另一种方法可以极大地方便我们：定义绝对位置。

事实上，非常简单的算法，比如我们例子中的算法，往往会以相同的方式使用所有线程索引。Numba 提供的函数大大方便了我们在内核中定义线程索引位置的工作。

通过 cuda.grid(ndim) 函数，我们可以知道当前线程在当前层次结构中的绝对位置。通过 ndim，我们可以定义内核所基于的维数。有了这个函数，内核函数代码就大大简化了：

@cuda.jit
def add_one(a):
pos = cuda.grid(1)
if pos < a.size:
a[pos] += 1
替换前面的代码并不会改变结果。

扩展到矩阵（二维数组）

我们可以做的另一个改动是将前面的示例扩展到二维数组的情况，因此我们将讨论矩阵而不是矢量。

要实现这种扩展，需要对代码进行一些修改。首先，我们可以定义一个矩阵 A，而不是向量 a。然后，我们进行适当的修改，将元素数减少到 4，以提高结果的可读性：

n = 4
A_host = np.random.random(n*n).reshape(n,n)
print(“Matrix A: \n %s” %A_host)
A_dev = cuda.to_device(A_host)
因此，我们在主机 A_host 和设备 A_dev 上都用一个 4×4 矩阵替换了之前的向量。然后，我们修改内核函数，使矩阵中的所有元素都能增加一个单位：

@cuda.jit
def add_one_2D(A):
x, y = cuda.grid(2)
if x < A.shape[0] and y < A.shape[1]:
A[x, y] += 1
有了绝对位置，定义二维核函数就变得非常容易了。至于线程层次结构，则需要一个新的几何结构，以便在处理二维向量时能做出最佳响应。因此，要为矩阵的每个元素分配一个线程，就必须定义如下新的层次结构：

threadsperblock = (128,128)
blockspergrid_x = math.ceil(A_host.shape[0] / threadsperblock[0])
blockspergrid_y = math.ceil(A_host.shape[1] / threadsperblock[1])
blockspergrid = (blockspergrid_x, blockspergrid_y)
add_one_2Dthreadsperblock, blockspergrid
最后，我们将把从设备获得的结果重新传输到主机，并在输出中显示出来：

运行代码后，我们将得到类似下面的结果：

Matrix A:
[[0.61817387 0.74884222 0.2364505 0.67749212]
[0.24190577 0.75021057 0.19960201 0.62161273]
[0.39070039 0.61669916 0.48010569 0.48047667]
[0.15454105 0.0797274 0.05984169 0.60568235]]
New Matrix A:
[[1.61817387 1.74884222 1.2364505 1.67749212]
[1.24190577 1.75021057 1.19960201 1.62161273]
[1.39070039 1.61669916 1.48010569 1.48047667]
[1.15454105 1.0797274 1.05984169 1.60568235]]
通过队列传输数据

另一项适用于所有示例（矢量和矩阵）的修改是使用队列将数据传输到 GPU 或从 GPU 传输到队列。编辑非常简单，只需创建一个 CUDA 流即可：stream = cuda.stream()

然后，将其添加到管理数据传输的两个函数中：

A_dev = cuda.to_device(A_host, stream=stream)
A_host = A_dev.copy_too_host(stream=stream)
这样，您就可以通过特定队列管理更高效的数据传输，在主机和设备之间正确、高效地管理数据传输。

两个矩阵之间的求和

我们进一步扩展了矩阵的情况，引入了多个二维向量之间的操作。在这种情况下，我们将把两个矩阵 A 和 B 的元素相加，得到第三个矩阵。这个例子展示了如何向内核引入多个输入向量，以及如何在代码的其余部分对它们进行管理。

现在，我们需要定义的不再是单一向量（矩阵），而是三个，其中一个是空的，因为它将存放计算结果：

n = 4
A_host = np.random.random(nn).reshape(n,n) B_host = np.random.random(nn).reshape(n,n)
C_host = np.zeros((n,n))
print(“Matrix A: \n %s” %A_host)
print(“Matrix B: \n %s” %B_host)
关于从主机到设备的数据传输，以及从设备到主机的数据传输，必须分别进行：

stream = cuda.stream()
A_dev = cuda.to_device(A_host, stream=stream)
B_dev = cuda.to_device(B_host, stream=stream)
C_dev = cuda.to_device(C_host, stream=stream)
至于从设备传输到主机的数据，在本例中，我们只需要传输 C_dev，因为只有它包含内核计算的结果。A_dev 和 B_dev 用于只读，并且与主机上的相同，因此不需要再次复制 A_dev 和 B_dev：

C_host = C_dev.copy_to_host(stream=stream)
print(“Matrix C: \n %s” %C_host)
既然我们已经管理好了二维矢量，那么接下来就需要对内核进行一些改动了，尽管一切仍然非常简单：

@cuda.jit
def add_one_2D(A,B,C):
x, y = cuda.grid(2)
if x < A.shape[0] and y < A.shape[1]:
C[x,y] = A[x, y] + B[x, y]
在上例中，线程的位置保持不变，因为两个矩阵 A 和 B 的元素之间的计算对应关系保持不变。只有 A 和 B 的对应元素，即具有相同 x 和 y 的元素，才会被相加，而结果将是 C 的元素始终位于相同的位置。在这种情况下，只需将不同矩阵之间的元素之和相加即可：C[x,y] = A[x, y] + B[x, y]

就线程层次结构而言，一切保持不变。

这次需要三个参数，也就是三个矩阵 A、B 和 C：

add_one_2D[threadsperblock, blockspergrid](A_dev, B_dev, C_dev)
让我们重写修改后的全部代码：

运行该代码将产生类似下面的结果：

矩阵间的乘法运算

对之前代码的最后一处修改是考虑矩阵之间的乘法运算，例如 A * B = C。相反，我们会发现这种运算完全可以在 GPU 上进行处理，与 CPU 相比，其性能确实非常出色：

正如我们所看到的，如果我们将内核排除在外，所产生的差异是微乎其微的。运行前面的代码，我们会得到类似下面的结果：

Matrix A:
[[0.87770781 0.72362124 0.96485584 0.66377971]
[0.60321934 0.51501588 0.81098102 0.9142394 ]
[0.89418345 0.43791226 0.75379705 0.93821043]
[0.55172303 0.07030213 0.97757898 0.5739226 ]]
Matrix B:
[[0.50570982 0.50927776 0.65589514 0.65859473]
[0.11216391 0.78912563 0.43172356 0.96117076]
[0.66527361 0.23458927 0.17906562 0.8516118 ]
[0.27048168 0.24474365 0.30308361 0.4918986 ]]
Matrix C:
[[1.38341763 1.232899 1.62075098 1.32237445]
[0.71538325 1.3041415 1.24270458 1.87541016]
[1.55945706 0.67250153 0.93286268 1.78982223]
[0.8222047 0.31504578 1.28066259 1.0658212 ]]
(base) [root@gczxagentx1 Chapter 6]# python numba_04.py
Matrix A:
[[0.89729956 0.41532718 0.0085074 0.14673637]
[0.1547752 0.28656798 0.13170739 0.56613537]
[0.42785439 0.01700211 0.84011128 0.05681843]
[0.71038642 0.27992718 0.32874299 0.14386456]]
Matrix B:
[[0.80083882 0.50297449 0.43292994 0.46281174]
[0.57941384 0.45732741 0.36802008 0.42373143]
[0.84677056 0.61206291 0.95717158 0.47114105]
[0.92708919 0.330491 0.63486587 0.84095001]]
Matrix C:
[[1.10248016 0.6949614 0.64261754 0.71867408]
[0.92637536 0.47661922 0.65795618 0.73120394]
[1.11655091 0.7559543 1.03169081 0.64881273]
[1.14284387 0.73408195 0.81656452 0.72325634]]

我们选择的示例只需要很少的执行时间，因为我们选择的是小尺寸矩阵。但随着矩阵大小的增加，情况很快就会不同。现在，我们将通过时间模块添加时间计量器。例如，我们将看到这种操作的执行时间是如何随着每块的线程数变化而变化的。让我们修改代码来记录执行时间：

从代码中我们可以看到，我们取消了矩阵值的打印，并在程序执行过程中从控制台插入了矩阵的大小，这样就不必每次都修改代码了。另一种解决方案是在程序启动时使用输入参数。

如果我们测量矩阵大小和区块大小（每个区块的线程数）变化所耗费的时间，就会得到类似下表的结果：

您可以看到计算矩阵乘法所用时间随矩阵大小增加而变化的趋势。报告显示了三种趋势，每种趋势每个区块的线程数都不同。可以看出，线程数越多，GPU 性能越好。

PyOpenCL

另一个允许我们通过 Python 开发 GPU 潜力的库是 pyOpenCL。该软件包是开放计算语言（OpenCL）的 Python 封装包。该框架基于 C99 语言的使用，允许您在不同的异构平台（包括 CPU 和 GPU）上工作。然后，使用 OpenCL 实现的程序将由适当的驱动程序编译，这些驱动程序使用针对每种技术的不同语言。因此，它的使用并不局限于英伟达™（NVIDIA®）显卡和 CUDA 的使用。因此，OpenCL 是将编程扩展到许多其他不同环境的绝佳解决方案，这些环境使用的是 IBM、AMD 或英特尔等不同制造商生产的 GPU 卡。

OpenCL 最初是由苹果公司创建的，但后来该项目转交给了非营利联盟 Khronos Group。由于 OpenCL 并不完全基于 CUDA，因此对于没有英伟达™（NVIDIA®）显卡但又想利用不同制造商生产的其他 GPU 以同样方式进行并行计算的人来说，这是一个极佳的替代方案。

安装 pyOpenCL

curl -L -O “https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-$(uname)-$(uname -m).sh”
bash ./Miniforge3-*.sh

(answer questions, pick install location)

source /WHERE/YOU/INSTALLED/MINIFORGE/bin/activate root
conda install pyopencl

CPU-based OpenCL driver

conda install pocl
conda install intel-opencl-rt
conda install intel-compute-runtime
conda install oclgrind
在系统上完成这些软件包的安装后，可以通过下面的示例检查是否一切顺利。在这方面，get_platforms() 是一个非常有用的函数，可以帮助我们调查正在运行的系统。该函数为我们提供了系统中能使用 OpenCL 的平台列表，我们还可以使用 get_devices()函数进一步获取其中设备（GPU）的信息。然后，我们使用这两个函数编写以下代码，以了解我们正在运行的系统的总体情况：

Vendor: Intel(R) Corporation
Version: OpenCL 3.0 LINUX
Device Intel(R) Xeon(R) Platinum 8352V CPU @ 2.10GHz
Max Clock Speed: 2100 MHz
Compute Units: 72
Local Memory: 32.0 Kb
Global Memory: 251.4 Gb
Platform Intel(R) FPGA Emulation Platform for OpenCL(TM)
Vendor: Intel(R) Corporation
Version: OpenCL 1.2 Intel(R) FPGA SDK for OpenCL(TM), Version 20.3
Device Intel(R) FPGA Emulation Device
Max Clock Speed: 2100 MHz
Compute Units: 72
Local Memory: 256.0 Kb
Global Memory: 251.4 Gb
Platform Portable Computing Language
Vendor: The pocl project
Version: OpenCL 3.0 PoCL 6.0 Linux, Release, RELOC, SPIR-V, LLVM 15.0.7, SLEEF, DISTRO, REMOTE, CUDA, POCL_DEBUG
Device cpu-skylake-avx512-Intel(R) Xeon(R) Platinum 8352V CPU @ 2.10GHz
Max Clock Speed: 3500 MHz
Compute Units: 72
Local Memory: 1280.0 Kb
Global Memory: 249.4 Gb
Device NVIDIA RTX A6000
Max Clock Speed: 1800 MHz
Compute Units: 84
Local Memory: 48.0 Kb
Global Memory: 47.32 Gb
我们可以看到，我们得到了系统中与 OpenCL 兼容的平台列表，以及这些平台中可以使用的设备。在前面的例子中，我们在主板上集成了英伟达（NVIDIA）和英特尔（Intel）两块显卡。

PyOpenCL 编程模型

使用 pyopencl 创建的程序在结构上必须非常精确，以便考虑到所使用系统的特定架构（如 GPU）和底层概念。因此，在代码中，我们必须始终提及主机（CPU）和设备（GPU）。

首先，必须开发程序，使其能够在主机上运行。然后，对程序进行编码和编译，确定要用于设备的不同部分。这将生成多个可执行程序包，每个程序包都是根据其适用设备的规则编译的（不同制造商生产的型号有不同的驱动程序），然后在运行时调度执行。这些特定的软件包由内核识别，内核是在设备上执行的代码部分。

由 OpenGL 框架管理的系统可以看作是分布式系统（我们在上一章中已经介绍过），其中异构处理器通过交换内核（而不是信息）相互连接。因此，与分布式系统一样，也需要一个系统来管理这些内核的交换和调度。因此，我们在 OpenCL 中使用了命令队列（Command Queue），这是一种数据结构，用于收集和积累设备接收到的不同内核。该结构的任务是管理各种内核的执行顺序。此外，作为 OpenCL 编程模型的一部分，还定义了一个上下文（Context），用于识别一组特定的设备，方便管理主机-设备系统内的内核和数据传输：

使用 pyOpenCL 开发程序

使用 pyOpenCL 实现程序必须遵循 OpenCL 编程模型的指导原则。由于 OpenCL 主要基于对 NumPy 数组的操作，我们可以从一个简单的例子开始，将其中一个向量的元素加倍。这一操作将在内核中实现，然后编译到 GPU（设备）上运行。

此时，在 GPU 编程模型中，第一步将是定义主机中的内存区域和其中包含的数据。在 OpenCL 中，这非常简单，只需使用 NumPy 数组即可。在我们这个简单的例子中，我们可以定义几个向量：a_host 是一个包含 4 个元素的向量，其中包含计算值；res_host 是一个向量，始终包含 4 个元素，其中必须包含计算结果：

a_host = np.array([0.1, 1.4, 2.3, 1.7])
a_host = a_host.astype(np.float32)
res_host = np.zeros((4), dtype=np.float32
下一步是定义内核调度管理对象：使用 Context() 构造函数创建上下文，使用 CommandQueue() 构造函数创建队列：

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)
由于我们要处理的是 GPU，因此必须明确管理设备上的内存分配以及从主机内存到设备内存的数据传输。OpenCL 通过 Buffer 类以简单的方式实现了这一点。然后，我们定义一个缓冲区，负责将 a_host 中包含的值传输到 GPU 设备的内存区域，该内存区域由变量 a_dev 标识。由于该内存区域只有读取功能，因此将其定义为 READ_ONLY。包含结果的向量也是如此。我们将创建一个与 a_dev 大小相同的 res_dev 向量，但这个用于包含结果的内存将是 WRITE_ONLY 内存：

mf = cl.mem_flags
a_dev = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a_host)
res_dev = cl.Buffer(ctx, mf.WRITE_ONLY, a_host.nbytes)
现在，我们已经完成了设备中的内存分配和接下来的数据传输，可以继续实现内核了。在 pyOpenCL 中，内核代码是作为 Program() 类构造函数的参数插入的：

prg = cl.Program(ctx, “””
__kernel void doubled(
__global const float *a_dev, __global float *res_dev)
{
int gid = get_global_id(0);
res_dev[gid] = a_dev[gid]*2;
}
“””).build()
我们可以看到，内核具有经典函数的结构，其中参数代表输入和返回值。事实上，内核不能有返回值，只能对作为参数传递的向量进行操作。一旦定义了内核，就必须调用它才能执行。一旦内核被声明和调用，它就会被分配给一个队列，该队列不仅负责将内核分派到设备，还负责使用 enqueue_copy()方法将包含结果的内存区域从设备传输到主机：

knl = prg.doubled
knl(queue, a_host.shape, None, a_dev, res_dev)
cl.enqueue_copy(queue, res_host, res_dev)
现在结果已经复制到主机的 res_host 数组中，可以输出得到的结果了：

print(“Vector a: %s” %a_host)
print(“Result: %s” %res_host)

运行前面的代码将得到类似下面的结果：

Choose platform:
[0]
[1]
[2]
Choice [0]:0
Choose device(s):
[0]
[1]
Choice, comma-separated [0]:0
Set the environment variable PYOPENCL_CTX=’0:0′ to avoid being asked again.
Vector a: [0.1 1.4 2.3 1.7]
Result: [0.2 2.8 4.6 3.4]
您可以立即看到，首先会询问您在哪个环境中工作（如果您的系统有多个可以使用 OpenCL 技术的平台，就会出现这种情况）。请输入您想要的值。如果按回车键，程序将在默认平台（值 0）上执行计算。为了稍后对这一选择进行补救，您可以在代码中插入环境变量，程序也会这样建议。因此，我们在前面代码的初始部分插入以下指令：

import os
os.environ[‘PYOPENCL_CTX’] = ‘0’
这样，我们就可以告诉编译器，我们希望在 NVIDIA 平台上运行程序。修改后重新运行程序，我们将直接得到想要的结果：

Vector a: [0.1 1.4 2.3 1.7]
Result: [0.2 2.8 4.6 3.4]
参考资料

软件测试精品书籍文档下载持续更新 https://github.com/china-testing/python-testing-examples 请点赞，谢谢！
本文涉及的python测试开发库
谢谢点赞！ https://github.com/china-testing/python_cn_resouce
python精品书籍下载
https://github.com/china-testing/python_cn_resouce/blob/main/python_good_books.md
Linux精品书籍下载 https://www.cnblogs.com/testing-/p/17438558.html
https://www.tutorialspoint.com/cuda/cuda_introduction_to_the_gpu.htm
矩阵间的乘法：示例

正如我们前面所看到的，使用 Numba 库，矩阵乘法是衡量 GPU 性能的一个很好的例子。这次让我们使用 pyOpenCL 库来实现它。例如，我们可以使用下面的代码将两个维度为 8 的矩阵 A * B 相乘，并得到包含运算结果的第三个矩阵 C：

import pyopencl as cl
import numpy as np

import os
os.environ[‘PYOPENCL_CTX’] = ‘0’

n = 8

a = np.random.randint(10, size=(nn)) b = np.random.randint(10, size=(nn))
c = np.zeros((n*n), dtype=np.float32)

a = a.astype(np.float32)
b = b.astype(np.float32)

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

mf = cl.mem_flags
a_buf = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a)
b_buf = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=b)
c_buf = cl.Buffer(ctx, mf.WRITE_ONLY, c.nbytes)

prg = cl.Program(ctx, “””
__kernel void multiply(ushort n,
ushort m, ushort p, __global float *a,
__global float *b, __global float *c)
{
int gid = get_global_id(0);
c[gid] = 0.0f;
int rowC = gid/p;
int colC = gid%p;
__global float *pA = &a[rowC*m];
__global float *pB = &b[colC];
for(int k=0; k<m; k++)
{
pB = &b[colC+kp]; c[gid] += ((pA++))(pB);
}
}
“””).build()

prg.multiply(queue, c.shape, None,
np.uint16(n), np.uint16(n), np.uint16(n),
a_buf, b_buf, c_buf)

cl.enqueue_copy(queue, c, c_buf)

print(“Matrix A”)
print(a.reshape(n,n))
print(“Matrix B”)
print(b.reshape(n,n))
print(“Matrix A*B”)
print(c.reshape(n, n))

运行这段代码将得到如下结果：

Choose device(s):
[0]
[1]
Choice, comma-separated [0]:0
Set the environment variable PYOPENCL_CTX=’0:0′ to avoid being asked again.
Matrix A
[[2. 1. 7. 9. 8. 5. 1. 2.]
[5. 0. 8. 9. 4. 7. 9. 9.]
[1. 6. 0. 3. 4. 2. 8. 1.]
[1. 8. 7. 3. 6. 2. 4. 9.]
[0. 2. 2. 7. 2. 7. 9. 3.]
[5. 4. 6. 8. 9. 8. 9. 1.]
[4. 7. 1. 1. 6. 0. 4. 2.]
[3. 3. 9. 3. 1. 7. 9. 3.]]
Matrix B
[[5. 3. 2. 1. 4. 6. 9. 4.]
[6. 9. 4. 4. 6. 4. 0. 1.]
[7. 6. 5. 1. 6. 5. 2. 9.]
[6. 5. 1. 6. 6. 4. 5. 4.]
[0. 8. 4. 1. 1. 0. 7. 9.]
[2. 7. 1. 7. 5. 4. 1. 7.]
[7. 3. 0. 4. 3. 4. 8. 0.]
[0. 8. 2. 6. 6. 5. 1. 6.]]
Matrix A*B
[[136. 220. 93. 126. 158. 121. 148. 227.]
[212. 288. 100. 210. 242. 215. 222. 267.]
[119. 150. 49. 99. 102. 87. 119. 78.]
[152. 278. 116. 148. 194. 154. 123. 209.]
[145. 181. 46. 157. 148. 125. 135. 133.]
[218. 290. 110. 182. 210. 181. 241. 253.]
[103. 162. 70. 73. 100. 87. 119. 102.]
[191. 213. 83. 146. 183. 166. 149. 184.]]

现在，让我们对之前的程序做一些修改，以引入对计算矩阵间乘积所用时间的测量：

import pyopencl as cl
import numpy as np
import time

import os
os.environ[‘PYOPENCL_CTX’] = ‘0’

n = 8

a = np.random.randint(10, size=(nn)) b = np.random.randint(10, size=(nn))
c = np.zeros((n*n), dtype=np.float32)

a = a.astype(np.float32)
b = b.astype(np.float32)

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

mf = cl.mem_flags
a_buf = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a)
b_buf = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=b)
c_buf = cl.Buffer(ctx, mf.WRITE_ONLY, c.nbytes)

t1 = time.perf_counter()
prg.multiply(queue, c.shape, None,
np.uint16(n), np.uint16(n), np.uint16(n),
a_buf, b_buf, c_buf)
t2 = time.perf_counter()

cl.enqueue_copy(queue, c, c_buf)
elapsed_time = t2 – t1
print(“Elapsed time %s” %elapsed_time)

我们取消了打印参与计算的矩阵的输出结果，取而代之的是插入时间模块的 perf_counter()函数，以计算矩阵间乘积所用的时间。

运行前面的代码，我们将得到类似下面的结果：

Elapsed time 0.007469899996067397
也就是说，我们获得了以秒为单位的时间输出。我们在计算机上进行一系列测量，每次执行时都修改 n 的值，例如将其赋值为 2 的幂次（4、8、16、32、64 等）。我们将得到类似的数值表：

从结果中我们可以看出，1024 以内所用时间几乎相同，然后随着 2 的幂次增加而增加。这让我们猜测这一耗时与 GPU 可用线程数有关。在达到一定值时，我们将有 n 个并行线程，它们计算矩阵（n x n）之间的乘积所需的时间相同。因此，根据这一机制，我们计算两个矩阵（4 x 4）之间的乘法所需的时间将与计算两个矩阵（512 x 512）之间的乘法所需的时间相同。

为了更好地理解 GPU 编程和计算所带来的潜力，让我们考虑一下两个矩阵之间按顺序进行的乘法运算。下面的代码可以作为一个例子，说明如何在不使用 GPU 的经典 Python 程序上实现两个矩阵 A 和 B 之间的乘积：

import numpy as np
import time

n = 8
a = np.random.randint(10, size=(nn)) b = np.random.randint(10, size=(nn))
c = np.zeros((n*n), dtype=np.float32)

a = a.astype(np.float32)
b = b.astype(np.float32)

def matmul(n, A, B, C):
# Iterate over the rows of A.
for i in range(n):
# Iterate over the columns of B
for j in range(n):
tmp = 0.0
# Iterate over the rows of B.
for k in range(n):
tmp += A[i * n + k] * B[k * n + j]
C[i * n + j] = tmp

matmul(n, a, b, c)
print(“Matrix A”)
print(a.reshape(n,n))
print(“Matrix B”)
print(b.reshape(n,n))
print(“Matrix A*B”)
print(c.reshape(n,n))
通过执行该代码，我们将得到与上一个使用 pyOpenCL 的示例相同的结果。数值是随机生成的，因此每次运行的结果都会有所不同，但计算方法和结果实际上是相同的。

此时，我们还修改了这段代码，以便能够测量矩阵之间的乘法运算随着矩阵大小的变化以及 n 值的变化所耗费的时间：

import numpy as np
import time

n = 8
a = np.random.randint(10, size=(nn)) b = np.random.randint(10, size=(nn))
c = np.zeros((n*n), dtype=np.float32)

a = a.astype(np.float32)
b = b.astype(np.float32)

t1 = time.perf_counter()
matmul(n, a, b, c)
t2 = time.perf_counter()
elapsed_time = t2 – t1
print(“Elapsed time %s” %elapsed_time)
通过对 pyOpenCL 代码中使用的 n 值执行一系列执行，我们将得到一系列时间：

我们可以看到，在这种情况下，时间的趋势与前一种情况完全不同。随着所用矩阵大小的增加，时间呈指数级增长，以至于超过一定大小后，时间变得非常紧迫。一个有趣的现象是，对于小数组（小于 32 个），通过 pyOpenCL 进行顺序计算比使用 GPU 计算要高效得多。下图（图 6.6）清楚地显示了两种编程方式的时间趋势：

从图中我们可以看到，对于小矩阵（最大阶数为 16），在 CPU 上进行的计算矩阵间乘积的顺序计算比在 GPU 上使用线程并行计算的效率更高。但在 CPU 上，计算时间的变化趋势是指数级的，而且对于更大尺寸的矩阵，计算时间会越来越长，令人望而却步。相反，我们在 GPU 上的表现却完全不同，尽管用于乘法运算的矩阵大小不断增加，但所需时间却保持不变。这主要归功于 GPU 上的大量线程，它们能够同时高效处理数千个线程的并行计算。

使用 pyopenCL 进行元素计算

PyOpenCL 提供了以隐式方式执行逐元素计算的可能性，并可在单个计算步骤中使用复杂的表达式。使用的函数是 ElementwiseKernel()。

让我们来看一个例子：

import numpy as np
import pyopencl as cl
import pyopencl.array
from pyopencl.elementwise import ElementwiseKernel

import os
os.environ[‘PYOPENCL_CTX’] = ‘0’

a_host = np.array([0.1 , 1.4, 2.3, 1.7])
a_host = a_host.astype(np.float32)
b_host = np.array([0.2 , 0.3, 1.0, 0.5])
b_host = b_host.astype(np.float32)

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

a_dev = cl.array.to_device(queue, a_host)
b_dev = cl.array.to_device(queue, b_host)

sum = ElementwiseKernel(ctx,
“float *a_dev, float *b_dev, float *res_dev”,
“res_dev[i] = a_dev[i] + b_dev[i]”,
“sum”)

res_dev = cl.array.empty_like(a_dev)
sum(a_dev, b_dev, res_dev)

print(“Vector A: %s” %a_host)
print(“Vector B: %s” %b_host)
print(“Vector Result: %s” %res_dev)
执行结果：

Vector A: [0.1 1.4 2.3 1.7]
Vector B: [0.2 0.3 1. 0.5]
Vector Result: [0.3 1.7 3.3 2.2]
使用 pyOpenCL 进行 MapReduce 计算

和前面的例子一样，还有一个构造允许隐式使用复杂表达式，比如针对 MapReduce 的表达式。PyOpenCL 为此提供了 ReductionKernel()：

import numpy as np
import pyopencl as cl
import pyopencl.array
import pyopencl.reduction

import os
os.environ[‘PYOPENCL_CTX’] = ‘0’

ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

ah = np.array([0.1 , 1.4, 2.3, 1.7])
ah = ah.astype(np.float32)
bh = np.array([0.2 , 0.3, 1.0, 0.5])
bh = bh.astype(np.float32)

a = cl.array.to_device(queue, ah)
b = cl.array.to_device(queue, bh)

krnl = pyopencl.reduction.ReductionKernel(ctx, np.float32, neutral=”0″,
reduce_expr=”a+b”, map_expr=”x[i]*y[i]”,
arguments=”__global float *x, __global float *y”)

res = krnl(a, b).get()

print(“Vector A: %s” %a)
print(“Vector B: %s” %b)
print(“Vector Result: %s” %res)
执行结果：

Vector A: [0.1 1.4 2.3 1.7]
Vector B: [0.2 0.3 1. 0.5]
Vector Result: [0.3 1.7 3.3 2.2]
结论
在本章中，我们了解了如何使用 GPU 编程来进一步增强并行计算。得益于 Numba 和 pyOpenCL 等库，我们现在能够利用 GPU 的计算潜力；得益于其特殊的内部架构，即使是复杂的计算也能并行执行。特别是，这些库都基于 CUDA 技术，该技术专门用于英伟达™（NVIDIA®）显卡，能最大限度地发挥其性能。然后，我们展示了作为 CUDA 基础的编程模型，这些模型允许我们在 Python 代码中插入内核，即编译成专门在 GPU 上执行的代码部分。通过本章，我们总结了所有并行编程模型和相关库。在下一章中，我们将分析前几章中涉及的这些技术和模型的所有可能应用领域。

声明：文中观点不代表本站立场。本文传送门：https://eyangzhen.com/424673.html

Python并行计算与高性能计算6使用 CUDA 通过 GPU 编程实现性能最大化

python numba_01.py

(answer questions, pick install location)

CPU-based OpenCL driver

作者专栏

pythontesting