1 GPU基础

1.1 GPU与CPU在数据处理中的对比

如今，各个领域——科学、商业、工程、媒体，甚至我们日常的网页浏览——都依赖于不断增长的数据量。单个桌面或服务器 CPU 就能处理整个工作流程的时代正在消逝。因此，随着数据量从兆字节（MB）到千兆字节（GB），再到兆兆字节（TB），我们会发现我们的脚本和应用程序正难以跟上这种增长的步伐。这不仅仅关乎数据量；还关乎数据的复杂性、速度以及快速获得结果的重要性。

我们的任务可能涉及处理医学影像、转换大型数据集进行分析、训练或部署神经网络、运行模拟，或向数百万人传输流媒体视频。即使是多核 CPU，其设计初衷也并非为了处理如此大量的并行重复性工作。您可能已经开始更多地了解 GPU 加速，以及它如何释放近乎神奇的性能，将原本需要数小时才能完成的任务缩短到几分钟甚至几秒钟。

在开始编写代码之前，我们应该了解为什么如此多的应用程序正在从 CPU 转向 GPU。本主题将帮助我们理解主要驱动力、硬件差异以及它们对我们日常数据挑战的实际影响。

1.1.1 CPU 和 GPU 如何处理数据？

计算机系统的中央处理器 (CPU) 通常被称为设备的“大脑”。它功能极其丰富，能够进行复杂的决策，运行整个操作系统，并同时管理数十个软件线程。CPU 内核功能强大且灵活。它配备了复杂的控制逻辑、巨大的缓存和极高的时钟频率。如今，台式机和服务器中的大多数 CPU 都拥有 4 到 64 个内核，每个内核都可以处理独立的指令流。

图片来源：https://developer.nvidia.com/blog/cuda-refresher-reviewing-the-origins-of-gpu-computing

然而，CPU 针对顺序处理进行了优化。您可能会发现我们的代码大量时间都花在循环或单线程例程上，而这些正是 CPU 的优势所在。但是，当处理海量数组时（例如过滤十亿个图像像素或对大型矩阵进行乘法运算），CPU 的设计就暴露出其局限性。虽然速度很快，但每个核心每秒只能处理有限的指令。添加更多核心固然有帮助，但很快就会达到收益递减的临界点。

另一方面，GPU（图形处理单元）采用了一种截然不同的方法。因此，它不再拥有几个“智能”核心，而是拥有数千个设计用于协同工作的简单、轻量级核心。虽然每个 GPU 核心的性能不如单个 CPU 核心强大，但它们擅长同时对海量数据块执行相同的操作。这被称为单指令多数据 (SIMD) 并行。

1.1.2 CPU 扩展的限制

有人可能会问：为什么我们不能直接添加更多 CPU 核心呢？毕竟，如今的服务器 CPU 都配备了数十个核心，而云服务提供商提供的服务器则配备了数百个虚拟 CPU (vCPU)。然而，问题远比简单的“越多越好”要复杂得多。

首先，CPU 的扩展成本很高。每个核心都需要复杂的控制单元和大型低延迟缓存，这会占用大量的硅片空间并消耗大量电力。其次，随着核心数量的增加，我们会遇到阿姆达尔定律——程序的速度取决于其最慢、最连续的部分。除非代码的每个部分都能独立运行，否则总会存在阻碍系统完美扩展的瓶颈。

例如，我们可能使用线程或多处理等库编写了多线程 Python 脚本。随着线程数量的增加，我们很快就会注意到全局锁、线程争用以及臭名昭著的全局解释器锁 (GIL) 限制了这些优势。在线程之间传输数据的过程会产生开销，从而减少可用的内存带宽。实际上，CPU 可以处理数十个并行线程，但无法处理数千个。

1.1.3 GPU 胜过 CPU

这正是 GPU 发挥作用的地方。想象一下，处理一幅 4,000 x 4,000 像素的图像。如果我们要对每个像素应用变换，则需要进行 1600 万次独立计算。这可能包括调整亮度或应用滤镜。即使拥有 16 或 32 个核心，CPU 也必须将这些工作负载划分到少数几个线程上。管理如此多的小任务并在核心之间移动数据的开销将是巨大的。
然而，GPU 正是为这项任务而生。凭借数千个计算核心，GPU 可以将每个像素的操作分配给各自的线程。突然之间，整个图像就可以并行处理，就像每个像素都有自己的处理器一样。数字上的区别显而易见：CPU 需要花费数秒才能完成的工作负载，在 GPU 上只需几分之一秒即可完成。这种架构不仅适用于图形处理。GPU 的对于任何遵循“相同操作，多个数据点”范式的工作负载来说，并行能力都是一大福音。这包括科学计算、神经网络推理、统计模拟和大规模数据分析。

1.1.4 深入了解硬件

