GPU 编程实战——GPU 基础入门

GPU 与 CPU 在数据处理上的比较

在这样一个时代，无论是科学、商业、工程、媒体，甚至我们日常的网页浏览，都依赖于日益增长的数据量。单个台式机或服务器的 CPU 足以应对整个工作流程的日子正在消逝。随着我们从兆字节迈向千兆字节，再到万亿字节的数据规模，你会发现脚本和应用程序确实难以跟上。这不仅仅是数据体量的问题，更在于其复杂度、处理速度，以及即时获取结果的重要性。

我们的任务可能包括处理医疗图像、转换大型数据集以进行分析、训练或部署神经网络、运行仿真，甚至向数百万用户推送视频流。即便拥有多核设计，CPU 也并非为大规模并行、重复性的计算而生。你或许已经听说过 GPU 加速技术，它能带来近乎“魔法”般的性能，让需要数小时完成的作业缩短到数分钟，甚至数秒。

在动手编写代码之前，我们应该弄清楚，为什么如今这么多应用都由 CPU 转向 GPU。本节将帮助我们理解这一趋势背后的主要驱动力、硬件架构的差异，以及这种转变对我们日常数据处理所带来的实际影响。

CPU 与 GPU 如何处理数据？

计算机系统的中央处理器（CPU）常被称为“计算机的大脑”。它功能极其通用，能够做出复杂决策、运行完整操作系统，并同时管理数十个软件线程。CPU 核心强大且灵活，配备了复杂的控制逻辑、大容量缓存和极高的主频。目前台式机和服务器里常见的 CPU 通常拥有 4 到 64 个核心，每个核心都能并行处理独立的指令流。

CPU 优化的是顺序处理。你可能会发现，我们的代码大部分时间都花在循环或单线程例程中——这些正是 CPU 最擅长的场景。但当需要处理海量数组（比如过滤十亿像素的图像，或对超大矩阵进行乘法运算）时，CPU 的设计局限就会显现。虽然每个核心速度很快，但每秒钟只能执行有限数量的指令。增加核心数量虽然能提升性能，却很快会遇到收益递减的瓶颈。

GPU（图形处理器）则采取了截然不同的策略：与少量“智能”核心不同，GPU 拥有成千上万个简单而轻量的核心，它们协同工作。单个 GPU 核心虽然不如 CPU 核心强大，但擅长在巨大数据块上同时执行同一操作，这就是所谓的“单指令多数据”（SIMD）并行。

CPU 扩展的天花板

有人会问：那我们为什么不继续增加 CPU 核心呢？毕竟，服务器 CPU 已经内置数十个核心，云服务商还提供数百 vCPU 的机器。但问题并非“核心越多越好”那么简单。

1. 成本高昂：每个 CPU 核心都需要复杂的控制单元和大容量、低延迟的缓存，这既占用硅片面积，也消耗大量功率。
2. 阿姆达尔定律：程序的执行速度受最慢的、最具顺序性的部分制约。除非代码的每一部分都能完全并行，否则总会存在瓶颈，无法实现完美扩展。

例如，即便在 Python 中使用 threading 或 multiprocessing 编写了多线程脚本，随着线程数增加，你很快便会遇到全局锁、线程争用以及著名的全局解释器锁（GIL）带来的限制。线程间传输数据的开销也会吞噬可用带宽。实际上，CPU 可以处理数十个并行线程，却远不能达到数千线程的规模。

GPU 如何胜出

这正是 GPU 大显身手的地方。想象一幅 4000×4000 像素的图像，我们要对每个像素都执行某种变换——比如调整亮度或应用滤镜，这就是 1600 万次独立计算。即便有 16 或 32 个核心，CPU 也只能将工作负载分割给少量线程，管理众多小任务和在核心间调度数据的开销非常可观。

而 GPU 正是为此而生。它将每个像素操作分配给不同的线程，同时并行处理整个图像，仿佛每个像素都拥有自己的处理器。结果非常明显：本需数秒或数十秒完成的任务，在 GPU 上只需短短几百毫秒。GPU 的并行架构不仅仅适用于图形渲染，对任何“相同操作，多数据点”范式的工作负载都大有裨益，如科学计算、神经网络推理、统计模拟和大规模数据分析。

