现代-CUDA-C---编程-2025-笔记-全-

现代 CUDA C++ 编程 2025 笔记（全）

001：使用并行算法加速应用

概述

在本节课中，我们将学习如何使用NVIDIA最新的编程工具在GPU上进行编程。我们将从理解CPU与GPU的基本区别开始，逐步学习如何通过简单的代码修改，利用并行算法和CUDA生态中的库（特别是Thrust）来显著加速应用程序。课程内容涵盖执行空间、并行算法、迭代器、内存空间等核心概念，旨在让初学者能够轻松上手GPU编程。

1. CPU与GPU的区别

上一节我们介绍了课程概述，本节中我们来看看CPU和GPU的核心区别。

我们可以用公交车和汽车来类比CPU和GPU。公交车一次可以运送很多人，但速度通常比汽车慢。汽车速度更快，但一次只能运送少数人。哪个“更快”取决于你的目标：是快速运送少数人，还是高效运送大批人。

CPU和GPU的关系与此类似。

CPU：延迟低（约100纳秒），但内存带宽有限（约100 GB/s）。它像汽车，处理少量数据时速度极快。
GPU：延迟较高（约500纳秒），但拥有极高的内存带宽（约1 TB/s）。它像公交车，虽然处理单个数据项稍慢，但能同时处理海量数据。

因此，选择CPU还是GPU取决于问题的规模：处理少量数据时CPU更高效；处理海量数据时，GPU的巨大带宽优势将带来显著的性能提升。

2. 一个简单的模拟问题

理解了硬件区别后，我们来看一个具体的例子：模拟多个杯子在室温下的冷却过程。

假设我们有多个杯子，每个都有不同的初始温度。我们想模拟它们随时间逐渐接近室温（20°C）的过程。每个杯子的温度变化遵循一个简单的独立公式：

T_next = T_prev + k * (T_room - T_prev)

其中，k是冷却系数。每个杯子的温度演化是独立的，互不影响。这种大量独立、可并行执行的任务，正是GPU擅长处理的。

以下是使用C++标准库在CPU上实现的代码：

#include <vector>
#include <algorithm>
#include <iostream>

int main() {
    float k = 0.5f; // 冷却系数
    std::vector<float> temps = {42.0f, 24.0f, 50.0f}; // 初始温度
    float room_temp = 20.0f;

    // 定义温度更新公式的Lambda函数
    auto update_temp = [k, room_temp](float t_prev) {
        return t_prev + k * (room_temp - t_prev);
    };

    // 模拟3个时间步长
    for (int step = 0; step < 3; ++step) {
        std::transform(temps.begin(), temps.end(), temps.begin(), update_temp);
        // 打印当前温度
        for (auto t : temps) std::cout << t << " ";
        std::cout << std::endl;
    }
    return 0;
}

std::transform算法将update_temp函数应用于temps向量的每个元素。目前，这一切都在CPU上串行执行。

3. 执行空间：指定代码运行位置

上一节我们看到了CPU上的串行实现，本节中我们来看看如何告诉编译器让代码在GPU上运行。

编译器（如g++用于CPU，nvcc用于NVIDIA GPU）负责将高级C++代码转换为特定硬件能执行的机器指令。然而，仅仅使用nvcc编译并不会自动让代码在GPU上运行。我们需要在代码中明确指定哪些部分应在GPU（设备）上执行，哪些部分应在CPU（主机）上执行。这就是执行空间的概念。

我们通常不从头编写所有GPU代码，而是使用库。NVIDIA提供了Thrust库，它提供了类似于C++标准模板库（STL）的接口，但算法能在GPU上执行。这让我们能够以熟悉的方式编写高性能GPU程序。

要将之前的冷却模拟代码移植到GPU，需要进行几处修改：

#include <thrust/device_vector.h>
#include <thrust/transform.h>
#include <thrust/execution_policy.h>

int main() {
    float k = 0.5f;
    thrust::device_vector<float> temps = {42.0f, 24.0f, 50.0f}; // 数据位于GPU
    float room_temp = 20.0f;

    // Lambda函数需要标明可在主机和设备上执行
    auto update_temp = [k, room_temp] __host__ __device__ (float t_prev) {
        return t_prev + k * (room_temp - t_prev);
    };

    for (int step = 0; step < 3; ++step) {
        // 使用thrust::transform并指定在设备上执行
        thrust::transform(thrust::device, temps.begin(), temps.end(), temps.begin(), update_temp);
        // ... 输出结果（需要将数据拷贝回CPU）
    }
    return 0;
}

关键修改点：

容器：std::vector 替换为 thrust::device_vector，确保数据分配在GPU内存中。
执行空间限定符：在Lambda函数前添加__host__ __device__，指示编译器为此函数生成CPU和GPU均可执行的代码。
算法与执行策略：std::transform 替换为 thrust::transform，并通过 thrust::device 参数明确要求算法在GPU设备上执行。

执行空间限定符（__host__, __device__）在编译时告诉编译器为哪个架构生成代码。执行策略（thrust::host, thrust::device）则在运行时告诉Thrust库在哪个硬件上执行算法。两者需要配合使用。

4. 并行算法与Thrust

我们已经学会了如何运行一个简单的变换操作，本节中我们来看看Thrust中更多的并行算法。

Thrust提供了丰富的并行算法，其接口与STL算法相似，使得移植代码变得简单。例如，计算一个温度序列的中位数，可以先排序再取中间值。在CPU上，我们使用std::sort。在GPU上，我们只需将其替换为thrust::sort。

// CPU版本 (串行)
#include <algorithm>
float median_cpu(std::vector<float>& data) {
    std::sort(data.begin(), data.end());
    return data[data.size() / 2];
}

// GPU版本 (并行)
#include <thrust/device_vector.h>
#include <thrust/sort.h>
float median_gpu(thrust::device_vector<float>& data) {
    thrust::sort(thrust::device, data.begin(), data.end());
    return data[data.size() / 2];
}

直接尝试在标记为__device__的函数中调用std::sort会失败，因为STL算法没有为GPU编译的版本。而thrust::sort是专为GPU设计的并行排序算法，能够充分利用GPU的众核优势。

5. 组合算法与性能陷阱

有时我们需要组合多个算法来实现复杂功能。例如，计算两个时间步之间温度的最大差值。朴素的方法是先计算差值数组，再求该数组的最大值。

// 朴素方法：先变换，后归约
thrust::device_vector<float> diff(N);
// 1. 变换：计算差值
thrust::transform(thrust::device, A.begin(), A.end(), B.begin(), diff.begin(),
                  [] __host__ __device__ (float a, float b) { return fabs(a - b); });
// 2. 归约：求最大值
float max_diff = thrust::reduce(thrust::device, diff.begin(), diff.end(),
                                0.0f, thrust::maximum<float>());

这种方法需要分配临时数组diff，导致 3N 次内存读取（两次读A/B，一次读diff）和 N 次内存写入（写diff）。如果手动编写循环，可以在单次遍历中同时计算差值和更新最大值，只需 2N 次读取和 1 次写入（存储最终结果）。这引出了一个问题：如何让Thrust以更高效的方式组合操作？

答案就是使用Fancy Iterators（花式迭代器）。

6. Fancy Iterators（花式迭代器）

上一节我们发现组合算法可能产生额外开销，本节中我们介绍一种强大的工具——Fancy Iterators来消除这种开销。

迭代器是指针的泛化，它提供了类似指针的接口（如解引用、递增），但行为可以自定义。Thrust提供的Fancy Iterators允许我们在数据“流动”的过程中进行即时转换，而无需物化（分配内存存储）中间结果。

以下是几种重要的Fancy Iterators：

计数迭代器 (`thrust::counting_iterator`)

它不存储任何序列，而是在被访问时动态生成连续的整数。

// 打印1到42。没有内存分配，数字是即时生成的。
thrust::for_each(thrust::device,
                 thrust::make_counting_iterator(1),
                 thrust::make_counting_iterator(43),
                 [] __host__ __device__ (int i) { printf("%d\n", i); });

变换迭代器 (`thrust::transform_iterator`)

它在迭代过程中，将给定的函数应用于底层迭代器的每个元素。

thrust::device_vector<int> vec = {1, 2, 3};
// 创建一个视图，使得所有元素看起来都乘以了2
auto doubled_view = thrust::make_transform_iterator(
    vec.begin(),
    [] __host__ __device__ (int x) { return x * 2; }
);
// 此时遍历doubled_view，会得到2, 4, 6，但原vec并未改变。

zip迭代器 (`thrust::zip_iterator`)

它将多个序列“压缩”成一个，迭代时同时返回多个序列中对应位置的元素组成的元组。

thrust::device_vector<int> A = {0, 1, 2};
thrust::device_vector<int> B = {5, 4, 2};
// 创建一个迭代器，同时遍历A和B
auto zipped = thrust::make_zip_iterator(thrust::make_tuple(A.begin(), B.begin()));
// 解引用zipped会得到一个元组，例如 (0, 5)

7. 使用Fancy Iterators优化