为了充分利用 GPU 编程，了解一些硬件内部结构会有所帮助。现代 CPU 可能拥有 8-32 个核心，每个核心的运行频率为 3-4 GHz，并拥有复杂的流水线和数兆字节的缓存。相比之下，典型的 NVIDIA GPU 可能拥有 5,000 到 10,000 个简单的 CUDA 核心，以较低的时钟速度运行，但被分组到流多处理器 (SM) 中。每个 SM 调度和管理 32 个线程组（称为 Warp），硬件可以同时保持多个 Warp 处于“运行”状态。

GPU 通过 PCI Express 连接到系统的其余部分，虽然速度很快，但仍然比片上内存带宽慢得多。您可以在 GPU 上访问多种类型的内存：全局内存（大容量、高延迟、高带宽）、共享内存（小容量、低延迟，位于 SM 内）以及寄存器（超高速，每个线程）。硬件旨在最大化吞吐量；即通过隐藏内存延迟并将计算与数据移动重叠来实现每秒处理的数据总量。

1.1.5 为什么 GPU 编程现在如此重要？

如果我们主要使用 Python 工作，我们就会知道使用 NumPy 或 Pandas 等库操作数据是多么容易。随着数据的增长，我们发现 NumPy 操作会变得缓慢，尤其是在执行大型矩阵乘法、约简或对大型数组进行元素级计算等操作时。GPU 编程曾经只为那些愿意深入研究 C++ 和复杂 CUDA API 的人所准备。值得庆幸的是，像 CuPy 和 PyCUDA 这样的库现在将 GPU 的强大功能直接引入我们的 Python 代码中，让我们能够使用熟悉的语法进行并行计算。您通常可以采用 NumPy 编写的代码，替换掉 import 语句，即可看到显著的加速——所有这些都无需离开 Python 的舒适区。

• 机器学习与深度学习推理
  假设我们正在部署一个神经网络，每秒处理数千张图像进行对象检测或分类。深度神经网络的前向传播包含许多大型矩阵乘法和元素激活。如果没有大规模集群，CPU 就无法在这种规模下进行实时推理。专为并行浮点运算而设计的 GPU 可以在极短的时间内提供结果，支持从云端到边缘的大规模 AI 服务。

• 高性能数据分析
  您可以处理 ETL 管道、处理日志或对数百万条记录进行统计分析。分组、聚合、过滤和直方图计算等任务与 GPU 的优势完美契合。RAPIDS 和 CuDF 等 Python 库利用 GPU 进行数据库式操作，使我们能够加速分析工作负载，使其远远超出 CPU 的能力。

• 实时可视化与渲染
  交互式可视化、3D 渲染和游戏都依赖于 GPU 每秒数十次更新数百万像素的能力。即使在图形处理领域之外，科学家和工程师也使用 GPU 来可视化模拟结果、绘制大型点云或制作分子动力学动画。CPU 根本无法应对这些高要求、高度并行的工作负载。

• 科学计算与模拟
  物理模拟——天气模型、流体动力学、蒙特卡洛模拟——需要更新数百万个粒子或网格点的状态。GPU 在这些高度并行的问题上表现出色，让我们能够运行更精确或更高分辨率的模型，而无需等待数小时才能获得结果。

在 GPU 上运行代码时，并非总能获得 1000 倍的加速。关键在于将工作负载与硬件的优势相匹配。对于足够大的数组，我们可以预期向量加法、矩阵乘法和约简等常见任务的加速可达 10 倍到 100 倍。对于较小的数据，将数据移入和移出 GPU 的开销可能会抵消任何优势。

为了亲眼见证这一点，我们将很快用 Python 运行一个基准测试，比较 CPU（NumPy）和 GPU（CuPy）上简单的数组操作。通过绘制运行时间与数组大小的关系图，我们发现 GPU 在处理大数据时绝对“胜出”，而对于传输开销占主导地位的小任务，CPU 有时速度更快。

1.1.6 GPU 编程有何不同？

GPU 擅长处理具有以下特点的工作负载：
● 多个数据点的操作相同（数据并行）
● 每个计算都是独立的（几乎没有线程间通信）
● 您需要处理大型数据集，而不仅仅是少数几个项目

GPU 不擅长的地方：本质上是顺序的、严重依赖复杂分支的或需要线程间持续通信的工作负载。对于这些工作负载，CPU 仍然至关重要。
随着我们转向 GPU 编程，我们的思维方式发生了转变。你不再考虑“一次一个操作”，而是问，“我能不能用一种方式来写，让所有元素都“一起用吗？”你开始关心内存布局、合并以及如何构建工作以最小化等待。在 Python 中，我们可以使用 CuPy 或 PyCUDA 自动完成大部分工作，但当我们了解如何将算法与硬件对齐时，我们将看到更好的结果。
我们经常会想到一些常见问题和真正的担忧，例如：
● GPU 编程总是有用吗？
● 学习 GPU 编程很难吗？
● 我的代码可以在任何地方运行吗？
● 我可以同时使用 GPU 和 CPU 吗？

