标签：Python 编程 线程 CUDA gpu GPU 我们

Python GPU 编程实用指南（一）

原文：zh.annas-archive.org/md5/ef7eb3c148e0cfdfe01c331f2f01557c

译者：飞龙

协议：CC BY-NC-SA 4.0

前言

问候和祝福！本文是关于使用 Python 和 CUDA 进行 GPU 编程的入门指南。GPU可能代表图形编程单元，但我们应该明确，这本书不是关于图形编程——它本质上是通用 GPU 编程的介绍，简称为GPGPU 编程。在过去的十年中，已经清楚地表明 GPU 除了渲染图形之外，也非常适合进行计算，特别是需要大量计算吞吐量的并行计算。为此，NVIDIA 发布了 CUDA 工具包，这使得几乎任何具有一些 C 编程知识的人都可以更轻松地进入 GPGPU 编程的世界。

《使用 Python 和 CUDA 进行 GPU 编程实践》的目标是尽快让您进入 GPGPU 编程的世界。我们努力为每一章设计了有趣和有趣的示例和练习；特别是，我们鼓励您在阅读时在您喜欢的 Python 环境中键入这些示例并运行它们（Spyder、Jupyter 和 PyCharm 都是合适的选择）。这样，您最终将学会所有必要的函数和命令，并获得编写 GPGPU 程序的直觉。

最初，GPGPU 并行编程似乎非常复杂和令人望而生畏，特别是如果您过去只做过 CPU 编程。有很多新概念和约定您必须学习，可能会让您觉得自己从零开始。在这些时候，您必须相信学习这个领域的努力不是徒劳的。通过一点主动性和纪律性，当您阅读完本文时，这个主题将对您来说变得轻而易举。

愉快的编程！

这本书是为谁准备的

这本书特别针对一个人，那就是我自己，2014 年，当时我正在尝试为数学博士研究开发基于 GPU 的模拟。我正在研究多本关于 GPU 编程的书籍和手册，试图对这个领域有一点了解；大多数文本似乎乐意在每一页向读者抛出无尽的硬件原理和术语，而实际的编程则被放在了次要位置。

这本书主要面向那些想要实际进行 GPU 编程，但不想陷入技术细节和硬件原理的人。在本书中，我们将使用适当的 C/C++（CUDA C）来编程 GPU，但我们将通过 PyCUDA 模块将其内联到 Python 代码中。PyCUDA 允许我们只编写我们需要的必要低级 GPU 代码，同时自动处理编译、链接和在 GPU 上启动代码的所有冗余工作。

本书涵盖的内容

第一章《为什么要进行 GPU 编程？》给出了一些我们应该学习这个领域的动机，以及如何应用阿姆达尔定律来估计将串行程序转换为利用 GPU 的潜在性能改进。

第二章《设置 GPU 编程环境》解释了如何在 Windows 和 Linux 下设置适当的 Python 和 C++开发环境以进行 CUDA 编程。

第三章《使用 PyCUDA 入门》展示了我们在使用 Python 编程 GPU 时最基本的技能。我们将特别看到如何使用 PyCUDA 的 gpuarray 类将数据传输到 GPU 和从 GPU 传输数据，以及如何使用 PyCUDA 的 ElementwiseKernel 函数编译简单的 CUDA 核函数。

[第四章]（5a5f4317-50c7-4ce6-9d04-ac3be4c6d28b.xhtml），核心，线程，块和网格，教授了编写有效的 CUDA 核心的基础知识，这些核心是在 GPU 上启动的并行函数。我们将看到如何编写 CUDA 设备函数（由 CUDA 核心直接调用的“串行”函数），并了解 CUDA 的抽象网格/块结构及其在启动核心中的作用。

[第五章]（ea648e20-8c72-44a9-880d-11469d0e291f.xhtml），流，事件，上下文和并发，涵盖了 CUDA 流的概念，这是一种允许我们在 GPU 上同时启动和同步许多内核的功能。我们将看到如何使用 CUDA 事件来计时内核启动，以及如何创建和使用 CUDA 上下文。

[第六章]（6d1c808f-1dc2-4454-b0b8-d0a36bc3c908.xhtml），调试和分析您的 CUDA 代码，填补了我们在纯 CUDA C 编程方面的一些空白，并向我们展示了如何使用 NVIDIA Nsight IDE 进行调试和开发，以及如何使用 NVIDIA 分析工具。

[第七章]（55146879-4b7e-4774-9a8b-cc5c80c04ed8.xhtml），使用 CUDA 库与 Scikit-CUDA，通过 Python Scikit-CUDA 模块简要介绍了一些重要的标准 CUDA 库，包括 cuBLAS，cuFFT 和 cuSOLVER。

[第八章]（d374ea77-f9e5-4d38-861d-5295ef3e3fbf.xhtml），CUDA 设备函数库和 Thrust，向我们展示了如何在我们的代码中使用 cuRAND 和 CUDA Math API 库，以及如何使用 CUDA Thrust C++容器。

[第九章]（3562f1e0-a53d-470f-9b4d-94fa41b1b2fa.xhtml），实现深度神经网络，作为一个巅峰，我们将学习如何从头开始构建整个深度神经网络，应用我们在文本中学到的许多想法。

[第十章]（5383b46f-8dc6-4e17-ab35-7f6bd35f059f.xhtml），使用已编译的 GPU 代码，向我们展示了如何将我们的 Python 代码与预编译的 GPU 代码进行接口，使用 PyCUDA 和 Ctypes。

[第十一章]（e853faad-3ee4-4df7-9cdb-98f74e435527.xhtml），CUDA 性能优化，教授了一些非常低级的性能优化技巧，特别是与 CUDA 相关的技巧，例如 warp shuffling，矢量化内存访问，使用内联 PTX 汇编和原子操作。

[第十二章]（2d464c61-de29-49fa-826a-a7437c368d6a.xhtml），从这里出发，概述了您将拥有的一些教育和职业道路，这些道路将建立在您现在扎实的 GPU 编程基础之上。

为了充分利用本书

这实际上是一个非常技术性的主题。为此，我们将不得不对读者的编程背景做一些假设。为此，我们将假设以下内容：

• 您在 Python 中具有中级的编程经验。

• 您熟悉标准的 Python 科学包，如 NumPy，SciPy 和 Matplotlib。

• 您具有任何基于 C 的编程语言的中级能力（C，C ++，Java，Rust，Go 等）。

• 您了解 C 中的动态内存分配概念（特别是如何使用 C 的malloc和free函数）。

GPU 编程主要适用于非常科学或数学性质的领域，因此许多（如果不是大多数）示例将利用一些数学知识。因此，我们假设读者对大学一年级或二年级的数学有一定了解，包括：

• 三角学（正弦函数：sin，cos，tan…）

• 微积分（积分，导数，梯度）

• 统计学（均匀分布和正态分布）

• 线性代数（向量，矩阵，向量空间，维度）。

如果您还没有学习过这些主题，或者已经有一段时间了，不要担心，因为我们将在学习过程中尝试回顾一些关键的编程和数学概念。

我们在这里将做另一个假设。请记住，本文中我们只会使用 CUDA，这是 NVIDIA 硬件的专有编程语言。因此，在开始之前，我们需要拥有一些特定的硬件。因此，我假设读者可以访问以下内容：

• 64 位 x86 英特尔/AMD PC

• 4 GB 或更多的 RAM

• 入门级 NVIDIA GTX 1050 GPU（Pascal 架构）或更高版本

读者应该知道，大多数旧的 GPU 可能会在本文中的大多数示例中正常工作，但本文中的示例仅在 Windows 10 下的 GTX 1050 和 Linux 下的 GTX 1070 上进行了测试。有关设置和配置的具体说明在第二章中给出，设置您的 GPU 编程环境。

下载示例代码文件

您可以从www.packt.com的帐户中下载本书的示例代码文件。如果您在其他地方购买了本书，可以访问www.packt.com/support并注册，文件将直接发送到您的邮箱。

您可以按照以下步骤下载代码文件：

1. 登录或注册www.packt.com。

2. 选择“支持”选项卡。

3. 单击“代码下载和勘误”。

4. 在搜索框中输入书名，然后按照屏幕上的说明操作。

下载文件后，请确保使用最新版本的解压缩或提取文件夹：

• WinRAR/7-Zip for Windows

• Zipeg/iZip/UnRarX for Mac

• 7-Zip/PeaZip for Linux

该书的代码包也托管在 GitHub 上，网址为github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA。如果代码有更新，将在现有的 GitHub 存储库上进行更新。

我们还有来自我们丰富书籍和视频目录的其他代码包，可在github.com/PacktPublishing/上找到。去看看吧！

下载彩色图像

我们还提供了一个 PDF 文件，其中包含本书中使用的屏幕截图/图表的彩色图像。您可以在此处下载：www.packtpub.com/sites/default/files/downloads/9781788993913_ColorImages.pdf。

使用的约定

本书中使用了许多文本约定。

CodeInText：表示文本中的代码单词、数据库表名、文件夹名、文件名、文件扩展名、路径名、虚拟 URL、用户输入和 Twitter 句柄。例如：“我们现在可以使用cublasSaxpy函数。”

代码块设置如下：

cublas.cublasDestroy(handle)
print 'cuBLAS returned the correct value: %s' % np.allclose(np.dot(A,x), y_gpu.get())

当我们希望引起您对代码块的特定部分的注意时，相关行或项目将以粗体显示：

def compute_gflops(precision='S'):

if precision=='S':
    float_type = 'float32'
elif precision=='D':
    float_type = 'float64'
else:
    return -1

任何命令行输入或输出都以以下方式编写：

$ run cublas_gemm_flops.py

粗体：表示新术语、重要单词或屏幕上看到的单词。例如，菜单或对话框中的单词会以这种方式出现在文本中。

警告或重要说明会出现在这样的地方。提示和技巧会出现在这样的地方。

第一章：为什么要进行 GPU 编程？

事实证明，除了能够为视频游戏渲染图形外，图形处理单元（GPU）还为普通消费者提供了一种便捷的方式进行大规模并行 计算——现在普通人可以在当地的电子商店购买一张价值 2000 美元的现代 GPU 卡，将其插入家中的个人电脑，几乎立即就可以用于计算能力，而这种计算能力在 5 年或 10 年前只能在顶级公司和大学的超级计算实验室中获得。近年来，GPU 的开放可访问性在许多方面已经显而易见，这可以通过简要观察新闻来揭示——加密货币挖矿者使用 GPU 生成比特币等数字货币，遗传学家和生物学家使用 GPU 进行 DNA 分析和研究，物理学家和数学家使用 GPU 进行大规模模拟，人工智能研究人员现在可以编程 GPU 来撰写剧本和作曲，而主要的互联网公司，如谷歌和 Facebook，使用带有 GPU 的服务器农场进行大规模机器学习任务……等等。

本书主要旨在让您迅速掌握 GPU 编程，以便您也可以尽快开始使用它们的强大功能，无论您的最终目标是什么。我们旨在涵盖如何编程 GPU 的核心要点，而不是提供 GPU 工作的复杂技术细节和原理图。在本书的末尾，我们将提供更多资源，以便您可以进一步专门化，并应用您对 GPU 的新知识。（有关特定所需的技术知识和硬件的进一步细节，请参阅本节后面的内容。）

在本书中，我们将使用CUDA，这是 NVIDIA 推出的通用 GPU（GPGPU）编程框架，最早发布于 2007 年。虽然 CUDA 专为 NVIDIA GPU 而设计，但它是一个成熟稳定的平台，相对容易使用，提供了一套无与伦比的第一方加速数学和人工智能相关库，并且在安装和集成方面几乎没有麻烦。此外，还有现成的标准化 Python 库，如 PyCUDA 和 Scikit-CUDA，使得渴望成为 GPU 程序员的人更容易接触 GPGPU 编程。出于这些原因，我们选择在本书中使用 CUDA。

CUDA 始终发音为 coo-duh，而不是缩写 C-U-D-A！CUDA 最初代表“计算统一设备架构”，但 Nvidia 已经放弃了这个缩写，现在将 CUDA 作为一个大写的专有名词。

我们现在将开始介绍 GPU 编程的旅程，并概述阿姆达尔定律。阿姆达尔定律是一种简单但有效的方法，用于估计将程序或算法转移到 GPU 上可以获得的潜在速度增益；这将帮助我们确定是否值得重新编写我们的代码以利用 GPU。然后，我们将简要回顾如何使用cProfile模块对我们的 Python 代码进行分析，以帮助我们找到代码中的瓶颈。

本章的学习成果如下：

• 了解阿姆达尔定律

• 在代码的上下文中应用阿姆达尔定律

• 使用cProfile模块对 Python 代码进行基本分析

技术要求

本章建议安装 Anaconda Python 2.7：

www.anaconda.com/download/

本章的代码也可以在 GitHub 上找到：

github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA

有关先决条件的更多信息，请查看本书的前言；有关软件和硬件要求，请查看github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA中的 README 部分。

并行化和阿姆达尔定律

在我们深入了解并发解锁 GPU 的潜力之前，我们首先要意识到它们的计算能力相对于现代英特尔/AMD 中央处理单元（CPU）的优势并不在于它的时钟速度比 CPU 更高，也不在于单个核心的复杂性或特定设计。一个单独的 GPU 核心实际上相当简单，并且与现代单个 CPU 核心相比处于劣势，后者使用了许多花哨的工程技巧，比如分支预测来减少计算的延迟。延迟指的是执行单个计算的开始到结束的持续时间。

GPU 的强大之处在于它的核心比 CPU 多得多，这意味着吞吐量有了巨大的提升。这里的吞吐量指的是可以同时执行的计算数量。让我们使用一个类比来更好地理解这意思。GPU 就像一条非常宽的城市道路，设计成可以同时处理许多行驶缓慢的汽车（高吞吐量，高延迟），而 CPU 就像一条狭窄的高速公路，一次只能容纳几辆车，但可以更快地将每辆车送到目的地（低吞吐量，低延迟）。

我们可以通过查看这些新 GPU 有多少核心来了解吞吐量的增加。举个例子，普通的英特尔或 AMD CPU 只有 2 到 8 个核心，而入门级的消费级 NVIDIA GTX 1050 GPU 有 640 个核心，而新的顶级 NVIDIA RTX 2080 Ti 有 4,352 个核心！我们可以利用这种大规模的吞吐量，只要我们知道如何正确地并行化任何我们希望加速的程序或算法。通过并行化，我们的意思是重写程序或算法，以便我们可以将工作负载并行地在多个处理器上同时运行。让我们从现实生活中想一个类比。

假设你正在建造一座房子，你已经准备好了所有的设计和材料。你雇了一个劳工，你估计需要 100 个小时来建造这座房子。假设这个特定的房子可以以这样的方式建造，即每增加一个劳工，工作就可以完美地分配给他们，也就是说，两个劳工需要 50 个小时，四个劳工需要 25 个小时，十个劳工需要 10 个小时来建造这座房子——建造你的房子所需的时间将是 100 除以你雇佣的劳工数量。这是一个可并行化的任务的例子。

我们注意到，这个任务对于两个劳工来说完成的速度是原来的两倍，对于十个劳工来说完成的速度是原来的十倍（也就是说，并行完成），而不是一个劳工独自建造房子（也就是说，串行完成）——也就是说，如果N是劳工的数量，那么速度将是N倍。在这种情况下，N被称为我们的任务并行化速度的加速比。

在我们开始编写给定算法的并行版本之前，我们经常首先估计一下并行化对我们任务可能带来的潜在 加速。这可以帮助我们确定是否值得花费资源和时间来编写我们程序的并行版本。因为现实生活比我们在这里给出的例子更复杂，很明显我们不可能始终完美地并行化每个程序——大多数情况下，我们的程序只有一部分可以很好地并行化，而其余部分将不得不串行运行。

使用阿姆达尔定律

我们现在将推导阿姆达尔定律，这是一个简单的算术公式，用于估计将一部分串行程序代码并行化到多个处理器上可能带来的潜在速度增益。我们将继续使用我们之前建造房子的类比来做这件事。

上次，我们只考虑了房子的实际物理建造作为整个时间持续时间，但现在，我们还将把设计房子所需的时间考虑在内。假设世界上只有一个人有能力设计你的房子——也就是你——并且你需要 100 小时来设计你的房子的计划。世界上没有其他人能够与你的建筑才华相比，因此这部分任务无法在其他建筑师之间分配，因此无论你拥有什么资源或可以雇佣多少人，设计你的房子都需要 100 小时。因此，如果你只有一名劳工来建造你的房子，建造你的房子所需的整个时间将是 200 小时——你设计它需要 100 小时，一名劳工建造它需要 100 小时。如果我们雇佣两名劳工，这将需要 150 小时——设计房子的时间仍然是 100 小时，而建造将需要 50 小时。很明显，建造房子所需的总时间将是 100 + 100 / N，其中N是我们雇佣的劳工数量。

现在，让我们退一步思考一下，如果我们只雇用一名劳工来建造房子需要多少时间——我们最终使用这个来确定我们雇用额外劳工时的加速度；也就是说，这个过程变得快了多少倍。如果我们只雇用一名劳工，我们会发现设计和建造房子需要相同的时间——100 小时。因此，我们可以说，设计所花费的时间是.5（50%），建造房子所花费的时间也是.5（50%）——当然，这两部分加起来是 1，也就是 100%。当我们增加劳工时，我们想要与这个进行比较——如果我们有两名劳工，建造的时间减半，因此与我们任务的原始串行版本相比，这将花费.5 + .5/2 = .75（75%）的时间，原始任务的.75 x 200 小时是 150 小时，因此我们可以看到这是有效的。此外，我们可以看到，如果我们有N名劳工，我们可以使用公式.5 + .5 / N 来计算我们并行化的建造所需的时间百分比。

现在，让我们确定通过增加额外的劳工我们获得的加速度。如果有两名劳工，建造一座房子只需要 75%的时间，我们可以取.75 的倒数来确定我们并行化的加速度——也就是说，加速度将是 1 / .75，比我们只有一名劳工时快大约 1.33 倍。在这种情况下，我们可以看到，如果有N名劳工，加速度将是 1 / (.5 + .5 / N)。

我们知道，随着我们增加越来越多的劳工，.5 / N 会缩小到接近 0，因此我们可以看到在并行化这个任务时，你可以获得的加速度总是有一个上限，即 1 / (.5 + 0) = 2。我们可以将原始串行时间除以估计的最大加速度，以确定此任务将花费的绝对最短时间——200 / 2 = 100 小时。

我们刚刚应用的用于确定并行编程中加速度的原则被称为阿姆达尔定律。它只需要知道原始串行程序中可并行执行时间的比例，即p，以及我们可用的处理器核心数N。

在这种情况下，不可并行化代码的执行时间比例始终为1-p，因此我们只需要知道p。

我们现在可以使用阿姆达尔定律来计算加速度，如下所示：

总之，Amdahl's Law 是一个简单的公式，允许我们粗略（非常粗略）地估计一个可以至少部分并行化的程序的潜在加速。这可以提供一个大致的想法，即是否值得编写特定串行程序的并行版本，前提是我们知道我们可以并行化代码的比例（p），以及我们可以在其上运行并行化代码的核心数（N）。

Mandelbrot 集

我们现在准备看一个非常标准的并行计算示例，我们将在本文中稍后重新讨论——一个生成Mandelbrot 集图像的算法。让我们首先确切地定义我们的意思。

对于给定的复数c，我们为定义一个递归序列，其中和对于。如果|z[n]|随着n增加到无穷大仍然受到 2 的限制，那么我们将说c是 Mandelbrot 集的成员。

回想一下，我们可以将复数可视化为驻留在二维笛卡尔平面上，其中x轴代表实部分量，y 轴代表虚部分量。因此，我们可以很容易地用一个非常吸引人（并且众所周知）的图形来可视化 Mandelbrot 集。在这里，我们将在复数笛卡尔平面上用浅色表示 Mandelbrot 集的成员，用深色表示非成员，如下所示：

现在，让我们考虑如何在 Python 中生成这个集合。首先，我们必须考虑一些事情——因为显然我们无法检查每一个复数是否在 Mandelbrot 集中，我们必须选择一个特定的范围进行检查；我们必须确定我们将考虑每个范围内的多少点（宽度，高度）；以及我们将检查的|z[n]|的最大值（max_iters）。我们现在可以准备实现一个生成 Mandelbrot 集图形的函数——在这里，我们通过在串行中迭代图中的每一个点来实现这一点。

我们将首先导入 NumPy 库，这是一个我们在本文中将大量使用的数值库。我们在simple_mandelbrot函数中实现这里。我们首先使用 NumPy 的linspace函数生成一个将充当离散复平面的格点（接下来的代码应该相当简单）：

import numpy as np