现在，我们可以使用Fancy Iterators来优化之前计算最大差值的问题。我们可以创建一个迭代器，它“看起来”像一个存储了绝对差值的序列，但实际上是在迭代时动态计算。

// 高效方法：使用zip和transform迭代器组合，在归约过程中即时计算差值
auto diff_iterator = thrust::make_transform_iterator(
    thrust::make_zip_iterator(thrust::make_tuple(A.begin(), B.begin())),
    [] __host__ __device__ (const thrust::tuple<float, float>& t) {
        return fabs(thrust::get<0>(t) - thrust::get<1>(t)); // 计算绝对值差
    }
);

float max_diff = thrust::reduce(thrust::device,
                                 diff_iterator,
                                 diff_iterator + N,
                                 0.0f,
                                 thrust::maximum<float>());

通过这种方式，thrust::reduce遍历的“数组”就是我们动态生成的差值序列。整个过程只需要读取A和B（共2N次），一次归约操作，并写入最终结果。完全避免了临时数组diff的分配和读写，性能显著提升，内存占用也更少。

8. 更复杂的案例：2D热扩散模拟

前面的例子都是独立操作，本节我们看一个更贴近实际、元素间有依赖关系的案例：二维热扩散模拟。

在2D热扩散中，一个单元格下一时刻的温度取决于它自身及其上下左右四个邻居的当前温度。我们使用一个5点模板（Stencil）来计算。虽然计算涉及邻居，但每个单元格的计算仍然可以并行进行，因为所有计算都基于同一时刻的数据。

我们需要将二维网格展平为一维数组进行处理。对于每个线程（处理一个单元格），我们需要根据其全局一维索引id计算出对应的二维坐标(row, col)。

int row = id / width;
int col = id % width;

然后，在Lambda函数中，我们根据row和col访问自身及邻居的数据，应用物理公式计算新温度，并注意处理边界条件（例如，边界温度保持不变）。

Thrust的transform或tabulate算法非常适合这种模式：每个线程处理一个独立的单元格。

thrust::tabulate(thrust::device,
                 output.begin(), output.end(),
                 [input_ptr, width, height] __host__ __device__ (int id) {
                     int row = id / width;
                     int col = id % width;
                     // 边界检查
                     if (row == 0 || row == height-1 || col == 0 || col == width-1) {
                         return input_ptr[id]; // 边界保持不变
                     }
                     // 内部点：使用5点模板计算新温度
                     float center = input_ptr[id];
                     float left   = input_ptr[id - 1];
                     float right  = input_ptr[id + 1];
                     float top    = input_ptr[id - width];
                     float bottom = input_ptr[id + width];
                     return some_physics_formula(center, left, right, top, bottom);
                 });

9. 使用CUDA C++标准库中的词汇类型

在编写GPU代码时，我们经常使用像std::pair, std::tuple这样的词汇类型。然而，STL中的这些类型通常没有标记为__device__，无法在GPU代码中直接使用。

NVIDIA提供了libcu++（CUDA C++ Standard Library），它包含了为GPU优化且带有__host__ __device__注解的对应类型，可以作为STL类型的直接替代品。

// 使用 libcu++ 中的 pair
#include <cuda/std/utility>
auto my_pair = cuda::std::make_pair(42, 3.14f); // 可在主机和设备代码中使用

此外，对于多维数组访问，手动计算索引容易出错且代码冗长。libcu++也提供了mdspan（多维视图），它允许我们以多维方式访问底层的一维数据，使代码更清晰、更安全。

#include <cuda/std/mdspan>
// 假设data_ptr指向一个height * width大小的float数组
cuda::std::mdspan<float, cuda::std::dextents<int, 2>> grid_2d(data_ptr, height, width);
// 现在可以像使用2D数组一样访问
float val = grid_2d(row, col); // 访问第row行，第col列
float neighbor_left = grid_2d(row, col - 1);

使用mdspan后，热扩散模拟中的邻居访问代码变得非常直观，极大地减少了索引计算错误。

10. 高级并行模式：按Key归约

在模拟中，我们可能想计算每行温度的总和。一种简单的方法是使用tabulate，让每个线程处理一行，并在线程内循环累加该行所有列。但这种方法在行数远少于列数时，GPU并行度很低，利用率不足。

更好的方法是使用按Key归约。我们将每个单元格的行号作为Key，温度值作为Value。thrust::reduce_by_key算法会自动将具有相同Key（即同一行）的所有Value归约起来（例如求和）。

// 1. 创建Key数组（行号）。可以使用计数迭代器和变换迭代器动态生成，避免分配。
auto row_id_iterator = thrust::make_transform_iterator(
    thrust::make_counting_iterator(0),
    [width] __host__ __device__ (int id) { return id / width; } // 动态计算行号
);

// 2. 执行按Key归约求和
thrust::device_vector<float> row_sums(height); // 存储每行的和
// 我们不需要输出的Key，使用discard_iterator丢弃
thrust::reduce_by_key(thrust::device,
                      row_id_iterator, row_id_iterator + total_cells, // Key范围
                      temperatures.begin(),                           // Value范围
                      thrust::make_discard_iterator(),                // 输出Key（丢弃）
                      row_sums.begin()                               // 输出Value（行和）
                     );

这种方法将工作均匀分布到所有处理单元格的线程上，并行度最大化，性能远高于tabulate内嵌循环的方法。

11. 内存空间与数据传输优化

到目前为止，我们主要关注计算优化。本节我们关注另一个关键性能因素：内存。

我们一直使用thrust::universal_vector，它位于统一内存中。统一内存让数据可以被CPU和GPU透明地访问，系统在背后自动进行数据迁移。然而，这种自动迁移是有成本的。如果CPU和GPU频繁交替访问数据，就会产生大量的隐式数据传输，严重拖慢程序。

例如，在模拟循环中：在GPU上计算一个时间步 -> 将结果保存到磁盘（需要CPU访问） -> 在GPU上计算下一个时间步。使用统一内存时，保存到磁盘会触发设备到主机的数据传输，而下一个时间步的计算又会触发主机到设备的数据传输，导致计算被传输延迟。

解决方案是使用显式的内存空间，并手动控制数据传输：

thrust::host_vector: 数据仅位于主机内存。
thrust::device_vector: 数据仅位于设备内存。
thrust::copy: 在主机和设备之间显式拷贝数据。

优化后的流程：

thrust::device_vector<float> d_data(N); // 设备数据
thrust::host_vector<float> h_data(N);   // 主机数据

for (int step = 0; step < steps; ++step) {
    // 1. 在GPU上模拟
    simulate_on_gpu(d_data);
    // 2. 显式拷贝到主机（仅在需要时）
    thrust::copy(d_data.begin(), d_data.end(), h_data.begin());
    // 3. 在CPU上保存到磁盘
    save_to_disk(h_data);
    // 4. 下一个循环，d_data仍在GPU上，无需传输即可直接计算
}

通过这种方式，我们确保了计算数据常驻GPU，只在必要时支付一次性的、显式的数据传输开销，从而避免了模拟循环中的性能抖动。

总结

本节课中我们一起学习了现代CUDA C++编程的基础知识，重点是如何使用Thrust库和并行算法加速应用。我们涵盖了以下核心内容：

硬件理解：CPU与GPU在延迟和带宽上的根本区别，决定了它们分别适合处理不同规模的问题。
执行模型：通过__host__/__device__限定符和Thrust的执行策略（thrust::device），明确控制代码在何处编译以及在何处执行。
Thrust库：作为GPU上的“STL”，提供了丰富的并行算法（如transform, sort, reduce），使得移植CPU代码到GPU变得简单。
Fancy Iterators：强大的工具，通过组合计数、变换、zip等迭代器，能够动态生成和转换数据视图，消除不必要的中间存储和内存访问，极大提升性能。
复杂模式实现：学习了如何利用tabulate处理映射问题，以及使用reduce_by_key实现高效的分组归约，以充分利用GPU并行性。
CUDA C++标准库：使用libcu++中的pair, mdspan等类型，编写更安全、更清晰且兼容GPU的代码。
内存管理：理解统一内存的便利与潜在开销，学会使用显式的host_vector/device_vector和thrust::copy来精确控制数据位置和传输，避免性能陷阱。

通过掌握这些概念，你已经能够开始使用GPU的强大并行能力来加速自己的C++应用程序。记住关键思路：将问题表达为数据并行的操作，利用高效的库和算法，并谨慎管理内存和数据传输。

002：异步性与CUDA流

在本节课中，我们将要学习如何通过异步编程和CUDA流来优化GPU程序的性能。我们将探讨如何让CPU和GPU的工作重叠，以及如何使用CUDA流来让GPU上的不同任务并发执行，从而更充分地利用系统资源。

概述

在上一节中，我们介绍了执行与内存空间、词汇数据类型和并行算法等核心概念，并创建了一个功能性的2D热方程模拟器。虽然它目前在GPU上运行，但我们的模拟器仍未充分利用系统的全部能力。