硬件深入剖析

要充分利用 GPU 编程，了解硬件内部结构很重要。现代 CPU 可能拥有 8–32 个核心，主频在 3–4 GHz 之间，拥有复杂的流水线和数 MB 缓存。相比之下，一块典型的 NVIDIA GPU 可能具备 5000–10000 个简单的 CUDA 核心，主频较低，但被划分为多个流式多处理器（SM）。每个 SM 会调度并管理一组 32 线程的“warp”，并能同时维持多个 warp“飞行”状态。

GPU 与系统其余部分通过 PCI Express 相连，这虽然速度很快，却仍比不上片上（on-chip）内存带宽。在 GPU 上，你会接触到多种内存类型：

• 全局内存：容量大、延迟高、带宽高
• 共享内存：SM 内部的小容量、低延迟内存
• 寄存器：每线程专属的超高速存储

硬件设计追求吞吐量最大化，即单位时间内处理的数据总量，通过隐藏内存延迟并在计算与数据传输间重叠来实现。

为什么现在要学 GPU 编程？

如果我们主要使用 Python，就知道借助 NumPy、pandas 操作数据有多方便，但当数据规模扩大时，NumPy 的性能瓶颈一览无余，尤其是大矩阵乘法、归约运算或大数组的逐元素计算。过去，GPU 编程往往意味着深入 C++ 和复杂的 CUDA API，但如今有了 CuPy、PyCUDA 等库，我们只需在 Python 中稍作调整，就能调用 GPU 加速，轻松获得显著提速。常见做法是将 import numpy as np 换成 import cupy as np，代码几乎不变，却能立刻体验到性能提升。

机器学习与深度学习推理

假设我们部署一个每秒处理数千张图像的神经网络，用于物体检测或分类。深度网络的前向推理充满大规模矩阵运算和逐元素激活函数。单靠 CPU，想实现实时推理需要庞大集群；而 GPU 专为并行浮点运算而设计，可在极短时间内完成，从云端到边缘都适用。

高性能数据分析

无论是 ETL 管道、日志处理，还是对数百万条记录的统计分析，分组、聚合、过滤或直方图计算都非常契合 GPU 的优势。RAPIDS、cuDF 等 Python 库让我们像操作 pandas 一样便捷地调用 GPU，大幅超越 CPU 性能。

实时可视化与渲染

交互式可视化、3D 渲染、游戏都依赖 GPU 每秒更新数百万像素的能力。但即便在非图形领域，科学家和工程师也用 GPU 可视化模拟结果、绘制大规模点云、动画分子动力学。CPU 难以承受如此高并行度的计算需求。

科学计算与仿真

天气模型、流体力学、蒙特卡洛模拟等物理仿真需更新数百万个粒子或网格点，属于“尴尬地并行”问题。GPU 在此类任务上表现尤为出色，让我们无需长时间等待即可运行更高精度或更高分辨率的模型。

并非所有场景都 1000× 提速

GPU 的优势在于匹配合适的工作负载。对于足够大的数组，常见向量加法、矩阵乘法和归约等操作通常能获得 10×–100× 的加速；而对于小规模数据，来回传输的开销可能抵消任何收益。

GPU 编程的特性

GPU 擅长的工作负载特点：

• 多数据点执行相同操作（数据并行）
• 各计算彼此独立（几乎无需线程间通信）
• 需处理大规模数据集

而对于高度顺序、分支复杂或频繁线程通信的任务，CPU 依然不可或缺。

思维模式的转变

进入 GPU 编程后，我们不再思考“每次一项操作”，而是问：“如何让所有元素同时执行？”我们开始关注内存布局、内存访问共同化（coalescing），以及如何安排任务以最小化等待。在 Python 中，CuPy、PyCUDA 帮我们自动完成大部分细节，但理解底层硬件能带来更佳性能。

常见疑问

• GPU 编程总是有用吗？ 要看数据大小与算法结构；若数组小或算法高度顺序化，CPU 可能更优。
• GPU 编程难学吗？ 借助高层库，不用写底层 CUDA C 代码也能上手。
• 代码能否跨平台运行？ 只要使用兼容库，通常可在支持 CUDA 的 NVIDIA GPU 上无缝运行。
• 能否同时用 GPU 和 CPU？ 是的，常见的混合工作流是：CPU 负责编排和数据准备，GPU 负责重负载计算。