def simple_mandelbrot(width, height, real_low, real_high, imag_low, imag_high, max_iters):

     real_vals = np.linspace(real_low, real_high, width)
     imag_vals = np.linspace(imag_low, imag_high, height)

     # we will represent members as 1, non-members as 0.

     mandelbrot_graph = np.ones((height,width), dtype=np.float32)

     for x in range(width):

         for y in range(height):

             c = np.complex64( real_vals[x] + imag_vals[y] * 1j  )           
             z = np.complex64(0)

             for i in range(max_iters):

                 z = z**2 + c

                 if(np.abs(z) > 2):
                     mandelbrot_graph[y,x] = 0
                     break

     return mandelbrot_graph

现在，我们想要添加一些代码来将 Mandelbrot 集的图像转储到 PNG 格式文件中，所以让我们在开头添加适当的标头：

from time import time
import matplotlib
# the following will prevent the figure from popping up
matplotlib.use('Agg')
from matplotlib import pyplot as plt

现在，让我们添加一些代码来生成 Mandelbrot 集并将其转储到文件中，并使用时间函数来计算这两个操作的时间：

if __name__ == '__main__':

     t1 = time()
     mandel = simple_mandelbrot(512,512,-2,2,-2,2,256, 2)
     t2 = time()
     mandel_time = t2 - t1

     t1 = time()
     fig = plt.figure(1)
     plt.imshow(mandel, extent=(-2, 2, -2, 2))
     plt.savefig('mandelbrot.png', dpi=fig.dpi)
     t2 = time()

     dump_time = t2 - t1

     print 'It took {} seconds to calculate the Mandelbrot graph.'.format(mandel_time)
     print 'It took {} seconds to dump the image.'.format(dump_time)

现在让我们运行这个程序（这也可以在 GitHub 存储库的文件夹1中的mandelbrot0.py文件中找到）：

生成 Mandelbrot 集大约需要 14.62 秒，转储图像大约需要 0.11 秒。正如我们所看到的，我们逐点生成 Mandelbrot 集；不同点的值之间没有相互依赖，因此，这是一个固有的可并行化函数。相比之下，转储图像的代码无法并行化。

现在，让我们从 Amdahl's Law 的角度来分析这个问题。如果我们在这里并行化我们的代码，我们可以得到什么样的加速？总的来说，程序的两部分共计大约需要 14.73 秒才能运行；因为我们可以并行化 Mandelbrot 集的生成，我们可以说可并行化代码的执行时间部分是 p = 14.62 / 14.73 = .99。这个程序有 99%的可并行性！

我们可能会得到什么样的加速？嗯，我目前正在使用一台配有 640 个核心的入门级 GTX 1050 GPU 的笔记本电脑；因此，当我们使用这个公式时，我们的N将是 640。我们计算速度提升如下：

这绝对非常好，这表明我们值得努力编程使我们的算法使用 GPU。请记住，阿姆达尔定律只是一个非常粗略的估计！当我们将计算卸载到 GPU 时，将会有其他考虑因素，比如 CPU 发送和接收数据到 GPU 的额外时间；或者卸载到 GPU 的算法只能部分并行化。

对代码进行性能分析

在前面的例子中，我们看到我们可以使用 Python 中的标准time函数来分别计时不同的函数和组件。虽然这种方法对我们的小例子程序效果很好，但对于调用许多不同函数的大型程序来说，这种方法并不总是可行，其中一些函数可能值得我们投入精力并行化，或者甚至在 CPU 上进行优化。我们的目标是找到程序的瓶颈和热点，即使我们在每个函数调用周围使用time，我们可能会错过一些东西，或者可能有一些系统或库调用我们甚至没有考虑到，这些调用可能会拖慢速度。在我们考虑重写代码在 GPU 上运行之前，我们必须始终遵循著名的美国计算机科学家唐纳德·克努斯的智慧话语：过早优化是万恶之源。

我们使用所谓的分析器来找到代码中的热点和瓶颈。分析器将方便地让我们看到程序花费最多时间的地方，并让我们相应地进行优化。

使用 cProfile 模块

我们将主要使用cProfile模块来检查我们的代码。这个模块是一个标准库函数，包含在每个现代 Python 安装中。我们可以从命令行运行分析器，使用-m cProfile，并指定我们要按每个函数花费的累积时间来组织结果，然后用>操作符将输出重定向到文本文件中。

这将在 Linux Bash 或 Windows PowerShell 命令行上都可以使用。

现在让我们试一试：

我们现在可以用我们喜欢的文本编辑器查看文本文件的内容。让我们记住，程序的输出将包含在文件的开头：

现在，由于我们没有删除原始示例中对time的引用，我们在开头的前两行看到了它们的输出。然后我们可以看到在程序中进行的函数调用总数，以及运行它所花费的累积时间。

随后，我们有一个程序中调用的函数列表，按照累积耗时最长的函数到最短的顺序排列；第一行是程序本身，而第二行是我们程序中的simple_mandelbrot函数。（请注意，这里的时间与我们用time命令测量的时间一致）。之后，我们可以看到许多与将 Mandelbrot 图形转储到文件相关的库和系统调用，所有这些调用所花费的时间相对较少。我们使用cProfile的输出来推断给定程序中的瓶颈在哪里。

总结

使用 GPU 而不是 CPU 的主要优势是其增加的吞吐量，这意味着我们可以在 GPU 上同时执行更多并行代码，而不是在 CPU 上；GPU 无法使递归算法或非并行化算法变得更快。我们看到一些任务，比如建房子的例子，只能部分并行化——在这个例子中，我们无法加快设计房子的过程（在这种情况下本质上是串行的），但我们可以通过雇佣更多的工人来加快施工的过程（在这种情况下是可并行化的）。

我们使用这个类比来推导阿姆达尔定律，这是一个公式，可以在我们知道可并行代码执行时间百分比以及我们将运行此代码所需的处理器数量时，给出程序潜在加速的粗略估计。然后，我们应用了阿姆达尔定律来分析一个生成曼德勃罗集并将其转储到图像文件的小程序，并且我们确定这将是一个很好的候选者，可以并行化到 GPU 上。最后，我们以简要概述使用cPython模块对代码进行性能分析结束；这使我们能够看到程序中的瓶颈在哪里，而无需显式计时函数调用。

现在我们已经有了一些基本概念，并且有了学习 GPU 编程的动力，我们将在下一章中设置基于 Linux 或 Windows 10 的 GPU 编程环境。然后我们将立即进入 GPU 编程的世界，在接下来的章节中，我们将实际编写一个基于 GPU 的曼德勃罗程序的版本，这是我们在本章中看到的。

问题

1. 本章的曼德勃罗示例中有三个for语句；然而，我们只能并行化前两个。为什么我们不能在所有的for循环上并行化呢？

2. 当我们将阿姆达尔定律应用于将串行 CPU 算法卸载到 GPU 时，有什么是阿姆达尔定律没有考虑到的吗？

3. 假设您独家获得了三个全新的绝密 GPU，它们在所有方面都相同，只是核心数量不同——第一个有 131,072 个核心，第二个有 262,144 个核心，第三个有 524,288 个核心。如果您将曼德勃罗示例并行化并卸载到这些 GPU 上（生成一个 512 x 512 像素的图像），第一个 GPU 和第二个 GPU 之间的计算时间会有区别吗？第二个 GPU 和第三个 GPU 之间呢？

4. 在阿姆达尔定律的背景下，您能想到指定某些算法或代码块为可并行化存在什么问题吗？

5. 为什么我们应该使用性能分析器而不只是使用 Python 的time函数？

第二章：设置 GPU 编程环境

现在我们将看到如何在 Windows 和 Linux 下设置适当的 GPU 编程环境。在这两种情况下，我们都需要采取几个步骤。我们将逐步进行这些步骤，注意 Linux 和 Windows 之间的任何差异。当然，您可以随意跳过或忽略不适用于您选择的操作系统的任何部分或注释。

读者应注意，本章仅涵盖 64 位 Intel/AMD PC 的两个平台——Ubuntu LTS（长期支持）版本和 Windows 10。请注意，任何基于 Ubuntu LTS 的 Linux 操作系统（如 Xubuntu，Kubuntu 或 Linux Mint）也同样适用于通用的 Unity/GNOME-based Ubuntu 版本。

我们建议使用 Python 2.7 而不是 Python 3.x。 Python 2.7 在本文中使用的所有库中都有稳定的支持，并且我们已经在 Windows 和 Linux 平台上使用 Python 2.7 测试了本书中给出的每个示例。 Python 3.x 用户可以使用本书，但应该注意 Python 2.7 和 Python 3.x 之间的区别。本书中的一些示例已经在 Python 3.7 上进行了测试，但需要标准更改，例如在 Python print函数中添加括号。

Packt 作者 Sebastian Raschka 博士在sebastianraschka.com/Articles/2014_python_2_3_key_diff.html提供了 Python 2.7 和 3.x 之间的关键区别列表。

我们特别建议 Windows 和 Linux 用户使用 Anaconda Python 2.7 版本，因为它可以在用户基础上安装，无需sudo或管理员权限，包含本文所需的所有数据科学和可视化模块，并使用快速预优化的 NumPy/SciPy 包，这些包利用了英特尔的数学核心库（MKL）。 （默认的 Linux /usr/bin/python安装对于本文也应该足够，但您可能需要手动安装一些包，如 NumPy 和 Matplotlib。）

Anaconda Python（2.7 和 3.x 版本）可以在www.anaconda.com/download/.下载到所有平台上。

其他受支持平台的用户（例如 macOS，Windows 7/8，Windows Server 2016，Red Hat/Fedora，OpenSUSE 和 CENTOS）应查阅官方的 NVIDIA CUDA 文档（docs.nvidia.com/cuda/）以获取更多详细信息。此外，还有其他硬件选择：对于对嵌入式系统或具有树莓派等开发板经验的读者，可能希望从基于 ARM 的 NVIDIA Jetson 开发板开始，而对于对云计算或 Web 编程感兴趣的读者，可能考虑远程使用适当的 Azure 或 AWS 实例。在这些情况下，鼓励读者阅读官方文档以设置其驱动程序，编译器和 CUDA 工具包。本章中的一些步骤可能适用，也可能不适用。

本章的学习目标是：

• 确保我们拥有适当的硬件

• 安装 NVIDIA GPU 驱动程序

• 设置适当的 C/C++编程环境

• 安装 NVIDIA CUDA 工具包

• 为 GPU 编程设置 Python 环境

技术要求

本章建议安装 Anaconda Python 2.7，网址为www.anaconda.com/download/.

本章的代码也可以在 GitHub 上找到，网址为github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA.

有关先决条件的更多信息，请查看本书的前言；有关软件和硬件要求，请查看github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA中的 README 部分。

确保我们拥有正确的硬件

对于本书，我们建议您至少具备以下硬件：

• 64 位英特尔/AMD PC

• 4GB RAM

• NVIDIA GeForce GTX 1050 GPU（或更高）

这种配置将确保您可以轻松学习 GPU 编程，在本书中运行所有示例，并且还可以运行一些其他新的有趣的基于 GPU 的软件，如 Google 的 TensorFlow（一种机器学习框架）或 Vulkan SDK（一种尖端的图形 API）。

请注意，您必须拥有 NVIDIA 品牌的 GPU 才能使用本书！ CUDA Toolkit 专为 NVIDIA 卡而设计，因此无法用于编程 Intel HD 或 Radeon GPU。

正如所述，我们将假设您使用的是 Windows 10 或 Ubuntu LTS（长期支持）版本。

Ubuntu LTS 发布通常具有 14.04、16.04、18.04 等形式的版本号。

Ubuntu LTS，大体上来说，是最主流的 Linux 版本，可以确保与新软件和工具包的最大兼容性。请记住，有许多基于 Ubuntu 的 Linux 变体，如 Linux Mint 或 Xubuntu，这些通常同样有效。（我个人发现 Linux Mint 在配备 GPU 的笔记本电脑上开箱即用效果相当不错。）

我们应该注意，我们假设您至少拥有一款入门级 GTX 1050（Pascal）GPU，或者在任何更新的架构中具有相当的性能。请注意，本书中的许多示例很可能在大多数旧 GPU 上运行，但作者只在 GTX 1050（在 Windows 10 下）和 GTX 1070（在 Linux 下）上进行了测试。虽然这些示例尚未在旧 GPU 上进行测试，但 2014 年的入门级 Maxwell 架构 GPU，如 GTX 750，也应足以满足本文的要求。

如果您使用的是台式 PC，请确保在继续之前已经按照所有包含的说明物理安装了 GPU。

检查您的硬件（Linux）

现在，我们将在 Linux 中进行一些基本检查，以确保我们拥有正确的硬件。首先让我们打开一个终端并切换到 bash 命令行——您可以通过在 Ubuntu 中快速按下组合键Ctrl + Alt + T来快速完成这一步。

现在，通过输入lscpu并按Enter来检查我们的处理器。会出现大量信息，但只需查看第一行，确保架构确实是 x86_64：

接下来，通过在 bash 提示符下输入free -g并再次按Enter来检查我们的内存容量。这将告诉我们在第一行的第一个条目中我们拥有的总内存量（以 GB 为单位），以及在接下来的行中交换空间中的内存量：

这绝对是足够的内存。

最后，让我们看看我们是否有适当的 GPU。NVIDIA GPU 通过 PCI 总线与我们的 PC 通信，因此我们可以使用lspci命令列出所有 PCI 硬件。通常会列出许多其他硬件，因此让我们使用grep命令仅过滤出 NVIDIA GPU，输入lspci | grep -e "NVIDIA"在 bash 提示符下：

这是一款 GTX 1070，幸运的是它超出了我们至少需要 GTX 1050 的要求。

检查您的硬件（Windows）

首先，我们必须打开 Windows 面板。我们可以通过按下Windows + R，然后在提示符处输入Control Panel来实现这一点，如下面的屏幕截图所示：

Windows 控制面板将弹出。现在点击系统和安全，然后选择以下屏幕上的系统。这将立即告诉我们我们拥有多少 RAM 以及我们是否拥有 64 位处理器：

要检查我们的 GPU，请点击此窗口左上角的设备管理器。然后 Windows 设备管理器将弹出；然后您可以选择显示适配器下拉框来检查您系统上的 GPU：

安装 GPU 驱动程序

如果您已经安装了 GPU 的驱动程序，您可能可以跳过此步骤；此外，一些版本的 CUDA 已经预先打包了最新的驱动程序。通常情况下，CUDA 对您安装的驱动程序非常挑剔，甚至可能无法与 CUDA Toolkit 驱动程序一起工作，因此您可能需要尝试几种不同的驱动程序，直到找到一个可用的。

一般来说，Windows 具有更好的 CUDA 驱动程序兼容性和更用户友好的安装比 Linux。Windows 用户可以考虑跳过此步骤，只使用与 CUDA Toolkit 捆绑的驱动程序，我们稍后将在本章中安装。然而，我们强烈建议 Linux 用户（特别是 Linux 笔记本用户）在继续之前，密切遵循本节中的所有步骤。

安装 GPU 驱动程序（Linux）

在 Ubuntu 中，NVIDIA GPU 的默认驱动程序是一个名为 Nouveau 的开源驱动程序；不幸的是，这在 CUDA 中根本不起作用，因此我们必须安装专有驱动程序。我们必须将特殊的graphics-drivers存储库添加到我们的软件包管理器中，以便能够将专有 NVIDIA 驱动程序下载到我们的 Ubuntu 系统中。我们通过在 bash 提示符中输入以下行来添加存储库：

sudo add-apt-repository ppa:graphics-drivers/ppa

由于这是一个sudo超级用户命令，您将需要输入您的密码。我们现在通过输入以下行来将系统与新的存储库同步：

sudo apt-get update

我们现在应该准备安装我们的驱动程序。从 Ubuntu 桌面，按下Windows + R，然后输入software and drivers：

软件和驱动程序设置菜单应该出现。从这里，点击标记为附加驱动程序的选项卡。您应该看到一系列可用的稳定专有驱动程序供您的 GPU 选择；选择您看到的最新的一个（在我的情况下，它是nvidia-driver-396，如下所示）：

选择最新的驱动程序后，点击应用更改。您将再次被提示输入您的sudo密码，然后驱动程序将安装；进度条应该出现。请注意，这个过程可能需要很长时间，而且可能会出现您的计算机“挂起”的情况；这个过程可能需要超过一个小时，所以请耐心等待。

最后，当过程完成时，重启您的计算机，返回到 Ubuntu 桌面。现在输入Windows + A，然后输入nvidia-settings（或者，从 bash 提示符中运行此程序）。NVIDIA X Server 设置管理器应该出现，并指示您正在使用适当的驱动程序版本：

安装 GPU 驱动程序（Windows）

重申一下-通常建议读者最初跳过此步骤，然后安装包含在 CUDA Toolkit 中的驱动程序。

Windows 的最新驱动程序可以直接从 NVIDIA 的www.nvidia.com/Download/下载。只需从下拉菜单中选择适用于您 GPU 的适当的 Windows 10 驱动程序，这些是可执行（.exe）文件。只需通过双击文件管理器中的文件来安装驱动程序。

建立 C++编程环境

现在我们已经安装了驱动程序，我们必须设置我们的 C/C++编程环境；Python 和 CUDA 都对它们可能集成的编译器和 IDE 有特殊要求，所以您可能需要小心。对于 Ubuntu Linux 用户，标准存储库的编译器和 IDE 通常可以完美地与 CUDA 工具包集成，而 Windows 用户可能需要更加小心。

设置 GCC、Eclipse IDE 和图形依赖项（Linux）

从 Ubuntu 桌面打开终端（Ctrl + Alt + T）。我们首先更新apt存储库如下：

sudo apt-get update

现在我们可以用一行额外的命令安装我们需要的 CUDA 一切：

sudo apt-get install build-essential binutils gdb eclipse-cdt

在这里，build-essential是带有gcc和g++编译器以及其他实用程序（如 make）的软件包；binutils有一些通用的实用程序，如 LD 链接器；gdb是调试器；Eclipse 是我们将要使用的 IDE。

让我们还安装一些额外的依赖项，这将允许我们使用以下命令运行 CUDA 工具包中包含的一些图形（OpenGL）演示：

sudo apt-get install freeglut3 freeglut3-dev libxi-dev libxmu-dev

现在您应该可以安装 CUDA 工具包了。

在 Windows 上设置 Visual Studio

在撰写本文时，只有一个版本的 Visual Studio 似乎完美地与 Python 和最新的 CUDA 工具包集成在一起——Visual Studio 2015；也就是说，Visual Studio 版本 14.0。

虽然可能可以在较新版本的 Visual Studio（例如 2017）下进行子安装，但我们建议读者直接在系统上安装带有 C/C++支持的 Visual Studio 2015。

Visual Studio Community 2015，这个软件的免费版本，可以在visualstudio.microsoft.com/vs/older-downloads/下载。

在这里，我们将进行最小化安装，只安装 CUDA 所需的组件。我们运行安装软件，并选择自定义安装：

点击下一步，然后点击编程语言的下拉框，然后选择 Visual C++（如果您需要其他包或编程语言，可以随意选择，但是对于 GPU 编程，我们只需要 Visual C++）：

这个安装过程可能需要一些时间。完成后，我们将准备安装 CUDA 工具包。

安装 CUDA 工具包

最后，我们开始接近我们的目标！现在我们通过访问developer.nvidia.com/cuda-downloads来下载我们的 CUDA 工具包。选择适当的操作系统，您将看到几个选项。对于 Windows 和 Linux，都有网络和本地安装选项。我倾向于在 Windows 和 Linux 下都使用本地安装选项，因为我更喜欢一次性下载整个软件包；如果有任何网络问题，那么您可以确保在安装 CUDA 工具包时不会发生问题。

安装 CUDA 工具包（Linux）

对于 Linux 用户，您将看到使用.deb包和.run文件的选择；对于大多数用户，我建议使用.deb文件，因为这将自动安装 CUDA 需要的任何缺少的软件包。.run文件安装在系统的高级软件包工具（APT）系统之外，它只是将适当的文件复制到系统的/usr二进制和库目录。如果您不想干扰系统的 APT 系统或存储库，并且对 Linux 有很好的理解，那么.run文件可能更合适。无论哪种情况，请仔细遵循网站上关于安装软件包的说明，这些说明可能会因版本而略有不同。

安装包完成后，您可能需要配置您的PATH和LD_SYSTEM_CONFIG环境变量，以便您的系统可以找到 CUDA 所需的适当的二进制可执行文件和库文件。我建议您通过将以下行附加到您用户目录中的.bashrc文件的末尾来完成这个步骤。使用您喜欢的文本编辑器，如gedit、nano、emacs或vim打开~/.bashrc文件，然后在文件的最底部添加以下行：

export PATH="/usr/local/cuda/bin:${PATH}
export LD_LIBRARY_PATH="/usr/local/cuda/lib64:${LD_LIBRARY_PATH}"

保存文件然后退出终端。您现在可以通过打开一个新的终端并输入nvcc --version然后按Enter来确保您已正确安装了工具包，这将给您工具包编译器的版本信息。（nvcc是命令行 CUDA C 编译器，类似于gcc编译器。）

安装 CUDA Toolkit（Windows）

对于 Windows 用户，您可以通过双击.exe文件并按照屏幕上的提示来安装包。