所有这些问题的答案都很简单。首先，我们必须考虑数据大小和算法结构。如果我们的数组很小或者我们的算法是高度顺序的，那么 CPU 可能是最佳选择。但是，一旦我们的数据增长，或者我们的计算可以被描述为“对每一行执行相同的操作”，GPU 就会脱颖而出。如果我们坚持使用现代 Python 库，除非我们愿意，否则我们不需要编写低级 CUDA C 代码。我们将获得熟悉的、感觉像 NumPy 的高级 API。本书将向您展示：逐步讲解如何从基本数据结构过渡到启动我们自己的 GPU 内核。

我们还需要兼容的 GPU（通常选择 NVIDIA 的 CUDA 库）和合适的驱动程序。对于许多云提供商来说，只需点击一下即可获得 GPU 实例。即使是消费级笔记本电脑也通常配备支持 CUDA 的硬件。一旦我们的环境搭建完成，工作流程就会像任何 Python 项目一样无缝衔接。我们必须记住，许多工作负载使用 CPU 进行编排、复杂逻辑和数据准备，而 GPU 则负责处理繁重的工作。我们通常会将数据复制到 GPU，运行内核，然后再将结果复制回来。这种混合工作流程在科学计算和机器学习中很常见。

1.1.7 让我们开始吧！

我们将通过对一个简单的操作进行 CPU 和 GPU 的基准测试，迈出实际编码的第一步。您将看到整个工作流程：用 Python 创建数据，在 CPU 上处理数据，然后使用 CuPy 或 PyCUDA 在 GPU 上运行相同的逻辑。您将比较运行时间，验证正确性，并学习如何解读我们的结果。整本书的基调是通过这种实践方法，我们将掌握所需的知识。每个概念都与一个特定的编码任务相关联，我们可以在自己的机器上运行和修改这些任务。随着学习的深入，我们将学习了解哪些操作需要外包，以及如何改进内存移动和内核配置。
您正处于一条实践之路的起点，这条道路将改变我们用 Python 实现高性能计算的方式。我们将学习如何识别数据并行的机会、分析代码，并在 CPU 和 GPU 之间轻松切换，从而使我们能够自信地应对现代数据挑战。如果您希望提高数据管道的速度、创建机器学习服务，或者只是想更快地分析数据，GPU 编程可以将我们的性能提升到一个全新的水平。我们将在实际代码中，观察这些想法如何转化为实际的加速效果，从而继续前进。

1.2 流式多处理器和核心概念

您已经看到了 GPU 加速的潜力，但要充分利用它，我们需要了解硬件内部的运作方式。与专注于高时钟速度和复杂指令逻辑的 CPU 不同，GPU 会组织其资源以实现并行吞吐量。每个现代 GPU 的核心都是一组流式多处理器 (SM：Streaming Multiprocessors)。每个 SM 都像一个强大的引擎，旨在并行执行数千个线程，并采用旨在隐藏延迟和保持数据传输的结构。

我们将深入剖析其结构。当我们观察 GPU 时，我们会发现数十个 SM 并排排列。每个 SM 包含一组简单的 CUDA 核心（基本计算单元）、特殊功能单元、寄存器、共享内存和 Warp 调度器。可以将 SM 视为一个独立的微型处理器，它管理着自己的线程池。SM 的数量因 GPU 型号而异——入门级显卡可能只有几个，而数据中心 GPU 则拥有 80 个或更多。

1.2.1 SM 如何工作？

当我们启动一个内核时，实际上是在告诉 GPU：“这是一个巨大的线程网格；请在所有线程上执行相同的代码。” GPU 会将我们的线程网格拆分成多个块，每个块最多包含数百个线程。然后，每个块会被分配给一个 SM。接下来，事情变得有趣起来：每个 SM 可以同时执行多个块，具体取决于资源可用性（寄存器、共享内存等）。

来源：https://developer.nvidia.com/blog/cuda-refresher-getting-started-with-cuda

在 SM 中，线程进一步分组为 Warp（线程束）——由 32 个线程组成的线程组，它们以同步方式运行。每个周期，SM 的 Warp 调度器都会选择一个或多个活动的 Warp，并将它们发送到 CUDA 核心进行执行。一个 Warp 中的所有线程同时执行同一条指令，但操作不同的数据。这是 SIMD（单指令多数据）编程的基础，也是 GPU 在处理大型数组或图像时如此高效的原因。

您可以想象一个大教室，每一行（warp）都遵循相同的学习计划，但每个学生都在自己的工作表上学习。SM 的 warp 调度程序可以同时管理多个这样的“行”，并在它们之间快速切换，以确保每个人都保持忙碌，即使某些线程需要等待内存或数据。

1.2.2 什么是占用率（occupancy）？

在 GPU 编程中，您会经常听到“占用率”这个术语。占用率衡量我们对 GPU 并行资源的利用效率，具体来说，每个 SM 上有多少个处于活动状态的 warp 相对于最大可能数量。更高的占用率意味着更多的 warp 准备运行，这有助于 SM 隐藏内存延迟。如果一个 warp 因等待数据而停滞，调度程序只需切换到另一个 warp 即可。
您可以通过在每个块中启动更多线程、减少每个线程的寄存器和共享内存使用量，或优化内核以避免瓶颈来提高占用率。然而，更高的占用率并不总是能保证更好的性能——有时，我们会找到一个平衡资源使用率和吞吐量的最佳点。