本书将手把手带你从基本数据结构，到编写并优化 GPU 内核，配合适当的硬件与驱动，一步步在 Python 中实现高性能并行计算。你将学习识别数据并行机会、使用分析工具剖析性能，并在 CPU 与 GPU 间自如切换，应对现代数据挑战。接下来，我们会通过实际代码，对比简单数组操作在 NumPy（CPU）与 CuPy（GPU）上的运行时，绘制性能随数组规模变化的曲线，直观感受 GPU 对大数据的压倒性优势。让我们立刻动手，开启高性能计算之旅！

流式多处理器与核心概念

你已经看到了 GPU 加速的潜力，但要充分利用它，我们需要了解硬件内部的运作机制。与注重高主频和复杂指令逻辑的 CPU 不同，GPU 将资源组织成并行吞吐的架构。每一代现代 GPU 的核心，都是一组流式多处理器（Streaming Multiprocessor，简称 SM）。每个 SM 都是一个强大的单元，能够并行执行数千个线程，其内部结构专门用于隐藏延迟并保持数据流动。

下面我们分解其“解剖结构”。从外部看，一块 GPU 上并排排列着数十个 SM。每个 SM 包含若干简单的 CUDA 核心（基本计算单元）、特殊功能单元、寄存器、共享内存和 warp 调度器。可以把 SM 想象成一个小型的独立处理器，管理着自己的线程池。不同型号的 GPU 拥有的 SM 数量不一：入门级显卡可能只有几个 SM，而数据中心级 GPU 则多达八十个或更多。

SM 如何工作？

当我们启动一个 kernel 时，相当于告诉 GPU：“这是一个庞大的线程网格，请在所有线程上执行相同的代码。”GPU 会将这个线程网格拆分成若干个 thread block，每个 block 含有最多数百个线程。然后，GPU 将这些 block 分配给各个 SM。值得注意的是，每个 SM 可以同时运行多个 block，具体数量取决于可用资源（如寄存器、共享内存等）。这样一来，GPU 就能高效地隐藏内存或执行延迟，在保持持续计算的同时最大化并行吞吐。

在 SM 内部的线程分组与 SIMD 并行
在 SM（流式多处理器）内部，线程会进一步被分成若干个 warp——每个 warp 包含 32 个同时执行的线程。每个时钟周期，SM 内的 warp 调度器会选择一个或多个活跃的 warp，并将它们派发到 CUDA 核心上执行。一个 warp 中的所有线程会同时执行相同的指令，但操作各自不同的数据片段。这正是 SIMD（单指令多数据）并行的基础，也因此 GPU 在处理大规模数组或图像时极为高效。

你可以将其想象成一个大型教室，每排同学（一个 warp）都按同一教学进度学习，但每位同学完成自己桌面上的习题。SM 的 warp 调度器可以同时管理多排同学，快速在它们之间切换，即便某些线程因等待内存或数据而暂停，调度器也能马上切换到下一个 warp，始终保持硬件忙碌。

什么是占用率（Occupancy）？

“占用率”是 GPU 编程中常见的术语，用来衡量我们对并行资源的利用效率，具体地说，就是每个 SM 上活跃 warp 的数量与其最大支持 warp 数量的比值。更高的占用率意味着更多的 warp 随时准备运行，有助于 SM 隐藏内存延迟——当一个 warp 因等待数据而停滞时，调度器只需切换到另一个活跃 warp 即可。

提高占用率的方法包括：

• 为每个 block 启动更多线程
• 减少每线程使用的寄存器和共享内存
• 优化 kernel 以避免资源瓶颈

但需要注意：更高的占用率并不总是带来更佳性能；往往需要在资源使用与吞吐量之间找到平衡的“甜蜜点”。

Warp 调度与性能模式

每个 SM 的 warp 调度器就像并行执行的“空中交通管制员”，它根据 warp 的可用性、就绪状态和资源约束来选择执行对象。当一个 warp 在等待内存或同步点时，调度器会迅速切换到另一个就绪 warp，确保硬件持续工作。这种快速切换开销极低，因为所有线程及其数据都常驻于 SM 内。