在本节中，我们将探索一些额外的技术来优化性能。首先，我们将讨论异步性，以及如何重叠通信和输入/输出操作以加速程序。然后，我们将学习CUDA流的概念，看看如何让计算与内存传输重叠。最后，我们将学习页锁定内存，了解如何使内存传输异步且更快。

2.2：异步性与重叠

上一节我们介绍了基本的GPU编程模型。本节中我们来看看如何通过异步操作来提升效率。

回顾我们模拟器的当前状态：首先将数据从GPU复制到主机，然后写入磁盘。一旦完成，我们使用Thrust在GPU上计算下一组温度值。我们可以看到，CPU启动了计算，然后GPU进行计算，而CPU在此期间空闲等待。

这是因为thrust::tabulate在底层调用CUDA时，会等待GPU完成计算后才将控制权交还给CPU。这是一个被浪费的机会。要编写高效的异构程序，我们需要充分利用所有资源，包括CPU。

那么，我们能否在GPU计算时，为CPU找到有意义的工作呢？答案是肯定的。首先，让我们看看数据依赖关系。写入磁盘的步骤读取的是上一步操作的主机内存结果，而计算步骤读取的是设备内存。由于这些操作使用不同的数据，它们彼此不依赖。因此，CPU可以在GPU计算下一个温度值时，将结果写入磁盘。

为了实现这种重叠，我们理想情况下需要将Thrust调用拆分为独立的启动和等待步骤。不幸的是，Thrust本身不支持这种拆分。但Thrust是基于另一个核心库CUB实现的，而CUB支持我们需要的异步性。

同步与异步操作的区别

首先，让我们更好地理解同步操作和异步操作之间的区别。

CUDA操作本质上是异步的。这意味着当你从主机启动一个在设备（GPU）上执行的操作时，它会立即将控制权交还给CPU。它启动GPU上的操作，然后直接返回控制权。

Thrust的做法则不同：它在底层使用CUDA启动计算，然后使用cudaDeviceSynchronize等待GPU完成，之后才将CPU控制权交还。这就是为什么thrust::tabulate操作是阻塞的。

以下是一个示例，展示了使用Thrust（同步）和CUB（异步）在CPU等待时间上的差异：

// 同步版本 (Thrust) - CPU等待时间随元素数量增加
auto start = std::chrono::high_resolution_clock::now();
thrust::tabulate(device_data.begin(), device_data.end(), computation_functor());
auto end = std::chrono::high_resolution_clock::now();
// CPU 等待了 (end - start) 时间

// 异步版本 (CUB) - CPU等待时间很短且恒定
auto start = std::chrono::high_resolution_clock::now();
cub::DeviceTransform::Transform(nullptr, temp_storage_size,
                                device_data_in, device_data_out,
                                computation_functor(), stream);
auto end = std::chrono::high_resolution_clock::now();
// CPU 几乎立即获得控制权，等待时间极短
// 如果需要等待结果，可以显式调用 cudaDeviceSynchronize()

使用CUB启动异步工作后，CPU被立即释放，因此它可以在GPU执行其他操作时并行工作。在我们的案例中，我们可以在GPU进行计算时，开始将数据写入磁盘。然后，我们只需要在GPU完成模拟步骤后进行同步。由于我们无论如何都需要等待GPU，我们并没有损失任何时间。相反，等待步骤变得更短，从而更高效地利用了系统。因此，我们可以预期这种方法会带来显著的性能提升。

之前我们的流程是：从设备复制到主机 -> 写入磁盘 -> 调用thrust::tabulate（进行计算并等待）-> 重复。
现在我们的目标是：从设备复制到主机 -> 调用CUB异步启动计算 -> 在GPU计算的同时，CPU写入磁盘 -> 在进入下一步复制前，等待计算完成。

练习：实现计算与I/O的重叠

在这个练习中，你需要将thrust::tabulate替换为使用cub::DeviceTransform::Transform，以便从需要等待数据的同步代码转变为异步代码，使得CPU可以在GPU工作时写入磁盘。你还需要在正确的位置调用cudaDeviceSynchronize，以确保在进行下一次数据复制时，数据确实已被处理。

以下是实现的关键步骤：

将 thrust::tabulate 替换为 cub::DeviceTransform::Transform。
确保在写入磁盘操作之后、下一次数据复制之前，添加 cudaDeviceSynchronize() 来等待GPU计算完成。

2.3：使用Nsight Systems进行性能分析

随着我们深入异步编程，理解和调试程序行为、追踪性能问题可能变得棘手。因为CPU和GPU上的事件同时发生，不容易理清。

为此，NVIDIA提供了Nsight Systems工具。它是一个软件，允许你可视化程序在CPU端和GPU端的执行情况。你可以在时间线上直观地看到CPU何时工作、何时启动异步工作、GPU何时工作。

练习：生成并分析Nsight Systems报告

在这个练习中，你将学习如何为你的代码生成Nsight Systems报告。

打开对应练习，执行命令生成Nsight Systems报告。
报告生成后，通过Nsight Systems UI打开它。
在UI中，你可以可视化主机和设备之间发生的事件。你可以选择时间区域来查看特定时间段内的详细信息。例如，你可以放大查看0.5秒到0.7秒之间发生的事件。
你的练习是：生成报告，使用Nsight Systems打开它，并尝试理解：GPU计算何时进行？CPU何时异步启动GPU计算？CPU何时写入磁盘？CPU何时等待？数据传输何时发生？请花时间熟悉这个工具。

分析报告后，我们可以看到重叠执行的效果：CPU启动GPU计算（耗时很短），然后GPU开始实际计算。在GPU进行长时间计算的同时，CPU时间线上显示正在执行写入磁盘的I/O操作（fopen, fwrite, fclose）。我们并不担心写入磁盘耗时较长，因为无论如何GPU计算都需要更长时间。通过重叠这两者，在任何时间点，CPU和GPU资源都得到了利用。

如果你对比早期同步版本的模拟器，你会看到清晰的性能差异。在同步版本中，每个CPU调用都精确对应一个GPU调用，并且CPU需要等待每个调用完成，最后才进行写入操作。而在异步版本中，CPU快速启动所有GPU调用（因为异步启动很快），然后在GPU异步计算的同时进行磁盘写入，最后再进行同步。这使得异步版本比同步版本快大约2倍。

2.4：使用NVTX标注代码

分析两个版本在时间线上的差异仍然需要手动将时间线上的事件与代码匹配，这可能有些繁琐。

幸运的是，我们有一个工具可以更轻松地将Nsight Systems中看到的内容与你的代码关联起来，这就是NVTX（NVIDIA Tools Extension）。当处理具有许多函数、GPU和CPU事件同时发生的复杂应用程序时，我们通常会使用NVTX。NVTX API允许你将自定义标记和范围直接插入到代码中。这样，当你在Nsight Systems报告中查看时，就会看到对应的标签。

例如，你可以在主循环中添加一个范围，命名为“Write Step”，并包含迭代ID。这样，在Nsight Systems中，你会看到“Write Step”实际上由fopen、fwrite和fclose组成。你不需要知道写入磁盘的具体底层调用，只需要知道在这个阶段你正在进行“写入磁盘”操作。

练习：使用NVTX标注代码步骤

你的练习是为所有不同的步骤添加NVTX范围标注：复制步骤、计算步骤、写入步骤以及等待步骤。通过添加这些标注，你将能够在Nsight Systems中更好地可视化所有事件。

添加标注后，你的Nsight Systems视图将变得更加清晰易懂。在CPU时间线上，你可以清楚地看到：首先进行复制，然后启动计算（注意，这只是启动命令，并非实际计算），接着是漫长的磁盘写入，最后是等待。在GPU部分，你可以看到复制步骤和占用大部分时间的计算步骤。现在，CPU和GPU之间发生的事情一目了然。

2.5：使用CUDA流重叠复制与计算

至此，我们已经学会了如何使用异步性来重叠CPU和GPU任务。通过将磁盘写入与GPU计算重叠，我们获得了巨大的性能提升。我们的代码现在的工作流程是：将数据从GPU复制到CPU，使用CUB启动异步计算，然后将结果写入磁盘。

然而，我们仍有改进空间。GPU计算实际上不必等待复制完成。我们可以将相同的重叠策略应用于复制和计算步骤。但这要求复制操作本身是异步的，而目前它并不是。

thrust::copy是通过一个更低级的API cudaMemcpyAsync实现的。这个函数可以在不阻塞（即异步地）的情况下在主机和设备内存之间复制数据。让我们仔细看看如何使用它来重叠复制和计算步骤。

首先，了解cudaMemcpyAsync如何被Thrust调用。cudaMemcpyAsync接收四个参数：目标内存地址指针、源内存地址指针、要复制的字节大小（注意是字节数，不是元素数量）以及复制方向（如主机到设备）。

异步操作的一个微妙之处是，错误并不总是在发生时立即浮现。因为操作是异步的，错误可能会在稍后才被捕获。例如，你可能启动了一个越界的内核，但GPU并没有立即开始执行。你可能在之后启动另一个内核时，才通过cudaMemcpyAsync的启动发现前一个内核的错误。因此，在每次CUDA调用后检查错误代码至关重要，并且要意识到错误可能比你预期的出现得晚。