安装完成后，重置您的系统。我们现在将通过检查nvcc编译器来确保 CUDA 已正确安装。在开始菜单下，点击Visual Studio 2015文件夹，然后点击 VS2015 x64 Native Tools Command Prompt。一个终端窗口将弹出；现在输入nvcc --version并按Enter，这应该会给您 NVIDIA 编译器的版本信息。

为 GPU 编程设置我们的 Python 环境

使用我们的编译器、集成开发环境和 CUDA 工具包正确安装在我们的系统上，我们现在可以为 GPU 编程设置一个合适的 Python 环境。这里有很多选择，但我们明确建议您使用 Anaconda Python Distribution。Anaconda Python 是一个独立且用户友好的分发版，可以直接安装在您的用户目录中，而且不需要任何管理员或sudo级别的系统访问权限来安装、使用或更新。

请记住，Anaconda Python 有两种版本——Python 2.7 和 Python 3。由于 Python 3 目前对我们将要使用的一些库的支持不是很好，我们将在本书中使用 Python 2.7，这仍然是广泛使用的。

您可以通过访问www.anaconda.com/download来安装 Anaconda Python，选择您的操作系统，然后选择下载分发版的 Python 2.7 版本。按照 Anaconda 网站上给出的说明安装分发版，这相对比较简单。现在我们可以为 GPU 编程设置本地 Python 安装。

我们现在将设置本书中可能是最重要的 Python 包：Andreas Kloeckner 的 PyCUDA 包。

安装 PyCUDA（Linux）

在 Linux 中打开一个命令行。通过在 bash 提示符下输入which python并按Enter来确保您的PATH变量正确设置为使用本地 Anaconda 安装的 Python（而不是系统范围的安装）（Anaconda 应该在安装过程中自动配置您的.bashrc）；这应该告诉您 Python 二进制文件在您的本地~/anaconda2/bin目录中，而不是在/usr/bin目录中。如果不是这种情况，请打开一个文本编辑器，并在您的~/.bashrc文件的末尾放置以下行export PATH="/home/${USER}/anaconda2/bin:${PATH}"，保存后，打开一个新的终端，然后再次检查。

有几种安装 PyCUDA 的选项。最简单的选项是从 PyPI 存储库安装最新稳定版本，方法是输入pip install pycuda。您还可以按照 PyCUDA 官方网站上的说明安装最新版本的 PyCUDA，网址为mathema.tician.de/software/pycuda/。请注意，如果您希望从不同的来源重新安装 PyCUDA，请确保首先使用pip uninstall pycuda卸载它。

创建一个环境启动脚本（Windows）

Windows 用户需要特别注意，他们的 Visual Studio 和 Anaconda Python 环境变量是否设置正确，以便使用 PyCUDA；否则，Python 将无法找到 NVIDIA 的nvcc CUDA 编译器或 Microsoft 的cl.exe C++编译器。幸运的是，包含了批处理脚本，可以自动为我们设置这些环境，但我们必须小心，每次想要进行 GPU 编程时都要执行这些脚本。

因此，我们将创建一个批处理脚本，通过连续调用这两个脚本来启动适当的 IDE 或命令行环境。 （此脚本也可在github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA/blob/master/2/launch-python-cuda-environment.bat上找到。）

请务必首先打开 Windows 记事本，并跟随操作：

首先找到您的 Visual Studio 的vcvars.bat文件的位置；对于 Visual Studio 2015，它位于C:\Program Files (x86)\Microsoft Visual Studio 14.0\VC\vcvarsall.bat。

在文本编辑器中输入以下行，然后按Enter：

call "C:\Program Files (x86)\Microsoft Visual Studio 14.0\VC\vcvarsall.bat" amd64

现在，我们需要调用 Anaconda 的activate.bat脚本来设置 Anaconda Python 环境变量；标准路径是Anaconda2\Scripts\activate.bat。我们还必须指示此脚本的参数是 Anaconda 库的位置。在我的情况下，我的启动脚本中的第二行将是call "C:\Users\%username%\Anaconda2\Scripts\activate.bat" C:\Users\%username%\Anaconda2。

最后，我们的批处理脚本的最后一行将启动您喜欢的任何环境——IDE 或命令行提示符，它将继承前两个脚本设置的所有必要环境和系统变量。如果您喜欢旧的标准 DOS 风格命令提示符，这行应该只是cmd。如果您喜欢从 PowerShell 工作，请将其更改为powershell。在某些情况下，特别是用于访问命令行pip和conda来更新 Python 库时，需要使用命令行。

最后，将此文件保存到桌面，并命名为launch-python-cuda-environment.bat。现在，您可以通过双击此文件来启动我们的 Python GPU 编程环境。

（请记住，如果您希望使用 Jupyter Notebook 或 Spyder Python IDE，您可以简单地通过jupyter-notebook或spyder从命令行启动它们，或者您可以制作一个批处理脚本，只需用适当的 IDE 启动命令替换cmd。）

安装 PyCUDA（Windows）

由于大多数 Python 库主要是由 Linux 用户编写和为 Linux 用户编写的，建议您从 Christoph Gohlke 的网站上安装预构建的 PyCUDA wheel 二进制文件，网址为：www.lfd.uci.edu/~gohlke/pythonlibs/#pycuda。下载一个文件，文件名为pycuda-2017.1.1+cuda(VERSION)-cp27-cp27m-win_amd64.whl，其中版本是您的 CUDA 版本号。现在，您可以通过在命令行中输入以下内容并用您的 PyCUDA wheel 的完整路径和文件名替换pycuda.whl来安装 PyCUDA：

pip install pycuda.whl

（或者，您可以尝试使用pip install pycuda从 PyPI 存储库安装 PyCUDA，或者按照 PyCUDA 网站上的说明操作。）

测试 PyCUDA

最后，我们到了可以看到我们的 GPU 编程环境是否真正起作用的时候。我们将运行下一章的一个小程序，该程序将查询我们的 GPU 并提供有关型号号码、内存、核心数量、架构等相关信息。从存储库中的目录3中获取 Python 文件（deviceQuery.py），也可以在github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA/blob/master/3/deviceQuery.py上找到。

如果您使用的是 Windows，请确保通过在桌面上启动我们在上一节中创建的.bat文件来启动 GPU 编程环境。否则，如果您使用的是 Linux，请打开一个 bash 终端。现在输入以下命令并按Enter键——python deviceQuery.py。

这将输出许多行数据，但前几行应该表明 PyCUDA 已经检测到您的 GPU，并且您应该在下一行看到型号号码：

恭喜，您现在已经准备好进入 GPU 编程的世界了！

总结

为 GPU 编程设置 Python 环境可能是一个非常微妙的过程。本文建议 Windows 和 Linux 用户都使用 Anaconda Python 2.7 发行版。首先，我们应该确保我们有正确的硬件进行 GPU 编程；一般来说，64 位 Windows 或 Linux PC，带有 4GB RAM 和 2016 年或之后的任何入门级 NVIDIA GPU 将足够满足我们的需求。Windows 用户应该注意使用一个既适用于 CUDA 工具包又适用于 Anaconda 的 Visual Studio 版本（如 VS 2015），而 Linux 用户在安装 GPU 驱动程序时应特别小心，并在其.bashrc文件中设置适当的环境变量。此外，Windows 用户应该创建一个适当的启动脚本，用于设置 GPU 编程环境，并应该使用预编译的 PyCUDA 库安装文件。

现在，我们的编程环境已经设置好了，接下来的一章我们将学习 GPU 编程的基础知识。我们将看到如何将数据写入 GPU 的内存，以及如何在 CUDA C 中编写一些非常简单的逐元素GPU 函数。（如果你看过经典的 1980 年代电影《功夫小子》，那么你可能会把下一章看作是 GPU 编程的“上蜡，下蜡”。）

问题

1. 我们可以在主处理器内置的英特尔 HD GPU 上运行 CUDA 吗？离散的 AMD Radeon GPU 呢？

2. 这本书的示例是使用 Python 2.7 还是 Python 3.7？

3. 在 Windows 中，我们使用什么程序来查看我们安装了什么 GPU 硬件？

4. 在 Linux 中，我们使用什么命令行程序来查看我们安装了什么 GPU 硬件？

5. 在 Linux 中，我们使用什么命令来确定系统有多少内存？

6. 如果我们不想改变我们的 Linux 系统的 APT 存储库，我们应该使用run还是deb安装程序来安装 CUDA？

第三章：开始使用 PyCUDA

在上一章中，我们设置了编程环境。现在，有了我们的驱动程序和编译器牢固地安装好，我们将开始实际的 GPU 编程！我们将首先学习如何使用 PyCUDA 进行一些基本和基础的操作。我们将首先看看如何查询我们的 GPU - 也就是说，我们将首先编写一个小的 Python 程序，告诉我们 GPU 的特性，如核心数量、架构和内存。然后，我们将花一些时间熟悉如何在 Python 和 GPU 之间传输内存，使用 PyCUDA 的gpuarray类以及如何使用这个类进行基本计算。本章的其余部分将花在展示如何编写一些基本函数（我们将称之为CUDA 内核），我们可以直接启动到 GPU 上。

本章的学习成果如下：

• 使用 PyCUDA 确定 GPU 特性，如内存容量或核心数量

• 理解主机（CPU）和设备（GPU）内存之间的区别，以及如何使用 PyCUDA 的gpuarray类在主机和设备之间传输数据

• 如何使用只有gpuarray对象进行基本计算

• 如何使用 PyCUDA 的ElementwiseKernel函数在 GPU 上执行基本的逐元素操作

• 理解函数式编程概念的 reduce/scan 操作，以及如何制作基本的缩减或扫描 CUDA 内核

技术要求

本章需要一台安装了现代 NVIDIA GPU（2016 年以后）的 Linux 或 Windows 10 PC，并安装了所有必要的 GPU 驱动程序和 CUDA Toolkit（9.0 以后）。还需要一个合适的 Python 2.7 安装（如 Anaconda Python 2.7），并安装了 PyCUDA 模块。

本章的代码也可以在 GitHub 上找到：github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA。

有关先决条件的更多信息，请查看本书的前言；有关软件和硬件要求，请查看github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA中的README部分。

查询您的 GPU

在我们开始编程 GPU 之前，我们应该真正了解一些关于其技术能力和限制的知识。我们可以通过进行所谓的GPU 查询来确定这一点。GPU 查询是一个非常基本的操作，它将告诉我们 GPU 的具体技术细节，如可用的 GPU 内存和核心数量。NVIDIA 在samples目录中（适用于 Windows 和 Linux）包含了一个纯 CUDA-C 编写的命令行示例deviceQuery，我们可以运行它来执行此操作。让我们看一下作者的 Windows 10 笔记本电脑（Microsoft Surface Book 2，配备了 GTX 1050 GPU）上产生的输出：

让我们来看看这里显示的所有技术信息的一些基本要点。首先，我们看到只安装了一个 GPU，设备 0 - 可能主机计算机安装了多个 GPU 并使用它们，因此 CUDA 将为每个GPU 设备指定一个独立的编号。有些情况下，我们可能需要明确指定设备编号，所以了解这一点总是很好的。我们还可以看到我们拥有的具体设备类型（这里是 GTX 1050），以及我们正在使用的 CUDA 版本。现在我们还要注意两件事：核心的总数（这里是 640），以及设备上的全局内存总量（在本例中为 2,048 兆字节，即 2 千兆字节）。

虽然您可以从deviceQuery中看到许多其他技术细节，但核心数量和内存量通常是您第一次在新 GPU 上运行时应该关注的前两件事，因为它们可以让您最直接地了解新设备的容量。

使用 PyCUDA 查询您的 GPU

现在，最后，我们将通过用 Python 编写我们自己的版本的deviceQuery来开始我们的 GPU 编程之旅。在这里，我们主要关注设备上可用内存的数量，计算能力，多处理器的数量和 CUDA 核心的总数。

我们将从以下方式初始化 CUDA：

import pycuda.driver as drv
drv.init()

请注意，我们将始终需要使用pycuda.driver.init()或通过导入 PyCUDA 的autoinit子模块import pycuda.autoinit来初始化 PyCUDA！

我们现在可以立即检查我们的主机计算机上有多少个 GPU 设备：

print 'Detected {} CUDA Capable device(s)'.format(drv.Device.count())

让我们在 IPython 中输入这个并看看会发生什么：

太好了！到目前为止，我已经验证了我的笔记本确实有一个 GPU。现在，让我们通过添加几行代码来迭代可以通过pycuda.driver.Device（按编号索引）单独访问的每个设备，以提取有关此 GPU（以及系统上的任何其他 GPU）的更多有趣信息。设备的名称（例如，GeForce GTX 1050）由name函数给出。然后我们使用compute_capability函数获取设备的计算能力和total_memory函数获取设备的总内存量。

计算能力可以被视为每个 NVIDIA GPU 架构的版本号；这将为我们提供一些关于设备的重要信息，否则我们无法查询，我们将在一分钟内看到。

我们将这样写：