这种架构的强大之处在于，它能在处理高延迟操作（如全局内存读取）时，通过重叠计算与数据传输最大化吞吐量。只要我们将代码结构化，暴露出足够多的并行任务，调度器便会自动替我们“隐藏”延迟。

运行简单的微基准测试

接下来，我们动手实践。可以使用 PyCUDA 或 CuPy 编写一个小型 kernel，测试不同 block 大小下 GPU 同时处理线程的能力及其执行时长变化。步骤示例：

1. 在 Python 中创建包含数百万元素的数组。
2. 编写一个 kernel，对每个元素执行简单运算（如加常数）。
3. 分别以 32、64、128、256、512、1024 线程/block 的配置多次启动 kernel，并测量执行时间。
4. 使用 CuPy/PyCUDA 提供的设备属性查询 GPU 上的 SM 数量及每个 SM 的最大线程数，将这些硬件参数与测试结果对照，找到“最契合”架构的启动配置。

• 线程过少：GPU 资源未充分利用。
• 线程过多：可能耗尽寄存器或共享内存，降低占用率，影响性能。

通过这些微基准测试，你会逐步摸索出高性能的三大模式：

1. 最大化每个 SM 的线程数，同时避免过度占用资源；
2. 启动足够多的 block，确保所有 SM 都持续繁忙；
3. 优化内存访问共线性（coalescing） ，让同一 warp 的线程读写相邻地址，提高带宽利用率。

理解了 SM、warp、占用率与调度器后，你就能像调音师一般“为硬件调校”启动参数与资源使用，让每个 kernel 都在最适合的节奏中运行，并实时观察优化效果。

线程块、网格与索引

想必你对 CUDA 已经有所耳闻。CUDA（Compute Unified Device Architecture，统一计算设备架构）让我们能够在自有代码中调用 NVIDIA GPU 的强大算力。它提供了一个框架，可在 GPU 上启动成千上万的轻量级线程。在 Python 中，我们可以借助 PyCUDA 或 CuPy 等库编写、编译并运行 CUDA kernel，同时仍保持在熟悉的 Python 环境中。为本书演示，我们将默认使用 Linux 系统，因为它对 CUDA 开发具有出色的支持。

线程层次结构
在 CUDA 中，平行计算通过两级层次结构来组织：网格（Grid）和线程块（Block）。每次启动 kernel 时，我们需要指定要运行多少个线程块，以及每个线程块中包含多少线程。网格中的每个线程都会执行相同的代码，但处理不同的数据部分。可以把数据想象成一张大型电子表格，CUDA 的思路就是将表格按区域划分，每个区域对应一个线程块（Block），然后在线程块内部再将每个单元格分配给对应线程（Thread）。借助 CUDA 架构，我们可以处理一维、二维，甚至三维的大规模数据集：将每个数据元素映射到网格中的某个线程。

一维数组映射

首先，我们需要一个支持 CUDA 的 Python 环境。若是初次使用，请在虚拟环境中安装 PyCUDA 或 CuPy，并确认 GPU 可见。可以通过导入库并打印设备属性来检查。例如，用 CuPy：

import cupy as cp
print(cp.cuda.runtime.getDeviceProperties(0))

接下来，使用一个简单的一维数组示例。假设有一个一百万长度的数组，要让每个线程对其对应元素求平方。启动 kernel 时，决定每个线程块中的线程数（如 256），以及所需线程块数：

blocks_per_grid = (num_elements + threads_per_block - 1) // threads_per_block

在 CUDA kernel 内部，每个线程计算自己在整个网格中的全局索引：

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 端，运行完 kernel 后，将输出数组从设备拷贝回主机并验证结果：

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))

# 将结果拷回主机并验证
result = cp.asnumpy(output_array)
assert (result == (cp.asnumpy(input_array) ** 2)).all()
print("1D kernel 执行成功！")

通过以上步骤，我们便将一维数组映射到 CUDA 的网格和线程块系统。

二维数组映射

假如要处理一张图像或矩阵这样二维数组，需要将行和列映射到网格的两个维度。CUDA 支持多维线程块和网格，让每个线程处理一个像素或矩阵元素。