那么，简单地将我们的示例中的thrust::copy替换为cudaMemcpyAsync是否足以重叠计算与复制？答案是否定的。因为默认情况下，GPU上的所有操作都是有序的。即使你调用了非阻塞的cudaMemcpyAsync，CPU会立即返回，但GPU仍然会按顺序执行它和之后的其他操作。这意味着GPU将等待复制完成，然后再进行下一个计算步骤。默认情况下，这两者不会重叠。

控制操作顺序的机制就是CUDA流。默认情况下，如果不指定流，则使用默认流，这意味着GPU上的一切都将按顺序发生。我们不必使用默认流，我们可以创建自己的CUDA流，实际上可以创建任意多个。当操作在不同的流中运行时，它们彼此之间不再有序，因此可以完全并发地执行。

这正是我们想要的。如果我们在一个流中发送异步复制，在另一个流中进行计算，那么这两个操作就可以重叠。

CUDA流的基本操作

以下是CUDA流的基本操作：

声明流对象：cudaStream_t copy_stream, compute_stream;
创建流：cudaStreamCreate(&copy_stream); cudaStreamCreate(&compute_stream);
在流中启动异步操作：将流作为参数传递给如cudaMemcpyAsync或CUB的函数。
流同步：cudaStreamSynchronize(stream) 使CPU等待指定流中的所有操作完成。
销毁流：cudaStreamDestroy(stream); 在使用完毕后销毁流。

通常，我们建议显式使用流和cudaStreamSynchronize，而不是使用默认流和cudaDeviceSynchronize。因为cudaDeviceSynchronize会导致CPU等待所有流上的所有GPU操作完成，而你可能只需要等待特定流中的操作。

引入数据竞争

让我们开始为计算和复制创建CUDA流。我们声明两个流，分配它们，然后调用cudaMemcpyAsync并传递复制流以启动异步复制。接着，我们通过向CUB传递计算流来启动计算。由于复制和计算发生在不同的流中，它们本应在GPU上并行执行。

然后，因为复制是异步的，我们必须确保它在读取主机数据之前确实完成了。我们需要在启动磁盘写入之前，在复制流上调用cudaStreamSynchronize。之后，我们还需要在进入下一步之前同步计算流。

然而，恭喜，我们刚刚引入了第一个数据竞争。让我们遍历前几次迭代来看看原因。第一次计算迭代从d_prev读取并写入d_temp，这没问题。但下一次迭代从d_temp读取并写回d_prev。由于这些操作在不同的流中运行，GPU可能会以意想不到的方式交错执行它们。这意味着当数据仍在从d_prev复制到主机时，它可能会被另一个流中的下一次计算迭代覆盖。我们创建了一个数据竞争，导致结果可能不一致。

解决方案：使用暂存缓冲区

像许多计算机科学问题一样，我们可以通过添加一个间接层来解决这个问题。我们可以分配一个额外的设备缓冲区作为暂存区。首先，将数据同步复制到这个暂存缓冲区，然后再从暂存缓冲区复制回CPU。由于没有其他操作写入这个暂存缓冲区，因此没有数据竞争的风险。

这是否违背了重叠计算与传输的初衷？乍一看，我们引入了额外的工作（一次复制和一个额外的缓冲区），这可能会使事情变慢。但如果我们仔细审视整个系统的带宽，这种方法就合理得多。GPU内部的数据复制（设备到设备）可以达到每秒数千GB的速度，而通过PCIe总线在CPU和GPU之间移动数据则要慢得多。设备到设备的复制速度非常快，与主机到设备的复制相比，其开销实际上可以忽略不计。

练习：实现流与异步复制

在这个任务中，你将用cudaMemcpyAsync替换thrust::copy，使设备到主机的传输异步化。你还需要将复制操作放在独立的流中，并按照幻灯片图所示精确地同步它们。最后，在Nsight Systems中分析代码，看其行为是否符合预期。

解决方案的关键步骤包括：

创建两个流（复制流和计算流）。
首先，使用thrust::copy（或cudaMemcpy）将设备数据同步复制到暂存缓冲区。这可以防止数据竞争。
然后，在复制流上调用cudaMemcpyAsync，将数据从暂存缓冲区异步复制到主机。使用thrust::raw_pointer_cast获取底层指针。
更新CUB变换调用，使其在我们创建的计算流上运行。
在写入磁盘前，同步复制流 (cudaStreamSynchronize(copy_stream))。
在进入下一次迭代前，同步计算流 (cudaStreamSynchronize(compute_stream))。

分析Nsight Systems的性能剖析图，有好消息也有坏消息。好消息是，我们对传输成本的直觉是正确的：设备到设备的复制（左侧）明显快于设备到主机的复制（右侧）。这证实了PCIe传输是真正的瓶颈，而添加设备到设备的复制开销很小。我们也不再存在数据竞争。

但坏消息是，即使我们使用了不同的流，复制和计算之间仍然没有重叠。计算似乎在等待复制完成后才启动。这与我们使用流的初衷相悖。我们遗漏了某些东西。

2.6：页锁定内存（Pinned Memory）

要理解为什么我们的复制和计算操作没有重叠，我们需要看看内存实际上是如何工作的。

当程序使用虚拟内存时，内存被划分为内存页。大多数时候，这些内存页由物理RAM支持。但如果RAM空间不足，操作系统可以将内存页从RAM交换到磁盘。这意味着在任何给定时刻，你都不能完全确定某些内存页是在RAM中还是在磁盘上。

为了确保特定的内存页位于RAM中并且不会被移动到磁盘，我们需要所谓的页锁定内存或固定内存。通过固定这些内存页，你可以确保操作系统不允许将这些内存页从RAM移回磁盘。

那么，这与GPU有什么关系呢？事实证明，GPU只能从固定内存中直接读取数据。那么，在我们没有为CPU和GPU之间的数据传输做任何特殊处理的情况下，事情是如何工作的呢？

在幕后，CUDA运行时使用一个小的固定缓冲区作为传输的暂存区。当你从常规的可分页内存复制数据时，驱动程序会将数据块移动到固定暂存缓冲区，然后将其发送到GPU，等待完成，然后重复这个过程。这个过程有效地使复制操作变成了同步的，即使我们使用了cudaMemcpyAsync，因为它需要反复执行这个过程，从而阻塞了其他操作。

好消息是，我们可以通过自己分配固定内存来避免这种情况。在Thrust中，有一个简单的接口：thrust::universal_host_pinned_vector。它允许我们分配保证是固定的内存，从而绕过这个隐藏的、导致复制同步化的暂存步骤。

练习：使用固定内存

在这个最后的练习中，你需要分配一个固定内存向量，而不是常规的主机向量。将容器从thrust::host_vector更改为thrust::universal_host_pinned_vector。完成更改后，启动Nsight Systems查看差异并进行可视化。

解决方案的更改非常简单：只需使用thrust::universal_host_pinned_vector替代thrust::host_vector。完成之后，设备到主机的传输现在可以与计算重叠了。在时间线上，我们可以看到GPU可以启动一个接一个的计算操作，而设备到主机的内存复制同时发生。这是因为现在数据位于固定区域，因此不需要先将数据分块复制到固定暂存区，它可以直接并行执行这两项操作。这样，我们通过使用流和固定内存，成功重叠了设备到主机的传输和计算，从而进一步提高了程序的效率。

总结

本节课中我们一起学习了如何通过异步编程和CUDA流优化GPU程序性能。关键要点如下：

优先使用异步性：当可以重叠主机和设备工作时，应使用异步操作。例如，将Thrust替换为CUB，利用其异步编程模型，使CPU任务与GPU任务重叠。调用CUB时，主机以异步方式启动设备计算，这意味着它不会等待GPU完成。无论GPU上要解决的问题规模多大，主机几乎都会立即返回，然后CPU可以继续执行其他工作，而GPU则在后台进行可能很耗时的计算。

使用Nsight Systems进行性能分析：我们学会了如何生成报告、打开报告，并查看CPU和GPU上的执行情况，识别时间浪费点或GPU利用率不足的情况。该工具也让我们能够可视化异步操作。

使用NVTX简化代码映射：NVTX允许我们通过手动在代码中有意义的地方添加范围标记，使性能可视化变得更加容易。我们可以直接将代码段映射到Nsight Systems报告中的事件。
利用CUDA流实现GPU内部并行：异步性通过CUB允许CPU和GPU工作并行化。而使用CUDA流则允许我们在GPU上并行化工作。例如，计算步骤和复制步骤可以并行执行。默认情况下，所有操作都在GPU的默认流上启动，因此工作将按顺序运行。为了并行运行可并行的操作，我们应该使用不同的流。cudaMemcpyAsync不仅允许我们异步启动主机到设备的内存复制，还允许我们添加流参数，从而能够在一个流上启动复制，在另一个流上启动计算。
使用固定内存加速传输：默认情况下，主机内存是可分页的，可能位于磁盘上。为了确保复制快速且非阻塞，我们应该使用正确的Thrust主机容器（如thrust::universal_host_pinned_vector）来确保内存是固定的，从而使复制操作非阻塞且快速。