1.2.3 Warp 调度和性能模式

每个 SM 中的 Warp 调度器就像一个用于并行执行的空中交通管制员。它会根据可用性、就绪情况和资源限制来选择运行哪些 Warp。当一个 Warp 正在等待内存或同步点时，调度器会快速切换入另一个 Warp，从而保持硬件繁忙。这种快速的上下文切换非常轻量，因为所有线程及其数据都驻留在 SM 上。
在处理高延迟操作（例如从全局内存读取）时，您就能体会到这种架构的真正威力。Warp 调度器会同时处理数十个 Warp、重叠计算和数据传输，因此 SM 永远不会空闲。如果我们构建代码以提供足够的并行工作，调度器就会为您处理延迟隐藏。

1.2.4 运行简单的微基准测试

我们将进行实践操作。您可以使用小型 PyCUDA 或 CuPy 内核探索 SM 和 Warp 调度概念。假设我们想要测量 GPU 可以并行处理的线程数，以及执行时间如何随不同的块大小而变化。这不仅仅是一个学术练习——它是调优代码以实现高性能的第一步。
首先，设置一个包含数百万个元素的数组。然后，编写一个内核，对每个元素执行一个简单的操作（例如，加一个常量）。启动内核，并使用不同的块大小（例如，每个块 32、64、128、256、512 和 1024 个线程），并测量每个块的执行时间。我们可能会观察到某些配置的运行速度比其他配置快得多，这取决于它们与我们 GPU 的 SM 架构的“匹配程度”。通过使用 Python 查询（CuPy 或 PyCUDA 提供设备属性）打印出 SM 的数量和每个 SM 的最大线程数，我们可以将启动配置与硬件的最佳性能相匹配。尝试用较少的线程运行相同的操作——你的 GPU 将无法充分利用。如果每个块的线程过多，我们可能会耗尽寄存器或共享内存的减少，导致占用率下降并损害性能。
在我们进行这些微基准测试的过程中，一些高性能模式变得清晰起来。首先，我们希望在不过度消耗资源的情况下最大化每个 SM 的线程数。其次，我们倾向于启动足够多的块，以便每个 SM 保持繁忙，尤其是在更大的 GPU 上。第三，我们注意到内存合并的重要性：构建数据访问，以便 Warp 中的线程从连续的内存地址读取数据，从而提高吞吐量。当我们将这些模式付诸实践时，即使是用于向量加法或直方图计算等任务的简单内核，其速度也会比简单的实现快几个数量级。
如果我们理解了 SM、Warp、占用率和调度程序，我们就可以自信地尝试启动参数和内核资源使用情况。我们不仅仅是编写代码，而是根据硬件的节奏调整程序。我们编写的每个内核都提供了一个实时观察决策效果的机会。

1.3 线程块、网格和索引

我相信您现在已经听说过 CUDA。这项技术使我们能够在自己的代码中充分利用 NVIDIA GPU 的强大功能。CUDA（统一计算设备架构：Compute Unified Device Architecture）提供了在 GPU 上启动数千个轻量级线程的框架。使用 Python，我们可以使用 PyCUDA 或 CuPy 等库来编写编译和运行 CUDA 内核的代码，同时又能保持我们惯用的 Python 环境。本书将使用 Linux 作为默认设置，因为它对 CUDA 开发提供了强大的支持。

（来源：https://developer.nvidia.com/blog/cuda-refresher-cuda-programming-model）

1.3.1 线程层次结构

在 CUDA 中，并行工作使用两级层次结构来组织：网格和块。每次启动内核时，我们都会指定要运行的块数以及每个块中要放置的线程数。网格中的每个线程执行相同的代码，但处理数据的不同部分。将我们的数据想象成一个大型电子表格。我们的想法是为电子表格的每个区域（或块）分配一个团队，并为该区域中的每个单元格（或线程）分配一个工作线程。CUDA 的架构使我们能够处理海量数据集：一维、二维甚至三维数组，方法是将每个数据元素映射到网格中其自己的线程。

1.3.2 映射一维数组

为此，我们需要一个支持 CUDA 的 Python 环境。如果这是我们第一次尝试，请继续在虚拟环境中安装 PyCUDA 或 CuPy，并确认我们的 GPU 可见。如果我们想要检查，只需导入我们的库并打印设备属性即可。例如，使用 CuPy：