在 kernel 内部，计算行和列索引：

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;
}

if 语句同样用于边界检查，以防网格大小无法整除数组维度时越界。

在 Python 端，调整启动配置：

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 执行并验证成功！")

三维数组扩展

同样模式可扩展到三维数据（如体积图像或仿真网格）。CUDA 允许指定三维线程块和网格，每个线程计算 (z, y, x) 索引来处理其对应的体素或数据点。关键始终是：将逻辑数据位置（一维、二维或三维）映射到 CUDA 的块和线程索引，做好边界检查，并确保每个线程仅在合法范围内工作。

掌握这些映射模式后，你就可以为任意形状的数据编写自定义 kernel。在你的高性能 Python 应用中，这些网格、线程块、索引映射与边界处理的基础模式会反复出现，为真实项目打下坚实基础。

主机–设备交互概述

接下来，你将亲身体验 Python 代码如何与 GPU 通信。这个流程包含几个核心步骤：内存分配、数据传输、启动 kernel，以及了解不同的内存空间。你已经安装好了 CuPy，接下来我们将通过一个最小但完整的脚本来涵盖所有要点，帮助你建立直观认识。

GPU 上的内存分配

在 CuPy 中，当我们分配数组时，实际上是在 GPU 端（设备内存）预留一块空间，而不是在主机的 RAM 上。与传统的 NumPy 数组不同，这些数组只能在 GPU 上使用。

import cupy as cp
# 在 GPU 上分配一个包含 10000 个 float 的数组
gpu_array = cp.zeros(10_000, dtype=cp.float32)
print("已在 GPU 上分配数组：", gpu_array)

主机到设备的数据传输

主 Python 脚本在 CPU（“主机”）上运行，但计算任务在 GPU（“设备”）上执行。要利用 GPU 加速，首先需要将数据从主机传输到设备：

import numpy as np
# 在主机上创建数据
host_data = np.arange(10_000, dtype=np.float32)
# 将数据传输到 GPU（设备）
device_data = cp.asarray(host_data)
print("数据已传输到 GPU。")

若要把结果带回主机，可使用：

# 从设备将数据传回主机
result_on_host = cp.asnumpy(device_data)
print("数据已传回主机，前五个元素：", result_on_host[:5])

启动一个简单的 Kernel