通过掌握这些概念和技术，你可以显著提高CUDA应用程序的性能和资源利用率。

003：通过CUDA内核实现新算法

在本节课中，我们将学习如何编写自定义的CUDA内核来实现特定算法。当现有的库（如Thrust、cuBLAS）无法满足需求时，这是必要的技能。我们将从理解CUDA内核的基础开始，逐步深入到原子操作、线程同步、共享内存以及协作算法等高级概念。

执行空间回顾

上一节我们探讨了异步并行算法如何提升性能。但如果你的特定用例没有现成的算法呢？这时，你需要使用CUDA内核编写自己的算法。

首先，我们回顾一下执行空间。到目前为止，我们讨论过：

__host__：指定函数在主机（CPU）上执行，编译器为其生成CPU指令。
__device__：指定函数在设备（GPU）上执行，编译器为其生成GPU指令。

还有一个我们尚未讨论的关键字：

__global__：与__device__类似，编译器会为函数生成GPU指令。但标记为__global__的函数可以从CPU调用，并在GPU上执行。这就是我们所说的CUDA内核。

从CPU调用CUDA内核的语法看起来有些特殊，使用了三重尖括号<<< >>>，我们稍后会讨论其中的参数含义。

另一个重要细节是，内核启动是异步的。CPU启动内核后，会继续执行后续代码，而不会等待GPU上的内核完成。

编写第一个CUDA内核

在开始编写自己的CUDA内核之前，我们先回顾一下模拟器当前的工作原理。我们使用mdspan获得温度网格的二维视图，然后使用thrust::device_transform将单元索引转换为新温度。thrust::device_transform内部使用CUDA内核实现，因此是异步运行的。

现在，让我们尝试重写模拟器，直接调用CUDA内核，而不是依赖Thrust。

以下是一个简单的CUDA内核示例，它遍历所有单元并将索引转换为新温度。

__global__ void naive_kernel(mdspan2d grid) {
    for (int i = 0; i < grid.extent(0); ++i) {
        for (int j = 0; j < grid.extent(1); ++j) {
            // 计算新温度...
        }
    }
}
// 启动内核
naive_kernel<<<1, 1, 0, stream>>>(grid_view);

注意，当我们启动内核时，使用了三重尖括号语法。最后一个参数是CUDA流，通过传递一个流，你可以在指定的流上异步运行计算。

GPU不会自动并行化你的代码。在这个例子中，单个GPU线程正在串行计算每个单元。总体而言，这使得我们的内核比使用cuBLAS慢大约10,000倍，因为我们只启动了一个线程来完成所有计算。

增加线程数量

我们可以轻松增加CUDA使用的线程数量。例如，要启动两个线程而不是一个，只需在第二个参数中传递2。在三重尖括号中，第二个参数代表你想要启动的线程数量。

使用两个线程几乎可以使内核速度翻倍，这是合理的，因为我们使用了两倍的线程。但仅仅启动更多线程是不够的。如果它们都在处理相同的单元，实际上不会有任何加速。

为了解决这个问题，我们需要修改内核，使每个线程能够区分自己并处理不同的单元子集。这时就需要用到threadIdx.x。

threadIdx.x是一个内置变量，在CUDA内核内部可用，它保存当前线程的索引。如果我们启动两个线程，第一个线程的threadIdx.x为0，第二个线程的为1。

使用这种方法，每个线程处理不同的单元子集。在内核中，我们首先读取threadIdx.x来识别当前线程索引，然后每个线程开始处理与其线程索引匹配的单元。

这意味着第一个线程处理单元0，第二个线程处理单元1。处理完一个单元后，我们通过将单元索引增加线程总数来前进。这样，第一个线程处理单元0、2、4...，第二个线程处理单元1、3、5...

这种方案允许我们并行计算所有单元，确保一个单元永远不会被多个线程处理。

另一个观察是，我们的代码中没有限制只能使用两个线程。我们可以使用更多线程。让我们尝试增加处理此计算内核的线程数量。

从2个线程增加到256个线程，性能显著提升，因为我们使用了更多线程。但可以看到，我们仍然远远达不到cuBLAS所能达到的性能。

线程块与网格

为什么不能直接添加更多线程呢？如果我们尝试用2048个线程启动内核，会得到一个神秘的错误：“invalid configure argument”。这是因为单个线程块中的线程数量是有限制的。

在CUDA中，线程被分组为线程块，所有线程块的集合称为网格。每个线程块可以包含32、64、128直到1024个线程，不能超过这个数量。

当你启动一个CUDA内核时，所有线程块将拥有完全相同数量的线程。当你使用threadIdx.x时，你得到的是线程在本地块内的索引。如果线程0在块0中，你会得到0；如果线程0在块1中，你也会得到0，因为threadIdx.x返回的是块内的索引，而不是跨块的索引。

你还可以使用blockDim.x获取块的宽度（即每个块的线程数），使用gridDim.x获取网格的长度（即线程块的数量），以及使用blockIdx.x获取当前块的ID。

通过组合这些信息，我们可以计算线程在整个网格中的全局索引。blockDim.x告诉我们每个块有多少线程，gridDim.x告诉我们有多少个块。将两者相乘，可以得到线程总数。

为了计算线程的全局索引，我们将块ID（blockIdx.x）乘以块的长度（blockDim.x），得到到达正确块的偏移量，然后加上本地线程ID（threadIdx.x）。

例如，最后一个块的第二个线程：块ID为2，乘以块长度2，得到偏移量4，然后加上线程ID 1，得到全局索引5。

选择块大小和网格大小

现在我们知道如何找到每个线程在网格中的唯一索引，也知道了线程总数。剩下的一个问题是：我应该使用什么样的块大小？不幸的是，没有一个适用于所有情况的神奇数字。通常，你必须根据问题微调这个数字。

一个有用的指导原则是始终使用32的倍数。例如，选择50或60没有太大意义。通常我们建议的默认值是256，但这取决于你的问题，你应该查看哪个值最好。

第二个问题是我们应该启动多少个块？网格大小通常取决于问题大小，因为我们希望最大化并行性。问题是如何将问题的元素分配给线程。

例如，如果我们希望每个线程只处理一个元素，使用简单的整数除法是行不通的。假设我们有6个单元要处理，选择块大小为4个线程。6除以4得到1，我们只会使用一个线程块。但这样，用一个4线程的块处理6个单元，每个线程处理一个元素，我们会遗漏最后两个元素。

在CUDA中，我们提供了cuda::std::ceil_div函数，它会向上取整。你可以使用它自动计算当你希望每个线程恰好处理一个元素时所需的正确线程块数量。

在我们的例子中，6和4，调用ceil_div将返回2而不是1，因此你将正确分配两个线程块，每个块有4个线程来处理所有6个单元。

要传递给内核你希望使用的每个块的线程数，它是第一个参数。你可以先声明块大小（例如256），然后计算需要多少个块。再次使用cuda::std::ceil_div，用总问题大小除以每个块的大小，得到所需的总块数。然后，在调用内核时，首先传递网格大小（即所需的线程块数量），然后是每个块的大小。

通过这些更改，我们现在启动了大约500万个线程，终于接近了cuBLAS的性能。当然，cuBLAS在底层有更高级的优化，你可能无法匹配其速度，但对于这样一个简单的内核，仅通过添加更多线程并正确计算ID和线程局部性就能实现这样的加速，已经非常不错了。

边界检查与调试

我们已经修改了模拟器，但尚未确认它是否仍然正确工作。因为我们的设置是对称的，网格顶部和底部的温度是相同的，所以我们可以通过检查对称性来进行快速冒烟测试。

想法很简单：对于给定的单元，我们通过从网格高度减去行索引来找到其镜像行。如果一切正确，对称单元的温度应该几乎相同。

在屏幕上，你会看到一个执行此检查的C++函数。作为本练习的一部分，你需要将此CPU C++函数转换为CUDA内核，并使用三重尖括号语法启动它，选择正确的块和线程数量。

要将对称检查函数转换为CUDA内核，我们必须用__global__说明符注释它。__global__允许你在GPU上启动代码。这个函数将被编译以在GPU上工作，并且可以从CPU以异步方式调用。

之后，我们需要使用我们讨论过的三重尖括号语法启动它。这里因为没有真正的并行工作，我们可以只启动一个包含单个线程的块。为了确保检查在正确的流上运行，我们只需指定最后一个参数为我们想要使用的流。

现在思考一下，整个行应该是对称的，而不仅仅是一组。在本练习中，你将修改内核，使每个线程恰好处理一列。

就像之前一样，如果我们使用块和线程索引，我们可以计算线程的全局ID。我们在顶部这样做是因为我们想知道我们在哪一列：我们使用块ID（我们的块索引）乘以块的宽度，然后加上我们在块内的当前线程ID。通过这样做，我们可以计算列索引。