for i in range(drv.Device.count()):

     gpu_device = drv.Device(i)
     print 'Device {}: {}'.format( i, gpu_device.name() )
     compute_capability = float( '%d.%d' % gpu_device.compute_capability() )
     print '\t Compute Capability: {}'.format(compute_capability)
     print '\t Total Memory: {} megabytes'.format(gpu_device.total_memory()//(1024**2))

现在，我们准备查看 PyCUDA 以 Python 字典类型形式提供给我们的 GPU 的一些剩余属性。我们将使用以下行将其转换为由字符串索引属性的字典：

    device_attributes_tuples = gpu_device.get_attributes().iteritems()
     device_attributes = {}

     for k, v in device_attributes_tuples:
         device_attributes[str(k)] = v

现在，我们可以使用以下内容确定设备上的多处理器数量：

    num_mp = device_attributes['MULTIPROCESSOR_COUNT']

GPU 将其各个核心划分为称为流处理器（SMs）的较大单元； GPU 设备将具有多个 SM，每个 SM 将根据设备的计算能力具有特定数量的 CUDA 核心。要明确：每个多处理器的核心数并不是由 GPU 直接指示的-这是由计算能力隐含给我们的。我们将不得不查阅 NVIDIA 的一些技术文件以确定每个多处理器的核心数（参见docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#compute-capabilities），然后创建一个查找表来给出每个多处理器的核心数。我们使用compute_capability变量来查找核心数：

    cuda_cores_per_mp = { 5.0 : 128, 5.1 : 128, 5.2 : 128, 6.0 : 64, 6.1 : 128, 6.2 : 128}[compute_capability]

现在我们可以通过将这两个数字相乘来最终确定设备上的总核心数：

    print '\t ({}) Multiprocessors, ({}) CUDA Cores / Multiprocessor: {} CUDA Cores'.format(num_mp, cuda_cores_per_mp, num_mp*cuda_cores_per_mp)

现在，我们可以通过迭代字典中剩余的键并打印相应的值来完成我们的程序：

    device_attributes.pop('MULTIPROCESSOR_COUNT')

     for k in device_attributes.keys():
         print '\t {}: {}'.format(k, device_attributes[k])

所以，现在我们终于完成了本文的第一个真正的 GPU 程序！（也可在github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA/blob/master/3/deviceQuery.py找到）。现在，我们可以按如下方式运行它：

现在，我们可以有点自豪，因为我们确实可以编写一个程序来查询我们的 GPU！现在，让我们真正开始学习使用我们的 GPU，而不仅仅是观察它。

使用 PyCUDA 的 gpuarray 类

就像 NumPy 的 array 类是 NumPy 环境中数值编程的基石一样，PyCUDA 的 gpuarray 类在 Python 中的 GPU 编程中扮演着类似的重要角色。它具有你从 NumPy 中熟悉和喜爱的所有功能——多维向量/矩阵/张量形状结构、数组切片、数组展开，以及用于逐点计算的重载运算符（例如 +、-、*、/ 和 **）。

gpuarray 对于任何新手 GPU 程序员来说都是一个不可或缺的工具。在我们继续之前，我们将花费这一部分时间来了解这种特定的数据结构，并对其有一个深入的理解。

使用 gpuarray 将数据传输到 GPU 和从 GPU 中传输数据

正如我们在 Python 中编写的先前的 deviceQuery 程序所示，GPU 有自己的内存，与主机计算机的内存分开，这被称为设备内存。（有时这更具体地被称为全局设备内存，以区分它与 GPU 上的其他缓存内存、共享内存和寄存器内存。）在大多数情况下，我们将 GPU 上的（全局）设备内存视为我们在 C 中动态分配的堆内存（使用 malloc 和 free 函数）或 C++（使用 new 和 delete 运算符）；在 CUDA C 中，这进一步复杂化，需要在 CPU 和 GPU 空间之间来回传输数据（使用诸如 cudaMemcpyHostToDevice 和 cudaMemcpyDeviceToHost 的命令），同时跟踪 CPU 和 GPU 空间中的多个指针，并执行适当的内存分配（cudaMalloc）和释放（cudaFree）。

幸运的是，PyCUDA 通过 gpuarray 类涵盖了所有的内存分配、释放和数据传输的开销。正如所述，这个类类似于 NumPy 数组，使用矢量/矩阵/张量形状结构信息来处理数据。gpuarray 对象甚至根据其生命周期自动执行清理，因此当我们完成后，我们不必担心释放存储在 gpuarray 对象中的任何 GPU 内存。

那么我们如何使用它将数据从主机传输到 GPU 呢？首先，我们必须将我们的主机数据包含在某种形式的 NumPy 数组中（我们称之为 host_data），然后使用 gpuarray.to_gpu(host_data) 命令将其传输到 GPU 并创建一个新的 GPU 数组。

现在让我们在 GPU 中执行一个简单的计算（在 GPU 上的常数点乘），然后使用 gpuarray.get 函数将 GPU 数据检索到一个新的数组中。让我们加载 IPython 并看看它是如何工作的（请注意，这里我们将使用 import pycuda.autoinit 初始化 PyCUDA）：

需要注意的一点是，当我们设置 NumPy 数组时，我们特别指定了主机上的数组类型为 NumPy float32 类型，并使用 dtype 选项；这与 C/C++ 中的浮点类型直接对应。一般来说，当我们发送数据到 GPU 时，最好使用 NumPy 明确设置数据类型。原因有两个：首先，由于我们使用 GPU 来提高应用程序的性能，我们不希望使用不必要的类型造成不必要的计算时间或内存开销；其次，由于我们很快将在内联 CUDA C 中编写代码的部分，我们必须非常具体地指定类型，否则我们的代码将无法正确工作，要记住 C 是一种静态类型语言。

记得为将要传输到 GPU 的 NumPy 数组明确设置数据类型。这可以在 numpy.array 类的构造函数中使用 dtype 选项来完成。

使用 gpuarray 进行基本的逐点算术运算

在最后一个示例中，我们看到我们可以使用（重载的）Python 乘法运算符（*）来将gpuarray对象中的每个元素乘以一个标量值（这里是 2）；请注意，逐点操作本质上是可并行化的，因此当我们在gpuarray对象上使用此操作时，PyCUDA 能够将每个乘法操作分配到单个线程上，而不是依次串行计算每个乘法（公平地说，一些版本的 NumPy 可以使用现代 x86 芯片中的高级 SSE 指令进行这些计算，因此在某些情况下性能将与 GPU 相当）。明确一点：在 GPU 上执行的这些逐点操作是并行的，因为一个元素的计算不依赖于任何其他元素的计算。

为了了解这些运算符的工作原理，我建议读者加载 IPython 并在 GPU 上创建一些gpuarray对象，然后玩几分钟，看看这些运算符是否与 NumPy 中的数组类似。以下是一些灵感：

现在，我们可以看到gpuarray对象的行为是可预测的，并且与 NumPy 数组的行为一致。（请注意，我们将不得不使用get函数从 GPU 中获取输出！）现在让我们比较一下 CPU 和 GPU 计算时间，看看在何时是否有任何优势进行这些操作。

速度测试

让我们编写一个小程序（time_calc0.py），对 CPU 上的标量乘法和 GPU 上的相同操作进行速度比较测试。然后，我们将使用 NumPy 的allclose函数比较两个输出值。我们将生成一个包含 5000 万个随机 32 位浮点值的数组（这将大约占用 48 兆字节的数据，因此在任何稍微现代的主机和 GPU 设备上都应该完全可行），然后我们将计算在两个设备上将数组乘以 2 所需的时间。最后，我们将比较输出值以确保它们相等。操作如下：

import numpy as np
import pycuda.autoinit
from pycuda import gpuarray
from time import time
host_data = np.float32( np.random.random(50000000) )

t1 = time()
host_data_2x =  host_data * np.float32(2)
t2 = time()

print 'total time to compute on CPU: %f' % (t2 - t1)
device_data = gpuarray.to_gpu(host_data)

t1 = time()
device_data_2x =  device_data * np.float32( 2 )
t2 = time()

from_device = device_data_2x.get()
print 'total time to compute on GPU: %f' % (t2 - t1)

print 'Is the host computation the same as the GPU computation? : {}'.format(np.allclose(from_device, host_data_2x) )

（您可以在之前提供给您的存储库中找到time_calc0.py文件。）

现在，让我们加载 IPython 并运行几次，以了解这些的一般速度，并查看是否有任何变化。（这里，这是在 2017 年的微软 Surface Book 2 上运行的，配备了 Kaby Lake i7 处理器和 GTX 1050 GPU。）：

我们首先注意到，每次计算的 CPU 计算时间大致相同（大约 0.08 秒）。然而，我们注意到，第一次运行时，GPU 计算时间比 CPU 计算时间慢得多（1.09 秒），并且在随后的运行中变得快得多，在每次后续运行中保持大致恒定（在 7 或 9 毫秒的范围内）。如果您退出 IPython，然后再次运行程序，将发生相同的情况。这种现象的原因是什么？好吧，让我们使用 IPython 内置的prun分析器进行一些调查工作。（这类似于第一章中介绍的cProfiler模块，为什么要进行 GPU 编程？）

首先，让我们将我们的程序作为文本加载到 IPython 中，然后通过 Python 的exec命令运行我们的分析器：

with open('time_calc0.py','r') as f:
     time_calc_code = f.read()

现在，我们在 IPython 控制台中键入%prun -s cumulative exec(time_calc_code)（带有前导%）并查看哪些操作花费了最多的时间：

现在，有一些可疑的对 Python 模块文件compiler.py的调用；这些调用总共大约需要一秒钟，比在这里进行 GPU 计算所需的时间略少。现在让我们再次运行一下，看看是否有任何差异：

请注意，这一次没有调用compiler.py。为什么呢？由于 PyCUDA 库的性质，GPU 代码通常在给定的 Python 会话中首次运行时使用 NVIDIA 的nvcc编译器进行编译和链接；然后它被缓存，如果再次调用代码，则不必重新编译。这甚至可能包括简单的操作，比如标量乘法！（我们最终会看到，通过使用第十章中的预编译代码或使用 NVIDIA 自己的线性代数库与 Scikit-CUDA 模块一起使用 CUDA 库，可以改善这一点，我们将在第七章中看到）。

在 PyCUDA 中，GPU 代码通常在运行时使用 NVIDIA 的nvcc编译器进行编译，然后从 PyCUDA 中调用。这可能会导致意外的减速，通常是在给定的 Python 会话中首次运行程序或 GPU 操作时。

使用 PyCUDA 的 ElementWiseKernel 执行逐点计算

现在，让我们看看如何使用 PyCUDA 的ElementWiseKernel函数直接在 GPU 上编写我们自己的逐点（或等效地，逐元素）操作。这就是我们之前对 C/C++编程的了解将变得有用的地方——我们将不得不在 CUDA C 中编写一点内联代码，这些代码是由 NVIDIA 的nvcc编译器在外部编译的，然后通过 PyCUDA 在运行时由我们的代码启动。

在本文中，我们经常使用术语kernel；通过kernel，我们总是指的是由 CUDA 直接启动到 GPU 上的函数。我们将使用 PyCUDA 的几个函数来生成不同类型的 kernel 的模板和设计模式，以便更轻松地过渡到 GPU 编程。

让我们直接开始；我们将从头开始重写代码，使用 CUDA-C 将gpuarray对象的每个元素乘以 2；我们将使用 PyCUDA 的ElementwiseKernel函数来生成我们的代码。您应该尝试直接在 IPython 控制台中输入以下代码。（不那么冒险的人可以从本文的 Git 存储库中下载，文件名为simple_element_kernel_example0.py）：

import numpy as np
import pycuda.autoinit
from pycuda import gpuarray
from time import time
from pycuda.elementwise import ElementwiseKernel
host_data = np.float32( np.random.random(50000000) )
gpu_2x_ker = ElementwiseKernel(
"float *in, float *out",
"out[i] = 2*in[i];",
"gpu_2x_ker")

让我们看看这是如何设置的；当然，这是几行内联 C。我们首先在第一行中设置输入和输出变量（"float *in, float *out"），这通常是指指向 GPU 上已分配内存的 C 指针的形式。在第二行中，我们使用"out[i] = 2*in[i];"定义了我们的逐元素操作，它将把in中的每个点乘以 2，并将其放在out的相应索引中。

请注意，PyCUDA 会自动为我们设置整数索引i。当我们使用i作为我们的索引时，ElementwiseKernel将自动在 GPU 的许多核心中并行化我们的计算。最后，我们给我们的代码片段起了一个内部 CUDA C kernel 的名称（"gpu_2x_ker"）。由于这是指 CUDA C 的命名空间而不是 Python 的，因此将其与 Python 中的名称相同是可以的（也很方便）。

现在，让我们进行速度比较：

def speedcomparison():
    t1 = time()
    host_data_2x =  host_data * np.float32(2)
    t2 = time()
    print 'total time to compute on CPU: %f' % (t2 - t1)
    device_data = gpuarray.to_gpu(host_data)
    # allocate memory for output
    device_data_2x = gpuarray.empty_like(device_data)
    t1 = time()
    gpu_2x_ker(device_data, device_data_2x)
    t2 = time()
    from_device = device_data_2x.get()
    print 'total time to compute on GPU: %f' % (t2 - t1)
    print 'Is the host computation the same as the GPU computation? : {}'.format(np.allclose(from_device, host_data_2x) )

if __name__ == '__main__':
    speedcomparison()

现在，让我们运行这个程序：

哇！看起来不太好。让我们从 IPython 中运行speedcomparison()函数几次：

正如我们所看到的，第一次使用给定的 GPU 函数后，速度显著增加。与前面的例子一样，这是因为 PyCUDA 在首次调用给定的 GPU kernel 函数时使用nvcc编译器编译我们的内联 CUDA C 代码。代码编译后，它将被缓存并在给定的 Python 会话的其余部分中重复使用。

现在，在我们继续之前，让我们再讨论一些重要的事情，这是非常微妙的。我们定义的小内核函数操作 C 浮点指针；这意味着我们将不得不在 GPU 上分配一些空的内存，该内存由out变量指向。再次看一下speedcomparison()函数中的这部分代码：

device_data = gpuarray.to_gpu(host_data)
# allocate memory for output
device_data_2x = gpuarray.empty_like(device_data)

与之前一样，我们通过gpuarray.to_gpu函数将一个 NumPy 数组（host_data）发送到 GPU，该函数会自动将数据分配到 GPU 并从 CPU 空间复制过来。我们将把这个数组插入到我们内核函数的in部分。在下一行，我们使用gpuarray.empty_like函数在 GPU 上分配空的内存。这类似于 C 中的普通malloc，分配一个与device_data大小和数据类型相同的数组，但不复制任何内容。现在我们可以将其用于内核函数的out部分。现在我们来看一下speedcomparison()中的下一行，看看如何将我们的内核函数启动到 GPU 上（忽略我们用于计时的行）：

gpu_2x_ker(device_data, device_data_2x)

再次，我们设置的变量直接对应于我们用ElementwiseKernel定义的第一行（这里是"float *in, float *out"）。

曼德勃罗重新审视

让我们再次从第一章“为什么使用 GPU 编程？”中生成曼德勃罗集的问题。原始代码可以在存储库的1文件夹中找到，文件名为mandelbrot0.py，在继续之前，您应该再次查看一下。我们看到该程序有两个主要组成部分：第一个是生成曼德勃罗集，第二个是将曼德勃罗集转储到 PNG 文件中。在第一章中，我们意识到我们只能并行生成曼德勃罗集，并且考虑到这占程序运行时间的大部分，这将是一个很好的候选算法，可以将其转移到 GPU 上。让我们看看如何做到这一点。（我们将避免重复定义曼德勃罗集，因此如果您需要更深入的复习，请重新阅读第一章“为什么使用 GPU 编程？”中的“曼德勃罗重新审视”部分）

首先，让我们基于原始程序中的simple_mandelbrot创建一个新的 Python 函数。我们将其称为gpu_mandelbrot，这将接受与之前完全相同的输入：

def gpu_mandelbrot(width, height, real_low, real_high, imag_low, imag_high, max_iters, upper_bound):

我们将从这里开始以稍微不同的方式进行。我们将首先构建一个复杂的晶格，其中包含我们将分析的复平面中的每个点。

在这里，我们将使用一些 NumPy 矩阵类型的技巧轻松生成晶格，然后将结果从 NumPy matrix类型转换为二维 NumPy array类型（因为 PyCUDA 只能处理 NumPy array类型，而不能处理matrix类型）。请注意我们非常小心地设置我们的 NumPy 类型：

    real_vals = np.matrix(np.linspace(real_low, real_high, width), dtype=np.complex64)
    imag_vals = np.matrix(np.linspace( imag_high, imag_low, height), dtype=np.complex64) * 1j
    mandelbrot_lattice = np.array(real_vals + imag_vals.transpose(), dtype=np.complex64)  

因此，我们现在有一个表示我们将生成曼德勃罗集的晶格的二维复杂数组；正如我们将看到的，我们可以在 GPU 内非常容易地操作这个数组。现在让我们将我们的晶格传输到 GPU，并分配一个数组来表示我们的曼德勃罗集：

    # copy complex lattice to the GPU
    mandelbrot_lattice_gpu = gpuarray.to_gpu(mandelbrot_lattice)    
    # allocate an empty array on the GPU
    mandelbrot_graph_gpu = gpuarray.empty(shape=mandelbrot_lattice.shape, dtype=np.float32)

重申一下——gpuarray.to_array函数只能操作 NumPy array类型，因此我们在将其发送到 GPU 之前一定要对其进行类型转换。接下来，我们必须使用gpuarray.empty函数在 GPU 上分配一些内存，指定数组的大小/形状和类型。同样，您可以将其视为类似于 C 中的malloc；请记住，由于gpuarray对象析构函数在作用域结束时自动处理内存清理，因此我们不必在以后释放或free这些内存。

当您使用 PyCUDA 函数gpuarray.empty或gpuarray.empty_like在 GPU 上分配内存时，由于gpuarray对象的析构函数管理所有内存清理，因此您不必在以后释放此内存。

现在我们准备启动内核；我们唯一需要做的更改是改变

我们还没有编写生成曼德勃罗集的内核函数，但让我们先写出这个函数的其余部分应该是怎样的：

    mandel_ker( mandelbrot_lattice_gpu, mandelbrot_graph_gpu, np.int32(max_iters), np.float32(upper_bound))

    mandelbrot_graph = mandelbrot_graph_gpu.get()

    return mandelbrot_graph

这就是我们希望我们的新内核行为的方式——第一个输入将是我们生成的复点阵（NumPy complex64类型），第二个将是指向二维浮点数组的指针（NumPy float32类型），它将指示哪些元素是曼德勃罗集的成员，第三个将是一个整数，表示每个点的最大迭代次数，最后一个输入将是用于确定曼德勃罗类成员资格的每个点的上限。请注意，我们在将所有输入传递给 GPU 时非常小心！

下一行将从 GPU 中检索我们生成的曼德勃罗集回到 CPU 空间，并返回结束值。（请注意，gpu_mandelbrot的输入和输出与simple_mandelbrot完全相同）。

现在让我们看看如何正确定义我们的 GPU 内核。首先，让我们在头部添加适当的include语句：

import pycuda.autoinit
from pycuda import gpuarray
from pycuda.elementwise import ElementwiseKernel

我们现在准备编写我们的 GPU 内核！我们将在这里展示它，然后逐行讨论：

mandel_ker = ElementwiseKernel(
"pycuda::complex<float> *lattice, float *mandelbrot_graph, int max_iters, float upper_bound",
"""
mandelbrot_graph[i] = 1;
pycuda::complex<float> c = lattice[i]; 
pycuda::complex<float> z(0,0);
for (int j = 0; j < max_iters; j++)
    {  
     z = z*z + c;     
     if(abs(z) > upper_bound)
         {
          mandelbrot_graph[i] = 0;
          break;
         }
    }         
""",
"mandel_ker")

首先，我们使用传递给ElementwiseKernel的第一个字符串设置我们的输入。我们必须意识到当我们在 CUDA-C 中工作时，特定的 C 数据类型将直接对应于特定的 Python NumPy 数据类型。再次注意，当数组被传递到 CUDA 内核时，它们被 CUDA 视为 C 指针。在这里，CUDA C int类型与 NumPy int32类型完全对应，而 CUDA C float类型对应于 NumPy float32类型。然后使用内部 PyCUDA 类模板进行复杂类型的转换——这里 PyCUDA ::complex<float>对应于 Numpy complex64。

让我们看看第二个字符串的内容，它用三个引号（"""）分隔。这使我们能够在字符串中使用多行；当我们在 Python 中编写更大的内联 CUDA 内核时，我们将使用这个。

虽然我们传入的数组在 Python 中是二维数组，但 CUDA 只会将它们视为一维数组，并由i索引。同样，ElementwiseKernel会自动为我们跨多个核心和线程索引i。我们将输出中的每个点初始化为 1，如mandelbrot_graph[i] = 1;，因为i将在曼德勃罗集的每个元素上进行索引；我们将假设每个点都是成员，除非证明相反。（再次说明，曼德勃罗集是在两个维度上的，实部和虚部，但ElementwiseKernel将自动将所有内容转换为一维集合。当我们再次在 Python 中与数据交互时，曼德勃罗集的二维结构将被保留。）

我们像在 Python 中一样为我们的c值设置适当的点阵点，如pycuda::complex<float> c = lattice[i];，并用pycuda::complex<float> z(0,0);将我们的z值初始化为0（第一个零对应于实部，而第二个对应于虚部）。然后我们使用一个新的迭代器j进行循环，如for(int j = 0; j < max_iters; j++)。（请注意，这个算法不会在j或任何其他索引上并行化——只有i！这个for循环将在j上串行运行——但整个代码片段将在i上并行化。）

然后，我们使用z = z*z + c;设置*z*的新值，按照曼德勃罗算法。如果这个元素的绝对值超过了上限（if(abs(z) > upper_bound)），我们将这个点设置为 0（mandelbrot_graph[i] = 0;），并用break关键字跳出循环。

在传递给ElementwiseKernel的最终字符串中，我们为内核赋予其内部 CUDA C 名称，这里是"mandel_ker"。

我们现在准备启动内核；我们唯一需要做的更改是将主函数中的引用从simple_mandelbrot更改为gpu_mandelbrot，然后我们就可以开始了。让我们从 IPython 中启动：

让我们检查转储的图像，以确保这是正确的：

这肯定是在第一章中生成的相同 Mandelbrot 图像，所以我们已经成功地将其实现到了 GPU 上！现在让我们看看我们得到的速度增加：在第一章中，我们花了 14.61 秒来生成这张图；而在这里，只花了 0.894 秒。请记住，PyCUDA 还必须在运行时编译和链接我们的 CUDA C 代码，并且需要花费时间来进行与 GPU 的内存传输。即使有了所有这些额外的开销，它仍然是一个非常值得的速度增加！（您可以在 Git 存储库中找到我们的 GPU Mandelbrot 的代码，文件名为gpu_mandelbrot0.py。）

对函数式编程的简要探讨

在我们继续之前，让我们简要回顾一下 Python 中用于函数式编程的两个函数——map和reduce。它们都被认为是函数式，因为它们都对函数进行操作。我们发现这些有趣，因为它们都对应于编程中的常见设计模式，所以我们可以替换输入中的不同函数，以获得多种不同（和有用的）操作。

让我们首先回顾 Python 中的lambda关键字。这允许我们定义一个匿名函数——在大多数情况下，这些可以被视为一次性函数，或者只希望使用一次的函数，或者可以在一行上定义的函数。让我们现在打开 IPython 并定义一个将数字平方的小函数，如pow2 = lambda x : x**2。让我们在一些数字上测试一下：

让我们回顾一下map作用于两个输入值：一个函数和给定函数可以作用的对象列表。map输出原始列表中每个元素的函数输出列表。现在让我们将我们的平方操作定义为一个匿名函数，然后将其输入到 map 中，并使用最后几个数字的列表进行检查，如map(lambda x : x**2, [2,3,4])：

我们看到map作为ElementwiseKernel！这实际上是函数式编程中的标准设计模式。现在，让我们看看reduce；它不是接收一个列表并直接输出相应列表，而是接收一个列表，在其上执行递归二进制操作，并输出一个单例。让我们通过键入reduce(lambda x, y : x + y, [1,2,3,4])来了解这种设计模式。当我们在 IPython 中键入这个时，我们将看到这将输出一个单个数字 10，这确实是1+2+3+4的和。您可以尝试用乘法替换上面的求和，并看到这确实适用于递归地将一长串数字相乘。一般来说，我们使用可结合的二进制操作进行缩减操作；这意味着，无论我们在列表的连续元素之间以何种顺序执行操作，都将始终得到相同的结果，前提是列表保持有序。（这与交换律不同。）

现在我们将看看 PyCUDA 如何处理类似于reduce的编程模式——使用并行扫描和归约内核。

并行扫描和归约内核基础

让我们看一下 PyCUDA 中一个复制 reduce 功能的基本函数——InclusiveScanKernel。（您可以在simple_scankernal0.py文件名下找到代码。）让我们执行一个在 GPU 上对一小组数字求和的基本示例：

import numpy as np
import pycuda.autoinit
from pycuda import gpuarray
from pycuda.scan import InclusiveScanKernel
seq = np.array([1,2,3,4],dtype=np.int32)
seq_gpu = gpuarray.to_gpu(seq)
sum_gpu = InclusiveScanKernel(np.int32, "a+b")
print sum_gpu(seq_gpu).get()
print np.cumsum(seq)

我们通过首先指定输入/输出类型（这里是 NumPy int32）和字符串"a+b"来构建我们的内核。在这里，InclusiveScanKernel自动在 GPU 空间中设置了名为a和b的元素，因此您可以将此字符串输入视为 Python 中的lambda a,b: a + b的类似物。我们实际上可以在这里放置任何（可结合的）二进制操作，只要我们记得用 C 语言编写它。

当我们运行sum_gpu时，我们会得到一个与输入数组大小相同的数组。数组中的每个元素表示计算中的每个步骤的值（我们可以看到，NumPy cumsum函数给出了相同的输出）。最后一个元素将是我们正在寻找的最终输出，对应于 reduce 的输出：

让我们尝试一些更具挑战性的东西；让我们找到一个float32数组中的最大值：

import numpy as np
import pycuda.autoinit
from pycuda import gpuarray
from pycuda.scan import InclusiveScanKernel
seq = np.array([1,100,-3,-10000, 4, 10000, 66, 14, 21],dtype=np.int32)
seq_gpu = gpuarray.to_gpu(seq)
max_gpu = InclusiveScanKernel(np.int32, "a > b ? a : b")
print max_gpu(seq_gpu).get()[-1]
print np.max(seq)

（您可以在名为simple_scankernal1.py的文件中找到完整的代码。）

在这里，我们所做的主要更改是用a > b ? a : b替换了a + b字符串。 （在 Python 中，这将在reduce语句中呈现为lambda a, b: max(a,b)）。在这里，我们使用了一个技巧，使用 C 语言的?运算符来给出a和b中的最大值。最后，我们显示了输出数组中结果元素的最后一个值，这将恰好是最后一个元素（我们总是可以用 Python 中的[-1]索引来检索）。

现在，让我们最后再看一个用于生成 GPU 内核的 PyCUDA 函数——ReductionKernel。实际上，ReductionKernel的作用类似于ElementwiseKernel函数，后面跟着一个并行扫描内核。哪种算法是使用ReductionKernel实现的一个好选择？首先想到的是线性代数中的点积。让我们记住计算两个向量的点积有两个步骤：

1. 将向量逐点相乘

2. 对结果的逐点乘积求和

这两个步骤也称为乘法和累加。现在让我们设置一个内核来执行这个计算：

dot_prod = ReductionKernel(np.float32, neutral="0", reduce_expr="a+b", map_expr="vec1[i]*vec2[i]", arguments="float *vec1, float *vec2") 

首先，注意我们为内核使用的数据类型（float32）。然后，我们使用arguments设置了我们的 CUDA C 内核的输入参数（这里是两个代表每个向量的浮点数组，用float *表示），并使用map_expr设置了逐点计算，这里是逐点乘法。与ElementwiseKernel一样，这是按i索引的。我们设置了reduce_expr，与InclusiveScanKernel一样。这将对数组执行元素操作的结果进行减少类型的操作。最后，我们使用neutral设置了中性元素。这是一个将作为reduce_expr的标识的元素；在这里，我们设置neutral=0，因为0在加法下始终是标识（在乘法下，1 是标识）。稍后在本书中更深入地讨论并行前缀时，我们将看到为什么我们必须设置这个。

摘要

我们首先学习了如何从 PyCUDA 查询我们的 GPU，并用此方法在 Python 中重新创建了 CUDA 的deviceQuery程序。然后我们学习了如何使用 PyCUDA 的gpuarray类及其to_gpu和get函数将 NumPy 数组传输到 GPU 的内存中。我们通过观察如何使用gpuarray对象来进行基本的 GPU 计算来感受了使用gpuarray对象的感觉，并且我们学会了使用 IPython 的prun分析器进行一些调查工作。我们发现，由于 PyCUDA 启动 NVIDIA 的nvcc编译器来编译内联 CUDA C 代码，有时在会话中首次运行 PyCUDA 的 GPU 函数时会出现一些任意的减速。然后我们学习了如何使用ElementwiseKernel函数来编译和启动逐元素操作，这些操作会自动并行化到 GPU 上。我们对 Python 中的函数式编程进行了简要回顾（特别是map和reduce函数），最后，我们介绍了如何使用InclusiveScanKernel和ReductionKernel函数在 GPU 上进行一些基本的 reduce/scan 类型计算。

现在我们已经掌握了编写和启动内核函数的绝对基础知识，我们应该意识到 PyCUDA 已经通过其模板为我们覆盖了编写内核的大部分开销。我们将在下一章学习 CUDA 内核执行的原则，以及 CUDA 如何将内核中的并发线程排列成抽象的网格和块。

问题

1. 在simple_element_kernel_example0.py中，我们在测量 GPU 计算时间时不考虑与 GPU 之间的内存传输。尝试使用 Python 时间命令测量gpuarray函数to_gpu和get的时间。考虑内存传输时间后，你会认为将这个特定函数卸载到 GPU 上值得吗？

2. 在第一章中，为什么进行 GPU 编程？，我们讨论了安德尔定律，这让我们对将程序的部分内容卸载到 GPU 上可能获得的收益有了一些了解。在本章中我们看到的两个问题，安德尔定律没有考虑到的是什么？

3. 修改gpu_mandel0.py以使用越来越小的复数格点，并将其与程序的 CPU 版本进行比较。我们可以选择足够小的格点，以至于 CPU 版本实际上比 GPU 版本更快吗？

4. 创建一个使用ReductionKernel的内核，该内核在 GPU 上获取两个相同长度的complex64数组，并返回两个数组中的绝对最大元素。

5. 如果一个gpuarray对象在 Python 中到达作用域的末尾会发生什么？

6. 你认为我们在使用ReductionKernel时为什么需要定义neutral？

7. 如果在ReductionKernel中我们设置reduce_expr ="a > b ? a : b"，并且我们正在操作 int32 类型，那么我们应该将"neutral"设置为什么？

第四章：核心、线程、块和网格

本章中，我们将看到如何编写有效的CUDA 核心。在 GPU 编程中，核心（我们可以互换使用术语CUDA 核心或核心函数）是一个可以直接从主机（CPU）启动到设备（GPU）的并行函数，而设备函数是一个只能从核心函数或另一个设备函数调用的函数。（一般来说，设备函数看起来和行为像普通的串行 C/C++函数，只是它们在 GPU 上运行，并且从核心函数并行调用。）

然后我们将了解 CUDA 如何使用线程、块和网格的概念来抽象 GPU 的一些基础技术细节（例如核心、warp 和流多处理器，我们将在本书的后面部分介绍），以及我们如何使用这些概念来减轻并行编程中的认知负担。我们将学习关于线程同步（块级和网格级），以及在 CUDA 中使用全局和共享****内存进行线程间通信。最后，我们将深入了解如何在 GPU 上实现我们自己的并行前缀类型算法（即我们在上一章中介绍的扫描/归约类型函数），这使我们能够将本章学到的所有原则付诸实践。

本章的学习成果如下：

• 理解核心和设备函数之间的区别

• 如何在 PyCUDA 中编译和启动核心，并在核心内使用设备函数

• 在启动核心的上下文中有效使用线程、块和网格，以及如何在核心内使用threadIdx和blockIdx

• 如何以及为什么在核心内同步线程，使用__syncthreads()来同步单个块中的所有线程，以及主机来同步整个块网格中的所有线程

• 如何使用设备全局和共享内存进行线程间通信

• 如何使用我们新获得的关于核心的所有知识来正确实现并行前缀和的 GPU 版本

技术要求

本章需要一台带有现代 NVIDIA GPU（2016 年以后）的 Linux 或 Windows 10 PC，并安装了所有必要的 GPU 驱动程序和 CUDA Toolkit（9.0 及以上）。还需要一个合适的 Python 2.7 安装（如 Anaconda Python 2.7），并安装了 PyCUDA 模块。

本章的代码也可以在 GitHub 上找到：

github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA

有关先决条件的更多信息，请查看本书的前言；有关软件和硬件要求，请查看github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA中的README部分。

核心

与上一章一样，我们将学习如何在 Python 代码中以内联 CUDA C 编写 CUDA 核心函数，并使用 PyCUDA 将它们启动到我们的 GPU 上。在上一章中，我们使用 PyCUDA 提供的模板来编写符合特定设计模式的核心，相比之下，我们现在将看到如何从头开始编写我们自己的核心，以便我们可以编写各种各样的核心，这些核心可能不属于 PyCUDA 涵盖的任何特定设计模式，并且我们可以更精细地控制我们的核心。当然，这些收益将以编程复杂性增加为代价；我们特别需要了解线程、块和网格及其在核心中的作用，以及如何同步我们的核心正在执行的线程，以及如何在线程之间交换数据。

让我们从简单开始，尝试重新创建我们在上一章中看到的一些逐元素操作，但这次不使用ElementwiseKernel函数；我们现在将使用SourceModule函数。这是 PyCUDA 中非常强大的函数，允许我们从头构建一个内核，所以通常最好从简单开始。

PyCUDA SourceModule 函数

我们将使用 PyCUDA 的SourceModule函数将原始内联 CUDA C 代码编译为可用的内核，我们可以从 Python 中启动。我们应该注意，SourceModule实际上将代码编译为CUDA 模块，这类似于 Python 模块或 Windows DLL，只是它包含一组编译的 CUDA 代码。这意味着我们必须使用 PyCUDA 的get_function“提取”我们想要使用的内核的引用，然后才能实际启动它。让我们从如何使用SourceModule的基本示例开始。

与以前一样，我们将从制作最简单的内核函数之一开始，即将向量乘以标量。我们将从导入开始：

import pycuda.autoinit
import pycuda.driver as drv
import numpy as np
from pycuda import gpuarray
from pycuda.compiler import SourceModule

现在我们可以立即开始编写我们的内核：

ker = SourceModule("""
__global__ void scalar_multiply_kernel(float *outvec, float scalar, float *vec)
{
 int i = threadIdx.x;
 outvec[i] = scalar*vec[i];
}
""")

因此，让我们停下来，对比一下在ElementwiseKernel中是如何完成的。首先，在 CUDA C 中声明内核函数时，我们要在前面加上__global__关键字。这将使编译器将该函数标识为内核。我们总是将其声明为void函数，因为我们总是通过传递指向一些空内存块的指针来获得输出值。我们可以像声明任何标准 C 函数一样声明参数：首先是outvec，这将是我们的输出缩放向量，当然是浮点数组指针。接下来是scalar，用一个简单的float表示；注意这不是一个指针！如果我们希望将简单的单例输入值传递给我们的内核，我们总是可以在不使用指针的情况下这样做。最后，我们有我们的输入向量vec，当然是另一个浮点数组指针。

单例输入参数可以直接从主机传递给内核函数，而无需使用指针或分配的设备内存。

让我们在继续测试之前先深入了解内核。我们记得ElementwiseKernel会自动并行化多个 GPU 线程，通过 PyCUDA 为我们设置的值i；每个单独线程的标识由threadIdx值给出，我们可以通过以下方式检索：int i = threadIdx.x;。

threadIdx用于告诉每个单独的线程其身份。这通常用于确定应在输入和输出数据数组上处理哪些值的索引。（这也可以用于使用标准 C 控制流语句（如if或switch）为特定线程分配不同的任务。）

现在，我们准备像以前一样并行执行标量乘法：outvec[i] = scalar*vec[i];。

现在，让我们测试这段代码：我们首先必须从我们刚刚使用SourceModule编译的 CUDA 模块中提取编译的内核函数的引用。我们可以使用 Python 的get_function来获取这个内核引用，如下所示：

scalar_multiply_gpu = ker.get_function("scalar_multiply_kernel")

现在，我们必须在 GPU 上放一些数据来实际测试我们的内核。让我们设置一个包含 512 个随机值的浮点数组，然后使用gpuarray.to_gpu函数将这些值复制到 GPU 的全局内存中的数组中。（我们将在 GPU 和 CPU 上将这个随机向量乘以一个标量，并查看输出是否匹配。）我们还将使用gpuarray.empty_like函数在 GPU 的全局内存中分配一块空内存块：

testvec = np.random.randn(512).astype(np.float32)
testvec_gpu = gpuarray.to_gpu(testvec)
outvec_gpu = gpuarray.empty_like(testvec_gpu)

现在，我们准备启动我们的内核。我们将标量值设置为2。（再次，由于标量是单例，我们不必将该值复制到 GPU，但是我们必须小心确保正确地进行类型转换。）在这里，我们必须使用block和grid参数明确设置线程数为512。我们现在准备启动：

scalar_multiply_gpu( outvec_gpu, np.float32(2), testvec_gpu, block=(512,1,1), grid=(1,1,1))

现在我们可以使用gpuarray输出对象中的get函数来检查输出是否与预期输出匹配，并将其与 NumPy 的allclose函数进行比较：

print "Does our kernel work correctly? : {}".format(np.allclose(outvec_gpu.get() , 2*testvec) )

（此示例的代码可在存储库中的simple_scalar_multiply_kernel.py文件中的4下找到。）

现在我们开始去掉前一章中学到的 PyCUDA 内核模板的训练轮——我们现在可以直接用纯 CUDA C 编写内核，并启动它在 GPU 上使用特定数量的线程。但是，在继续使用内核之前，我们必须更多地了解 CUDA 如何将线程结构化为抽象单位块和网格。

线程、块和网格

到目前为止，在本书中，我们一直认为线程这个术语是理所当然的。让我们退后一步，看看这究竟意味着——线程是在 GPU 的单个核心上执行的一系列指令—核心和线程不应被视为同义词！事实上，可以启动使用的线程数量比 GPU 上的核心数量多得多。这是因为，类似于英特尔芯片可能只有四个核心，但在 Linux 或 Windows 中运行数百个进程和数千个线程，操作系统的调度程序可以快速在这些任务之间切换，使它们看起来是同时运行的。GPU 以类似的方式处理线程，允许在成千上万的线程上进行无缝计算。

多个线程在 GPU 上以抽象单位块执行。您应该回忆一下我们如何从标量乘法内核中的threadIdx.x获得线程 ID；末尾有一个x，因为还有threadIdx.y和threadIdx.z。这是因为您可以在三个维度上对块进行索引，而不仅仅是一个维度。为什么我们要这样做？让我们回忆一下有关从第一章中计算 Mandelbrot 集的示例，为什么使用 GPU 编程？和第三章，使用 PyCUDA 入门。这是在二维平面上逐点计算的。因此，对于这样的算法，我们可能更倾向于在两个维度上对线程进行索引。同样，在某些情况下，使用三个维度可能是有意义的——在物理模拟中，我们可能需要在 3D 网格内计算移动粒子的位置。

块进一步以称为网格的抽象批次执行，最好将其视为块的块。与块中的线程一样，我们可以使用blockIdx.x、blockIdx.y和blockIdx.z给出的常量值在网格中的最多三个维度上对每个块进行索引。让我们看一个示例来帮助我们理解这些概念；为了简单起见，我们这里只使用两个维度。

康威的生命游戏

《生命游戏》（通常简称为 LIFE）是一种细胞自动机模拟，由英国数学家约翰·康威于 1970 年发明。听起来很复杂，但实际上非常简单——LIFE 是一个零玩家的“游戏”，由一个二维二进制格子组成，其中的“细胞”被认为是“活着的”或“死了的”。这个格子通过以下一组规则进行迭代更新：

• 任何活细胞周围少于两个活邻居的细胞会死亡

• 任何活细胞周围有两个或三个邻居的细胞会存活

• 任何活细胞周围有三个以上的邻居的细胞会死亡

• 任何死细胞周围恰好有三个邻居的细胞会复活

这四条简单的规则产生了一个复杂的模拟，具有有趣的数学特性，而且在动画时也非常美观。然而，在晶格中有大量的细胞时，它可能运行得很慢，并且通常在纯串行 Python 中编程时会导致动画不流畅。然而，这是可以并行化的，因为很明显晶格中的每个细胞可以由一个单独的 CUDA 线程管理。

我们现在将 LIFE 实现为一个 CUDA 核函数，并使用 matplotlib.animation 模块来进行动画。这对我们来说现在很有趣，因为我们将能够在这里应用我们对块和网格的新知识。

我们将首先包括适当的模块，如下所示：

import pycuda.autoinit
import pycuda.driver as drv
from pycuda import gpuarray
from pycuda.compiler import SourceModule
import numpy as np
import matplotlib.pyplot as plt 
import matplotlib.animation as animation

现在，让我们通过 SourceModule 来编写我们的核函数。我们将首先使用 C 语言的 #define 指令来设置一些我们将在整个核函数中使用的常量和宏。让我们看看我们将设置的前两个，_X 和 _Y：

ker = SourceModule("""
#define _X  ( threadIdx.x + blockIdx.x * blockDim.x )
#define _Y  ( threadIdx.y + blockIdx.y * blockDim.y )

首先让我们记住这里 #define 的工作原理——它会在编译时用定义的值（在括号中）直接替换 _X 或 _Y 的任何文本，也就是说，它为我们创建了宏。（作为个人风格的问题，我通常会在所有的 C 宏之前加上下划线。）

在 C 和 C++中，#define 用于创建宏。这意味着 #define 不会创建任何函数或设置正确的常量变量——它只允许我们在编译之前通过交换文本来在我们的代码中以简写方式编写东西。

现在，让我们具体讨论一下 _X 和 _Y 的含义——这将是单个 CUDA 线程在我们用于 LIFE 的二维晶格上的笛卡尔 x 和 y 值。我们将在一个二维网格上启动核函数，由二维块组成，这些块将对应整个细胞晶格。我们将使用线程和块常量来找到晶格上的笛卡尔点。让我们看一些图表来说明这一点。驻留在二维 CUDA 块中的线程可以被可视化如下：

此时，你可能会想知道为什么我们不在一个单独的块上启动我们的核函数，这样我们就可以将 _X 设置为 threadIdx.x，将 _Y 设置为 threadIdx.y，然后就完成了。这是由于 CUDA 对我们施加了块大小的限制——目前只支持由最多 1024 个线程组成的块。这意味着我们只能将我们的细胞晶格的尺寸最大设为 32 x 32，这将导致一个相当无聊的模拟，最好在 CPU 上完成，所以我们将在网格上启动多个块。（我们当前块的尺寸将由 blockDim.x 和 blockDim.y 给出，这将帮助我们确定目标 x 和 y 坐标，正如我们将看到的。）

同样，和之前一样，我们可以确定我们在二维网格中的块是哪个，使用 blockIdx.x 和 blockIdx.y：

在我们稍微思考一下数学之后，应该很清楚 _X 应该被定义为 (threadIdx.x + blockIdx.x * blockDim.x)，而 _Y 应该被定义为 (threadIdx.y + blockIdx.y * blockDim.y)。（添加括号是为了不干扰宏插入代码时的运算顺序。）现在，让我们继续定义剩下的宏：

#define _WIDTH  ( blockDim.x * gridDim.x )
#define _HEIGHT ( blockDim.y * gridDim.y  )

#define _XM(x)  ( (x + _WIDTH) % _WIDTH )
#define _YM(y)  ( (y + _HEIGHT) % _HEIGHT )

_WIDTH和_HEIGHT宏将分别给出我们单元格格子的宽度和高度，这应该从图表中清楚地看出。让我们讨论_XM和_YM宏。在我们的 LIFE 实现中，我们将端点“环绕”到格子的另一侧 - 例如，我们将考虑-1的x值为_WIDTH - 1，y值为-1为_HEIGHT - 1，同样地，我们将考虑_WIDTH的x值为0，y值为_HEIGHT为0。我们为什么需要这个？当我们计算给定单元格的存活邻居数时，我们可能处于某个边缘，邻居可能是外部点 - 定义这些宏来调制我们的点将自动为我们覆盖这一点。请注意，在使用 C 的模运算符之前，我们必须添加宽度或高度 - 这是因为，与 Python 不同，C 中的模运算符对于整数可以返回负值。

我们现在有一个最终的宏要定义。我们记得 PyCUDA 将二维数组作为一维指针传递到 CUDA C 中；二维数组从 Python 以按行的方式传递到一维 C 指针中。这意味着我们必须将格子上给定的笛卡尔（x，y）点转换为指向格子对应的指针中的一维点。在这里，我们可以这样做：

#define _INDEX(x,y)  ( _XM(x)  + _YM(y) * _WIDTH )

由于我们的单元格格子是按行存储的，我们必须将y值乘以宽度以偏移到对应行的点。我们现在终于可以开始我们的 LIFE 实现了。让我们从 LIFE 最重要的部分开始 - 计算给定单元格的存活邻居数。我们将使用 CUDA 设备函数来实现这一点，如下所示：

__device__ int nbrs(int x, int y, int * in)
{
     return ( in[ _INDEX(x -1, y+1) ] + in[ _INDEX(x-1, y) ] + in[ _INDEX(x-1, y-1) ] \
                   + in[ _INDEX(x, y+1)] + in[_INDEX(x, y - 1)] \
                   + in[ _INDEX(x+1, y+1) ] + in[ _INDEX(x+1, y) ] + in[ _INDEX(x+1, y-1) ] );
}

设备函数是以串行方式编写的 C 函数，由内核中的单个 CUDA 线程调用。也就是说，这个小函数将由我们内核中的多个线程并行调用。我们将把我们的单元格格子表示为 32 位整数的集合（1 表示活细胞，0 表示死细胞），所以这对我们的目的是有效的；我们只需要添加周围当前单元格的值。

CUDA 设备函数是由内核中的单个 CUDA 线程调用的串行 C 函数。虽然这些函数本身是串行的，但它们可以由多个 GPU 线程并行运行。设备函数本身不能由主机计算机启动到 GPU 上，只能由内核启动。

我们现在准备编写 LIFE 的内核实现。实际上，我们已经完成了大部分的艰苦工作 - 我们检查当前线程单元格的邻居数量，检查当前单元格是生还是死，然后使用适当的 switch-case 语句来根据 LIFE 的规则确定下一次迭代的状态。我们将使用两个整数指针数组作为内核 - 一个将用作输入参考上一次迭代（lattice），另一个将用作输出我们将计算的迭代（lattice_out）的参考。

__global__ void conway_ker(int * lattice_out, int * lattice  )
{
   // x, y are the appropriate values for the cell covered by this thread
   int x = _X, y = _Y;

   // count the number of neighbors around the current cell
   int n = nbrs(x, y, lattice);

    // if the current cell is alive, then determine if it lives or dies for the next generation.
    if ( lattice[_INDEX(x,y)] == 1)
       switch(n)
       {
          // if the cell is alive: it remains alive only if it has 2 or 3 neighbors.
          case 2:
          case 3: lattice_out[_INDEX(x,y)] = 1;
                  break;
          default: lattice_out[_INDEX(x,y)] = 0;                   
       }
    else if( lattice[_INDEX(x,y)] == 0 )
         switch(n)
         {
            // a dead cell comes to life only if it has 3 neighbors that are alive.
            case 3: lattice_out[_INDEX(x,y)] = 1;
                    break;
            default: lattice_out[_INDEX(x,y)] = 0;         
         }

}
""")

conway_ker = ker.get_function("conway_ker")

我们记得用三重括号关闭内联 CUDA C 段落，然后用get_function获取对我们的 CUDA C 内核的引用。由于内核只会一次更新格子，我们将在 Python 中设置一个简短的函数，它将涵盖更新格子的所有开销以用于动画：

def update_gpu(frameNum, img, newLattice_gpu, lattice_gpu, N):    

frameNum参数只是 Matplotlib 动画模块对于我们可以忽略的更新函数所需的一个值，而img将是我们单元格格子的代表图像，这是动画模块所需的，将被迭代显示。

让我们关注另外三个参数—newLattice_gpu和lattice_gpu将是我们保持持久的 PyCUDA 数组，因为我们希望避免在 GPU 上重新分配内存块。lattice_gpu将是细胞数组的当前一代，对应于内核中的lattice参数，而newLattice_gpu将是下一代晶格。N将指示晶格的高度和宽度（换句话说，我们将使用N x N晶格）。

我们使用适当的参数启动内核，并设置块和网格大小如下：

    conway_ker(newLattice_gpu, lattice_gpu, grid=(N/32,N/32,1), block=(32,32,1) )    

我们将块大小设置为 32 x 32，使用(32, 32, 1)；因为我们只使用两个维度来表示我们的细胞晶格，所以我们可以将z维度设置为 1。请记住，块的线程数限制为 1,024 个线程—32 x 32 = 1024，所以这样可以工作。（请记住，32 x 32 没有什么特别之处；如果需要，我们可以使用 16 x 64 或 10 x 10 等值，只要总线程数不超过 1,024。）

CUDA 块中的线程数最多为 1,024。

现在我们来看一下网格值—在这里，因为我们使用 32 的维度，很明显N（在这种情况下）应该是 32 的倍数。这意味着在这种情况下，我们只能使用 64 x 64、96 x 96、128 x 128 和 1024 x 1024 等晶格。同样，如果我们想使用不同大小的晶格，那么我们将不得不改变块的维度。（如果这不太清楚，请查看之前的图表，并回顾一下我们如何在内核中定义宽度和高度宏。）

现在我们可以使用get()函数从 GPU 的内存中获取最新生成的晶格，并为我们的动画设置图像数据。最后，我们使用 PyCUDA 切片操作[:]将新的晶格数据复制到当前数据中，这将复制 GPU 上先前分配的内存，这样我们就不必重新分配了：

    img.set_data(newLattice_gpu.get() )    
    lattice_gpu[:] = newLattice_gpu[:]

    return img

让我们设置一个大小为 256 x 256 的晶格。现在我们将使用numpy.random模块中的 choice 函数为我们的晶格设置初始状态。我们将随机用 1 和 0 填充一个N x N的整数图表；通常，如果大约 25%的点是 1，其余的是 0，我们可以生成一些有趣的晶格动画，所以我们就这样做吧：

if __name__ == '__main__':
    # set lattice size
    N = 256

    lattice = np.int32( np.random.choice([1,0], N*N, p=[0.25, 0.75]).reshape(N, N) )
    lattice_gpu = gpuarray.to_gpu(lattice)

最后，我们可以使用适当的gpuarray函数在 GPU 上设置晶格，并相应地设置 Matplotlib 动画，如下所示：

lattice_gpu = gpuarray.to_gpu(lattice)
    lattice_gpu = gpuarray.to_gpu(lattice)
    newLattice_gpu = gpuarray.empty_like(lattice_gpu) 

    fig, ax = plt.subplots()
    img = ax.imshow(lattice_gpu.get(), interpolation='nearest')
    ani = animation.FuncAnimation(fig, update_gpu, fargs=(img, newLattice_gpu, lattice_gpu, N, ) , interval=0, frames=1000, save_count=1000) 

    plt.show()

现在我们可以运行我们的程序并享受展示（代码也可以在 GitHub 存储库的4目录下的conway_gpu.py文件中找到）：

线程同步和互通

现在我们将讨论 GPU 编程中的两个重要概念—线程同步和线程互通。有时，我们需要确保每个线程在继续任何进一步的计算之前都已经到达了代码中完全相同的行；我们称之为线程同步。同步与线程互通相辅相成，也就是说，不同的线程之间传递和读取输入；在这种情况下，我们通常希望确保所有线程在传递数据之前都处于计算的相同步骤。我们将从学习 CUDA __syncthreads设备函数开始，该函数用于同步内核中的单个块。

使用 __syncthreads()设备函数

在我们之前的康威生命游戏的示例中，我们的内核每次由主机启动时只更新了晶格一次。在这种情况下，同步所有在启动的内核中的线程没有问题，因为我们只需要处理已经准备好的晶格的上一个迭代。

现在假设我们想做一些稍微不同的事情——我们想重新编写我们的内核，以便在给定的细胞点阵上执行一定数量的迭代，而不是由主机一遍又一遍地重新启动。这一开始可能看起来很琐碎——一个天真的解决方案将是只需在内联conway_ker内核中放置一个整数参数来指示迭代次数和一个for循环，进行一些额外的琐碎更改，然后就完成了。

然而，这引发了竞争条件的问题；这是多个线程读取和写入相同内存地址以及由此可能产生的问题。我们的旧conway_ker内核通过使用两个内存数组来避免这个问题，一个严格用于读取，一个严格用于每次迭代写入。此外，由于内核只执行单次迭代，我们实际上是在使用主机来同步线程。

我们希望在 GPU 上进行多次完全同步的 LIFE 迭代；我们还希望使用单个内存数组来存储点阵。我们可以通过使用 CUDA 设备函数__syncthreads()来避免竞争条件。这个函数是一个块级同步屏障——这意味着在一个块内执行的每个线程在到达__syncthreads()实例时都会停止，并等待直到同一块内的每个其他线程都到达__syncthreads()的同一调用，然后线程才会继续执行后续的代码行。

__syncthreads()只能同步单个 CUDA 块内的线程，而不能同步 CUDA 网格内的所有线程！

让我们现在创建我们的新内核；这将是之前 LIFE 内核的修改，它将执行一定数量的迭代然后停止。这意味着我们不会将其表示为动画，而是作为静态图像，因此我们将在开始时加载适当的 Python 模块。（此代码也可在 GitHub 存储库的conway_gpu_syncthreads.py文件中找到）：

import pycuda.autoinit
import pycuda.driver as drv
from pycuda import gpuarray
from pycuda.compiler import SourceModule
import numpy as np
import matplotlib.pyplot as plt 

现在，让我们再次设置计算 LIFE 的内核：

ker = SourceModule("""

当然，我们的 CUDA C 代码将放在这里，这将大致与以前相同。我们只需要对内核进行一些更改。当然，我们可以保留设备函数nbrs。在我们的声明中，我们将只使用一个数组来表示细胞点阵。我们可以这样做，因为我们将使用适当的线程同步。我们还必须用一个整数表示迭代次数。我们设置参数如下：

__global__ void conway_ker(int * lattice, int iters) {

我们将继续与以前类似，只是使用for循环进行迭代：

 int x = _X, y = _Y; 
 for (int i = 0; i < iters; i++)
 {
     int n = nbrs(x, y, lattice); 
     int cell_value;

让我们回想一下以前，我们直接在数组中设置新的细胞点阵值。在这里，我们将在cell_value变量中保存值，直到块内的所有线程都同步。我们以前也是类似地进行，使用__syncthreads阻止执行，直到确定了当前迭代的所有新细胞值，然后才在点阵数组中设置值：

 if ( lattice[_INDEX(x,y)] == 1)
 switch(n)
 {
 // if the cell is alive: it remains alive only if it has 2 or 3 neighbors.
 case 2:
 case 3: cell_value = 1;
 break;
 default: cell_value = 0; 
 }
 else if( lattice[_INDEX(x,y)] == 0 )
 switch(n)
 {
 // a dead cell comes to life only if it has 3 neighbors that are alive.
 case 3: cell_value = 1;
 break;
 default: cell_value = 0; 
 } 
 __syncthreads();
 lattice[_INDEX(x,y)] = cell_value; 
 __syncthreads();
 } 
}
""")

我们现在将像以前一样启动内核并显示输出，迭代点阵 100 万次。请注意，由于每个块的线程限制为 1,024 个，我们在网格中只使用一个块，大小为 32 x 32。（再次强调，__syncthreads仅在块中的所有线程上工作，而不是在网格中的所有线程上工作，这就是为什么我们在这里限制自己使用单个块的原因）：

conway_ker = ker.get_function("conway_ker")
if __name__ == '__main__':
 # set lattice size
 N = 32
 lattice = np.int32( np.random.choice([1,0], N*N, p=[0.25, 0.75]).reshape(N, N) )
 lattice_gpu = gpuarray.to_gpu(lattice)
 conway_ker(lattice_gpu, np.int32(1000000), grid=(1,1,1), block=(32,32,1))
 fig = plt.figure(1)
 plt.imshow(lattice_gpu.get())

当我们运行程序时，我们将得到以下所需的输出（这是随机 LIFE 点阵在一百万次迭代后会收敛到的结果！）：

使用共享内存

从先前的例子中，我们可以看到内核中的线程可以使用 GPU 全局内存中的数组进行相互通信；虽然可以使用全局内存进行大多数操作，但使用共享内存可以加快速度。这是一种专门用于单个 CUDA 块内线程相互通信的内存类型；与全局内存相比，使用共享内存的优势在于纯线程间通信速度更快。不过，与全局内存相反，存储在共享内存中的内存不能直接被主机访问——共享内存必须首先由内核自己复制回全局内存。

在继续之前，让我们先退一步思考一下我们的意思。让我们看看我们刚刚看到的迭代 LIFE 内核中声明的一些变量。首先看看x和y，这两个整数保存着特定线程单元的笛卡尔坐标。请记住，我们正在使用_X和_Y宏设置它们的值。（尽管编译器优化，我们希望将这些值存储在变量中以减少计算，因为直接使用_X和_Y将在我们的代码中引用这些宏时每次重新计算x和y的值）：

 int x = _X, y = _Y; 

我们注意到，对于每个单个线程，点阵中将对应于x和y的唯一笛卡尔点。同样，我们使用一个变量n，它声明为int n = nbrs(x, y, lattice);，来表示特定单元周围的存活邻居的数量。这是因为，当我们通常在 CUDA 中声明变量时，它们默认是每个单独线程的本地变量。请注意，即使我们在线程内部声明数组如int a[10];，也将有一个大小为 10 的数组，它是每个线程的本地数组。

本地线程数组（例如，在内核内部声明int a[10];）和指向全局 GPU 内存的指针（例如，以int * b形式作为内核参数传递的值）可能看起来和行为类似，但实际上非常不同。对于内核中的每个线程，将有一个单独的a数组，其他线程无法读取，但将有一个单独的b，它将保存相同的值，并且对所有线程都是同样可访问的。

我们准备使用共享内存。这使我们能够声明在单个 CUDA 块内的线程之间共享的变量和数组。这种内存比使用全局内存指针（我们到目前为止一直在使用的）要快得多，同时减少了指针分配内存的开销。

假设我们想要一个大小为 10 的共享整数数组。我们声明如下——__shared__ int a[10] 。请注意，我们不必局限于数组；我们可以按如下方式创建共享的单例变量：__shared__ int x。

让我们重新编写 LIFE 的迭代版本的一些行，以利用共享内存。首先，让我们将输入指针重命名为p_lattice，这样我们可以在我们的共享数组上使用这个变量名，并在我们的代码中懒惰地保留所有对lattice的引用。由于我们将坚持使用 32 x 32 个单元的点阵，我们设置新的共享lattice数组如下：

__global__ void conway_ker_shared(int * p_lattice, int iters)
{
 int x = _X, y = _Y;
 __shared__ int lattice[32*32];

现在我们必须将全局内存p_lattice数组中的所有值复制到lattice中。我们将以完全相同的方式索引我们的共享数组，因此我们可以在这里使用我们旧的_INDEX宏。请注意，在复制后我们确保在我们继续 LIFE 算法之前放置__syncthreads()，以确保所有对 lattice 的内存访问完全完成：

 lattice[_INDEX(x,y)] = p_lattice[_INDEX(x,y)];
 __syncthreads();

内核的其余部分与以前完全相同，只是我们必须将共享的点阵复制回 GPU 数组。我们这样做，然后关闭内联代码：

 __syncthreads();
 p_lattice[_INDEX(x,y)] = lattice[_INDEX(x,y)];
 __syncthreads();
} """)