CuPy 支持在 Python 中直接编写自定义 CUDA kernel。例如，将数组中每个元素乘以 2：

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;
    }
}
'''
module = cp.RawModule(code=kernel_code)
multiply_by_two = module.get_function('multiply_by_two')

# 准备参数
n = 10_000
threads_per_block = 256
blocks_per_grid = (n + threads_per_block - 1) // threads_per_block

# 启动 kernel
multiply_by_two((blocks_per_grid,), (threads_per_block,), (gpu_array, n))
print("已启动 kernel，将元素乘以二。")

# 将结果取回主机并验证
cpu_result = cp.asnumpy(gpu_array)
print("Kernel 执行后前五个元素：", cpu_result[:5])

理解 GPU 的内存空间

到目前为止，我们接触了两种主要的“内存空间”：

• 主机内存（Host memory）：CPU 系统 RAM，存放 NumPy 数组并执行 Python 脚本。
• 设备内存（Device memory）：GPU 上的内存，存放 CuPy 数组并供 kernel 读写。

GPU 编程的核心在于合理安排何时以及在何处进行数据传输，最大化在设备端完成计算，减少不必要的主机–设备往返。

此外，更高级的 kernel 还可使用：

• 共享内存：线程块内部共享的小容量低延迟存储，用于线程协作算法。
• 寄存器：每个线程私有的超高速存储。
• 常量内存与纹理内存：只读或具有空间局部性的数据专用存储，后续章节将详细探讨。

完整工作流回顾

1. 在 GPU 上分配数组
2. 在主机与设备之间传输数据
3. 编写并启动简单的 kernel
4. 观察并管理 GPU 内存使用

几乎所有支持 GPU 加速的项目，都遵循上述步骤。随着熟练度提升，你会尽量让更多计算留在设备端，仅在输入与最终输出时才进行主机–设备传输。掌握这些之后，就可以将 GPU 工作流提升到新的水平。

简单向量加法 Kernel

我们将构建第一个自定义 GPU kernel 来完成基本任务，比如向量加法。每位新手 GPU 程序员都要做这个练习，才能算入门。下面我们分别用 CuPy 和 PyCUDA 介绍如何创建、启动并验证一个将两个向量逐元素相加的 kernel，还会将 GPU 结果与 CPU 结果对比，确保一切正常。

安装 PyCUDA

我们已经准备好了 CuPy，用于大多数高层次的 GPU 数组操作。如果你想获得更细粒度的控制、直接使用 CUDA C，可以安装并使用 PyCUDA。它让我们在 Python 中编写原生 CUDA kernel，并灵活管理内存、启动方式和设备上下文。

pip install pycuda

在主机和设备上准备向量

首先，用 NumPy（CPU）和 CuPy（GPU）准备输入数据：

import numpy as np
import cupy as cp

# 向量长度
N = 1_000_000

# 在主机上创建两个随机向量
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)

编写向量加法 Kernel

我们的 kernel 接收两个输入数组 a、b，并将它们逐元素相加，结果写入第三个数组 c。以下是 CUDA C 版本的 kernel 代码：

kernel_code = r'''
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];
    }
}
'''

使用 CuPy 运行 Kernel

CuPy 提供了 RawModule 接口，可直接在 Python 中编译并启动 CUDA kernel。首先分配输出数组：

c_gpu = cp.empty_like(a_gpu)
threads_per_block = 256
blocks_per_grid = (N + threads_per_block - 1) // threads_per_block

module = cp.RawModule(code=kernel_code)
vector_add = module.get_function('vector_add')

# 启动 kernel
vector_add((blocks_per_grid,), (threads_per_block,), (a_gpu, b_gpu, c_gpu, N))

# 将结果拷回主机并验证
c_result = cp.asnumpy(c_gpu)
print("GPU 和 CPU 结果是否一致？", np.allclose(c_result, c_cpu))

输出应为 True，表示我们的第一个自定义 kernel 在 CuPy 下运行成功！

使用 PyCUDA 编写并运行相同 Kernel

接下来，用 PyCUDA 更深入地掌控 CUDA 框架。首先导入模块并准备 GPU 数组：

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

# 将输入数据传输到 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)

编译 kernel 并获取函数句柄：

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

使用 PyCUDA 的启动语法调用 kernel：

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)
)

最后把结果取回主机并验证：

c_result_py = c_gpu_py.get()
print("PyCUDA GPU 结果与 CPU 是否匹配？", np.allclose(c_result_py, c_cpu))

同样会得到 True。

至此，我们已经使用 CuPy 和 PyCUDA 分别编写、编译并启动了自定义 CUDA kernel，完成了向量加法，并验证了 GPU 计算结果与 CPU 等价。这一模式是后续所有复杂 GPU 操作的基础。

总结

总体而言，本章为实用 GPU 编程奠定了坚实基础。我们首先了解了高性能数据处理的现实需求，认识到在现代工作负载下，传统 CPU 经常遇到性能瓶颈；同时明白了为什么拥有数千轻量级核心的 GPU 能显著加速诸如机器学习推理和大规模数据分析等任务。接着，我们深入剖析了流式多处理器（SM）的设计，学习了 SM、warp 以及占用率等概念如何高效地实现并行执行，还掌握了 GPU 隐藏内存延迟、持续保持硬件繁忙的各种技巧。

随后，我们着手将多维数组映射到 CUDA 的网格和线程块，熟练掌握了线程层次结构、索引计算与健壮的边界检查方法。我们在 Python 中完整演练了 GPU 数据处理工作流：从使用 CuPy 在设备上分配内存、在主机与设备之间传输数据，到启动自定义 kernel，再将结果拷回主机进行验证。我们还分别在 CuPy 和 PyCUDA 中编写并运行了第一个自定义 CUDA C kernel，体验了从 CPU 上准备数据到验证 GPU 输出的全流程。

凭借这些技能与实践习惯，我们已经做好准备，在本书后续章节中深入探索更高级、更具创造性的 GPU 编程应用。