然后，当我们想要调用实际的内核时，就像之前一样，我们可以声明块宽度（我们使用默认值256），计算网格大小，检查问题的总大小（在我们的例子中是宽度），然后使用cuda::std::ceil_div除以块大小。

当我们调用内核时，我们只需传递网格参数和块参数。这样，每个线程恰好处理一列，并行检查整个行。

但这样做，我们刚刚导致了第一次越界访问。因为我们向上取整以确保每个线程至少有一个元素，我们可能会启动比实际问题大小更多的线程。在这种情况下，我们启动了两个线程块，每个块有4个线程，但我们总共只有6个元素。第一个块将用其4个线程处理前4个元素，第二个块也有4个线程，但只有两个元素要处理，所以最后两个线程将尝试访问实际上不存在的元素。

这两个线程发出的地址将越界。我们需要找到一种方法来修复这个问题。一个简单的方法是在内核中添加所谓的边界检查。在计算全局线程ID之后，我们只需要将此值与总列数进行比较，如果没有工作要做，就跳过它。

在我们计算了列ID之后，我们可以检查列ID是否小于宽度（我们可以通过mdspan的extent(1)找到），只有当我的列ID小于这个宽度时，我才进入这里进行检查。

如果我们不知道这些边界检查，我们如何自己发现这个问题呢？在NVIDIA，我们提供了一个名为Compute Sanitizer的特殊工具，帮助检测CUDA内核中的越界访问和其他内存问题。

如果你使用幻灯片上显示的标志编译此模拟器代码，Compute Sanitizer将准确指出非法访问发生在哪一行。例如，你可以看到我们拥有所有需要的信息：我们知道它进行了无效的全局读取（大小4字节，因为我们使用浮点数），我们知道它来自这个特定的内核symmetry_check_kernel，我们甚至知道是第7行。我们知道它来自哪个线程和哪个块，我们甚至知道最近的有效分配是什么。

通过这种方式，Compute Sanitizer确实为我们提供了理解问题所需的所有信息：我们缺少的确实是一个边界检查条件。

实际上，mdspan也可以自行检测越界访问。默认情况下它是关闭的，但如果你在编译时定义了_CCCL_ENABLE_ASSERTIONS，就像我们在顶部展示的那样，这将使mdspan在访问时检测你是否越界。如果发生越界，它将触发一个断言。同样，你将知道是哪个线程在哪个块中进行了越界访问。这是通过检查当我们通过mdspan访问数组时，我们尝试访问的索引是否大于mdspan的大小来实现的。

通过添加这些边界检查，我们刚刚修复了这个错误。

直方图计算与数据竞争

现在很自然地会问：为什么我们需要线程层次结构的所有这些复杂性，比如线程、线程块、网格？为了理解为什么需要所有这些，我们首先需要稍微修改一下问题，看看为什么线程块层次结构实际上很有意义。

现在让我们为模拟器产生的温度创建一个直方图。直方图帮助我们查看某些温度范围出现的频率。

首先，我们将整个温度范围划分为不同的区间。在我们的例子中，我们将有10个不同的区间。对于模拟中的每个单元，我们确定其温度将落入哪个区间。将温度划分到特定区间的简单方法是向下取整。

例如，如果单元温度为4度，我们除以10得到0，因此我们知道这个4度应该进入区间0。如果单元温度为15度，我们想将15放入第一个区间，同样只需将15除以10得到1，因此我们知道这个15度必须放入区间1。

一旦我们将每个单元分配到正确的区间，我们只需计算有多少个单元落入每个区间，这就可以成为直方图中条形的高度。

现在让我们思考构建我们自己的内核签名来完成这个直方图。因为它是直方图，只需要一个维度，所以我们不需要mdspan，因为mdspan是多维的跨度，这里我们只有一个维度。这就是为什么这里我们使用cuda::std::span而不是cuda::std::mdspan。

在CUDA中，就像在C++中一样，span是访问数据的首选方式，而不是仅仅使用原始指针，因为它更安全（可以检测越界），你也可以避免对指针进行奇怪的操作，还可以访问大小等。

要构造一个span，就像mdspan一样，只需传递你想要视为跨度的底层容器的数据和大小即可。然后，要访问它，就像任何常规容器（如向量）一样，只需使用方括号来访问。

将温度和直方图传递给内核后，我们再次需要计算我们的全局线程ID，使用与之前相同的公式：块ID乘以块维度加上线程ID。然后每个线程加载其对应的温度值，并将其除以区间宽度（在我们的例子中是10），以再次找到正确的区间索引，我们将基于温度值在该区间上加1。

在生产环境中，你还需要确保这些内存访问在边界内，但为了简单起见，在幻灯片上我们将跳过边界检查，并假设问题大小确实是块大小的倍数。

接下来，每个线程从其分配的区间读取当前计数，增加值，然后将其写回直方图。我们首先加载整个直方图值，然后递增。

例如，我发现我的温度是4，我想让它进入第一个区间，通过除以区间宽度的计算将返回0，因为4在0到10之间，所以我得到区间0。然后，我将通过直方图获取区间0的值（假设是13），我只想增加1，因为我有一个新值要添加到此区间0，然后我可以将结果写回内存。

但是，如果我这样做，我在直方图中看不到任何东西，所以看起来出了问题。我们在网格中大约有400万个单元，条形应该在百万范围内，但我们什么也看不到。发生了什么？

问题在于我们的内核在代码的高亮行中存在数据竞争。所有数百万个线程都在到处读取和写入相同的内存位置。让我们看看这里到底发生了什么。

假设我们只有两个线程。线程0想从区间0读取，线程1也想从区间0读取。它们都发出相同的内存读取请求。线程1首先发出对直方图的内存读取，然后内存子系统将响应存储在直方图区间中的当前值（在我们的例子中是0）。然后，线程1将1加到该变量并将结果存储在new中。但同时，当线程1这样做时，线程2也在向直方图发出读取。由于此加法尚未在内存中发生，线程2也从直方图读取值0。然后线程2也将1加到旧值并想再次存储它，但同时，线程1正在将值1写入直方图。所以这里1被写入直方图，现在直方图值是1。但问题是线程2也有值1，它也将开始将1写入直方图值。这个线程不知道另一个线程当前正在写入值1，因此应该使用这个1值来增加自己的计数器，但为时已晚，值已经被读取，所以它也会将1写入此直方图。

最终，直方图内部的值将是1，而两个不同的线程都试图向此直方图值加1。这就是两个线程时发生的情况。现在想象一下，当数百万个线程试图这样做时会发生什么。

这就是为什么我们的直方图看起来是空的，因为当我们递增时，没有考虑到其他线程也在尝试做同样的事情，这就是我们所说的数据竞争。

原子操作

为了修复这个竞争，我们需要确保“读取-修改-写入”这个序列被当作一个单一的、不可分割的内存操作来处理。我们不希望先读取，然后修改，再写入；我们希望整个操作在一个块中完成。

在C++中，我们可以使用原子操作来实现这一点。你可以将原子操作视为存储指令本身，而不仅仅是字节；我们真的希望整个操作一步完成。在右边的例子中，我们将“加一”存储到内存中，当操作到达内存系统时，它读取当前值，递增它，然后在单个步骤中将结果写回。

在左边，之前我们是先读取，然后更新局部变量，再写入。现在，因为我们执行的是fetch_add（在C++中也可以这样做），所有操作都在一个步骤中完成。通过这样做，我们防止了任何线程同时读取或修改，每个线程将发出一个原子操作，而整个“读取-修改-写入”将在一个序列中发生。

幸运的是，CUDA也通过cuda::std::atomic_ref提供原子操作，它允许你将任何现有的内存位置视为原子变量。你获取任何当前可用的内存位置，可以将其视为原子引用。

在这个例子中，我们通过将count0包装到atomic_ref中来创建对count0的原子引用。这意味着，然后我们可以通过这个引用（它是一个原子引用）对count0执行原子操作。我们现在可以执行原子加、原子减、原子与等操作，每次我们这样做时，这些操作都将是原子的。它们都将在单个操作中执行，不会有先读取、再修改、后写入的情况。

这将允许我们并行修改多个线程的计数，而没有任何竞争条件。让我们看看如何将其放入我们自己的内核中。

在我们的内核中，现在会发生的是：线程1将在原子引用上使用fetch_add向直方图加1，线程2也将想在直方图上执行fetch_add加1。

现在可能发生的情况是：线程1执行fetch_add，这将到达内存，执行读取、修改和写入，所有操作都在一个序列中完成，结果将存储回内存。所以现在当线程2也想执行其操作时，当它最终到达内存时，它将正确地读取新值（即1），因为第二个fetch_add不能在第一个fetch_add之前执行。当这个fetch_add到达时，第一个已经完成，这就是为什么我们可以正确地看到1，然后正确地加1，并正确地将最终的直方图值（现在是2）更新回内存。现在我们知道结果是正确的，并且在此过程中没有丢失任何增量。