现在我们可以像以前一样运行，使用完全相同的测试代码。（此示例可以在 GitHub 存储库中的conway_gpu_syncthreads_shared.py中找到。）

并行前缀算法

现在，我们将利用我们对 CUDA 核心的新知识来实现并行前缀算法，也称为扫描设计模式。我们已经在上一章中以 PyCUDA 的InclusiveScanKernel和ReductionKernel函数的形式看到了这种简单的例子。现在让我们更详细地了解这个想法。

这个设计模式的核心动机是，我们有一个二元运算符，也就是说，一个作用于两个输入值并给出一个输出值的函数（比如—+，，（最大值），（最小值）），和元素的集合，我们希望从中高效地计算。此外，我们假设我们的二元运算符是可结合的—这意味着，对于任意三个元素x、y和z，我们总是有:。

我们希望保留部分结果，也就是n-1个子计算—。并行前缀算法的目的是高效地产生这些n个和。在串行操作中，通常需要O(n)的时间来产生这些n个和，我们希望降低时间复杂度。

当使用术语“并行前缀”或“扫描”时，通常意味着一个产生所有这些n个结果的算法，而“减少”/“归约”通常意味着只产生单个最终结果，。（这是 PyCUDA 的情况。）

实际上，并行前缀算法有几种变体，我们将首先从最简单（也是最古老的）版本开始，这就是所谓的天真并行前缀算法。