>>> import cupy as cp
>>> print(cp.cuda.runtime.getDeviceProperties(0))
{'name': b'NVIDIA GeForce RTX 4090', 'totalGlobalMem': 25241190400, 'sharedMemPerBlock': 49152, 'regsPerBlock': 65536, 'warpSize': 32, 'maxThreadsPerBlock': 1024, 'maxThreadsDim': (1024, 1024, 64), 'maxGridSize': (2147483647, 65535, 65535), 'clockRate': 2520000, 'totalConstMem': 65536, 'major': 8, 'minor': 9, 'textureAlignment': 512, 'texturePitchAlignment': 32, 'multiProcessorCount': 128, 'kernelExecTimeoutEnabled': 0, 'integrated': 0, 'canMapHostMemory': 1, 'computeMode': 0, 'maxTexture1D': 131072, 'maxTexture2D': (131072, 65536), 'maxTexture3D': (16384, 16384, 16384), 'concurrentKernels': 1, 'ECCEnabled': 0, 'pciBusID': 172, 'pciDeviceID': 0, 'pciDomainID': 0, 'tccDriver': 0, 'memoryClockRate': 10501000, 'memoryBusWidth': 384, 'l2CacheSize': 75497472, 'maxThreadsPerMultiProcessor': 1536, 'isMultiGpuBoard': 0, 'cooperativeLaunch': 1, 'cooperativeMultiDeviceLaunch': 1, 'deviceOverlap': 1, 'maxTexture1DMipmap': 32768, 'maxTexture1DLinear': 268435456, 'maxTexture1DLayered': (32768, 2048), 'maxTexture2DMipmap': (32768, 32768), 'maxTexture2DLinear': (131072, 65000, 2097120), 'maxTexture2DLayered': (32768, 32768, 2048), 'maxTexture2DGather': (32768, 32768), 'maxTexture3DAlt': (8192, 8192, 32768), 'maxTextureCubemap': 32768, 'maxTextureCubemapLayered': (32768, 2046), 'maxSurface1D': 32768, 'maxSurface1DLayered': (32768, 2048), 'maxSurface2D': (131072, 65536), 'maxSurface2DLayered': (32768, 32768, 2048), 'maxSurface3D': (16384, 16384, 16384), 'maxSurfaceCubemap': 32768, 'maxSurfaceCubemapLayered': (32768, 2046), 'surfaceAlignment': 512, 'asyncEngineCount': 2, 'unifiedAddressing': 1, 'streamPrioritiesSupported': 1, 'globalL1CacheSupported': 1, 'localL1CacheSupported': 1, 'sharedMemPerMultiprocessor': 102400, 'regsPerMultiprocessor': 65536, 'managedMemory': 1, 'multiGpuBoardGroupID': 0, 'hostNativeAtomicSupported': 0, 'singleToDoublePrecisionPerfRatio': 64, 'pageableMemoryAccess': 1, 'concurrentManagedAccess': 1, 'computePreemptionSupported': 1, 'canUseHostPointerForRegisteredMem': 1, 'sharedMemPerBlockOptin': 101376, 'pageableMemoryAccessUsesHostPageTables': 0, 'directManagedMemAccessFromHost': 0, 'uuid': b'\xea\xce)\xa0\xbe\xf8\x9d\xfe\x06w\x96\xdb\x89\x84,\xc4', 'luid': b'', 'luidDeviceNodeMask': 0, 'persistingL2CacheMaxSize': 51904512, 'maxBlocksPerMultiProcessor': 24, 'accessPolicyMaxWindowSize': 134213632, 'reservedSharedMemPerBlock': 1024}

之后，我们将从一个简单的一维数组开始。假设我们有一个包含一百万个元素的数组，并且我们希望每个线程对其分配的元素进行平方。启动内核时，我们会确定每个块的线程数（例如 256），以及所需的块数（num_elements //threads_per_block + 1）。
在 CUDA 内核中，每个线程都会计算出其唯一的索引：

int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < num_elements) {
   output[idx] = input[idx] * input[idx];
}

blockIdx.x * blockDim.x + threadIdx.x 的计算为每个线程在整个网格中提供了唯一的全局索引。 if 语句是我们的边界检查——某些线程的索引可能超出了数组的末尾，因此我们阻止它们越界写入。
在 Python 端，运行内核后，我们将输出数组从设备复制回主机并检查结果：

import cupy as cp
num_elements = 1_000_000
input_array = cp.arange(num_elements, dtype=cp.float32)
output_array = cp.empty_like(input_array)
threads_per_block = 256
blocks_per_grid = (num_elements + threads_per_block - 1) // threads_per_block
kernel_code = '''
extern "C" __global__
void square(const float *input, float *output, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
   if (idx < n) {
       output[idx] = input[idx] * input[idx];
   }
}
'''
module = cp.RawModule(code=kernel_code)
square_kernel = module.get_function('square')
square_kernel((blocks_per_grid,), (threads_per_block,),
             (input_array, output_array, num_elements))
# Copy results to host and verify
result = cp.asnumpy(output_array)
assert (result == (cp.asnumpy(input_array) ** 2)).all()
print("1D kernel executed successfully!")

我们刚刚将一维数组映射到 CUDA 的网格和块系统上。

1.3.3 映射二维数组

现在，假设我们要处理一个二维数组，例如图像或矩阵。您需要将行和列映射到两个网格维度。 CUDA 支持多维块和网格，允许每个线程处理单个像素或矩阵条目。
现在，在我们的内核中，我们需要计算行和列的索引：