在本练习中，你需要修改我们现有的直方图代码，不再使用三行代码进行“读取-修改-写入”（这会导致数据竞争），而是使用我们讨论过的atomic_ref来使代码工作。

修改实际上只有几行代码。我们需要做的就像我们在示例中展示的那样，将我们要写入的直方图区间包装到atomic_ref中。现在，当我们通过这个原子引用写入时，我们正在写入直方图区间，但我们是原子地写入，所以我们确信当我们在这里发出fetch_add时，这个fetch_add操作是一个原子操作，因此它是一个执行所有“读取-修改-写入”的单一序列，而不是我们自己操作（这会导致数据竞争）。

修复之后，我们的结果终于正确了，我们可以看到有意义的直方图。注意初始温度接近零，意味着第一个区间拥有最多的单元，随着热量在网格中扩散，我们看到更多区间开始填充。

性能问题与私有化

但现在我们还有另一个问题：性能非常差，我们只达到了6 GB/s，这远远低于GPU的最大内存带宽。我们真的想要比6 GB/s更好的性能。

上次我们看到如此差的性能是因为序列化，当时我们只使用了几个线程，而没有使用所有线程。但这次我们的代码中没有任何循环，我们使用了所有线程，并行地使用线程和线程块，这是怎么回事？

但这次，序列化来自原子操作本身。是的，它们解决了数据竞争问题，但由于设计原因，它们不能并行运行。所有线程都针对相同的内存位置，现在我们有许多许多原子操作一个接一个地排队。当我们启动大约256个线程时，总共有大约16,000个块，我们发出了400万个原子操作，所有这些操作都在排队。

这显然太多了，我们不想排队那么多原子操作，因为这会使代码非常慢。我们有什么选择可以解决这个问题？

一种解决方案是所谓的私有化。我们不是让所有线程更新同一个全局直方图，并首先在所有区间上发出太多原子操作，而是可以为每个线程块分配一个小的私有直方图。

每个线程块将有一个私有直方图，可以在其中安全地使用原子操作。这不会减少原子操作的总数，但现在它分布在更多的内存位置上。不是所有线程都写入相同的内存位置，现在线程块首先写入自己的内存位置，然后再传播回全局直方图。完成后，块聚合其私有直方图的计数，并使用一个内存原子操作将它们更新到主直方图。

在我们简单的两个块的例子中，每个块递增自己的本地块直方图，然后才将结果传播回完整的最终直方图。虽然在这个例子中，显然不是很有帮助，因为我们只有两个块，但想象一下真实情况，我们有16,000个块。在这个场景中，我们将有16,000个私有直方图并行更新，随后只有16,000个原子操作。之前我们在全局内存中进行了400万个原子操作，它们都在排队；现在因为每个线程块只做一个，我们在全局内存中只有16,000个原子操作。

为了实现这种私有化，我们需要为内核添加一个新参数，我们称之为块直方图。接下来，我们为当前线程块的私有化直方图部分创建一个子跨度。subspan函数接受一个偏移量和一个大小，我们通过将块索引乘以直方图大小来计算偏移量。第一个块的偏移量为0，第二个块的偏移量为10，等等。这确保每个块写入自己的内存切片，每个块将拥有自己的私有化直方图切片。

然后就像之前一样，每个线程计算其区间索引，但是，我们现在递增私有化的块直方图，而不是全局直方图，因为我们想递增每个块的私有化直方图，而不是全局的。之后，每个块将其部分计数从其私有化直方图添加回全局直方图内部。在CUDA代码片段中，块中的前10个线程每个处理一个区间，用于块的私有直方图，并并行化对全局直方图的写入，它们再次使用原子引用来安全地将其本地计数添加到最终直方图。

但不幸的是，这样做我们引入了一个新的错误。

线程同步

让我们仔细看看发生了什么以及为什么这里有一个错误。在代码中，我们假设所有线程在下一段代码开始之前都完成了对本地块直方图的更新。我们假设所有线程在写入全局块直方图之前都写入了本地块直方图。

然而，CUDA线程在一个块中并发运行，并且没有任何顺序保证。这意味着你块中的一些线程可能甚至还没有开始，而其他线程可能已经开始退出你的内核。所以你无法保证线程块中所有线程的调度顺序，一个线程可能处于代码的最开始，而一些线程可能处于代码的最后。

一种可能的情况是，一些线程可能在所有块中的线程完成更新之前读取了块直方图，然后才写入全局直方图。

为了解决这个问题，我们必须以某种方式确保所有线程在读取之前都完成了对本地块直方图的更新。为此，CUDA提供了__syncthreads()。

__syncthreads()是一个特殊的CUDA函数，它充当块中所有线程的屏障。当一个线程到达此屏障时，它会等待，直到线程块中的所有其他线程也到达此屏障。

在这个意义上，它类似于标准C++库中的std::barrier，但这里有一个重要的区别：块中的每个线程必须从相同的逻辑路径调用__syncthreads()。

__syncthreads()将解决我们的功能问题，但仍然存在性能问题。

线程作用域

如果我们回顾CUDA概述，你会注意到两组不同的原子引用：一组在cuda::std命名空间中，另一组在cuda命名空间中。这自然引出一个问题：两者之间有什么区别？

cuda::std::atomic_ref和cuda::atomic_ref提供了基本相同的接口。然而，cuda::atomic_ref通过添加一个thread_scope参数扩展了cuda::std::atomic_ref。选择正确的线程作用域可以显著影响性能。

让我们试着更好地理解什么是线程作用域。线程作用域表示可以使用给定原子进行同步的线程集合。作用域可以是整个系统、仅设备，甚至只是线程块。

例如，当所有线程彼此相关时，对于原子操作，你会使用cuda::thread_scope_system。这意味着来自任何GPU的线程都可以与GPU中的任何其他线程甚至CPU线程同步，它是整个系统。

cuda::std::atomic_ref实际上与没有std的cuda::atomic_ref相同，但作用域是系统。下一个作用域，当它不是整个系统时，是设备。device作用域意味着单个GPU内的任何线程。使用此作用域，多个GPU无法向同一内存位置发出原子操作，GPU和CPU也不能写入同一内存位置，但GPU内的任何线程都可以自动写入一个位置，而同一GPU的任何其他线程也正在写入此相同内存位置。

最后一个线程作用域是thread_scope_block，它将原子操作限制在单个块内的线程。只要它们在同一线程块内，多个线程可以对同一内存位置执行原子操作，但如果线程来自另一个线程块，则不行。

在我们的例子中，我们可以看到，因为每个线程块都在更新自己的本地直方图，这非常适合使用thread_scope_block。

这引出了下一个练习。这次你必须使用线程块同步来确保所有线程在读取和写入全局直方图之前都完成了对块直方图的更新。如果时间允许，你还可以通过添加线程作用域来提高原子操作的效率。

首先，我们将cuda::std::atomic_ref替换为cuda::atomic_ref。然后，我们将块直方图上的原子操作作用域限制为thread_scope_block，因为我们知道，当我们写入这个对每个块私有的直方图时，我们只由同一线程块中的不同线程执行原子操作。这就是为什么我们可以使用cuda::thread_scope_block。

然后我们还添加了__syncthreads()，因为在完成直方图更新后，我们想确保所有线程在更新全局直方图之前都完成了对此块直方图的更新。

同样，在更新此全局直方图时，我们可以使用scope_device而不是默认的系统原子，因为我们知道我们只从设备写入直方图。但我们必须使用device而不是block，因为我们知道这个直方图是全局的，由所有块中的所有线程共享。

仅通过这些更改，我们就从6 GB/s一路提升到100 GB/s，这是一个巨大的改进，仅仅通过改变我们执行的原子操作的作用域。

但我们还可以做得更好。

共享内存

在进一步优化直方图内核之前，让我们回顾一下到目前为止学到的东西。我们刚刚介绍了原子操作、线程作用域和同步。但我们最初的问题呢？最初的问题是：为什么我们首先需要线程层次结构，为什么需要线程块？让我们再看看线程块同步和线程块作用域的原子操作，它们仅在给定线程块的线程之间可用。

跨线程组同步要昂贵得多。线程层次结构在这里有好处。仅同步块内的线程远比同步所有块内的所有线程成本低。

但是，在某些线程组内是否有其他可访问的设施？为了找出答案，我们来谈谈GPU架构。

GPU由统一的构建块组成，我们称之为SM（流式多处理器），这就是你在屏幕上看到的绿色部分。单个GPU可以包含数百个SM。

每个SM有很多核心，但也有一个本地L1缓存。全局内存，即主要的GPU内存，位于这些SM之外，这就是为什么它被称为全局内存。访问全局内存的延迟比访问L1缓存中的内容要大得多，因为L1缓存是SM本地的，而且访问L2也比访问全局内存快。

因为GPU与CPU相比，通常延迟不高但带宽很好，这就是为什么我们希望通过利用L2和L1缓存来使延迟尽可能好。

像线程块同步这样的功能是专门构建在SM硬件中的。因为当我们进行线程同步时，只是针对一个线程块，我们不是同步整个网格，所以它只需要在SM级别运行，这就是为什么它更快。