天真并行前缀算法

这个算法是原始版本的天真并行前缀算法；这个算法是“天真的”，因为它假设给定n个输入元素，，进一步假设n是二进制的（也就是说，对于某个正整数k），我们可以在n个处理器（或n个线程）上并行运行算法。显然，这将对我们可以处理的集合的基数n施加严格的限制。然而，只要满足这些条件，我们就有一个很好的结果，即其计算时间复杂度仅为O(log n)。我们可以从算法的伪代码中看到这一点。在这里，我们将用表示输入值，用表示输出值：

input: x0, ..., xn-1 initialize:
for k=0 to n-1:
    yk := xk begin:
parfor i=0 to n-1 :
    for j=0 to log2(n):
        if i >= 2j :
            yi := yi  yi - 2^j else:
            continue
        end if
    end for
end parfor
end
output: y0, ..., yn-1

现在，我们可以清楚地看到这将花费O(log n)的渐近时间，因为外部循环是在parfor上并行化的，内部循环需要log2。经过几分钟的思考，很容易看出y[i]值将产生我们期望的输出。

现在让我们开始我们的实现；在这里，我们的二元运算符将简单地是加法。由于这个例子是说明性的，这个核心代码将严格地针对 1,024 个线程。

让我们先设置好头文件，然后开始编写我们的核心代码：

import pycuda.autoinit
import pycuda.driver as drv
import numpy as np
from pycuda import gpuarray
from pycuda.compiler import SourceModule
from time import time

naive_ker = SourceModule("""
__global__ void naive_prefix(double *vec, double *out)
{
     __shared__ double sum_buf[1024]; 
     int tid = threadIdx.x; 
     sum_buf[tid] = vec[tid];

因此，让我们看看我们有什么：我们将我们的输入元素表示为双精度 GPU 数组，即double *vec，并用double *out表示输出值。我们声明一个共享内存sum_buf数组，我们将用它来计算我们的输出。现在，让我们看看算法本身的实现：

 int iter = 1;
 for (int i=0; i < 10; i++)
 {
     __syncthreads();
     if (tid >= iter )
     {
         sum_buf[tid] = sum_buf[tid] + sum_buf[tid - iter]; 
     } 
     iter *= 2;
 }
 __syncthreads();

当然，这里没有parfor，它是隐式的，通过tid变量表示线程编号。我们还可以通过从初始化为 1 的变量开始，然后在每次 i 的迭代中乘以 2 来省略log[2]和2^i的使用。（请注意，如果我们想更加技术化，我们可以使用位移运算符来实现这一点。）我们将i的迭代限制在 10 次，因为2¹⁰ = 1024。现在我们将结束我们的新内核如下：

 __syncthreads();
 out[tid] = sum_buf[tid];
 __syncthreads();

}
""")
naive_gpu = naive_ker.get_function("naive_prefix")

现在让我们看一下内核后面的测试代码：

if __name__ == '__main__':
 testvec = np.random.randn(1024).astype(np.float64)
 testvec_gpu = gpuarray.to_gpu(testvec)

 outvec_gpu = gpuarray.empty_like(testvec_gpu)
 naive_gpu( testvec_gpu , outvec_gpu, block=(1024,1,1), grid=(1,1,1))

 total_sum = sum( testvec)
 total_sum_gpu = outvec_gpu[-1].get()

 print "Does our kernel work correctly? : {}".format(np.allclose(total_sum_gpu , total_sum) )

我们只关心输出中的最终总和，我们使用outvec_gpu[-1].get()来检索它，需要回想一下，"-1"索引给出了 Python 数组中的最后一个成员。这将是vec中每个元素的总和；部分总和在outvec_gpu的先前值中。（此示例可以在 GitHub 存储库中的naive_prefix.py文件中看到。）

由于其性质，并行前缀算法必须在n个线程上运行，对应于一个大小为 n 的数组，其中n是二进制的（这意味着n是 2 的某个幂）。然而，我们可以将这个算法扩展到任意非二进制大小，假设我们的运算符具有单位元素（或等效地，中性元素）——也就是说，存在某个值e，使得对于任何x值，我们有 。在运算符为+的情况下，单位元素是 0；在运算符为  的情况下，它是 1；然后我们只需用一系列e值填充  的元素，以便我们有新集合的二进制基数 。

包含与排他前缀

让我们停下来，做一个非常微妙但非常重要的区分。到目前为止，我们一直关心接受形式为  的输入，并产生形式为  的总和数组作为输出。产生这种输出的前缀算法称为包含；在包含前缀算法的情况下，每个索引处的相应元素包含在输出数组的相同索引处的总和中。这与排他前缀算法形成对比。排他前缀算法不同之处在于，它同样接受形式为  的n个输入值，并产生长度为n的输出数组 。

这很重要，因为一些高效的前缀算法的变体天生就是排他的。我们将在下一个小节中看到一个例子。

请注意，排他算法产生的输出与包含算法几乎相同，只是右移并省略了最后一个值。因此，我们可以从任一算法中轻松获得等效输出，只要我们保留  的副本。

一个高效的并行前缀算法

在我们继续使用新算法之前，我们将从两个角度看待朴素算法。在理想情况下，计算时间复杂度为O(log n)，但这仅当我们有足够数量的处理器用于我们的数据集时成立；当数据集的基数（元素数量）n远大于处理器数量时，这将成为O(n log n)时间复杂度的算法。

让我们定义一个与我们的二进制运算符 相关的新概念——这里并行算法的工作是执行期间所有线程对此运算符的调用次数。同样，跨度是线程在内核执行期间进行调用的次数；而整个算法的跨度与每个单独线程的最长跨度相同，这将告诉我们总执行时间。

我们寻求特别减少算法在所有线程中执行的工作量，而不仅仅是关注跨度。在简单前缀的情况下，所需的额外工作在可用处理器数量不足时会花费更多时间；这些额外工作将溢出到有限数量的可用处理器中。

我们将介绍一种新的算法，它是工作高效的，因此更适合有限数量的处理器。这包括两个独立的部分——向上扫描（或减少）阶段和向下扫描阶段。我们还应该注意，我们将看到的算法是一种独占前缀算法。

向上扫描阶段类似于单个减少操作，以产生由减少算法给出的值，即  ；在这种情况下，我们保留所需的部分和（）以实现最终结果。然后，向下扫描阶段将对这些部分和进行操作，并给出最终结果。让我们看一些伪代码，从向上扫描阶段开始。（接下来的小节将立即深入到伪代码的实现中。）

工作高效的并行前缀（向上扫描阶段）

这是向上扫描的伪代码。（注意parfor覆盖j变量，这意味着此代码块可以并行化，由j索引的线程）：

input: x0, ..., xn-1initialize:
    for i = 0 to n - 1:
        yi := xi
begin:
for k=0 to log2(n) - 1:
    parfor j=0 to n - 1: 
        if j is divisible by 2k+1:
            yj+2^(k+1)-1 = yj+2^k-1  yj +2^(k+1) -1else:
            continue
end
output: y0, ..., yn-1

工作高效的并行前缀（向下扫描阶段）

现在让我们继续向下扫描，它将操作向上扫描的输出：

input: x0, ..., xn-1 initialize:
    for i = 0 to n - 2:
        yi := xi
    yn-1 := 0
begin:
for k = log2(n) - 1 to 0:
    parfor j = 0 to n - 1: 
        if j is divisible by 2k+1:
            temp := yj+2^k-1
            yj+2^k-1 := yj+2^(k+1)-1
            yj+2^(k+1)-1 := yj+2^(k+1)-1  temp
        else:
            continue
end
output: y0 , y1 , ..., yn-1

工作高效的并行前缀 — 实现

作为本章的压轴，我们将编写该算法的实现，该算法可以在大于 1,024 的任意大小的数组上操作。这意味着这将在网格和块上操作；因此，我们将不得不使用主机进行同步；此外，这将要求我们为向上扫描和向下扫描阶段实现两个单独的内核，这将作为两个阶段的parfor循环，以及作为向上和向下扫描的外部for循环的 Python 函数。

让我们从向上扫描内核开始。由于我们将从主机迭代重新启动此内核，我们还需要一个指示当前迭代（k）的参数。我们将使用两个数组进行计算，以避免竞争条件——x（用于当前迭代）和x_old（用于之前的迭代）。我们声明内核如下：