int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < height && col < width) {
    output[row * width + col] = input[row * width + col] + 1.0f;
}

在 Python 中，我们可以调整启动配置：

import cupy as cp

height, width = 1024, 1024
input_array = cp.random.rand(height, width).astype(cp.float32)
output_array = cp.empty_like(input_array)
block = (16, 16)
grid = ((width + block[0] - 1) // block[0],
       (height + block[1] - 1) // block[1])
kernel_code = '''
extern "C" __global__
void increment(const float *input, float *output, int width, int height) {
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    if (row < height && col < width) {
       output[row * width + col] = input[row * width + col] + 1.0f;
   }
}
'''
module = cp.RawModule(code=kernel_code)
increment_kernel = module.get_function('increment')
increment_kernel(grid, block, (input_array, output_array, width, height))
result = cp.asnumpy(output_array)
assert (result == cp.asnumpy(input_array) + 1).all()
print("2D kernel executed and verified successfully!")

之后，每个线程计算其行和列。边界检查确保如果网格大小与数组尺寸不完全匹配，我们不会越过边界。

您可以将此模式扩展到三维数据，例如体积图像或模拟网格。 CUDA 允许我们指定三维块和网格，每个线程计算 (z, y, x) 索引来处理其自身的体素或数据点。关键始终如一：将逻辑数据位置（无论是一维、二维还是三维）映射到 CUDA 的线程和块索引，检查边界，并确保每个线程仅在有效数据范围内工作。

了解了这一点后，我们就可以开始对任何数据形状使用自定义内核了。将 Python 数组映射到 CUDA 网格和块、处理边界以及验证输出的基础知识已经包含在我们的工具包中。这些模式几乎会体现在我们构建的每个 GPU 应用程序中，为实际的高性能 Python 代码奠定基础。

1.4 主机与设备交互概述

您将亲身体验我们的 Python 代码如何与 GPU 通信。此过程包含几个核心步骤：内存分配、数据传输、启动内核以及了解所涉及的不同内存空间。您已经安装了 CuPy，因此可以开始将理论付诸实践。我们将通过编写一个涵盖所有基本知识的简洁而完整的脚本来构建我们的直觉。

1.4.1 GPU 上的内存分配

在 CuPy 中，当我们分配数组时，实际上是直接在 GPU 上预留内存，而不是在计算机的主内存上。这是一个巨大的思维转变。如果我们使用 CuPy，我们创建的数组将保存在设备上。
这里有一个简单的方法：

import cupy as cp
# Allocate an array of 10,000 floats on the GPU
gpu_array = cp.zeros(10_000, dtype=cp.float32)
print("Array allocated on GPU:", gpu_array)

这一行代码为我们分配了一块 GPU 内存。与 NumPy 数组不同，这些数组不能直接用于仅使用 CPU 的函数。

1.4.2 主机到设备，以及从设备到主机

主要的 Python 脚本在 CPU（“主机”）上运行，但繁重的工作在 GPU（“设备”）上进行。要使用 GPU 加速，我们首先需要将数据从主机移动到设备。
因此，我们首先创建一个 NumPy 数组，然后将其传输到 GPU。

import numpy as np
import cupy as cp

# Create data on the host
host_data = np.arange(10_000, dtype=np.float32)
# Transfer data to the GPU (device)
device_data = cp.asarray(host_data)
print("Data transferred to GPU.")

# Bring the data back from device to host
result_on_host = cp.asnumpy(device_data)
print("Data transferred back to host. First five elements:", result_on_host[:5])

您经常会发现自己需要来回移动数据，尤其是在我们需要验证输出、保存结果或与其他需要 NumPy 数组的库交互时。

1.4.3 启动简单内核

CuPy 的一大亮点在于它允许我们直接在 Python 中编写自定义 CUDA 内核。假设有一个数组，我们想将每个元素乘以 2。你可以使用 CuPy 的 RawKernel 接口来编写这个内核。

import cupy as cp

# 1. 内核代码定义
kernel_code = r'''
extern "C" __global__
void multiply_by_two(float* data, int n) {
    int idx = blockDim.x * blockIdx.x + threadIdx.x;
    if (idx < n) {
        data[idx] *= 2.0f;
    }
}
'''

# 2. 编译内核
module = cp.RawModule(code=kernel_code)
multiply_by_two = module.get_function('multiply_by_two')

# 3. 设置参数和数据
n = 10_000
threads_per_block = 256
blocks_per_grid = (n + threads_per_block - 1) // threads_per_block

# 4. 创建并初始化GPU数组
gpu_array = cp.random.rand(n, dtype=cp.float32)

print("Original array (first 5 elements):", gpu_array[:5].get()) # .get() 将数据从GPU传回CPU

# 5. 启动内核
multiply_by_two((blocks_per_grid,), (threads_per_block,), (gpu_array, n))

# 6. 等待内核执行完毕并验证结果
cp.cuda.stream.get_current_stream().synchronize()
print("Kernel launched and synchronized.")
print("Modified array (first 5 elements):", gpu_array[:5].get()) # 再次将数据传回CPU，查看结果

1.4.4 理解 GPU 内存空间

所以到目前为止，您已经使用了两个主要的“内存空间”：
● 主机内存（您的主系统 RAM），NumPy 数组存放于此，Python 脚本也从这里开始。
● 另一个是设备内存（GPU），CuPy 数组在此分配，内核也在此读写数据。
GPU 编程的魅力在于了解何时移动数据以及每个操作在何处进行。CuPy 数组上的大多数计算都完全在设备上进行，从而保持快速运行并避免不必要的传输。
此外，高级内核还可以访问：
● 共享内存，这是一个在块中的线程之间共享的小型、快速的内存区域，对于需要协作的算法非常有用。
● 寄存器，每个线程私有的超高速存储空间。
● 常量内存和纹理内存，它们是用于只读或空间相干数据的专用空间，我们将在后面的章节中探讨。
有了这些，我们现在有了一个完整的工作流程，如下所示：
● 在 GPU 上分配数组
● 在主机和设备之间移动数据
● 编写并启动简单的内核
● 观察和管理 GPU 内存使用情况
每个支持 GPU 加速的项目都会遵循这些步骤。随着我们水平的提升，我们希望尽可能多地在设备上进行计算，并且仅将主机-设备传输用于输入和最终输出。掌握这些知识后，您就可以将我们的 GPU 工作流程提升到一个新的水平。

参考资料

• 软件测试精品书籍文档下载持续更新 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
• python八字排盘 https://github.com/china-testing/bazi
• 联系方式：钉ding或V信: pythontesting

1.5 简单向量加法内核

我们即将构建第一个自定义 GPU 内核，用于执行向量加法等基本任务。每个新加入团队的 GPU 程序员都必须完成这项练习才能加入团队。我们将学习如何使用 CuPy 和 PyCUDA 创建、启动和验证一个用于将两个向量相加的内核。我们还会将 GPU 上的结果与 CPU 上的结果进行比较，以确保一切正常。

1.5.1 安装 PyCUDA

我们已经准备好了 CuPy，它非常适合大多数高级 GPU 阵列工作。如果您想要更多控制权并直接访问 CUDA C 语言，可以使用 PyCUDA。它允许我们编写原始 CUDA 内核，并赋予我们对内存管理、内核启动和设备上下文的强大控制权，所有这些都可以通过 Python 实现。
首先，让我们在 Python 环境中安装 PyCUDA。

pip install pycuda
``

安装完成后，我们可以导入 PyCUDA 模块并立即开始使用其功能。

### 1.5.2 主机和设备上的向量

- CuPy

我们首先使用 NumPy（CPU）和 CuPy（GPU）设置输入数据：
```python
import numpy as np
import cupy as cp

# 设置向量大小
N = 1_000_000
# 在主机（CPU）上创建两个随机输入向量
a_host = np.random.rand(N).astype(np.float32)
b_host = np.random.rand(N).astype(np.float32)
# 计算CPU上的结果，用于后续验证
c_cpu = a_host + b_host
# 使用 CuPy 将输入向量传输到 GPU
a_gpu = cp.asarray(a_host)
b_gpu = cp.asarray(b_host)

# CUDA C++ 内核代码
kernel_code = r'''
// 声明为全局函数，可以在GPU上执行
extern "C" __global__
// 定义向量相加函数，接受三个数组指针和数组大小
void vector_add(const float* a, const float* b, float* c, int n) {
    // 计算当前线程在整个网格中的全局唯一索引
    int idx = blockDim.x * blockIdx.x + threadIdx.x;
    // 边界检查，确保线程索引不超过数组大小
    if (idx < n) {
        // 执行向量加法操作
        c[idx] = a[idx] + b[idx];
    }
}
'''

# 在GPU上创建用于存放结果的空数组
c_gpu = cp.empty_like(a_gpu)
# 定义每个块的线程数
threads_per_block = 256
# 计算需要的块数，确保所有元素都能被处理到
blocks_per_grid = (N + threads_per_block - 1) // threads_per_block
# 编译CUDA内核代码
module = cp.RawModule(code=kernel_code)
# 从编译模块中获取内核函数
vector_add = module.get_function('vector_add')
# 启动内核，并传入网格、块的维度以及函数参数
vector_add((blocks_per_grid,), (threads_per_block,), (a_gpu, b_gpu, c_gpu, N))

# 将GPU上的结果数组传输回主机（CPU）
c_result = cp.asnumpy(c_gpu)
# 比较GPU和CPU的结果，检查它们是否几乎相等
print("Are GPU and CPU results equal?", np.allclose(c_result, c_cpu))

• 使用 CuPy 运行内核

import numpy as np
import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule

# 设置向量大小
N = 1_000_000
# 在主机（CPU）上创建两个随机输入向量
a_host = np.random.rand(N).astype(np.float32)
b_host = np.random.rand(N).astype(np.float32)

# 计算CPU上的结果，用于后续验证
c_cpu = a_host + b_host

# 在GPU上分配内存
# pycuda.mem_alloc 分配的内存需要手动管理
a_gpu = cuda.mem_alloc(a_host.nbytes)
b_gpu = cuda.mem_alloc(b_host.nbytes)
c_gpu = cuda.mem_alloc(a_host.nbytes)

# 将数据从主机（CPU）传输到设备（GPU）
cuda.memcpy_htod(a_gpu, a_host)
cuda.memcpy_htod(b_gpu, b_host)

# CUDA C++ 内核代码
kernel_code = """
__global__ void vector_add(const float* a, const float* b, float* c, int n) {
    int idx = blockDim.x * blockIdx.x + threadIdx.x;
    if (idx < n) {
        c[idx] = a[idx] + b[idx];
    }
}
"""

# 编译内核代码
# pycuda.compiler.SourceModule 负责编译
mod = SourceModule(kernel_code)
vector_add = mod.get_function("vector_add")

# 定义每个块的线程数和需要的块数
threads_per_block = 256
blocks_per_grid = (N + threads_per_block - 1) // threads_per_block

# 启动内核
# 注意参数的传递方式：(grid_dim, block_dim, shared_mem, stream, *args)
# 这里我们没有使用共享内存和流，所以为0和None
vector_add(a_gpu, b_gpu, c_gpu, np.int32(N),
           grid=(blocks_per_grid, 1, 1),
           block=(threads_per_block, 1, 1))

# 从设备（GPU）将结果传输回主机（CPU）
c_result = np.empty_like(a_host)
cuda.memcpy_dtoh(c_result, c_gpu)

# 比较GPU和CPU的结果，检查它们是否几乎相等
print("Are GPU and CPU results equal?", np.allclose(c_result, c_cpu))

# 释放GPU内存
a_gpu.free()
b_gpu.free()
c_gpu.free()

另外一种更简洁的实现：

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

# --- 1. 定义缺失的变量和代码 ---
kernel_code = """
__global__ void vector_add(const float* a, const float* b, float* c, int n) {
    int idx = blockDim.x * blockIdx.x + threadIdx.x;
    if (idx < n) {
        c[idx] = a[idx] + b[idx];
    }
}
"""

# 设置向量大小
N = 1_000_000
# 在主机（CPU）上创建两个随机输入向量
a_host = np.random.rand(N).astype(np.float32)
b_host = np.random.rand(N).astype(np.float32)

# 计算CPU结果用于验证
c_cpu = a_host + b_host

threads_per_block = 256
blocks_per_grid = (N + threads_per_block - 1) // threads_per_block
# ------------------------------------

# Transfer input data to GPU
a_gpu_py = gpuarray.to_gpu(a_host)
b_gpu_py = gpuarray.to_gpu(b_host)
c_gpu_py = gpuarray.empty_like(a_gpu_py)

mod = SourceModule(kernel_code)
vector_add_func = mod.get_function("vector_add")

# --- 2. 修正内核启动参数 ---
vector_add_func(
    a_gpu_py, b_gpu_py, c_gpu_py, np.int32(N),
    block=(threads_per_block, 1, 1), grid=(blocks_per_grid, 1, 1)
)

# 获取GPU结果
c_result_py = c_gpu_py.get()
print("PyCUDA GPU result matches CPU?", np.allclose(c_result_py, c_cpu))

现在，我们应该再次看到 True。两个库都给出了相同的结果，证实我们的内核可以正常工作！这样，我们就可以在 CuPy 和 PyCUDA 中编写、编译并启动一个自定义 CUDA 内核了。我们以闪电般的速度执行了元素向量加法，并验证了结果与我们的 CPU 计算结果一致。这是所有后续 GPU 编程的基础，因为每个复杂的操作都是基于这些模式构建的。

1.6 小结

本章讲解了所有必要的基础模块，为实际的 GPU 编程奠定了坚实的基础。我们首先了解了现实世界对高性能数据处理的需求，了解了传统 CPU 在现代工作负载下经常遇到的瓶颈，并理解了拥有数千个轻量级核心的 GPU 为何能够加速机器学习推理和大规模分析等任务。我们研究了流多处理器的设计，了解了 SM、warp 和占用率指标如何让您高效地并行运行。我们还学习了 GPU 隐藏内存延迟并保持所有任务繁忙的技巧。

然后，我们深入研究了将多维数组映射到 CUDA 网格和块。我们掌握了线程层次结构、索引计算和稳健边界处理的概念。我们用 Python 完成了整个 GPU 数据处理工作流程，从使用 CuPy 设置设备内存，到在主机和设备之间移动数据、启动自定义内核，再到将结果复制回来以验证所有内容。我们尝试使用 CuPy 和 PyCUDA 编写并运行了我们的第一个自定义 CUDA C 内核，并经历了从在 CPU 上准备数据到验证 GPU 输出的整个过程。凭借这些技能和习惯，我们将在本书的其余部分深入探索 GPU 编程的更高级、更具创造性的用法。