现在让我们思考一下我们的直方图内核是如何映射到硬件上的。每个线程块在单个SM上运行，块直方图存储在全局内存中。然后每个线程向我们的私有块直方图发出原子操作，这很可能由L2缓存支持。这效率不高，因为我们实际上不需要在内核完成后保留此块直方图，它们只是临时存储，而L1缓存正好位于SM上，目前我们没有使用它。

如果我们可以直接在L1缓存中分配此块直方图，那将对我们非常高效，首先不会浪费全局内存，而且速度会更快，因为我们将直接从L1缓存读取和写入。

幸运的是，我们可以在CUDA中做到这一点。CUDA提供了一个显式的、可编程的线程块本地内存，我们称之为共享内存。

共享内存位于L1缓存所在的位置，它具有与L1缓存相似的延迟和带宽。共享内存的唯一缺点是它有限制，你不能拥有无限量的共享内存，因为每个SM的L1缓存数量有限。

如果你想分配共享内存，只需将__shared__说明符添加到变量中。在这个例子中，我们只是声明了一个大小为4的数组，即共享内存中的4个整数。

但不幸的是，共享内存中的变量不能有非平凡的构造函数，只能是像int、float等平凡类型。

在线程块内，你可以像任何其他数组一样读取和写入这些数组。在这个例子中，每个线程会将其自己的线程索引写入数组的当前扩展元素中。线程0将0写入位置0，线程1将1写入位置1，等等。

为了不忘记任何数据竞争，我们需要在完成对此共享内存的任何写入之后、想要进行任何读取之前，不要忘记调用__syncthreads()。

为了强调这个共享数组确实在线程块内共享，你可以看到这里我让第一个线程打印出存储的每个值，所以所有线程在此时都会看到所有其他线程的更新，并且它将打印值0、1、2、3，因为我们使用了__syncthreads()。

__syncthreads()不仅强制我们这里的四个线程完成，而且还使更改对其他线程可见。这就是为什么在我们完成三次写入并进行同步之后，线程0能够看到，并且我们知道它能够看到，因为当我们打印结果时，我们正确地看到了0、1、2、3。这是因为__syncthreads()允许线程块中的所有线程看到同一线程块中其他线程所做的所有修改。

现在让我们进入下一个练习。这次，你必须将块直方图分配在共享内存中，而不是全局内存中，这应该有望提高我们代码的性能。

首先，我们需要更改函数签名，因为我们不再需要块直方图跨度，现在我们不再使用全局内存。为了分配这个共享内存，我们想要的是一个块直方图，它将再次对每个线程块是局部的。

这就是为什么我们使用区间数量（这里是num_bins）的大小来分配它。然后我们想确保因为它是一个直方图，但初始化为零，因为然后我们将逐步向其中加一，所以它需要从零开始。这就是为什么我们使用我们的线程将前num_bins个元素分配为0。然后我们不应该忘记使用__syncthreads()，因为我们刚刚写入共享内存，所以在真正读取内存之前，我们需要同步。

最后，我们再次将这个共享块直方图包装到原子引用中，就像之前一样。然后其余的代码可以保持不变，我们做相同的fetch_add，相同的__syncthreads()，相同的一切，只是现在我们在共享内存中而不是全局内存中做。

通过进行这个小小的更改，我们再次从100 GB/s一路提升到400 GB/s，即4倍加速。现在，我们的内核功能完善且性能良好。

协作算法

问题是，我们在这里犯了一个更高层次的错误。当我们第一次开始使用GPU时，我们依赖加速库，但现在因为我们编写自己的内核，我们试图从头开始重新实现一切，而我们本可以使用库。

编写CUDA内核并不意味着你必须重新发明一切。是的，也许你的特定用例非常特殊，所以你无法使用cuBLAS、Thrust或cuFFT，因此你需要一个CUDA内核。但即使在你自己的CUDA内核内部，我们仍然有一些工具可以提供给你。

CUDA提供了一系列我们称之为协作库的工具，可以加速你的内核并缩短开发时间。

例如，除了从主机启动的并行算法（如Thrust和cuBLAS）之外，我们还有设备端的版本。我们有像cuBLAS DX和cuFFT DX这样的东西，可以直接从设备启动BLAS线性代数操作或傅里叶变换操作。

cuRand也是如此。但“协作”在这种情况下意味着什么，为什么我们称之为协作库？

我们说协作是因为我们在CUDA中有不同类型的算法。让我们以排序为例。当我们谈论串行算法时，它由单个线程调用和执行。如果两个不同的线程各自调用串行算法，一个线程的输入不会影响另一个线程的输出。

另一方面，协作算法由多个线程调用，并由所有线程同时集体执行。你可以将此协作算法的输入想象为一个恰好被分割到不同线程的大型虚拟数组，输出也是如此。

在幻灯片上的例子中，协作排序从第一个线程获取元素1、3，从第二个线程获取2、4，排序后，第一个线程的输入变为1、2，第二个线程变为3、4，即使这两个值最初来自线程1。这就是为什么我们说它是协作的，输入和输出由所有线程共享。

我们还使用了并行算法，它们由单个线程调用，但在幕后由许多线程执行。从用户的角度来看，这些算法看起来是串行的，因为你只是像调用任何其他函数一样调用它们。然而，在幕后，它们分布在多个线程上以使任务更快。

到目前为止我们使用过的一些并行算法的例子有thrust::transform或cuda::device::transform。

现在让我们仔细看看这些协作算法。在幻灯片上，你可以看到协作归约可能如何工作。这里，所有线程都向协作归约提供一个值。线程1贡献值0，线程2贡献值1，等等。

在协作算法内部，这些值被写入共享内存。在线程块同步之后，一组较小的线程计算部分和并将其写回，再进行一次同步后，总和我们例子中的6被写回第一个线程。

虽然这是一个非常基本的模型，但它说明了协作算法的两个关键点：它们依赖共享内存在线程之间交换数据，并且可能包含同步（这里我们看到有两个__syncthreads()使归约工作）。

这意味着块中的任何线程如果未能调用协作算法，整个内核可能会出现问题。我们需要所有线程在调用此协作算法时都被包含在内。

现在让我们看看这种直觉如何应用于实际代码。在屏幕上，你可以看到由cuRand提供的协作归约。与传统的面向函数的接口不同，cuRand将其协作算法公开为模板化结构。

模板参数用于针对手头的问题专门化算法。例如，这里我们可以指定我们想要归约特定类型（这里是整数），并且我们需要指定我们的块宽度将是256个线程。

协作算法有一个嵌套的临时存储标签，它指定每种类型以及协作算法进行内部通信所需的临时存储量。

我们在共享内存中分配此类型的实例，然后在构造协作算法时传递对该共享内存实例的引用。最后，成员函数提供特定协作算法的不同变体。

如果我们把所有这些东西放在一起，协作算法的使用看起来像这样。我们首先通过声明我们想要归约整数并且线程块大小为4来实例化块归约。然后我们在共享内存中分配临时存储。接下来，我们构造此协作算法结构的一个实例。最后，当此对象可用时，我们最终可以从块中的所有线程调用协作算法。

cuRand提供了许多线程块级别的通用算法，你可以看到我们有交换、直方图、加载/存储、归约、扫描、洗牌等。我们可以看到，我们这里需要的是块直方图，而不是自己实现块直方图并可能产生未高度优化的代码，我们可以使用cuRand块直方图。

这个块直方图算法有一些模板参数。首先，我们必须指定直方图的类型，在我们的例子中是整数。然后我们需要指定线程块大小，在我们的例子中是256。然后我们需要指定每个线程将贡献多少个区间。最后，我们必须指定我们的直方图中有多少个区间。

此协作原语的使用类似于我们在归约案例中看到的。这引出了本节的最终练习。这次，我们必须在直方图内核中使用cuRand块直方图，而不是自己实现它。

我们首先从网格加载单元，就像之前一样。这次，我们在栈上本地分配区间，使用C数组而不是跨度。然后我们实例化块直方图类型，在共享内存中分配其临时存储。然后我们可以调用直方图方法来使直方图发生。最后，我们不应该忘记，在之后进行任何内存读写之前，我们仍然需要同步。

总结

在本节课中，我们一起学习了如何编写自定义CUDA内核来实现特定算法。我们涵盖了以下关键概念：

执行空间：回顾了__host__、__device__和__global__函数说明符，理解了CUDA内核的启动机制和异步特性。
线程层次结构：深入理解了线程、线程块和网格的概念，学会了使用threadIdx、blockIdx、blockDim和gridDim计算线程的全局索引。
内核编写与优化：从编写一个简单的串行内核开始，逐步通过增加线程、正确计算索引来并行化，并学习了如何选择块大小和网格大小。
边界检查与调试：认识到越界访问的风险，学会了在内核中添加边界检查，并了解了使用Compute Sanitizer和

posted @ 2026-03-29 09:47 布客飞龙III 阅读(65) 评论(0) 收藏举报

刷新页面返回顶部

龙哥盟