up_ker = SourceModule("""
__global__ void up_ker(double *x, double *x_old, int k)
{

现在让我们设置tid变量，它将是网格中所有块中所有线程中当前线程的标识。我们使用与我们之前看到的康威生命游戏的原始网格级实现相同的技巧：

int tid =  blockIdx.x*blockDim.x + threadIdx.x;

我们现在将使用 C 位移运算符直接从k生成 2^k 和 2^(k+1 )。我们现在将j设置为tid乘以_2k1—这将使我们能够删除“如果j可被 2^(k+1)整除”的部分，如伪代码中所示，从而使我们只启动所需数量的线程：

 int _2k = 1 << k;
 int _2k1 = 1 << (k+1);

 int j = tid* _2k1;

我们可以使用 CUDA C 中的左位移运算符（<<）轻松生成二进制（2 的幂）整数。请记住，整数 1（即 2⁰）表示为 0001，2（2¹）表示为 0010，4（2²）表示为 0100，依此类推。因此，我们可以使用1 << k操作计算 2^k。

现在我们可以运行上扫描阶段的单行代码，注意j确实可以被 2^(k+1)整除，因为它的构造方式是这样的：

 x[j + _2k1 - 1] = x_old[j + _2k -1 ] + x_old[j + _2k1 - 1];
}
""")

我们已经编写完我们的核心代码了！但这当然不是完整的上扫描实现。我们还需要在 Python 中完成其余部分。让我们拿到我们的核心代码并开始实现。这基本上是按照伪代码进行的，我们应该记住，我们通过使用[:]从x_gpu复制到x_old_gpu来更新x_old_gpu，这将保留内存分配，并仅复制新数据而不是重新分配。还要注意，我们根据要启动的线程数量设置我们的块和网格大小 - 我们尝试保持我们的块大小为 32 的倍数（这是本文中的经验法则，我们在第十一章中详细介绍为什么我们特别使用 32，CUDA 性能优化）。我们应该在文件开头加上from __future__ import division，因为我们将使用 Python 3 风格的除法来计算我们的块和核心大小。

有一点需要提到的是，我们假设x的长度是二进制长度 32 或更大 - 如果您希望将其操作在其他大小的数组上，可以通过用零填充我们的数组来轻松修改，然而：

up_gpu = up_ker.get_function("up_ker")

def up_sweep(x):
    x = np.float64(x)
    x_gpu = gpuarray.to_gpu(np.float64(x) )
    x_old_gpu = x_gpu.copy()
    for k in range( int(np.log2(x.size) ) ) : 
        num_threads = int(np.ceil( x.size / 2**(k+1)))
        grid_size = int(np.ceil(num_threads / 32))

        if grid_size > 1:
            block_size = 32
        else:
            block_size = num_threads

        up_gpu(x_gpu, x_old_gpu, np.int32(k) , block=(block_size,1,1), grid=(grid_size,1,1))
        x_old_gpu[:] = x_gpu[:]

    x_out = x_gpu.get()
    return(x_out)

现在我们将开始编写下扫描。同样，让我们从核心开始，它将具有伪代码中内部parfor循环的功能。它与之前类似 - 再次使用两个数组，因此在这里使用temp变量与伪代码中是不必要的，我们再次使用位移运算符来获得 2^k 和 2^(k+1)的值。我们计算j与之前类似：

down_ker = SourceModule("""
__global__ void down_ker(double *y, double *y_old, int k)
{
 int j = blockIdx.x*blockDim.x + threadIdx.x;

 int _2k = 1 << k;
 int _2k1 = 1 << (k+1);

 int j = tid*_2k1;

 y[j + _2k - 1 ] = y_old[j + _2k1 - 1];
 y[j + _2k1 - 1] = y_old[j + _2k1 - 1] + y_old[j + _2k - 1];
}
""")

down_gpu = down_ker.get_function("down_ker")

现在我们可以编写一个 Python 函数，该函数将迭代地启动核心，这对应于下扫描阶段的外部for循环。这类似于上扫描阶段的 Python 函数。从伪代码中看到的一个重要区别是，我们必须从外部for循环中的最大值迭代到最小值；我们可以使用 Python 的reversed函数来做到这一点。现在我们可以实现下扫描阶段：

def down_sweep(y):
    y = np.float64(y)
    y[-1] = 0
    y_gpu = gpuarray.to_gpu(y)
    y_old_gpu = y_gpu.copy()
    for k in reversed(range(int(np.log2(y.size)))):
        num_threads = int(np.ceil( y.size / 2**(k+1)))
        grid_size = int(np.ceil(num_threads / 32))

        if grid_size > 1:
            block_size = 32
        else:
            block_size = num_threads

        down_gpu(y_gpu, y_old_gpu, np.int32(k), block=(block_size,1,1), grid=(grid_size,1,1))
        y_old_gpu[:] = y_gpu[:]
    y_out = y_gpu.get()
    return(y_out)

在实现了上扫描和下扫描阶段之后，我们的最后任务是轻而易举地完成：

def efficient_prefix(x):
        return(down_sweep(up_sweep(x)))

我们现在已经完全实现了一个与主机同步的高效并行前缀算法！（这个实现可以在存储库中的work-efficient_prefix.py文件中找到，还有一些测试代码。）

总结

我们从康威的生命游戏实现开始，这给了我们一个关于 CUDA 核心的许多线程是如何在块-网格张量类型结构中组织的想法。然后，我们通过 CUDA 函数__syncthreads()深入研究了块级同步，以及通过使用共享内存进行块级线程互通；我们还看到单个块有一定数量的线程，我们可以操作，所以在创建将使用多个块跨越更大网格的核心时，我们必须小心使用这些功能。

我们概述了并行前缀算法的理论，并最后实现了一个天真的并行前缀算法，作为一个单个核心，可以操作大小受限的数组，该数组与___syncthreads同步，并在内部执行for和parfor循环，并且实现了一个高效的并行前缀算法，该算法跨两个核心和三个 Python 函数实现，可以操作任意大小的数组，核心充当算法的内部parfor循环，而 Python 函数有效地充当外部for循环并同步核心启动。

问题

1. 更改simple_scalar_multiply_kernel.py中的随机向量，使其长度为 10,000，并修改内核定义中的i索引，以便它可以在网格形式的多个块中使用。看看现在是否可以通过将块和网格参数设置为block=(100,1,1)和grid=(100,1,1)来启动这个内核超过 10,000 个线程。

2. 在上一个问题中，我们启动了一个同时使用 10,000 个线程的内核；截至 2018 年，没有 NVIDIA GPU 具有超过 5,000 个核心。为什么这仍然有效并产生预期的结果？

3. 天真的并行前缀算法在数据集大小为n的情况下，时间复杂度为 O(log n)，假设我们有n个或更多处理器。假设我们在具有 640 个核心的 GTX 1050 GPU 上使用天真的并行前缀算法。在n >> 640的情况下，渐近时间复杂度会变成什么？

4. 修改naive_prefix.py以在大小为 1,024 的数组上运行（可能是非二进制的）。

5. __syncthreads() CUDA 设备函数只在单个块中同步线程。我们如何在整个网格的所有块中同步所有线程？

6. 您可以通过这个练习说服自己，第二个前缀和算法确实比天真的前缀和算法更有效。假设我们有一个大小为 32 的数据集。在这种情况下，第一个和第二个算法所需的“加法”操作的确切数量是多少？

7. 在高效并行前缀的实现中，我们使用 Python 函数来迭代我们的内核并同步结果。为什么我们不能在内核中使用for循环，同时小心地使用__syncthreads()？

8. 为什么在 CUDA C 中实现天真的并行前缀在一个单独的内核中处理自己的同步比使用两个内核和 Python 函数实现高效并行前缀更有意义，并且让主机处理同步更有意义？

第五章：流、事件、上下文和并发性

在之前的章节中，我们看到在与 GPU 交互时，主机执行的两个主要操作是：

• 将内存数据从 GPU 复制到 GPU，以及从 GPU 复制到内存

• 启动内核函数

我们知道，在单个内核中，其许多线程之间存在一定程度的并发性；然而，在多个内核和 GPU 内存操作之间，还存在另一种并发性。这意味着我们可以同时启动多个内存和内核操作，而无需等待每个操作完成。然而，另一方面，我们必须有一定的组织能力，以确保所有相互依赖的操作都得到同步；这意味着我们不应该在其输入数据完全复制到设备内存之前启动特定的内核，或者在内核执行完成之前，不应该将已启动内核的输出数据复制到主机。

为此，我们有所谓的CUDA 流—流是在 GPU 上按顺序运行的一系列操作。单独一个流是没有用的—重点是通过使用多个流在主机上发出的 GPU 操作来获得并发性。这意味着我们应该交错启动与不同流对应的 GPU 操作，以利用这个概念。

我们将在本章中广泛涵盖流的概念。此外，我们还将研究事件，这是流的一个特性，用于精确计时内核，并指示主机在给定流中已完成哪些操作。

最后，我们将简要介绍 CUDA 上下文。上下文可以被视为类似于操作系统中的进程，因为 GPU 将每个上下文的数据和内核代码隔离并封装起来，使其与当前存在于 GPU 上的其他上下文相分离。我们将在本章末尾看到这方面的基础知识。

本章的学习成果如下：

• 理解设备和流同步的概念

• 学习如何有效地使用流来组织并发的 GPU 操作

• 学习如何有效地使用 CUDA 事件

• 理解 CUDA 上下文

• 学习如何在给定上下文中显式同步

• 学习如何显式创建和销毁 CUDA 上下文

• 学习如何使用上下文允许主机上的多个进程和线程共享 GPU 使用

技术要求

本章需要一台带有现代 NVIDIA GPU（2016 年以后）的 Linux 或 Windows 10 PC，并安装了所有必要的 GPU 驱动程序和 CUDA Toolkit（9.0 及以上）。还需要一个合适的 Python 2.7 安装（如 Anaconda Python 2.7），并安装了 PyCUDA 模块。

本章的代码也可以在 GitHub 上找到：

github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA

有关先决条件的更多信息，请查看本书的前言，有关软件和硬件要求，请查看github.com/PacktPublishing/Hands-On-GPU-Programming-with-Python-and-CUDA中的 README。

CUDA 设备同步

在我们可以使用 CUDA 流之前，我们需要了解设备同步的概念。这是一种操作，其中主机阻塞任何进一步的执行，直到发给 GPU 的所有操作（内存传输和内核执行）都已完成。这是为了确保依赖于先前操作的操作不会被无序执行—例如，确保 CUDA 内核启动在主机尝试读取其输出之前已完成。

在 CUDA C 中，设备同步是通过cudaDeviceSynchronize函数执行的。这个函数有效地阻止了主机上的进一步执行，直到所有 GPU 操作完成。cudaDeviceSynchronize是如此基本，以至于它通常是 CUDA C 大多数书籍中最早涉及的主题之一-我们还没有看到这一点，因为 PyCUDA 已经在需要时自动为我们调用了这个函数。让我们看一个 CUDA C 代码的例子，看看如何手动完成这个操作：

// Copy an array of floats from the host to the device.
cudaMemcpy(device_array, host_array, size_of_array*sizeof(float), cudaMemcpyHostToDevice);
// Block execution until memory transfer to device is complete.
cudaDeviceSynchronize();
// Launch CUDA kernel.
Some_CUDA_Kernel <<< block_size, grid_size >>> (device_array, size_of_array);
// Block execution until GPU kernel function returns.
cudaDeviceSynchronize();
// Copy output of kernel to host.
cudaMemcpy(host_array,  device_array, size_of_array*sizeof(float), cudaMemcpyDeviceToHost); // Block execution until memory transfer to host is complete.
cudaDeviceSynchronize();

在这段代码块中，我们看到我们必须在每个单独的 GPU 操作之后直接与设备进行同步。如果我们只需要一次调用单个 CUDA 内核，就像这里看到的那样，这是可以的。但是，如果我们想要同时启动多个独立的内核和操作不同数据数组的内存操作，跨整个设备进行同步将是低效的。在这种情况下，我们应该跨多个流进行同步。我们现在就来看看如何做到这一点。

使用 PyCUDA 流类

我们将从一个简单的 PyCUDA 程序开始；它只是生成一系列随机的 GPU 数组，对每个数组进行简单的内核处理，然后将数组复制回主机。然后我们将修改它以使用流。请记住，这个程序将完全没有任何意义，只是为了说明如何使用流以及一些基本的性能提升。（这个程序可以在 GitHub 存储库的5目录下的multi-kernel.py文件中看到。）

当然，我们将首先导入适当的 Python 模块，以及time函数：

import pycuda.autoinit
import pycuda.driver as drv
from pycuda import gpuarray
from pycuda.compiler import SourceModule
import numpy as np
from time import time

我们现在将指定要处理多少个数组-在这里，每个数组将由不同的内核启动进行处理。我们还指定要生成的随机数组的长度，如下所示：

num_arrays = 200
array_len = 1024**2

我们现在有一个在每个数组上操作的内核；它只是迭代数组中的每个点，并将其乘以 2 再除以 2，重复 50 次，最终保持数组不变。我们希望限制每个内核启动将使用的线程数，这将帮助我们在 GPU 上的许多内核启动之间获得并发性，以便我们可以让每个线程通过for循环迭代数组的不同部分。（再次提醒，这个内核函数除了用来学习流和同步之外，完全没有任何用处！）如果每个内核启动使用的线程太多，以后获得并发性将更加困难：

ker = SourceModule(""" 
__global__ void mult_ker(float * array, int array_len)
{
     int thd = blockIdx.x*blockDim.x + threadIdx.x;
     int num_iters = array_len / blockDim.x;

     for(int j=0; j < num_iters; j++)
     {
         int i = j * blockDim.x + thd;

         for(int k = 0; k < 50; k++)
         {
              array[i] *= 2.0;
              array[i] /= 2.0;
         }
     }
}
""")

mult_ker = ker.get_function('mult_ker')

现在，我们将生成一些随机数据数组，将这些数组复制到 GPU，迭代地在每个数组上启动我们的内核，然后将输出数据复制回主机，并使用 NumPy 的allclose函数来断言它们是否相同。我们将使用 Python 的time函数来计算从开始到结束的所有操作的持续时间：

data = []
data_gpu = []
gpu_out = []

# generate random arrays.
for _ in range(num_arrays):
    data.append(np.random.randn(array_len).astype('float32'))

t_start = time()

# copy arrays to GPU.
for k in range(num_arrays):
    data_gpu.append(gpuarray.to_gpu(data[k]))

# process arrays.
for k in range(num_arrays):
    mult_ker(data_gpu[k], np.int32(array_len), block=(64,1,1), grid=(1,1,1))

# copy arrays from GPU.
for k in range(num_arrays):
    gpu_out.append(data_gpu[k].get())

t_end = time()

for k in range(num_arrays):
    assert (np.allclose(gpu_out[k], data[k]))

print 'Total time: %f' % (t_end - t_start)

我们现在准备运行这个程序。我现在就要运行它：

所以，这个程序完成需要了将近三秒的时间。我们将进行一些简单的修改，使我们的程序可以使用流，然后看看我们是否可以获得任何性能提升（这可以在存储库中的multi-kernel_streams.py文件中看到）。

首先，我们注意到对于每个内核启动，我们有一个单独的数据数组进行处理，这些数组存储在 Python 列表中。我们将不得不为每个单独的数组/内核启动对创建一个单独的流对象，所以让我们首先添加一个空列表，名为streams，用来保存我们的流对象：

data = []
data_gpu = []
gpu_out = []
streams = []

我们现在可以生成一系列流，我们将使用它们来组织内核启动。我们可以使用pycuda.driver子模块和Stream类来获取流对象。由于我们已经导入了这个子模块并将其别名为drv，我们可以用以下方式用新的流对象填充我们的列表：

for _ in range(num_arrays):
    streams.append(drv.Stream())

现在，我们将首先修改将数据传输到 GPU 的内存操作。考虑以下步骤：

1. 查找第一个循环，使用gpuarray.to_gpu函数将数组复制到 GPU。我们将希望切换到此函数的异步和适用于流的版本gpu_array.to_gpu_async。（现在我们还必须使用stream参数指定每个内存操作应该使用哪个流）：

for k in range(num_arrays):
    data_gpu.append(gpuarray.to_gpu_async(data[k], stream=streams[k]))

1. 我们现在可以启动我们的内核。这与以前完全相同，只是我们必须通过使用stream参数来指定要使用哪个流：

for k in range(num_arrays):
    mult_ker(data_gpu[k], np.int32(array_len), block=(64,1,1), grid=(1,1,1), stream=streams[k])

1. 最后，我们需要从 GPU 中取出数据。我们可以通过将gpuarray get函数切换为get_async来实现这一点，并再次使用stream参数，如下所示：

for k in range(num_arrays):
    gpu_out.append(data_gpu[k].get_async(stream=streams[k]))

我们现在准备运行我们适用于流的修改后的程序：

在这种情况下，我们获得了三倍的性能提升，考虑到我们需要进行的修改非常少，这并不算太糟糕。但在我们继续之前，让我们尝试更深入地理解为什么这样做有效。

让我们考虑两个 CUDA 内核启动的情况。在我们启动内核之前和之后，我们还将执行与每个内核相对应的 GPU 内存操作，总共有六个操作。我们可以通过图表来可视化在 GPU 上随时间发生的操作，如此—在x轴上向右移动对应于时间持续时间，而y轴对应于在特定时间执行的 GPU 上的操作。这可以用以下图表来描述：

很容易想象为什么流在性能提高方面效果如此好——因为单个流中的操作直到所有必要的先前操作完成之后才会被阻塞，我们将获得不同 GPU 操作之间的并发性，并充分利用我们的设备。这可以通过并发操作的大量重叠来看出。我们可以将基于流的并发性随时间的变化可视化如下：

使用 CUDA 流进行并发康威生命游戏

我们现在将看到一个更有趣的应用——我们将修改上一章的 LIFE（康威的生命游戏）模拟，以便我们将同时显示四个独立的动画窗口。（建议您查看上一章的这个示例，如果您还没有的话。）

让我们从上一章的存储库中获取旧的 LIFE 模拟，应该在4目录中的conway_gpu.py下。现在我们将把这个修改为基于新的 CUDA 流的并发 LIFE 模拟。（我们将在稍后看到的这种基于流的新模拟也可以在本章目录5中的conway_gpu_streams.py文件中找到。）

转到文件末尾的主函数。我们将设置一个新变量，用num_concurrent表示我们将同时显示多少个并发动画（其中N表示模拟栅格的高度/宽度，与以前一样）。我们在这里将其设置为4，但您可以随意尝试其他值：

if __name__ == '__main__':

    N = 128
    num_concurrent = 4

我们现在需要一组num_concurrent流对象，并且还需要在 GPU 上分配一组输入和输出栅格。当然，我们只需将这些存储在列表中，并像以前一样初始化栅格。我们将设置一些空列表，并通过循环填充每个适当的对象，如下所示（请注意我们如何在每次迭代中设置一个新的初始状态栅格，将其发送到 GPU，并将其连接到lattices_gpu）：

streams = []
lattices_gpu = []
newLattices_gpu = []

for k in range(num_concurrent):
    streams.append(drv.Stream())
    lattice = np.int32( np.random.choice([1,0], N*N, p=[0.25, 0.75]).reshape(N, N) )
    lattices_gpu.append(gpuarray.to_gpu(lattice)) 
    newLattices_gpu.append(gpuarray.empty_like(lattices_gpu[k])) 

由于我们只在程序启动期间执行此循环一次，并且几乎所有的计算工作都将在动画循环中进行，因此我们实际上不必担心实际使用我们刚刚生成的流。

我们现在将使用 Matplotlib 使用子图函数设置环境；请注意，我们可以通过设置ncols参数来设置多个动画图。我们将有另一个列表结构，它将对应于动画更新中所需的图像imgs。请注意，我们现在可以使用get_async和相应的流来设置这个：

fig, ax = plt.subplots(nrows=1, ncols=num_concurrent)
imgs = []

for k in range(num_concurrent):
    imgs.append( ax[k].imshow(lattices_gpu[k].get_async(stream=streams[k]), interpolation='nearest') )

在主函数中要更改的最后一件事是以ani = animation.FuncAnimation开头的倒数第二行。让我们修改update_gpu函数的参数，以反映我们正在使用的新列表，并添加两个参数，一个用于传递我们的streams列表，另一个用于指示应该有多少并发动画的参数：

ani = animation.FuncAnimation(fig, update_gpu, fargs=(imgs, newLattices_gpu, lattices_gpu, N, streams, num_concurrent) , interval=0, frames=1000, save_count=1000)    

现在我们需要对update_gpu函数进行必要的修改以接受这些额外的参数。在文件中向上滚动一点，并按照以下方式修改参数：

def update_gpu(frameNum, imgs, newLattices_gpu, lattices_gpu, N, streams, num_concurrent):

现在，我们需要修改这个函数，使其迭代num_concurrent次，并像以前一样设置imgs的每个元素，最后返回整个imgs列表：

for k in range(num_concurrent):
    conway_ker( newLattices_gpu[k], lattices_gpu[k], grid=(N/32,N/32,1), block=(32,32,1), stream=streams[k] )
     imgs[k].set_data(newLattices_gpu[k].get_async(stream=streams[k]) )
     lattices_gpu[k].set_async(newLattices_gpu[k], stream=streams[k])

 return imgs

注意我们所做的更改——每个内核都在适当的流中启动，而get已经切换为与相同流同步的get_async。

最后，在循环中的最后一行将 GPU 数据从一个设备数组复制到另一个设备数组，而不进行任何重新分配。在此之前，我们可以使用简写切片操作符[:]直接在数组之间复制元素，而不在 GPU 上重新分配任何内存；在这种情况下，切片操作符符号充当 PyCUDA set函数的别名，用于 GPU 数组。 （当然，set是将一个 GPU 数组复制到另一个相同大小的数组，而不进行任何重新分配的函数。）幸运的是，确实有一个流同步版本的这个函数，set_async，但我们需要明确地使用这个函数来调用这个函数，明确指定要复制的数组和要使用的流。

我们现在已经完成并准备运行。转到终端并在命令行输入python conway_gpu_streams.py来观看展示：

事件

事件是存在于 GPU 上的对象，其目的是作为操作流的里程碑或进度标记。事件通常用于在设备端精确计时操作的时间持续时间；到目前为止，我们所做的测量都是使用基于主机的 Python 分析器和标准 Python 库函数，如time。此外，事件还可以用于向主机提供有关流状态和已完成的操作的状态更新，以及用于显式基于流的同步。

让我们从一个不使用显式流并使用事件来测量仅一个单个内核启动的示例开始。（如果我们的代码中没有显式使用流，CUDA 实际上会隐式定义一个默认流，所有操作都将放入其中）。

在这里，我们将使用与本章开头相同的无用的乘法/除法循环内核和标题，并修改大部分以下内容。我们希望为此示例运行一个长时间的单个内核实例，因此我们将生成一个巨大的随机数数组供内核处理，如下所示：

array_len = 100*1024**2
data = np.random.randn(array_len).astype('float32')
data_gpu = gpuarray.to_gpu(data)

现在，我们使用pycuda.driver.Event构造函数构造我们的事件（当然，pycuda.driver已经被我们之前的导入语句别名为drv）。

我们将在这里创建两个事件对象，一个用于内核启动的开始，另一个用于内核启动的结束（我们将始终需要两个事件对象来测量任何单个 GPU 操作，很快我们将看到）：

start_event = drv.Event()
end_event = drv.Event()

现在，我们准备启动我们的内核，但首先，我们必须使用事件记录函数标记start_event实例在执行流中的位置。我们启动内核，然后标记end_event在执行流中的位置，并且使用record：

start_event.record()
mult_ker(data_gpu, np.int32(array_len), block=(64,1,1), grid=(1,1,1))
end_event.record()

事件具有一个二进制值，指示它们是否已经到达，这是由函数 query 给出的。让我们在内核启动后立即为两个事件打印状态更新：

print 'Has the kernel started yet? {}'.format(start_event.query())
 print 'Has the kernel ended yet? {}'.format(end_event.query())

现在让我们运行这个，看看会发生什么：

我们的目标是最终测量内核执行的时间持续，但内核似乎还没有启动。在 PyCUDA 中，内核是异步启动的（无论它们是否存在于特定流中），因此我们必须确保我们的主机代码与 GPU 正确同步。

由于end_event是最后一个，我们可以通过此事件对象的synchronize函数阻止进一步的主机代码执行，直到内核完成；这将确保在执行任何进一步的主机代码之前内核已经完成。让我们在适当的位置添加一行代码来做到这一点：

end_event.synchronize()

print 'Has the kernel started yet?  {}'.format(start_event.query())

print 'Has the kernel ended yet? {}'.format(end_event.query())

最后，我们准备测量内核的执行时间；我们可以使用事件对象的time_till或time_since操作来与另一个事件对象进行比较，以获取这两个事件之间的时间（以毫秒为单位）。让我们使用start_event的time_till操作在end_event上：

print 'Kernel execution time in milliseconds: %f ' % start_event.time_till(end_event)

时间持续可以通过 GPU 上已经发生的两个事件之间的time_till和time_since函数来测量。请注意，这些函数始终以毫秒为单位返回值！

现在让我们再次运行我们的程序：

（此示例也可在存储库中的simple_event_example.py文件中找到。）

事件和流

我们现在将看到如何在流方面使用事件对象；这将使我们对各种 GPU 操作的流程具有高度复杂的控制，使我们能够准确了解每个单独流的进展情况，甚至允许我们在忽略其他流的情况下与主机同步特定流。

首先，我们必须意识到这一点——每个流必须有自己专用的事件对象集合；多个流不能共享一个事件对象。让我们通过修改之前的示例multi_kernel_streams.py来确切了解这意味着什么。在内核定义之后，让我们添加两个额外的空列表——start_events和end_events。我们将用事件对象填充这些列表，这些事件对象将对应于我们拥有的每个流。这将允许我们在每个流中计时一个 GPU 操作，因为每个 GPU 操作都需要两个事件：

data = []
data_gpu = []
gpu_out = []
streams = []
start_events = []
end_events = []

for _ in range(num_arrays):
    streams.append(drv.Stream())
    start_events.append(drv.Event())
    end_events.append(drv.Event())

现在，我们可以通过修改第二个循环来逐个计时每个内核启动，使用事件的开始和结束记录。请注意，由于这里有多个流，我们必须将适当的流作为参数输入到每个事件对象的record函数中。还要注意，我们可以在第二个循环中捕获结束事件；这仍然可以完美地捕获内核执行持续时间，而不会延迟启动后续内核。现在考虑以下代码：

for k in range(num_arrays):
    start_events[k].record(streams[k])
    mult_ker(data_gpu[k], np.int32(array_len), block=(64,1,1), grid=(1,1,1), stream=streams[k])

for k in range(num_arrays):
    end_events[k].record(streams[k])

现在我们将提取每个单独内核启动的持续时间。让我们在迭代断言检查之后添加一个新的空列表，并通过time_till函数将其填充为持续时间：

kernel_times = []
for k in range(num_arrays):
   kernel_times.append(start_events[k].time_till(end_events[k]))

现在让我们在最后添加两个print语句，告诉我们内核执行时间的平均值和标准偏差：

print 'Mean kernel duration (milliseconds): %f' % np.mean(kernel_times)
print 'Mean kernel standard deviation (milliseconds): %f' % np.std(kernel_times)

我们现在可以运行这个：

（此示例也可在存储库中的multi-kernel_events.py中找到。）

我们看到内核持续时间的标准偏差相对较低，这是好事，考虑到每个内核在相同的块和网格大小上处理相同数量的数据——如果存在较高的偏差，那将意味着我们在内核执行中使用 GPU 的不均匀程度很高，我们将不得不重新调整参数以获得更高的并发水平。

上下文

CUDA 上下文通常被描述为类似于操作系统中的进程。让我们来看看这意味着什么——进程是计算机上运行的单个程序的实例；操作系统内核之外的所有程序都在一个进程中运行。每个进程都有自己的一组指令、变量和分配的内存，一般来说，它对其他进程的操作和内存是盲目的。当一个进程结束时，操作系统内核会进行清理，确保进程分配的所有内存都已被释放，并关闭进程使用的任何文件、网络连接或其他资源。（好奇的 Linux 用户可以使用命令行top命令查看计算机上运行的进程，而 Windows 用户可以使用 Windows 任务管理器查看）。

类似于进程，上下文与正在使用 GPU 的单个主机程序相关联。上下文在内存中保存了所有正在使用的 CUDA 内核和分配的内存，并且对其他当前存在的上下文的内核和内存是盲目的。当上下文被销毁（例如，在基于 GPU 的程序结束时），GPU 会清理上下文中的所有代码和分配的内存，为其他当前和未来的上下文释放资源。到目前为止，我们编写的程序都存在于一个单一的上下文中，因此这些操作和概念对我们来说是不可见的。

让我们也记住，一个单独的程序开始时是一个单独的进程，但它可以分叉自身以在多个进程或线程中运行。类似地，一个单独的 CUDA 主机程序可以在 GPU 上生成和使用多个 CUDA 上下文。通常，当我们分叉主机进程的新进程或线程时，我们将创建一个新的上下文以获得主机端的并发性。（然而，应该强调的是，主机进程和 CUDA 上下文之间没有确切的一对一关系）。

与生活中的许多其他领域一样，我们将从一个简单的例子开始。我们将首先看看如何访问程序的默认上下文并在其间进行同步。

同步当前上下文

我们将看到如何在 Python 中显式同步我们的设备上下文，就像在 CUDA C 中一样；这实际上是 CUDA C 中最基本的技能之一，在其他大多数关于这个主题的书籍的第一章或第二章中都有涵盖。到目前为止，我们已经能够避免这个话题，因为 PyCUDA 自动使用pycuda.gpuarray函数（如to_gpu或get）执行大多数同步；否则，在本章开头我们看到的to_gpu_async或get_async函数的情况下，同步是由流处理的。

我们将谦卑地开始修改我们在第三章中编写的程序，使用 PyCUDA 入门，它使用显式上下文同步生成 Mandelbrot 集的图像。（这在存储库的3目录下的文件gpu_mandelbrot0.py中可用。）

我们在这里不会获得任何性能提升，与我们最初的 Mandelbrot 程序相比；这个练习的唯一目的是帮助我们理解 CUDA 上下文和 GPU 同步。

看一下头部，我们当然看到了import pycuda.autoinit一行。我们可以使用pycuda.autoinit.context访问当前上下文对象，并且可以通过调用pycuda.autoinit.context.synchronize()函数在当前上下文中进行同步。

现在让我们修改gpu_mandelbrot函数以处理显式同步。我们看到的第一行与 GPU 相关的代码是这样的：

mandelbrot_lattice_gpu = gpuarray.to_gpu(mandelbrot_lattice)

我们现在可以将其更改为显式同步。我们可以使用to_gpu_async异步将数据复制到 GPU，然后进行同步，如下所示：

mandelbrot_lattice_gpu = gpuarray.to_gpu_async(mandelbrot_lattice) pycuda.autoinit.context.synchronize()

然后我们看到下一行使用gpuarray.empty函数在 GPU 上分配内存。由于 GPU 架构的性质，CUDA 中的内存分配始终是自动同步的；这里没有异步内存分配的等价物。因此，我们保持这行与之前一样。

CUDA 中的内存分配始终是同步的！

现在我们看到接下来的两行 - 我们的 Mandelbrot 内核通过调用mandel_ker启动，并且我们通过调用get复制了 Mandelbrot gpuarray对象的内容。在内核启动后我们进行同步，将get切换为get_async，最后进行最后一行同步：

mandel_ker( mandelbrot_lattice_gpu, mandelbrot_graph_gpu, np.int32(max_iters), np.float32(upper_bound))
pycuda.autoinit.context.synchronize()
mandelbrot_graph = mandelbrot_graph_gpu.get_async()
pycuda.autoinit.context.synchronize()

现在我们可以运行这个程序，它将像第三章中的Getting Started with PyCUDA一样在磁盘上生成一个 Mandelbrot 图像。

（此示例也可在存储库中的gpu_mandelbrot_context_sync.py中找到。）

手动上下文创建

到目前为止，我们一直在所有 PyCUDA 程序的开头导入pycuda.autoinit；这实际上在程序开始时创建一个上下文，并在结束时销毁它。

让我们尝试手动操作。我们将制作一个小程序，只是将一个小数组复制到 GPU，然后将其复制回主机，打印数组，然后退出。

我们从导入开始：

import numpy as np
from pycuda import gpuarray
import pycuda.driver as drv

首先，我们使用pycuda.driver.init函数初始化 CUDA，这里被别名为drv：

drv.init()

现在我们选择要使用的 GPU；在拥有多个 GPU 的情况下是必要的。我们可以使用pycuda.driver.Device选择特定的 GPU；如果您只有一个 GPU，就像我一样，可以使用pycuda.driver.Device(0)访问它，如下所示：

dev = drv.Device(0)

现在我们可以使用make_context在此设备上创建一个新上下文，如下所示：

ctx = dev.make_context()

现在我们有了一个新的上下文，这将自动成为默认上下文。让我们将一个数组复制到 GPU，然后将其复制回主机并打印出来：

x = gpuarray.to_gpu(np.float32([1,2,3]))
print x.get()

现在我们完成了。我们可以通过调用pop函数销毁上下文：

ctx.pop()

就是这样！我们应该始终记得在程序退出之前销毁我们显式创建的上下文。

（此示例可以在存储库中的本章目录下的simple_context_create.py文件中找到。）

主机端多进程和多线程

当然，我们可以通过在主机的 CPU 上使用多个进程或线程来实现并发。让我们现在明确区分主机端操作系统进程和线程。

每个存在于操作系统内核之外的主机端程序都作为一个进程执行，并且也可以存在于多个进程中。进程有自己的地址空间，因为它与所有其他进程同时运行并独立运行。进程通常对其他进程的操作视而不见，尽管多个进程可以通过套接字或管道进行通信。在 Linux 和 Unix 中，使用 fork 系统调用生成新进程。

相比之下，主机端线程存在于单个进程中，多个线程也可以存在于单个进程中。单个进程中的多个线程并发运行。同一进程中的所有线程共享进程内的相同地址空间，并且可以访问相同的共享变量和数据。通常，资源锁用于在多个线程之间访问数据，以避免竞争条件。在编译语言（如 C、C++或 Fortran）中，多个进程线程通常使用 Pthreads 或 OpenMP API 进行管理。

线程比进程更轻量级，操作系统内核在单个进程中的多个线程之间切换任务要比在多个进程之间切换任务快得多。通常，操作系统内核会自动在不同的 CPU 核心上执行不同的线程和进程，以建立真正的并发性。

Python 的一个特点是，虽然它通过threading模块支持多线程，但所有线程都将在同一个 CPU 核心上执行。这是由于 Python 是一种解释脚本语言的技术细节所致，与 Python 的全局标识符锁（GIL）有关。要通过 Python 在主机上实现真正的多核并发，不幸的是，我们必须使用multiprocessing模块生成多个进程。（不幸的是，由于 Windows 处理进程的方式，multiprocessing模块目前在 Windows 下并不完全可用。Windows 用户将不幸地不得不坚持单核多线程，如果他们想要在主机端实现任何形式的并发。）

我们现在将看到如何在 Python 中使用两个线程来使用基于 GPU 的操作；Linux 用户应该注意，通过将threading的引用切换到multiprocessing，将Thread的引用切换到Process，可以很容易地将其扩展到进程，因为这两个模块看起来和行为都很相似。然而，由于 PyCUDA 的性质，我们将不得不为每个将使用 GPU 的线程或进程创建一个新的 CUDA 上下文。让我们立即看看如何做到这一点。

主机并发的多个上下文

让我们首先简要回顾一下如何在 Python 中创建一个可以通过简单示例返回值给主机的单个主机线程。（这个例子也可以在存储库中的single_thread_example.py文件的5下看到。）我们将使用threading模块中的Thread类来创建Thread的子类，如下所示：

import threading
class PointlessExampleThread(threading.Thread):

现在我们设置我们的构造函数。我们调用父类的构造函数，并在对象中设置一个空变量，这将是线程返回的值：

def __init__(self):
    threading.Thread.__init__(self)
    self.return_value = None

现在我们在我们的线程类中设置运行函数，这是在启动线程时将被执行的内容。我们只需要让它打印一行并设置返回值：

def run(self):
    print 'Hello from the thread you just spawned!'
    self.return_value = 123

最后，我们需要设置 join 函数。这将允许我们从线程中接收一个返回值：

def join(self):
    threading.Thread.join(self)
    return self.return_value

现在我们已经设置好了我们的线程类。让我们将这个类的一个实例作为NewThread对象启动，通过调用start方法来生成新的线程，然后通过调用join来阻塞执行并从主机线程获取输出：

NewThread = PointlessExampleThread()
NewThread.start()
thread_output = NewThread.join()
print 'The thread completed and returned this value: %s' % thread_output

现在让我们运行这个：

现在，我们可以在主机上的多个并发线程之间扩展这个想法，通过多个上下文和线程来启动并发的 CUDA 操作。现在我们来看最后一个例子。让我们重新使用本章开头的无意义的乘法/除法内核，并在我们生成的每个线程中启动它。

首先，让我们看一下导入部分。由于我们正在创建显式上下文，请记住删除pycuda.autoinit，并在最后添加一个threading导入：

import pycuda
import pycuda.driver as drv
from pycuda import gpuarray
from pycuda.compiler import SourceModule
import numpy as np
from time import time
import threading 

我们将使用与之前相同的数组大小，但这次线程的数量和数组的数量将直接对应。通常情况下，我们不希望在主机上生成超过 20 个左右的线程，所以我们只会使用10个数组。因此，现在考虑以下代码：

num_arrays = 10
array_len = 1024**2

现在，我们将把旧的内核存储为一个字符串对象；由于这只能在一个上下文中编译，我们将不得不在每个线程中单独编译这个内核：

kernel_code = """ 
__global__ void mult_ker(float * array, int array_len)
{
     int thd = blockIdx.x*blockDim.x + threadIdx.x;
     int num_iters = array_len / blockDim.x;
    for(int j=0; j < num_iters; j++)
     {
     int i = j * blockDim.x + thd;
     for(int k = 0; k < 50; k++)
     {
         array[i] *= 2.0;
         array[i] /= 2.0;
     }
 }
}
"""

现在我们可以开始设置我们的类。我们将像以前一样创建threading.Thread的另一个子类，并设置构造函数以接受一个输入数组作为参数。我们将用None初始化一个输出变量，就像以前一样：

class KernelLauncherThread(threading.Thread):
    def __init__(self, input_array):
        threading.Thread.__init__(self)
        self.input_array = input_array
        self.output_array = None

现在我们可以编写run函数。我们选择我们的设备，在该设备上创建一个上下文，编译我们的内核，并提取内核函数引用。注意self对象的使用：

def run(self):
    self.dev = drv.Device(0)
    self.context = self.dev.make_context()
    self.ker = SourceModule(kernel_code)
    self.mult_ker = self.ker.get_function('mult_ker')

现在我们将数组复制到 GPU，启动内核，然后将输出复制回主机。然后我们销毁上下文：

self.array_gpu = gpuarray.to_gpu(self.input_array)
self.mult_ker(self.array_gpu, np.int32(array_len), block=(64,1,1), grid=(1,1,1))
self.output_array = self.array_gpu.get()
self.context.pop()

最后，我们设置了 join 函数。这将把output_array返回到主机：

 def join(self):
     threading.Thread.join(self)
     return self.output_array

我们现在已经完成了我们的子类。我们将设置一些空列表来保存我们的随机测试数据、线程对象和线程输出值，与之前类似。然后我们将生成一些随机数组进行处理，并设置一个核发射器线程列表，这些线程将分别操作每个数组：

data = []
gpu_out = []
threads = []
for _ in range(num_arrays):
    data.append(np.random.randn(array_len).astype('float32'))
for k in range(num_arrays):
 threads.append(KernelLauncherThread(data[k]))

现在我们将启动每个线程对象，并通过使用join将其输出提取到gpu_out列表中：

for k in range(num_arrays):
    threads[k].start()

for k in range(num_arrays):
    gpu_out.append(threads[k].join())

最后，我们只需对输出数组进行简单的断言，以确保它们与输入相同：

for k in range(num_arrays):
    assert (np.allclose(gpu_out[k], data[k]))

这个示例可以在存储库中的multi-kernel_multi-thread.py文件中看到。

总结

我们通过学习设备同步和从主机上同步 GPU 操作的重要性开始了本章；这允许依赖操作在进行之前完成。到目前为止，这个概念一直被 PyCUDA 自动处理，然后我们学习了 CUDA 流，它允许在 GPU 上同时执行独立的操作序列而无需在整个 GPU 上进行同步，这可以大大提高性能；然后我们学习了 CUDA 事件，它允许我们在给定流中计时单个 CUDA 内核，并确定流中的特定操作是否已发生。接下来，我们学习了上下文，它类似于主机操作系统中的进程。我们学习了如何在整个 CUDA 上下文中显式同步，然后看到了如何创建和销毁上下文。最后，我们看到了如何在 GPU 上生成多个上下文，以允许主机上的多个线程或进程使用 GPU。

问题

1. 在第一个示例中内核的启动参数中，我们的内核每个都是在 64 个线程上启动的。如果我们将线程数增加到并超过 GPU 中的核心数，这会如何影响原始版本和流版本的性能？

2. 考虑本章开头给出的 CUDA C 示例，它演示了cudaDeviceSynchronize的使用。您认为在不使用流，只使用cudaDeviceSynchronize的情况下，是否可能在多个内核之间获得一定程度的并发？

3. 如果您是 Linux 用户，请修改上一个示例，使其在进程上运行而不是线程。

4. 考虑multi-kernel_events.py程序；我们说内核执行持续时间的标准偏差很低是件好事。如果标准偏差很高会有什么坏处？

5. 在上一个示例中，我们只使用了 10 个主机端线程。列出两个原因，说明为什么我们必须使用相对较少数量的线程或进程来在主机上启动并发的 GPU 操作。

标签：Python,编程,线程,CUDA,gpu,GPU,我们
From： https://www.cnblogs.com/apachecn/p/18140316