贝鲁特-GPU-编程笔记-全-

贝鲁特 GPU 编程笔记（全）

GPU计算：10：归约（Reduction） 🧮

概述

在本节课中，我们将要学习一种新的并行模式——归约（Reduction）。归约是一种将一组输入值合并为单个输出值的操作，例如求和、求积、求最大值或最小值。我们将以求和为例，探讨如何在GPU上高效地实现归约操作。

回顾：模板计算（Stencil）

上一节我们介绍了模板计算，这是一种在网格上计算某点值，该值依赖于该点及其邻居点的计算模式。我们实现了3D模板计算，并首次使用了三维线程网格。

以下是模板计算的核心概念：

数据重用：在计算输出元素时，多个线程会重复使用相同的输入数据。
共享内存优化：通过将输入数据块加载到共享内存中，减少对全局内存的重复访问，从而提高计算与内存访问的比率。
寄存器平铺（Register Tiling）：我们进一步优化，将某些仅被单个线程使用的数据平面存储在寄存器中，而非共享内存，以减少共享内存的使用并提升性能。

什么是归约？ 🤔

归约是一种将一组输入值通过一个关联且可交换的操作，合并为单个输出值的操作。常见的归约操作包括求和、求积、求最小值和最大值。

每个归约操作都有一个定义明确的单位元：

求和的单位元是 0。
求积的单位元是 1。
求最小值的单位元是 INT_MAX（或 FLT_MAX）。
求最大值的单位元是 INT_MIN（或 FLT_MIN）。

串行归约

首先，我们看看如何实现串行归约。以求和为例，伪代码如下：

float sum = 0; // 初始化为单位元
for (int i = 0; i < n; i++) {
    sum += input[i]; // 应用归约操作
}

如你所见，归约本质上是顺序的，循环迭代之间存在依赖关系，这给并行化带来了挑战。

并行归约：归约树 🌳

为了并行化归约，我们引入归约树的概念。其核心思想是将数据分成对，让多个线程并行计算部分和，然后递归地合并这些部分和。

以下是归约树的步骤说明：

第一步：假设有8个元素，我们创建4个线程，每个线程并行计算两个相邻元素的和。结果得到4个部分和。
第二步：使用2个线程，并行计算上一步得到的4个部分和中的两对。结果得到2个部分和。
第三步：使用1个线程，计算最后两个部分和，得到最终的总和。

对于一个包含 n 个元素的数组，归约树需要 log₂(n) 步完成。每一步之后，活跃的线程数量减半。

GPU上的归约挑战与策略

在GPU上实现归约树面临一个关键限制：线程只能在同一个线程块内进行同步，不同线程块间的线程无法同步。

因此，我们采用分段归约策略：

将整个输入数组分割成多个段。
每个线程块独立地对其负责的数据段执行归约树计算（因为块内线程可以同步），产生一个部分和。
所有线程块将其部分和写入一个全局数组。
最后，我们需要将这些部分和合并为最终结果。合并方法有多种：
- 启动一个新的归约内核来处理部分和数组。
- 使用原子操作（Atomic Operations）将部分和累加到一个全局变量中。
- 如果部分和数量很少，直接在CPU上相加。

本节课，我们将重点探讨如何在单个线程块内实现高效的归约树。

基础归约树实现（低效版本）

我们首先实现一个直观但低效的版本。在这个版本中，每个线程负责处理数组中相隔较远的元素。

核心代码逻辑如下：

// stride 从1开始，每次翻倍，直到等于blockDim.x
for (int stride = 1; stride <= blockDim.x; stride *= 2) {
    // 只有线程ID是stride的倍数的线程参与计算
    if (threadIdx.x % stride == 0) {
        input[global_idx] += input[global_idx + stride];
    }
    __syncthreads(); // 每一步都需要同步
}

这个实现存在几个明显问题：

内存访问未合并：线程访问的内存地址不连续，无法利用GPU的合并内存访问特性，导致带宽利用率低。
控制流发散（Control Divergence）：随着步长增加，越来越多的线程不参与计算，但由于它们属于同一个Warp（GPU执行单元），这些空闲线程仍然会消耗执行资源，拖慢整个Warp的速度。
线程利用率低：计算后期，大部分线程处于闲置状态。

优化归约树实现（高效版本）

为了解决上述问题，我们重新组织线程的工作方式。让线程块中的线程连续地处理输入数据段的开头部分。

优化后的核心代码逻辑如下：

// 每个线程加载相邻的两个元素到共享内存（或直接相加）
unsigned int i = segment_start + threadIdx.x;
input_shared[threadIdx.x] = input[i] + input[i + blockDim.x];
__syncthreads();

// stride 从blockDim.x/2开始，每次减半，直到1
for (int stride = blockDim.x / 2; stride >= 1; stride >>= 1) {
    // 只有ID小于stride的线程参与计算
    if (threadIdx.x < stride) {
        input_shared[threadIdx.x] += input_shared[threadIdx.x + stride];
    }
    __syncthreads();
}

这个优化带来了以下好处：

内存访问合并：在初始加载和后续每一步中，活跃线程访问的共享内存地址都是连续的，符合合并访问条件。
减少控制流发散：在非最终步骤中，整个Warp要么全部活跃，要么全部闲置，避免了Warp内部的分支发散。
更高的线程利用率：在计算初期充分利用了所有线程。

进一步优化：使用共享内存

在基础版本中，我们直接修改全局内存中的输入数组。一个更好的做法是使用共享内存作为中间存储。

使用共享内存的优点：

避免污染输入：原始输入数据保持不变。
提升速度：共享内存的访问速度远快于全局内存。
减少全局内存访问：只有在初始加载和最终存储时才访问全局内存。

实现方法：

每个线程将其负责的初始元素对加载并相加，结果存入共享内存数组。
后续的归约树步骤完全在共享内存数组上进行。
最后，由线程0将共享内存中的最终结果写回全局内存。

高级优化：线程粗化（Thread Coarsening） ⚙️

归约的并行化成本主要来自同步和控制流发散。线程粗化 是一种通过让每个线程处理更多数据来摊销这些成本的技术。

其核心思想是：与其让硬件在资源不足时串行执行多个线程块（每个块都产生同步和发散开销），不如我们主动让一个线程块处理更大的数据段。

实现方式：

定义一个粗化因子（例如4）。
每个线程不再只处理1对元素，而是循环处理 粗化因子 * 2 个元素，并在寄存器中累加一个局部和。
这个循环阶段没有同步，也没有控制流发散，因为所有线程执行相同次数的循环。
循环结束后，每个线程将其局部和写入共享内存。
然后再执行共享内存中的归约树（此时才有同步和发散）。

伪代码示意：

float thread_sum = 0;
// 每个线程处理 coarsening_factor * 2 个元素
for (int tile = 0; tile < coarsening_factor * 2; ++tile) {
    int load_idx = global_idx + tile * blockDim.x;
    if (load_idx < array_size) { // 处理边界条件
        thread_sum += input[load_idx];
    }
}
// 将局部和存入共享内存
shared_mem[threadIdx.x] = thread_sum;
__syncthreads();
// ... 后续进行共享内存上的归约树 ...

线程粗化的权衡：

优点：显著减少了同步次数和整体控制流发散，通常能提升性能。
缺点：牺牲了透明可扩展性。优化效果依赖于对硬件并行能力的估计（如每个流多处理器能同时驻留的线程块数量），可能需要针对特定硬件调整粗化因子。

边界条件处理

当输入数组大小不是线程块大小的整数倍时，最后一个线程块可能包含超出数组边界的元素。我们需要在从全局内存加载数据时进行边界检查。

简单的处理方法是：

float value = 0;
if (load_index < total_array_size) {
    value = input[load_index];
}
// 然后使用value进行计算

确保在初始加载阶段处理边界，后续在共享内存中的计算就无需再检查边界。

总结

本节课我们一起学习了GPU上的归约操作。我们从串行归约出发，引入了并行归约树的概念。我们探讨了在GPU上实现归约的挑战，特别是线程块间同步的限制，从而采用了分段归约的策略。

我们重点实现了单个线程块内的归约树，并逐步进行了优化：

基础实现：直观但存在内存访问未合并和控制流发散的问题。
优化实现：重组线程工作方式，实现了内存访问合并并减少了控制流发散。
共享内存优化：利用更快的共享内存作为工作区，减少全局内存访问。
线程粗化：通过让每个线程处理更多数据，来摊销同步和发散的开销，进一步提升性能。

最后，我们还讨论了边界条件的处理方法。归约是一种重要的并行模式，掌握其优化技巧对于编写高性能GPU程序至关重要。

GPU计算：第11讲：扫描（Kogge-Stone）🔍

概述

在本节课中，我们将要学习一种新的并行模式——扫描（Scan）。我们将重点介绍一种特定的并行扫描实现方法：Kogge-Stone算法。我们还将学习一种名为“双缓冲”的优化技术。最后，我们会讨论并行算法的工作效率概念。

回顾：归约（Reduction）

上一节我们介绍了归约操作。归约是一种将一组输入值通过某种运算符（如求和、求积、求最小值等）合并为单个值的操作。这个运算符需要满足结合律和交换律，并且需要一个单位元。

顺序归约（以求和为例）的伪代码如下：

sum = 0; // 单位元
for (int i = 0; i < n; i++) {
    sum = sum + input[i]; // 运算符
}

并行归约通常使用归约树来实现。对于一个包含 n 个元素的数组，我们可以在 log n 步内完成归约。然而，每一步都需要线程块内的全局同步。由于无法在不同线程块之间进行全局同步，我们通常采用分段归约的方法：每个线程块先进行本地归约，产生部分和，然后再对这些部分和进行归约。

我们还讨论了通过线程粗化来优化并行归约，以减少同步开销和步骤数。

扫描（Scan）操作介绍

本节中我们来看看扫描操作。扫描操作接收一个输入数组 X 和一个满足结合律的运算符，并返回一个输出数组 Y。扫描有两种形式：包含式扫描和排除式扫描。

包含式扫描：输出数组 Y 中的每个元素 Y[i] 是输入数组 X 中从 X[0] 到 X[i]（包含）所有元素的组合结果。
- 公式：Y[i] = X[0] op X[1] op ... op X[i]
排除式扫描：输出数组 Y 中的每个元素 Y[i] 是输入数组 X 中从 X[0] 到 X[i-1]（不包含 X[i]）所有元素的组合结果。Y[0] 是单位元。
- 公式：Y[0] = identity; Y[i] = X[0] op X[1] op ... op X[i-1] (对于 i > 0)

示例（运算符为加法）：

输入: X = [3, 6, 7, 4, ...]
包含式扫描输出: Y_inc = [3, 9, 16, 20, ...]
排除式扫描输出: Y_exc = [0, 3, 9, 16, ...]

顺序扫描的伪代码如下：

// 包含式扫描
Y[0] = X[0];
for (int i = 1; i < n; i++) {
    Y[i] = Y[i-1] + X[i]; // 使用加法运算符
}

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/14e88fa3bece47823cac167515ef8a0e_7.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/14e88fa3bece47823cac167515ef8a0e_8.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/14e88fa3bece47823cac167515ef8a0e_10.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/14e88fa3bece47823cac167515ef8a0e_12.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/14e88fa3bece47823cac167515ef8a0e_14.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/14e88fa3bece47823cac167515ef8a0e_15.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/14e88fa3bece47823cac167515ef8a0e_17.png)

// 排除式扫描
Y[0] = 0; // 单位元
for (int i = 1; i < n; i++) {
    Y[i] = Y[i-1] + X[i-1]; // 使用加法运算符
}

分段扫描（Segmented Scan）

与归约类似，并行扫描也需要线程间的同步。由于无法跨线程块同步，解决方案是进行分段（或分层）扫描。

以下是分段扫描的步骤：

本地扫描：每个线程块负责输入数组的一个片段，并在其内部进行并行扫描。
收集部分和：每个线程块计算其片段的“部分和”（即该片段所有元素的组合结果），并将其存储到一个全局数组中。
扫描部分和：对存储部分和的全局数组执行一次扫描操作。
更新结果：每个线程块将其扫描结果与“前驱块的部分和扫描结果”相加，得到最终的全局扫描结果。

这种方法允许不同线程块并行处理各自的片段，然后通过一次额外的扫描和加法阶段来合并结果。扫描部分和数组通常通过递归启动新的内核来实现。

Kogge-Stone 并行扫描算法

现在，我们深入探讨如何在单个线程块内实现并行扫描。我们将使用 Kogge-Stone 算法来实现包含式扫描。

该算法的核心思想是通过多轮迭代，逐步扩大每个元素所“看到”的前驱元素范围。

算法步骤（以包含式扫描、8个元素、加法运算符为例）：

初始化：每个线程负责一个输入元素，并将其值加载到共享内存缓冲区。
第1轮（步长=1）：对于索引 i >= 1 的线程，执行 output[i] += output[i-1]。完成后，output[1] 包含了 X[0]+X[1]。
第2轮（步长=2）：对于索引 i >= 2 的线程，执行 output[i] += output[i-2]。完成后，output[3] 包含了 X[0] 到 X[3] 的和。
第3轮（步长=4）：对于索引 i >= 4 的线程，执行 output[i] += output[i-4]。完成后，所有 output[i] 都包含了从 X[0] 到 X[i] 的和。

对于一个包含 n 个元素的块，该算法需要 log n 轮迭代。在每一轮中，大约有一半的线程是活跃的。

基础实现的关键点：

需要将数据从全局内存加载到共享内存以提升性能。
在读取和写入之间必须进行同步（__syncthreads()），以防止数据竞争。一个简单的实现需要在每轮迭代中进行两次同步：一次确保所有读取完成，另一次确保所有写入完成。

优化：双缓冲（Double Buffering）

上一节的基础实现中，我们在同一块共享内存缓冲区上进行读写，因此需要严格的同步来避免竞争。双缓冲技术可以消除其中一次同步。

双缓冲原理：
我们创建两个共享内存缓冲区：bufferA 和 bufferB。在每一轮迭代中：

其中一个作为输入缓冲区（in_buf），另一个作为输出缓冲区（out_buf）。
活跃线程从 in_buf 读取数据，进行计算，并将结果写入 out_buf。
不活跃的线程（即本轮不需要进行加法操作的线程）则简单地将 in_buf 中的对应值复制到 out_buf。
一轮迭代结束后，交换 in_buf 和 out_buf 的指针，进行下一轮。

优势：由于读写发生在不同的物理缓冲区，我们无需在单次迭代内区分读后同步和写前同步，只需在每轮迭代结束后进行一次 __syncthreads() 即可。这减少了同步开销，提升了性能。

实现排除式扫描

基于包含式扫描的实现，可以很容易地修改为排除式扫描。

方法：在将数据加载到共享内存时，将整个片段向右移动一位。

对于线程块内的第一个线程（索引0），将其共享内存位置初始化为单位元（如0）。
对于其他线程（索引 i > 0），加载输入数组中索引为 (global_i - 1) 的元素。
然后，对这个“移位后”的数组运行包含式扫描的 Kogge-Stone 算法。
最后，在计算部分和以及后续的添加阶段需要做细微调整，以确保所有元素都被正确计入。

另一种思路是直接在全局内存层面将输入数组移位，然后调用包含式扫描内核。

工作效率（Work Efficiency）分析

一个并行算法如果其执行的总操作量（如浮点加法次数）与对应的顺序算法相同，则称其为工作高效的。

顺序扫描：对于 n 个元素，需要 n-1 次加法操作。时间复杂度为 O(n)。
Kogge-Stone 并行扫描：分析其操作总数。对于 n 个元素，有 log n 轮。总操作数约为 n * log n - n。时间复杂度为 O(n log n)。

因此，Kogge-Stone 算法不是工作高效的。它通过更多的并行操作来换取更少的执行步数（log n 步）。在计算资源充足、并行度完美的情况下，这能带来加速。但如果硬件资源有限（例如，线程块被序列化执行），由于其总工作量更大，实际执行时间可能反而比高效的低并行度算法更长。

总结

本节课我们一起学习了：

扫描操作的定义、包含式与排除式的区别及其顺序实现。
如何在GPU上通过分段扫描的策略处理大规模数组。
Kogge-Stone 并行扫描算法在线程块内的具体步骤和实现。
使用共享内存和双缓冲技术来优化扫描内核的性能。
如何基于包含式扫描修改代码来实现排除式扫描。
并行算法的工作效率概念，并分析了Kogge-Stone算法在工作效率上的不足。

下一节课，我们将学习另一种工作更高效的并行扫描算法，并对比它们的性能。

GPU计算：第12讲：扫描算法（Brent-Kung）🔬

概述

在本节课中，我们将学习一种不同于上次讨论的并行扫描算法——Brent-Kung并行算法。我们将比较它与Kogge-Stone算法的异同，分析其工作效率，并探讨线程粗化（Thread Coarsening）优化在扫描算法中的应用。

回顾：扫描模式与Kogge-Stone算法

上一节我们介绍了扫描模式以及使用Kogge-Stone方法进行并行扫描。扫描模式是指，给定一个输入数组，输出数组的每个元素将是其前面所有元素的组合（对于包含性扫描），或者仅是其前面元素的组合（对于排他性扫描）。我们以加法为例，但该模式适用于其他结合性操作符，如乘法、最小值、最大值等。

并行执行扫描需要同步，而不同线程块之间无法同步。因此，我们采用分段处理的方法：

将输入数组分段。
每个线程块扫描其负责的段，并将部分和存储在一个子数组中。
扫描这个部分和数组。
每个线程块将排在其前面的块的部分和扫描结果加到其内部元素上，最终得到全局扫描结果。

Kogge-Stone算法在每个线程块内执行扫描，其核心操作是让每个线程迭代地加上与其相距特定步长（stride）的元素，步长从1开始，每次翻倍，直到覆盖整个数组。由于存在读写冲突（同一迭代中，一个线程读取的位置可能被另一个线程写入），我们需要在读写之间进行同步。此外，内存位置被重用，因此我们使用共享内存来提升性能。

为了减少同步开销，我们引入了双缓冲（Double Buffering）优化。线程不再每次都读写同一个缓冲区，而是使用一个输入缓冲区和一个输出缓冲区，每次迭代后交换它们，这样每次迭代只需要一次同步。

Kogge-Stone算法的工作效率分析显示，其时间复杂度为O(n log n)，而顺序扫描算法为O(n)。这表明并行算法做了更多冗余工作，如果由于资源限制导致线程块串行执行，这种低工作效率可能会影响性能。

Brent-Kung并行扫描算法

本节中，我们来看看另一种并行扫描算法——Brent-Kung方法。与Kogge-Stone相比，它旨在提高工作效率。

Brent-Kung并行包含性扫描算法分为两个阶段：

1. 归约阶段（Reduction Step）

这个阶段类似于归约树。

第一步，线程将相邻的两个元素相加。
第二步，线程将相距两个步长的元素相加。
以此类推，步长每次翻倍，直到完成对整个数组的归约。
在此过程中，部分扫描值已经就绪，但并非全部。

2. 后归约阶段（Post-Reduction Step）

这个阶段用于完成剩余未就绪的扫描值。

从最大的步长开始，反向逐步减小步长。
在每个步长下，线程将其负责的元素值加到其右侧相距该步长的元素上。
通过这种方式，最终补全所有扫描值。

算法图解示意（以8个元素为例）：

归约阶段 (步长s=1,2,4):
x0 x1 x2 x3 x4 x5 x6 x7
  \/    \/    \/    \/   (s=1)
   a     b     c     d
     \   /       \   /   (s=2)
       e           f
           \     /       (s=4)
               g

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/5cad7391a7b1791668885f3b803b3542_1.png)

后归约阶段 (步长s=4,2,1):
               g (设为0用于排他性扫描)
           /     \       (s=4)
       e           f
     /   \       /   \   (s=2)
   a     b     c     d
  / \   / \   / \   / \  (s=1)
x0 x1 x2 x3 x4 x5 x6 x7

算法比较与分析

以下是两种方法的对比：

特性	Kogge-Stone	Brent-Kung
总加法操作数	多 (O(n log n))	少 (O(n))
步骤数	少 (log n)	多 (~2 log n)
工作效率	较低	较高
是否需要双缓冲	是	否

Brent-Kung工作量分析：

归约阶段：log n 步，总操作数约为 n-1。
后归约阶段：log n - 1 步，总操作数约为 n - 2log n。
总计：约 2log n - 1 步，2n - 2log n - 2 次操作，时间复杂度为 O(n)。

因此，Brent-Kung算法用更多的并行步骤换取了更高的工作效率（更少的总体操作）。

Brent-Kung算法实现要点

接下来，我们探讨实现Brent-Kung算法的关键优化和注意事项。

共享内存的使用

与Kogge-Stone类似，使用共享内存有两个好处：

数据复用：减少访问全局内存的次数。
促进合并访问：在Brent-Kung中，直接使用全局内存的访问模式是非合并的。先将数据以合并方式加载到共享内存，可以大幅提升内存访问效率。

无需双缓冲

在Brent-Kung算法中，同一迭代内没有任何数据元素会被不同的线程同时读写，因此不需要双缓冲优化，使用单个缓冲区即可。

控制流分歧与线程重索引

在Kogge-Stone中，每个线程固定负责一个数据元素，控制流分歧不严重。但在Brent-Kung中，如果固定线程到元素，每个迭代都会有大量线程不工作，导致严重的控制流分歧，浪费GPU计算资源。

为了最小化控制流分歧，我们采用线程重索引策略：

不再让线程固定负责某个数据元素。
在每次迭代中，根据当前步长，动态计算每个活动线程应负责处理哪个元素。
目标是让活动线程的ID总是从0开始连续排列，从而尽可能填满前几个线程束（Warp），减少线程束内部分线程空闲的情况。

索引计算：
经过推导，对于给定的线程索引 t 和当前步长 s，线程在归约阶段和后归约阶段应操作的元素索引 i 可以用同一个公式计算：
i = (t + 1) * 2 * s - 1
在代码中，我们需要检查计算出的索引 i 及其目标索引（如 i - s 或 i + s）是否在数组边界内。

线程重索引是避免控制流分歧的通用策略，其核心思想是根据计算阶段动态分配工作，以保持线程的活跃度。

从包含性扫描到排他性扫描

对于排他性扫描，有两种实现方式：

基于包含性扫描转换：先进行包含性扫描，然后将结果数组整体右移一位，并在开头插入0。
修改Brent-Kung的后归约阶段：可以设计一个专门用于排他性扫描的后归约步骤。其归约阶段与包含性扫描相同，但在后归约阶段，通过不同的数据移动和相加顺序，直接得到排他性扫描结果。这种方法更精巧，也是练习线程索引计算的良好案例。

线程粗化优化

我们曾提到，并行化可能带来工作效率下降的开销。如果硬件资源不足导致线程块串行执行，这种开销将得不偿失。线程粗化 就是用来减少因过度并行化导致工作效率损失的优化技术。

在线程块内部应用线程粗化进行扫描：

块内分段：每个线程块将其负责的段进一步划分为更小的子段。
线程顺序扫描：每个线程顺序地扫描其分配到的子段。顺序扫描是工作高效的。
收集部分和：每个线程将其子段的总和（部分和）放入一个数组中。
并行扫描部分和：线程块使用之前讨论的并行扫描算法（如Kogge-Stone）扫描这个部分和数组。
添加偏移量：每个线程将扫描得到的、对应其之前所有线程的部分和结果，加到其子段的所有元素上。

这实质上是将“分段扫描”的思想应用到了线程级别。通过让每个线程处理多个元素（粗化），减少了并行扫描的工作量，从而提高了整体工作效率。

性能提升：在示例中，应用线程粗化（因子为8）后，Kogge-Stone扫描算法的运行时间从约4.5毫秒减少到约1.5毫秒，获得了近3倍的加速。

总结

本节课我们一起学习了：

Brent-Kung并行扫描算法：它通过归约和后归约两个阶段，以更多的步骤换取了O(n)的更高工作效率。
算法实现优化：包括使用共享内存促进合并访问、利用其无读写冲突的特性避免双缓冲，以及通过线程重索引来最小化控制流分歧。
排他性扫描的实现：可通过转换或修改后归约阶段实现。
线程粗化技术：通过让每个线程顺序处理多个数据，减少并行开销，显著提升扫描算法在实际GPU上的性能。这体现了在工作效率和并行度之间进行权衡的设计思想。

通过对比Kogge-Stone和Brent-Kung，我们深入理解了并行算法设计中工作量、步骤数和实际硬件执行效率之间的复杂关系。

GPU计算：13：直方图与原子操作

在本节课中，我们将学习一个新的并行模式——直方图。我们将探讨其基本概念、并行实现时遇到的挑战，以及如何使用原子操作来解决数据竞争问题。最后，我们将介绍两种关键的优化技术：私有化和线程粗化，以提升直方图计算的性能。

回顾：扫描操作

上一节我们介绍了并行扫描操作，特别是Brent-Kung方法。本节中，我们来看看直方图。

我们讨论了两种并行扫描方法：Kogge-Stone和Brent-Kung。

Kogge-Stone方法需要 log n 步完成，总操作数为 O(n log n)。
Brent-Kung方法需要 2 log n - 1 步，但只执行 O(n) 次操作，因此具有更高的工作效率。

我们还对扫描操作应用了多种优化，包括使用共享内存（这同时实现了访存合并）以及通过重新索引线程来最小化控制流分化。最后，我们探讨了线程粗化如何通过让每个线程顺序处理一段数据，再并行扫描部分和，来提升算法的工作效率。

什么是直方图？

直方图用于近似表示数据集的分布。其工作原理是将输入数据集中可能出现的值范围划分为若干个“区间”（或称为“桶”），然后统计落入每个区间的数据数量。

一个常见的应用是颜色直方图。例如，一张图像的像素值通常用8位整数表示，范围是0到255。颜色直方图可以统计每个像素值出现的次数，从而概括图像的色彩分布、明暗程度等信息。

一个顺序实现的直方图伪代码如下：

for (int i = 0; i < num_pixels; i++) {
    unsigned char pixel_value = image[i];
    bins[pixel_value]++; // 递增对应区间的计数
}

并行直方图的挑战

如何并行化直方图操作？一个直观的想法是将图像数据分配给多个线程，每个线程负责更新自己对应的区间计数。我们可以为每个像素分配一个线程。

以下是初步的CUDA内核代码框架：

__global__ void histogram_kernel(unsigned char* image, int* bins, int num_pixels) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < num_pixels) {
        unsigned char pixel_value = image[i];
        bins[pixel_value]++; // 存在数据竞争！
    }
}

然而，这个实现存在严重问题。当多个线程同时读取、修改并写入同一个内存位置（例如bins[10]）时，会发生数据竞争。因为bins[pixel_value]++并非原子操作，它包含加载、加一、存储多个步骤，这些步骤的交错执行可能导致最终计数错误。

原子操作：解决数据竞争

为了避免数据竞争，对同一内存位置的并发“读-改-写”操作必须是互斥的。在CPU上，我们通常使用锁（如互斥锁）来实现互斥。但在GPU上，使用锁可能导致死锁（例如，同一线程束中的线程相互等待）。

GPU提供了更优的解决方案：原子操作。原子操作能保证“读-改-写”作为一个不可分割的单一硬件指令执行。在此期间，硬件会确保没有其他线程能访问该内存位置。

CUDA提供了多种原子操作，对于直方图，我们使用atomicAdd：

int atomicAdd(int* address, int val); // 将val加到address指向的值，返回旧值

修改我们的内核，使用原子操作：

__global__ void histogram_kernel_atomic(unsigned char* image, int* bins, int num_pixels) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < num_pixels) {
        unsigned char pixel_value = image[i];
        atomicAdd(&bins[pixel_value], 1); // 使用原子加操作
    }
}

现在，代码能正确运行并得到准确结果。

优化一：私有化

虽然原子操作解决了正确性问题，但性能可能不佳。全局内存的原子操作延迟很高，且当大量线程更新同一个区间时，会发生严重竞争，导致序列化执行。

私有化是一种优化技术。其核心思想是：不让所有线程块直接更新全局直方图，而是让每个线程块维护一个私有的直方图副本。

每个线程块在共享内存中分配一个私有直方图数组，并初始化为0。
线程块内的线程更新这个私有副本（可以使用更快的共享内存原子操作）。
当线程块处理完所有分配的数据后，将私有副本中的非零值累加到全局直方图中（此时仍需使用全局原子操作，但次数大大减少）。

这种方法将线程块间的竞争转移到了线程块内部，而共享内存的原子操作比全局内存快得多，从而显著提升了性能。

优化二：线程粗化

线程粗化是另一种优化手段。在最初的实现中，我们为每个像素分配一个线程，这会产生大量线程块和私有副本，导致最终向全局内存提交的次数较多。

通过线程粗化，我们减少线程块数量，让每个线程处理多个输入元素。

每个线程块负责更大的一段输入数据。
每个线程循环处理该段数据中的多个像素，并更新其所在线程块的私有直方图。
确保以合并访问的方式加载输入数据以保持内存访问效率。

这样做的好处是：

减少了线程块数量，从而减少了需要创建和提交的私有副本数量。
进一步降低了全局原子操作的总次数。
更多的更新发生在快速的共享内存私有副本中。

结合私有化和线程粗化，可以最大限度地减少昂贵的全局内存原子操作，从而大幅提升直方图计算的并行性能。

总结

本节课我们一起学习了直方图这一并行模式。我们首先了解了直方图的基本概念和顺序实现。然后，我们尝试并行化时遇到了数据竞争问题，并引入了GPU的原子操作作为解决方案。接着，为了提升原子操作带来的性能瓶颈，我们探讨了两种关键优化：私有化（通过创建本地副本来减少全局竞争）和线程粗化（通过让每个线程处理更多工作来减少私有副本数量）。这些技术不仅适用于直方图，也是优化许多存在输出竞争的并行算法的通用策略。

GPU计算：第14讲：合并模式 🧩

在本节课中，我们将要学习一个新的并行计算模式：合并。我们将探讨如何将两个已排序的数组合并成一个更大的有序数组，并学习如何在GPU上高效地并行实现这一操作。

概述

上一节课我们介绍了直方图模式，并引入了原子操作和私有化优化。直方图用于统计数据集中落入各个区间的值的数量。由于多个线程可能同时更新同一个区间，我们遇到了数据竞争问题。在CPU上，我们通常使用互斥锁来解决，但在GPU上，锁可能导致死锁。因此，CUDA提供了原子操作，如 atomicAdd，来确保内存操作的原子性。为了减少全局内存上的原子操作竞争，我们使用了私有化优化，即每个线程块拥有直方图的私有副本，最后再合并到全局副本中。

本节课，我们将学习合并模式。合并操作接收两个已排序的列表，并将它们组合成一个单一的、更大的有序列表。我们将从顺序合并算法开始，然后探讨如何将其并行化，并解决并行化过程中的关键挑战。

顺序合并算法

在讨论并行合并之前，我们先来实现一个顺序合并算法。这个算法将作为我们后续并行实现的基准。

顺序合并的思路是使用三个指针：i 指向数组A，j 指向数组B，k 指向输出数组C。我们比较 A[i] 和 B[j]，将较小的元素放入 C[k]，并移动相应的指针。当其中一个数组耗尽后，我们将另一个数组的剩余元素直接复制到C中。

以下是顺序合并的代码实现：

void mergeSequential(int* A, int* B, int* C, unsigned int M, unsigned int N) {
    unsigned int i = 0; // 索引数组A
    unsigned int j = 0; // 索引数组B
    unsigned int k = 0; // 索引数组C

    // 合并两个数组，直到其中一个耗尽
    while (i < M && j < N) {
        if (A[i] <= B[j]) {
            C[k++] = A[i++];
        } else {
            C[k++] = B[j++];
        }
    }

    // 如果数组A还有剩余元素，复制到C
    while (i < M) {
        C[k++] = A[i++];
    }

    // 如果数组B还有剩余元素，复制到C
    while (j < N) {
        C[k++] = B[j++];
    }
}

这个算法的时间复杂度是 O(M + N)，其中M和N分别是两个输入数组的大小。

并行合并的挑战与思路

顺序合并算法本质上是串行的，因为每一步的输出都依赖于上一步的比较结果。那么，如何在GPU上并行化这个操作呢？

一个直观的想法是对输出进行分区，而不是对输入进行分区。我们将输出数组C划分为多个大小相等的段，并让每个GPU线程负责合并生成其中一个段。这样，每个线程的工作量大致相同，有利于负载均衡。

然而，关键挑战在于：每个线程如何找到其负责的输出段所对应的输入段（即在数组A和B中的起始位置）？

给定一个输出索引 k（线程负责段的起始位置），我们需要找到对应的输入索引 i（在A中）和 j（在B中）。我们称 i 和 j 为 k 的协同秩。

寻找协同秩

我们的目标是：给定 k，找到 i 和 j。

首先，我们观察到输出数组C中 k 之前的元素数量，等于输入数组A中 i 之前的元素数量加上输入数组B中 j 之前的元素数量。因此，我们有基本关系：

k = i + j

所以，如果我们能找到 i，就可以通过 j = k - i 计算出 j。

接下来，我们需要确定 i 的取值范围。i 是数组A的索引，因此：

i 的下界：i >= 0 且 i >= k - N（因为 j = k - i 必须 <= N）。所以下界是 max(0, k - N)。
i 的上界：i <= M 且 i <= k（因为 j = k - i 必须 >= 0）。所以上界是 min(M, k)。

因此，i 的取值范围是：[max(0, k - N), min(M, k)]。

在这个范围内，我们需要找到正确的 i。由于数组A和B都是有序的，我们可以在这个区间内进行二分查找。

判断一个候选值 i 是否正确（以及对应的 j = k - i）的条件是：

A[i-1] <= B[j]
B[j-1] <= A[i]

如果 i 的猜测值太高，我们会发现 A[i-1] > B[j]。
如果 i 的猜测值太低，我们会发现 B[j-1] > A[i]。

基于这些规则，我们可以实现一个二分查找函数来找到正确的协同秩 i。

并行合并内核实现

现在，我们可以实现并行的合并内核了。每个线程的执行步骤如下：

计算自己负责的输出段的起始索引 k。
调用 coRank 函数，通过二分查找找到对应的输入索引 i。
计算 j = k - i。
为了进行合并，线程还需要知道其输入段的大小。这可以通过计算下一个线程的起始索引 k_next，并再次调用 coRank 找到 i_next 和 j_next 来获得。那么，该线程负责的A段大小就是 i_next - i，B段大小是 j_next - j。
最后，线程调用顺序合并函数，将其负责的A和B的子段合并到输出数组C的对应位置。

以下是并行合并内核的伪代码框架：

__global__ void mergeParallel(int* A, int* B, int* C, unsigned int M, unsigned int N, unsigned int elementsPerThread) {
    unsigned int tid = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned int k = tid * elementsPerThread; // 该线程负责段的起始位置

    if (k < M + N) { // 边界检查
        // 1. 找到本段起始的协同秩
        unsigned int i = coRank(A, B, M, N, k);
        unsigned int j = k - i;

        // 2. 找到下一段起始的协同秩，以确定本段大小
        unsigned int k_next = min(k + elementsPerThread, M + N);
        unsigned int i_next = coRank(A, B, M, N, k_next);
        unsigned int j_next = k_next - i_next;

        // 3. 执行顺序合并
        mergeSequential(&A[i], &B[j], &C[k], i_next - i, j_next - j);
    }
}

性能优化：使用共享内存

我们实现的初始并行版本存在一个主要的性能问题：内存访问不连续（非合并访问）。

在 coRank 函数中，每个线程独立进行二分查找，会随机访问全局内存中的A和B数组。
在顺序合并阶段，不同线程访问的输入子段位置也各不相同。

这导致了大量低效的全局内存访问。优化方法是利用共享内存：

以线程块为单位查找协同秩：每个线程块只计算一次其负责的整个输出大段的协同秩（i_block, j_block），而不是每个线程都算。这减少了全局内存上的二分查找次数。
协作加载到共享内存：线程块内的所有线程协作，以合并访问的方式将A和B的对应大段从全局内存加载到共享内存中。
在共享内存中执行操作：每个线程在共享内存的数组副本上进行二分查找（寻找其子段的协同秩）和顺序合并。共享内存的延迟远低于全局内存，且随机访问代价小。
协作写回全局内存：所有线程协作，将共享内存中合并好的结果以合并访问的方式写回全局内存。

这种优化策略将耗时的、非合并的全局内存访问，转换为了快速的共享内存访问，并确保了进出全局内存的访问是高效的合并访问。

总结

本节课我们一起学习了合并模式。我们从基础的顺序合并算法出发，理解了其工作原理。然后，我们探讨了在GPU上并行化合并操作的核心思想：对输出进行分区，并让每个线程负责一个输出段。

实现并行合并的关键在于解决如何为每个输出段找到对应的输入段这一问题，我们引入了协同秩的概念，并利用二分查找在有序数组中高效地找到它。

最后，我们分析了初始并行实现的性能瓶颈——非合并的全局内存访问，并提出了使用共享内存进行优化的策略：以线程块为单位组织计算，协作加载数据，在共享内存中完成核心操作，再协作写回结果。

通过本节课的学习，你应该掌握了合并模式的基本原理、并行化方法以及针对GPU架构的性能优化思路。

GPU计算｜CMPS 297S396AA：15：排序算法

在本节课中，我们将学习两种重要的并行排序算法：基数排序和归并排序。我们将重点探讨基数排序的原理、并行化策略及其优化方法，并简要介绍如何将之前学习的并行归并模式应用于归并排序。

概述

上一节课我们介绍了并行归并模式。我们提到，有序归并是将两个有序列表合并成一个有序列表的操作。为了并行化这个操作，一种方法是将要生成的输出列表划分为相等的段，每个线程负责一段。每个线程需要在两个输入数组A和B中找到对应的输入段，然后顺序合并这两个段。并行归并的关键挑战在于每个线程如何找到其在A和B中对应的输入段。

我们定义了“协同秩”的概念：对于输出数组C中的索引K，其在数组A中的协同秩I，以及在数组B中的协同秩J。我们发现J = K - I，因此核心任务是找到I。我们通过设定I的边界并进行二分搜索来确定I，判断条件基于一个观察：A中I之前的元素应小于B中J之后的元素，反之亦然。

每个线程通过二分搜索找到I和J，然后顺序合并输入段到输出段。这个操作的内存访问不是合并的，因为二分搜索是随机访问，且后续的顺序合并中，每个线程访问的输入段也不同。为了优化，我们让一个线程块中的一个线程确定整个块的输入段范围（I_block, J_block），然后块内线程协作，以合并的方式将这些段加载到共享内存中。接着，所有线程在共享内存中进行二分搜索和合并，最后再以合并的方式将结果从共享内存写回全局内存。

我们还提到了线程聚合。并行化归并的代价是每个线程都需要进行二分搜索。通过为每个线程分配多个输出元素（即一个输出段），我们分摊了二分搜索的成本，这本身就是一种线程聚合。

本节课，我们将开始讨论排序算法，主要关注基数排序，并简要介绍归并排序。

基数排序 🧮

基数排序是一种非比较型排序算法。它的工作原理是，根据基数（或进制）将待排序的关键字分配到多个“桶”中。如果输入关键字采用某种位置计数系统表示，我们就使用该基数，并迭代处理关键字的每一位数字。每次迭代，都根据当前处理的数字位将关键字分配到对应的桶中。

一个重要的原则是，每次迭代时，需要保留前一次迭代在每个桶内建立的顺序。通常，在计算机中实现基数排序时，我们喜欢使用2的幂作为基数，因为这简化了对二进制数的处理。每次迭代处理关键字中固定的一组比特位。

为了简化说明，我们将从基数为2（即每次处理1个比特位）的情况开始，然后再扩展到更大的基数。

单比特基数排序演示

假设我们有一个数组需要排序。使用1比特基数排序时，我们从最低有效位开始处理。

第一步：处理最低有效位
我们将所有关键字根据其最低有效位是0还是1分配到两个桶中。将所有最低有效位为0的元素放在一起，为1的元素放在其后。此时，数组仅按最低有效位排序。

第二步：处理次低有效位
接下来，我们处理次低有效位。根据该位的值（0或1）再次分配元素。关键点在于，在分配时，我们必须保持每个桶内元素在上一步中已建立的相对顺序。这样，数组现在就按最低的两个比特位排序了。

后续步骤
我们重复这个过程，从低到高依次处理每一位。每次迭代都根据当前比特位分配元素，并保持前一次迭代在每个桶内建立的顺序。处理完所有比特位后，数组就完全排序了。

基数排序是非比较型排序算法，通常要求关键字具有固定大小，以便可以划分为数字位。如果关键字不满足此条件，则需要使用比较型排序算法，如我们稍后将看到的归并排序。

并行化基数排序的关键步骤

基数排序的核心操作可以归结为：在每次迭代中，根据当前关注的比特位，将数组元素分离，使所有该位为0的元素在前，为1的元素在后，同时保持原有顺序。

问题转化为：如何为输入数组中的每个元素找到其在输出（排序后）数组中的目标索引？

分析目标索引：

对于一个0元素，其目标索引等于其左侧0元素的数量。这又等于该元素的索引 减去 其左侧1元素的数量。
对于一个1元素，其目标索引等于数组中0元素的总数 加上 其左侧1元素的数量。而0元素的总数等于数组总大小 减去 1元素的总数。

因此，无论是0还是1元素，其目标索引的计算都依赖于一个关键信息：每个元素左侧的1元素数量。

解决方案：使用独占扫描

提取比特位：为每个元素提取当前要处理的比特位，生成一个由0和1组成的数组 B。
执行独占扫描：对数组 B 执行独占扫描操作。扫描结果数组 S 中，S[i] 就表示原数组中第 i 个元素左侧的1元素数量。同时，我们也能得到整个数组 B 的总和，即1元素的总数 total_ones。
计算目标索引：对于索引为 i 的元素：
- 如果它是 0 (B[i] == 0)，其目标索引 dest = i - S[i]。
- 如果它是 1 (B[i] == 1)，其目标索引 dest = (n - total_ones) + S[i]，其中 n 是数组大小。
分散写入：每个线程根据计算出的目标索引，将其负责的元素写入输出数组的对应位置。

并行化：上述步骤很容易并行化。我们可以为每个输入元素分配一个线程。提取比特位是局部操作。扫描操作我们已经知道如何高效并行化。计算目标索引和分散写入也是简单的逐元素操作。

优化：改善合并访问

然而，上述朴素方法存在一个性能问题：分散写入步骤的内存访问不是合并的。相邻的线程可能将元素写入输出数组中相距很远的位置，导致内存访问效率低下。

优化策略：局部排序与分桶写入
我们可以利用共享内存来改善合并访问。思路是，先在线程块内部进行局部排序和分桶，然后再以合并的方式将整个桶写入全局内存。

局部分离：每个线程块将其负责的输入元素加载到共享内存中。块内线程协作，在共享内存中执行上述的扫描和分离操作。这样，在该块的共享内存中，所有当前位为0的元素被聚集在一起，所有为1的元素被聚集在另一处。
计算全局桶偏移：每个线程块需要知道它的0桶和1桶应该写入全局输出数组的哪个起始位置。
- 每个块计算并输出它包含的0元素数量 (count0) 和1元素数量 (count1)。
- 将所有块的 (count0, count1) 数组扁平化，得到一个一维数组。
- 对这个一维数组执行独占扫描。扫描结果中，对应每个块 count0 的位置给出了该块0桶的全局起始索引；对应每个块 count1 的位置给出了该块1桶的全局起始索引（需要加上所有块0元素的总数）。
合并写入：每个线程块根据得到的全局起始索引，将其共享内存中的0桶和1桶分别连续地写入全局内存。由于一个桶内的元素是连续的，且由块内相邻的线程写入，因此这些写入操作是高度合并的。

这种优化将非合并的随机写入，转换为了两次合并的顺序写入（先写所有0桶，再写所有1桶），显著提升了内存带宽利用率。

基数选择与线程聚合

基数大小的权衡
我们之前使用了1比特基数。对于一个32位的关键字，这需要32次迭代。为了减少迭代次数，我们可以使用更大的基数，例如2比特基数。这样只需要16次迭代。在局部，可以通过连续执行两次1比特分离操作来实现一次2比特分离。

然而，更大的基数意味着更多的桶（2比特对应4个桶）。虽然迭代次数减少，但每个线程块需要向更多不同的全局内存地址范围写入数据，这可能损害合并访问的效果。因此，基数的选择需要在迭代次数和合并访问效率之间取得平衡。

线程聚合以改善合并访问
另一种优化手段是线程聚合，即为每个线程分配多个输入元素。这样，每个线程块处理的数据量更大，产生的每个桶也更大。当块将大桶写入全局内存时，由于桶内元素更多，写入操作的连续性更好，从而进一步改善了合并访问。

线程聚合增加了每个线程的工作量，但通过提高内存访问效率，整体性能通常能得到提升。

归并排序 🔀

基数排序并非适用于所有类型的数据（例如可变长度关键字）。我们仍然需要一种可并行化的比较型排序算法，归并排序是一个流行的选择。

归并排序是分治算法：

分解：将列表递归地分成两半，直到子列表足够小（例如，只有一个元素）。
解决：排序这些子列表（对于很小的列表，排序很简单）。
合并：将已排序的子列表两两合并，形成更大的有序列表，直到最终合并成完整的有序列表。

并行化归并排序

我们可以利用上一课学习的并行归并模式来并行化归并排序。

并行策略：

底层排序：最初，将输入数组划分为许多小的片段。可以使用任何排序方法（甚至是另一个内核）对这些小片段进行排序。这些排序操作可以并行进行。
逐层合并：排序后，我们得到许多小的有序数组。然后开始合并阶段：
- 在第一次合并中，将相邻的两个有序数组合并为一个。每个合并操作可以分配一个或多个线程块来并行执行（使用我们学过的并行归并算法）。
- 由于存在多个独立的合并对，这些合并操作之间也可以并行执行。
递归合并：重复此过程。在每一层，合并后的数组变大，但需要合并的对数减少。因此，并行性从跨多个合并操作之间的并行，逐渐转向单个大型合并操作内部的并行（使用更多线程块协作完成一次合并）。

通过这种方式，归并排序的整个合并过程可以被有效地映射到GPU的并行架构上。

总结

本节课我们一起学习了两种重要的并行排序算法。

我们深入探讨了基数排序，这是一种非比较型排序。我们理解了其按位分配、保持顺序的原理，并分析了其并行化方法。重点是，我们讨论了如何通过局部扫描分离和全局桶偏移计算来优化分散写入时的内存合并访问。同时，我们也探讨了基数大小选择与线程聚合对性能的影响。

随后，我们简要介绍了归并排序，这是一种比较型排序。我们看到了如何将之前学习的并行归并模式应用于归并排序的合并阶段，通过在不同层级和同一合并操作内部利用并行性，来实现高效的排序。

理解这些排序算法的并行化策略及其优化，对于在GPU上实现高性能计算至关重要。

GPU计算：16：稀疏矩阵计算（COO与CSR格式）

在本节课中，我们将要学习稀疏矩阵计算，并重点探讨两种经典的稀疏矩阵存储格式：坐标格式（COO）和压缩稀疏行格式（CSR）。我们将以稀疏矩阵-向量乘法（SPMV）作为案例研究，分析不同格式在并行化、内存访问和性能方面的权衡。

什么是稀疏矩阵？

在深入探讨存储格式之前，我们首先需要理解什么是稀疏矩阵。

一个稠密矩阵是指矩阵中大多数元素都是非零值的矩阵。相反，一个稀疏矩阵则是指矩阵中包含大量零元素的矩阵。在实际应用中，许多系统产生的矩阵都是稀疏的，这意味着矩阵中绝大多数元素的值都是零。

利用矩阵的稀疏性，我们可以获得多方面的优势：

节省内存/存储：通过压缩矩阵，只存储非零值，无需为零值分配空间。
节省内存带宽：计算时无需加载零值。
节省计算时间：计算时无需处理零值。

因此，稀疏矩阵计算的核心挑战之一，就是设计高效的存储格式，以便在节省内存的同时，还能优化计算性能。

稀疏矩阵存储格式概览

存在多种稀疏矩阵存储格式，每种格式都有其设计目标和适用场景。在本课程中，我们将重点介绍四种格式，它们足以展示稀疏矩阵处理的核心概念：

坐标格式
压缩稀疏行格式
ELLPACK格式
锯齿状对角线格式

今天我们将专注于前两种格式：COO和CSR。

设计存储格式时，通常需要考虑以下几个关键因素：

空间效率：格式能节省多少内存。
灵活性：是否便于添加、重排或删除矩阵中的元素。
可访问性：是否便于以特定方式（如按行或按列）访问数据。
内存访问模式：在GPU上下文中，格式是否支持内存合并访问。
负载均衡：在GPU上下文中，格式是否有助于最小化控制流发散。

我们将以稀疏矩阵-向量乘法作为具体应用场景，来评估COO和CSR格式在这些因素上的表现。SPMV的运算可以表示为：
y = A * x
其中 A 是稀疏矩阵，x 是稠密输入向量，y 是稠密输出向量。

坐标格式

首先，让我们来看看坐标格式。

COO格式的原理

COO格式的存储方式非常直观：它为矩阵中的每一个非零元素存储三个值：元素值、行索引和列索引。

假设我们有一个稀疏矩阵，其非零元素分布如下（白色方格代表零值）：

在COO格式中，我们会创建三个数组：

values: 按顺序存储所有非零值，例如 [1, 7, 5, 3, 9, 2, 8]。
row_indices: 存储每个非零值对应的行索引，例如 [0, 0, 1, 1, 1, 2, 2]。
col_indices: 存储每个非零值对应的列索引，例如 [0, 1, 0, 2, 3, 1, 2]。

这样，数组中的第 i 个元素就完整地描述了矩阵中的一个非零值：(values[i], row_indices[i], col_indices[i])。

基于COO格式的SPMV实现

在COO格式下，实现SPMV最自然的并行化策略是：为每一个非零元素分配一个线程。

每个线程的执行步骤如下：

根据线程ID i，从 values、row_indices 和 col_indices 数组中读取其负责的非零值 val、行号 row 和列号 col。
从输入向量 x 中读取 x[col]。
计算部分结果 partial = val * x[col]。
将这个部分结果累加到输出向量 y 的 y[row] 位置。

由于多个线程可能同时更新同一个输出元素 y[row]，因此第4步需要使用原子操作来避免竞态条件。

以下是该内核函数的伪代码示例：

__global__ void spmv_coo_kernel(COOMatrix coo, float* x, float* y) {
    int i = blockIdx.x * blockDim.x + threadIdx.x; // 线程索引，对应非零元素索引
    if (i >= coo.num_nonzeros) return;

    int row = coo.row_indices[i];
    int col = coo.col_indices[i];
    float val = coo.values[i];

    float partial = val * x[col];
    atomicAdd(&y[row], partial); // 必须使用原子操作
}

COO格式的优缺点分析

基于SPMV应用，我们来总结COO格式的优缺点。

优点：

灵活性高：非零元素可以任意顺序存储，添加新元素非常容易（只需追加到数组末尾）。
可访问性：给定一个非零元素，可以立即获取其行和列索引。这使得跨非零元素并行化非常容易。
内存合并访问：线程按顺序访问 values、row_indices、col_indices 数组，访问模式是连续的，具有良好的内存合并性。
负载均衡：每个线程只处理一个非零元素，工作量相同，没有控制流发散。

缺点：

空间效率较低：需要为每个非零元素存储行和列索引，如果矩阵非常稀疏，这仍比稠密存储好，但相比其他格式有额外开销。
特定访问模式困难：给定一个行号，很难快速找到该行所有的非零元素（需要搜索整个 row_indices 数组）。因此，难以实现“每个输出元素一个线程”的并行策略。
需要原子操作：在SPMV中，由于多个线程可能写入同一内存位置，必须使用开销较大的原子操作。

压缩稀疏行格式

为了克服COO格式在SPMV中需要原子操作的缺点，我们接下来看压缩稀疏行格式。

CSR格式的原理

CSR格式按行组织非零元素。它存储三个数组：

values: 按行优先顺序存储所有非零值。同一行的非零值在数组中连续存放。
col_indices: 存储每个非零值对应的列索引，与 values 数组一一对应。
row_ptr: 这是一个长度等于 行数 + 1 的数组。row_ptr[i] 给出了第 i 行的第一个非零元素在 values 和 col_indices 数组中的起始索引。第 i 行的非零元素范围是 row_ptr[i] 到 row_ptr[i+1]-1。

使用之前同一个矩阵的例子，在CSR格式下：

values: [1, 7, 5, 3, 9, 2, 8] (注意：非零值已按行分组)
col_indices: [0, 1, 0, 2, 3, 1, 2]
row_ptr: [0, 2, 5, 7]。解读：第0行非零元素索引从0开始（到1结束）；第1行从索引2开始（到4结束）；第2行从索引5开始（到6结束）。

基于CSR格式的SPMV实现

在CSR格式下，实现SPMV最自然的并行化策略是：为输出向量的每一个元素（即矩阵的每一行）分配一个线程。

每个线程的执行步骤如下：

线程ID row 直接对应其负责的输出行号。
通过 row_ptr[row] 和 row_ptr[row+1] 确定该行非零元素在 values 和 col_indices 数组中的范围。
遍历该范围内的所有非零元素：对于每个元素，读取其值 val 和列号 col，计算 val * x[col]，并累加到一个局部变量 sum 中。
遍历结束后，将 sum 写入输出向量 y[row]。

由于每个输出元素 y[row] 只被一个线程写入，因此不需要原子操作。

以下是该内核函数的伪代码示例：

__global__ void spmv_csr_kernel(CSRMatrix csr, float* x, float* y) {
    int row = blockIdx.x * blockDim.x + threadIdx.x; // 线程索引，对应行号
    if (row >= csr.num_rows) return;

    float sum = 0.0f;
    int row_start = csr.row_ptr[row];
    int row_end = csr.row_ptr[row + 1];

    for (int j = row_start; j < row_end; ++j) {
        int col = csr.col_indices[j];
        float val = csr.values[j];
        sum += val * x[col];
    }
    y[row] = sum; // 无需原子操作
}

CSR格式的优缺点分析

现在，我们来分析CSR格式在SPMV应用中的表现。

优点：

空间效率更高：相比COO，它用更小的 row_ptr 数组替代了 row_indices 数组，节省了存储空间。
可访问性：给定一个行号，可以快速定位到该行所有非零元素。这使得按行并行化（如SPMV）非常高效。
无需原子操作：在SPMV中，每个输出元素由唯一线程计算和写入，消除了同步开销。
输出写入合并：线程按行号顺序写入输出向量 y，写入操作是合并的。

缺点：

灵活性差：添加或删除非零元素可能涉及移动大量数据和更新 row_ptr 数组，成本较高。
特定访问模式困难：给定一个非零元素，很难找到其行号（需在 row_ptr 中搜索）。给定一个列号，也很难找到该列所有非零元素。因此，难以实现“每个非零元素一个线程”的并行策略。
内存访问未合并：线程在遍历其负责的行时，对 values 和 col_indices 数组的访问是跳跃式的，不同线程访问的位置不相邻，导致内存访问无法合并。
负载不均衡/控制流发散：不同行的非零元素数量可能差异很大，导致线程间工作量不均。有些线程很快完成，有些则要处理很多元素，造成GPU利用率下降和控制流发散。

扩展：CSC格式
与CSR对应的是压缩稀疏列格式。CSC按列存储非零元素，并有一个 col_ptr 数组。它使得按列访问非零元素变得容易，适用于需要列遍历的运算。

性能对比与总结

在本节课中，我们一起学习了两种基础的稀疏矩阵存储格式：COO和CSR，并基于SPMV运算分析了它们的实现与特性。

通过实际代码测试，我们可能会发现，对于某些矩阵，即使COO需要原子操作，其性能也可能优于CSR。这是因为COO格式对矩阵数据的访问是完美合并的，而CSR格式的未合并访问和控制流发散带来的性能损失，有时会超过原子操作的开销。此外，CSR格式在空间效率上通常优于COO。

两种格式的核心权衡在于：

COO：灵活性和合并访问更优，但需要原子操作。
CSR：空间效率更高且无需原子操作，但存在未合并访问和控制流发散。

显然，我们希望能有一种格式，既能像CSR一样避免原子操作，又能像COO一样实现内存合并访问，同时还能改善负载均衡。这正是我们下节课将要探讨的内容。下一次，我们将介绍ELLPACK和锯齿状对角线格式，它们旨在解决CSR格式在内存访问模式和负载均衡方面的缺陷。

GPU计算：第17讲：稀疏矩阵计算（ELL与JDS格式）🚀

在本节课中，我们将继续学习稀疏矩阵计算，重点介绍两种新的稀疏矩阵存储格式：ELL（ELLPACK）格式和JDS（Jagged Diagonal Storage，锯齿状对角线存储）格式。我们将探讨它们的设计原理、优缺点，并分析它们在稀疏矩阵-向量乘法（SpMV）中的性能表现。

上一讲我们介绍了稀疏矩阵的基本概念以及COO和CSR存储格式。本节中，我们将深入探讨旨在优化内存访问和控制流一致性的ELL和JDS格式。

稀疏矩阵计算回顾

稀疏矩阵是指矩阵中绝大多数元素为零的矩阵。通过不存储零元素，我们可以节省内存；通过不加载和计算零元素，我们可以节省内存带宽和计算时间。选择不同的存储格式会影响计算性能，尤其是在GPU上执行SpMV运算时。

设计存储格式时，我们主要考虑以下几个因素：

空间效率：节省存储空间的程度。
灵活性：添加、删除或重排元素的难易程度。
可访问性：查找特定行、列或非零元素的便利性。
内存访问模式：是否支持合并内存访问。
负载均衡：在GPU上下文中，是否能最小化控制流发散。

我们以SpMV运算作为评估这些格式性能的基准。SpMV的公式如下：
y = A * x
其中 A 是稀疏矩阵，x 是稠密输入向量，y 是稠密输出向量。

ELL（ELLPACK）格式

ELL格式旨在解决CSR格式中内存访问不合并的问题。其核心思想是通过填充（Padding）使每一行具有相同数量的元素，然后按列主序存储，从而实现线程间的合并访问。

ELL格式的构造

以下是ELL格式的构造步骤：

按行分组：像CSR一样，将每行的非零元素及其列索引分组。
填充行：对每行进行填充，使所有行具有相同数量的“元素”（包括实际非零元和填充的无效元素）。填充数量由矩阵中最大非零元数/行决定。
列主序存储：将填充后的二维数组（行 x 最大非零元数）按列主序展开成一维数组进行存储。这意味着所有行的第一个元素连续存储，接着是所有行的第二个元素，依此类推。

以下是一个矩阵转换为ELL格式的示例：

原矩阵：
[1, 0, 0, 7]
[0, 5, 3, 0]
[9, 0, 0, 2]
[0, 0, 8, 0]

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/212d6a91d61aa14c0d3f4dbf0b2d76a9_35.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/212d6a91d61aa14c0d3f4dbf0b2d76a9_36.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/212d6a91d61aa14c0d3f4dbf0b2d76a9_38.png)

最大非零元数/行 = 3
填充并转置后（按列主序存储）：
值数组:    [1, 5, 9, 8,  7, 3, 2, 0,  0, 0, 0, 0]
列索引数组: [0, 1, 0, 2,  3, 2, 3, -1, -1,-1,-1,-1]

（注：-1或特定值表示填充位置）

SpMV在ELL格式上的并行化

在ELL格式上实现SpMV的典型方法是为每一行分配一个线程。每个线程负责计算一个输出元素 y[row]。

每个线程的执行流程如下：

线程根据其ID确定负责的行 row。
线程循环遍历该行的元素（循环上限为该行的实际非零元数，以避免计算填充元素）。
在每次迭代 i 中，计算在值/列索引数组中的位置：index = i * num_rows + row。
读取 values[index] 和 col_idx[index]。
执行计算 y[row] += values[index] * x[col_idx[index]]。

由于 index 的计算依赖于连续的 row 值，且 row 直接映射到连续的线程ID，因此对 values 和 col_idx 数组的访问是合并的。

ELL格式的优缺点分析

上一节我们介绍了ELL格式的原理和实现，本节中我们来看看它的优缺点。

优点：

内存合并访问：SpMV实现了合并内存访问，优于CSR。
灵活性适中：对于未满的行（即有填充位的行），可以较容易地添加新非零元（替换填充位）。
可访问性较好：给定行，容易找到其所有非零元；给定非零元在数组中的索引，也容易反推出其所在行和列。

缺点：

空间效率低：填充可能导致大量空间浪费，特别是当某行非零元远多于其他行时。
控制流发散：与CSR类似，不同行的非零元数量不同，线程循环次数不同，导致控制流发散。
转换开销：从其他格式转换为ELL格式需要成本。

我们之前实现的COO和CSR内核性能对比显示，对于特定矩阵，COO可能因更好的合并性而优于CSR。ELL格式的目标是结合两者的优点。

混合ELL/COO格式

ELL格式的一个主要问题是当矩阵中有一行或几行特别长（非零元特别多）时，会导致其他行出现大量填充，严重浪费空间。

为了解决这个问题，可以采用混合ELL/COO格式：

设定一个阈值K（例如，大多数行的非零元数）。
对每行的前K个非零元使用ELL格式存储。
将超过K个的非零元存储在额外的COO格式结构中。

这样既保留了ELL格式的合并访问优势，又通过COO格式收纳了“超长”部分，减少了填充浪费，也增加了添加元素的灵活性（超长部分可加入COO列表）。

JDS（锯齿状对角线存储）格式

JDS格式的主要设计目标是减少SpMV中的控制流发散。它通过按行非零元数对行进行排序，并调整存储布局来实现。

JDS格式的构造

以下是JDS格式的构造步骤：

按行非零元数排序：根据每行非零元的数量对矩阵行进行降序排序。行号信息保存在一个perm排列数组中。
按“锯齿对角线”存储：将排序后的矩阵，将其非零元按“对角线”方式提取并连续存储。具体来说，将所有行的第一个非零元连续存储，接着是所有行的第二个非零元，依此类推。这形成了一个锯齿状的布局。
迭代指针数组：需要一个iter_ptr数组，来指示第i个“对角线”（即所有行的第i个非零元）在值数组中的起始位置。

示例：对之前矩阵的行按非零元数排序后，再按JDS规则存储。

原矩阵行（已排序）：
行1: [9, 0, 0, 2]  (2个非零元)
行3: [0, 5, 3, 0]  (2个非零元)
行0: [1, 0, 0, 7]  (2个非零元)
行2: [0, 0, 8, 0]  (1个非零元)

JDS存储：
值数组:      [9, 5, 1, 8,  2, 3, 7]
列索引数组:   [0, 1, 0, 2,  3, 2, 3]
排列数组 perm: [1, 3, 0, 2] // 排序后第0行是原第1行...
迭代指针 iter_ptr: [0, 4, 7] // 第0个对角线起始于0，第1个起始于4...

SpMV在JDS格式上的并行化

在JDS格式上实现SpMV，同样为排序后的每一行分配一个线程。

每个线程的执行流程如下：

线程ID对应排序后的行 sorted_row。
通过 perm[sorted_row] 找到原始行号 orig_row，该线程负责输出 y[orig_row]。
线程循环遍历“对角线”。在迭代 i，所有活跃线程访问第 i 个对角线。
线程检查 sorted_row 是否小于 iter_ptr[i+1] - iter_ptr[i]（即该行在第 i 个对角线上是否有元素）。如果有，则计算索引 index = iter_ptr[i] + sorted_row。
读取 values[index] 和 col_idx[index] 并进行累加计算。
随着迭代 i 增加，非零元数少的行（排序在后面的行）会提前完成计算，线程变为不活跃。由于行已排序，不活跃的线程是连续位于末端的，这极大地减少了控制流发散。

JDS格式的优缺点分析

优点：

优秀的控制流一致性：行按长度排序，使得线程组（如warp）同时变得不活跃，最小化了控制流发散。
良好的内存合并访问：每个迭代中，活跃线程访问连续的内存位置。
无填充开销：相比ELL，空间效率更高。

缺点：

灵活性差：添加/删除非零元会改变行长度和排序，更新数据结构成本高。
可访问性复杂：给定原矩阵行号，需要查找perm数组才能定位数据。给定列号，很难找到所有相关非零元。
额外存储开销：需要存储perm和iter_ptr数组。

性能对比显示，对于示例矩阵，JDS内核的性能优于ELL、CSR和COO内核，这主要归功于其出色的内存访问合并性和控制流一致性。

总结

本节课中我们一起学习了两种重要的稀疏矩阵存储格式：ELL和JDS。

ELL格式通过填充和列主序存储，优化了SpMV中的内存合并访问，但可能因填充导致空间浪费。
JDS格式通过按行非零元数排序和锯齿对角线存储，同时优化了内存合并访问和控制流一致性，通常能获得最佳性能，但数据结构更复杂且灵活性差。

需要强调的是，没有一种存储格式在所有情况下都是最优的。格式的选择取决于具体的稀疏矩阵模式、目标计算内核（不仅是SpMV）以及硬件特性。ELL和JDS格式为我们提供了在追求高性能SpMV计算时的有力工具。

GPU计算：第18讲：图处理

概述

在本节课中，我们将学习如何在GPU上进行图处理。我们将首先回顾稀疏矩阵计算，然后探讨如何利用稀疏矩阵的存储格式来表示图。接着，我们将介绍两种主要的图处理并行化方法：顶点中心法和边中心法，并以广度优先搜索为例，详细分析这两种方法的实现、性能差异及其适用场景。

回顾：稀疏矩阵存储格式

上一节我们介绍了稀疏矩阵计算，本节我们来看看如何将这些知识应用于图处理。我们之前讨论过四种稀疏矩阵存储格式：COO、CSR、ELL和JDS。

COO格式：存储非零元素的行索引、列索引和值。
CSR格式：存储行指针、列索引和值，便于按行访问。
ELL格式：通过填充使所有行具有相同数量的非零元素，并按列主序存储，以实现内存访问合并。
JDS格式：对行按长度排序并存储排序信息，然后按列主序存储，以减少控制流分歧并实现内存访问合并。

图的表示

图由顶点和边组成。我们可以用邻接矩阵来逻辑上表示一个图。邻接矩阵通常非常稀疏，因此可以使用我们学过的稀疏矩阵存储格式来表示图。

为了简化讨论，我们专注于无权图（所有边权重为1）和无向图（邻接矩阵对称）。这意味着在存储时，我们通常不需要存储值数组，只需关注非零元素的位置。

以下是使用不同稀疏格式表示图的方法：

COO格式表示图：相当于存储边的列表。包含两个数组：源顶点数组和目的顶点数组。对于加权图，还需要第三个权重数组。
- 源顶点数组：[0, 0, 1, 2, 3, 3]
- 目的顶点数组：[1, 2, 0, 0, 1, 2]
CSR格式表示图：便于查找某个顶点的所有邻居。包含源指针数组和目的顶点数组。
- 源指针数组：[0, 2, 5, 6] （表示顶点i的邻居在目的数组中的起始位置）
- 目的顶点数组：[1, 2, 0, 2, 3, 0]

对于无向图，CSR和CSC格式是等价的。图的性质（如顶点度数的分布）会影响并行化策略和存储格式的选择。

图处理的并行化方法

主要有两种并行化图处理的方法：顶点中心法和边中心法。

顶点中心法：为图中的每个顶点分配一个线程。该线程负责执行与该顶点相关的操作，通常需要访问该顶点的邻居列表。因此，CSR或CSC格式更适合这种方法，因为它们便于快速查找某个顶点的所有邻居。
边中心法：为图中的每条边分配一个线程。该线程负责执行与该边相关的操作，通常需要访问该边的源顶点和目的顶点。因此，COO格式更适合这种方法，因为它便于快速查找每条边的端点。

在某些算法中（例如三角形计数），可能需要同时查找边的端点以及端点的邻居，这时就需要混合表示，同时维护COO和CSR两种格式。

案例研究：广度优先搜索

我们将以广度优先搜索为例，具体分析顶点中心法和边中心法的实现。BFS的目标是找到从源顶点到图中所有其他顶点的距离（层级）。

BFS算法从源顶点（层级0）开始，逐层向外探索：

访问当前层级所有顶点的未访问邻居。
将这些邻居标记为下一层级。
重复此过程，直到没有新的顶点被访问。

顶点中心法 - 自顶向下方法

在自顶向下方法中，每次迭代为每个顶点启动一个线程。线程检查其负责的顶点是否属于上一层级。如果是，则该线程遍历该顶点的所有邻居，并将未访问的邻居标记为当前层级。

这种方法类似于在BFS树中，父顶点寻找其子顶点。

以下是该方法的CUDA内核伪代码核心逻辑：

unsigned int vertex = threadIdx.x + blockIdx.x * blockDim.x;
if (vertex < num_vertices) {
    if (level[vertex] == current_level - 1) { // 顶点属于上一层级
        for (int edge = csr.src_ptr[vertex]; edge < csr.src_ptr[vertex+1]; edge++) {
            unsigned int neighbor = csr.dst[edge];
            if (level[neighbor] == INF) { // 邻居未访问
                level[neighbor] = current_level; // 标记为当前层级
                visited_flag = 1; // 设置已访问新顶点标志
            }
        }
    }
}

缺点：存在控制流分歧。只有少数顶点（上一层级的顶点）的线程会执行大量工作（遍历邻居），而其他线程立即退出。此外，对于高度数顶点，其线程需要遍历很长的邻居列表。

顶点中心法 - 自底向上方法

在自底向上方法中，每次迭代同样为每个顶点启动一个线程。但线程检查其负责的顶点是否尚未被访问。如果是，则该线程遍历该顶点的所有邻居，检查是否有邻居属于上一层级。一旦找到，就将自身顶点标记为当前层级并提前退出循环。

这种方法类似于潜在的子顶点通过检查其父顶点来确定自己是否属于当前层级。

以下是该方法的CUDA内核伪代码核心逻辑：

unsigned int vertex = threadIdx.x + blockIdx.x * blockDim.x;
if (vertex < num_vertices) {
    if (level[vertex] == INF) { // 顶点未访问
        for (int edge = csr.src_ptr[vertex]; edge < csr.src_ptr[vertex+1]; edge++) {
            unsigned int neighbor = csr.dst[edge];
            if (level[neighbor] == current_level - 1) { // 邻居属于上一层级
                level[vertex] = current_level; // 将自身标记为当前层级
                visited_flag = 1; // 设置已访问新顶点标志
                break; // 提前退出
            }
        }
    }
}

优点：对于高度数顶点，一旦找到一个已访问的邻居，线程就可以提前退出，减少了不必要的遍历。在后期迭代中，当大多数顶点已被访问时，只有少数未访问顶点的线程需要工作，效率较高。
缺点：在早期迭代中，由于已访问顶点很少，大多数未访问顶点的线程会遍历所有邻居却一无所获，造成大量冗余工作。

方向优化BFS

结合自顶向下和自底向上方法的优点：

早期迭代：使用自顶向下方法。此时已访问顶点（前沿）很小，工作量大。
后期迭代：切换为自底向上方法。此时未访问顶点集合变小，自底向上方法更高效。
这种自适应策略通常能获得最佳整体性能。

边中心法

在边中心法中，每次迭代为图中的每条边启动一个线程。线程检查其负责的边的源顶点是否属于上一层级，并且目的顶点是否未被访问。如果条件满足，则将目的顶点标记为当前层级。

以下是该方法的CUDA内核伪代码核心逻辑：

unsigned int edge = threadIdx.x + blockIdx.x * blockDim.x;
if (edge < num_edges) {
    unsigned int src = coo.src[edge];
    unsigned int dst = coo.dst[edge];
    if (level[src] == current_level - 1 && level[dst] == INF) {
        level[dst] = current_level;
        visited_flag = 1;
    }
}

优点：没有循环，每个线程只处理一条边，控制流分歧小，并行度非常高。
缺点：每次迭代都启动所有边的线程，即使大部分边不参与当前迭代，存在冗余工作。

数据结构的影响

最佳BFS并行化策略高度依赖于图的结构：

高度数图：顶点度数差异大。
- 自底向上和边中心法表现更好，因为它们能更好地处理负载不平衡（高度数顶点）。
- 自顶向下方法表现较差，因为处理高度数顶点的线程会成为瓶颈。
低度数图：顶点度数相对均匀。
- 自顶向下方法表现更好，因为每次迭代的前沿规模可控，工作负载均衡。
- 典型的例子是道路网络图。

因此，没有一种算法在所有数据集上都是最优的，需要根据图的特点进行选择。

BFS与SpMV的相似性

图处理与稀疏线性代数之间存在紧密联系。BFS操作，特别是自底向上的顶点中心法，与稀疏矩阵-向量乘法在计算模式上非常相似：

SpMV：for each row i: for each nonzero (i, j): output[i] += matrix[i][j] * input[j]
BFS（自底向上）：for each vertex v: for each neighbor u of v: if (level[u] == L-1) then level[v] = L

许多图算法可以表述为稀疏线性代数操作。这样做的优势是可以利用成熟的高性能稀疏线性代数库。缺点是，这种表述有时并非解决特定图问题的最直接或最高效的方法。

冗余工作与前沿优化

本节课讨论的所有方法都有一个共同点：每次迭代都检查所有顶点或边，无论它们是否与当前迭代相关。

优点：实现简单，并行度高，线程间无需同步。
缺点：产生大量冗余工作，许多线程空转。

下次课程将介绍另一种策略：前沿优化。该策略旨在只处理与当前迭代相关的顶点或边（即“前沿”）。这需要动态构建和维护这些前沿集合，虽然减少了冗余工作，但引入了线程间的同步和通信开销。这是一种在并行度、冗余工作和同步开销之间的权衡。

总结

本节课我们一起学习了GPU上的图处理。我们回顾了稀疏矩阵格式如何用于表示图，并深入探讨了两种主要的并行化范式：顶点中心法和边中心法。通过以广度优先搜索为例，我们分析了自顶向下、自底向上、方向优化以及边中心等具体实现方法的原理、代码结构和性能特点。我们了解到，图的结构（如度数分布）显著影响算法性能，因此需要根据数据特征选择合适的方法。最后，我们指出了当前方法存在的冗余工作问题，并预告了下节课将探讨的前沿优化策略。

GPU计算：19：图处理（第二部分）🚀

在本节课中，我们将继续学习图处理算法。我们将重点探讨如何通过前沿队列和私有化等技术来减少冗余计算，并优化BFS（广度优先搜索）在GPU上的实现性能。

概述

上一节我们介绍了图处理的基础知识，包括图的表示方法（如COO、CSR格式）以及BFS的几种并行实现（顶点中心化、边中心化）。我们观察到，之前的实现方法在每次迭代中都会检查所有顶点或边，这导致了大量的冗余工作。

本节中，我们将学习一种更高效的方法：前沿队列。这种方法的核心思想是，在每次迭代中，我们只处理那些在上一层级中被访问过的顶点（即“前沿”），从而避免为不相关的顶点启动线程。虽然这引入了更多的同步开销，但能显著减少计算冗余。

前沿队列方法

上一节我们介绍了BFS的几种实现，它们都存在检查冗余的问题。本节中我们来看看如何通过前沿队列来优化。

我们的目标是：在每一层级，只检查属于上一层级的顶点。为此，我们将在每个层级构建一个队列（称为“前沿”），其中包含在该层级被访问的顶点。在下一层级，我们只启动线程来处理这个前沿队列中的顶点。

以下是该方法的步骤：

初始化前沿队列，仅包含源顶点。
对于当前前沿队列中的每个顶点，检查其所有邻居。
如果邻居未被访问，则将其标记为已访问（属于下一层级），并加入下一个前沿队列。
将下一个前沿队列设置为当前前沿，重复步骤2-3，直到没有新的顶点可访问。

这种方法减少了线程的浪费启动，但需要同步操作来安全地向共享队列中添加顶点。

基础实现与原子操作

让我们开始实现基础的前沿队列方法。首先，我们需要在每次迭代中启动与前沿队列中顶点数量相等的线程。

以下是核心步骤的伪代码描述：

// 假设有以下参数：
// prev_frontier: 存储上一层级顶点的数组
// num_prev: 上一层级顶点数量
// curr_frontier: 用于存储当前层级（下一层级）顶点的数组
// num_curr: 当前层级顶点数量的指针（初始为0）
// level: 存储每个顶点所属层级的数组
// current_level: 当前正在处理的层级号

// 每个线程处理prev_frontier中的一个顶点
int tid = blockIdx.x * blockDim.x + threadIdx.x;
if (tid < num_prev) {
    int vertex = prev_frontier[tid];
    // 遍历该顶点的所有出边（邻居）
    for (int edge = graph.src_ptr[vertex]; edge < graph.src_ptr[vertex+1]; edge++) {
        int neighbor = graph.dst[edge];
        // 关键步骤：使用原子比较交换操作来检查并标记邻居
        // 如果neighbor的层级还是未访问状态（例如UINT_MAX），则将其设置为current_level
        unsigned int old_level = atomicCAS(&level[neighbor], UINT_MAX, current_level);
        if (old_level == UINT_MAX) {
            // 成功标记该邻居为本层级，将其加入当前前沿队列
            int pos = atomicAdd(num_curr, 1); // 原子地获取队列中的写入位置
            curr_frontier[pos] = neighbor;
        }
    }
}

代码解释：

atomicCAS (Compare And Swap): 这是一个原子操作。它比较 level[neighbor] 的当前值是否等于 UINT_MAX（表示未访问）。如果是，则将其原子地设置为 current_level 并返回旧值 UINT_MAX；如果不是，则不做修改，直接返回当前值。这确保了即使多个线程同时发现同一个未访问邻居，也只有一个线程能成功将其标记并加入队列，从而避免了重复添加。
atomicAdd: 这是一个原子加法操作。多个线程需要安全地向 curr_frontier 队列末尾添加元素。atomicAdd(num_curr, 1) 会原子地将 num_curr 的值加1，并返回加1之前的值，这个返回值就是该线程可以安全写入的队列索引。

通过使用这些原子操作，我们解决了向共享队列添加元素时的竞态条件问题。

优化：队列私有化

基础实现中，所有线程都通过原子操作竞争同一个全局队列计数器 (num_curr)，这在高并发时会导致严重的序列化瓶颈和访问延迟。

为了解决这个问题，我们可以应用私有化策略。思路是让每个线程块维护自己的私有队列（放在共享内存中），线程只向本线程块的私有队列添加元素。待所有线程处理完后，再由每个线程块一次性将其私有队列的内容合并到全局队列中。

以下是优化步骤：

在共享内存中为每个线程块声明一个私有队列 local_queue 及其计数器 local_count。
线程在处理顶点时，使用原子操作向 local_count 和 local_queue 添加元素（共享内存原子操作速度更快）。
设置一个私有队列的大小上限。如果写满，则溢出的元素直接写入全局队列（回退机制）。
处理结束后，通过线程块内同步 (__syncthreads())。
由线程块内的一个线程（如threadIdx.x == 0）使用一次原子操作，向全局计数器申请一块连续空间，大小为该线程块私有队列的元素数量。
该线程将申请到的起始索引存入共享变量，再次同步后，线程块内所有线程协作将私有队列中的数据拷贝到全局队列的指定位置。

这种方法的优点是：

将全局原子操作冲突分散到各个线程块内的共享内存原子操作，冲突大大减少。
从私有队列向全局队列写入数据时，内存访问模式是合并的，效率更高。

进一步优化：减少内核启动开销

在前沿队列方法中，我们为图的每一层级都启动了一次内核。内核启动本身有一定开销，并且需要在CPU和GPU之间传输前沿队列的大小等信息。

我们可以观察到，在图搜索的开始和结束阶段，前沿队列通常很小（可能只需要一个线程块就能处理）。对于这些连续的小层级，我们可以尝试在单个内核启动中处理多个层级。

优化思路：

在内核中，不仅处理当前前沿队列来生成下一队列，还可以检查下一队列的大小。
如果下一队列的大小也足够小（例如，小于一个线程块能处理的最大顶点数），那么我们可以直接在同一个线程块内，使用 __syncthreads() 进行同步，然后继续处理下一层级的顶点。
这样就避免了为这个小层级再次启动内核的开销。

这种优化将多个小的网格启动合并为一次，特别适用于社交网络等图中常见的长尾层级分布，可以进一步减少总体运行时间。

总结

本节课中我们一起学习了如何优化GPU上的图处理算法，特别是BFS。

前沿队列方法：通过只处理活跃顶点（前沿）来消除冗余计算，核心是使用原子操作 (atomicCAS, atomicAdd) 来安全地维护队列。
队列私有化：通过让每个线程块维护私有队列，减少了全局原子操作的冲突，并改善了内存访问模式，从而提升了性能。
减少内核启动：通过在内核中处理连续的小型前沿队列，合并多个层级的计算，减少了内核启动和CPU-GPU通信的开销。

这些优化技术体现了GPU编程中常见的权衡：通过增加同步和更复杂的逻辑来减少总计算量和开销，从而在合适的场景下获得显著的性能提升。理解图的特性和硬件行为对于选择和应用这些优化至关重要。

GPU计算：01：课程介绍与动机

在本节课中，我们将要学习GPU计算的基本概念、其出现的历史背景，以及它与传统CPU在设计哲学上的根本区别。我们将从计算机技术的发展历史讲起，理解为什么并行计算和GPU变得如此重要。

计算机技术的发展与摩尔定律

首先，我们来回顾一下计算机技术在过去半个世纪的发展历程。计算机科学和电子工程专业的学生应该对此非常熟悉。你们还记得摩尔定律吗？

摩尔定律并非物理定律，而是英特尔联合创始人戈登·摩尔提出的一个预测。他预测，单位面积集成电路上可容纳的晶体管数量，大约每18到24个月便会增加一倍。

这个预测在过去的半个世纪里基本成立，半导体行业也以此为目标进行发展。然而，正如一位同学指出的，这个趋势正在走向终结。因为晶体管尺寸已经小到难以继续缩小的物理极限。

晶体管数量的持续翻倍带来了一个好处：晶体管越小，开关速度就越快。这意味着我们可以提高处理器的时钟频率。因此，在2005年之前，处理器频率也遵循着类似的翻倍趋势。如果你在2005年之前购买了一台新电脑，即使运行相同的程序，速度也会比两年前的旧电脑快一倍，这主要归功于频率的提升。

频率停滞与并行计算的兴起

上一节我们介绍了频率提升带来的“免费午餐”，但这种情况在2005年左右发生了变化。我们继续获得更多的晶体管，但无法再让它们更快地开关了。有谁知道为什么吗？

原因在于功耗和散热。让晶体管更快地开关会产生更多热量。如果散热技术无法及时带走这些热量，CPU就会过热损坏。因此，在2005年左右，处理器频率的提升停滞了。

这意味着，仅仅依靠购买更新的电脑，不再能自动获得程序运行速度的翻倍提升。当然，单线程性能并未完全停滞，它仍在缓慢提升，主要得益于编译器技术的进步，以及利用额外晶体管实现的更复杂的硬件技术，例如分支预测和乱序执行。

由于单线程性能提升放缓，我们获得了更多晶体管却无法显著提升性能，那么人们开始用这些额外的晶体管做什么呢？

答案就是增加核心数量。硬件设计者开始在处理器中集成更多核心。这标志着并行计算开始成为主流。软件也必须跟上硬件的步伐，开发者需要重写程序以利用并行能力。免费的“性能午餐”就此结束。

系统设计的两种取向：延迟导向与吞吐量导向

在深入探讨CPU和GPU的区别之前，我们先来了解两种通用的系统设计思路。这不仅适用于计算机系统，也适用于任何旨在完成任务的系统设计。

这两种设计思路分别是延迟导向设计和吞吐量导向设计。

延迟导向设计旨在最小化完成单个任务所需的时间。一个典型的例子是汽车。汽车的设计目标是尽快将一个人从A点送到B点。

吞吐量导向设计则旨在最大化在给定时间内完成的任务数量，即使每个独立任务耗时更长。一个典型的例子是公共汽车。公交车一次可以运送很多人，虽然每个人到达目的地的时间可能比开车要长。

选择哪种设计取决于任务规模。运送1到4个人，汽车是高效的选择。但若要运送20个人，公交车只需一趟，而汽车需要往返多次，此时公交车就成为了更具吸引力的解决方案。

CPU与GPU：延迟与吞吐量的设计哲学

理解了通用设计思路后，我们将其应用到处理器设计上。CPU是延迟导向设计的典型代表，而GPU则是吞吐量导向设计的典型代表。

换句话说，如果你希望以最快速度完成一个任务，应该使用CPU。但如果你希望在特定时间内完成许多小任务，GPU就成为了更佳选择。

以下是这两种设计哲学如何具体影响CPU和GPU设计的几个方面：

1. 算术逻辑单元

CPU：拥有少量但非常强大的ALU。这些ALU占用大量芯片面积，经过高度优化，旨在将单个算术运算（如加法、乘法）的延迟降到最低。
GPU：拥有大量小型、简单的ALU。每个ALU执行运算的延迟较长，但通过大量ALU并行工作以及流水线技术，可以实现极高的整体运算吞吐量。

2. 缓存

CPU：配备大容量缓存。目的是将访问慢速DRAM内存的高延迟操作，转化为访问快速缓存的低延迟操作，从而降低单线程的访存延迟。
GPU：缓存容量较小。这意味着缓存命中率较低，更多访问需要去到延迟更高的内存。但这样做的好处是可以将更多的芯片面积用于计算单元（ALU），而不是缓存。

3. 控制逻辑

CPU：采用复杂昂贵的控制逻辑来降低延迟，例如分支预测、数据前递、乱序执行等。这些技术需要大量额外的硬件支持。
GPU：控制逻辑相对简单。运算分支、数据冲突等带来的延迟更高，但这是为了节省芯片面积，以容纳更多计算单元，从而换取高吞吐量。

4. 多线程与延迟隐藏
由于CPU操作延迟较低，流水线中的“气泡”（空闲周期）相对较少。为了进一步隐藏这些短暂的延迟，CPU采用适度规模的多线程（例如，每个核心同时运行2个线程，即超线程技术）。当一个线程的指令因等待而无法执行时，可以执行另一个线程的指令来填充流水线。

GPU则不同。其ALU运算、内存访问和控制决策的延迟都非常高，流水线中容易出现大量气泡。为了隐藏这些高延迟，GPU采用了大规模的多线程。单个GPU核心上会同时调度并准备执行数十甚至上百个线程。当某个线程需要等待时，硬件可以立刻从众多就绪线程中挑选另一个线程的指令来执行，从而保持计算单元的忙碌，实现高吞吐量。

5. 时钟频率

CPU：通常具有较高的时钟频率。
GPU：由于集成了海量的计算单元，功耗已经很大，因此通常运行在相对较低的时钟频率以控制总功耗和散热。

GPU的起源与发展

上一节我们对比了CPU和GPU的设计哲学，那么GPU最初是从何而来的呢？

GPU的全称是图形处理单元。它最初是为图形渲染而设计的。图形工作负载有一个独特的特点：高度并行。例如，渲染屏幕上的成千上万个像素，这些像素的计算在很大程度上是相互独立的。因此，GPU从诞生之初就被设计成一种拥有大规模并行能力和高吞吐量的架构。

后来，人们意识到，尽管GPU是为图形设计的，但它同样适用于其他需要高吞吐量的计算任务。这催生了通用GPU计算的趋势。

在2007年之前，编程GPU只能通过特定的图形API（如OpenGL、Direct3D）。想要进行通用计算的研究人员，必须将自己的计算问题“伪装”成图形渲染问题，过程非常繁琐。

2007年，英伟达发布了CUDA平台。这是一个革命性的变化，它提供了一个直接的编程接口，让开发者能够像编写通用程序一样为GPU编写代码。同时，GPU硬件架构也进行了相应扩展，以更好地支持通用计算。这标志着GPU计算时代的真正开始。

此后，GPU在科学计算、深度学习等领域被广泛采用。根据2020年11月的全球超级计算机TOP500榜单，前10名中有6台使用了GPU加速。在能效榜单（Green500）上，前10名中有7台是GPU加速的超级计算机。虽然单个GPU功耗可能更高，但由于其强大的并行计算能力能极大缩短任务完成时间，因此从完成单位计算量所消耗的总能量来看，GPU往往更具能效优势。

为什么是GPU取得了成功？

并行处理器有多种设计可能，为什么最终是GPU取得了巨大成功？这背后有技术和市场的双重原因。

芯片设计制造成本极高，需要巨大的销量来分摊成本。当并行计算成为主流时，GPU已经拥有了一个庞大的现有市场——游戏产业。这使得GPU厂商在向通用计算领域拓展时，拥有其他潜在竞争对手所不具备的规模优势和成本分摊能力。

根据2019年的市场划分，GPU收入中58%来自游戏，数据中心占24%（其中包含科学计算和深度学习）。如果一家公司只为科学计算市场设计一款专用并行芯片，其销量将难以支撑高昂的研发制造成本。而GPU厂商则可以依托游戏市场的稳定收入，来发展和推广其通用计算能力。

这也带来一个有趣的考量：GPU架构的改进必须兼顾所有主要市场。一项可能大幅提升科学计算性能但会损害游戏性能的改进，很可能不会被采纳，因为游戏仍是最大的客户基础。

课程总结与教材

本节课中，我们一起学习了：

摩尔定律的终结如何促使计算范式从提升单核频率转向增加核心数量，即并行计算。
系统设计的两种基本哲学：延迟导向（如CPU）与吞吐量导向（如GPU），以及它们在ALU、缓存、控制逻辑和多线程策略上的具体体现。
GPU的起源（图形处理）及其向通用计算演进的关键里程碑（CUDA的发布）。
GPU在当今高性能与高效能计算中的重要地位，以及其取得成功的市场因素。

本课程将主要关注如何使用GPU进行通用目的编程，重点案例将来自科学计算和深度学习领域。

推荐教材：《Programming Massively Parallel Processors: A Hands-on Approach》 by David Kirk and Wen-mei Hwu。课程内容将与教材章节对应，建议课后阅读以加深理解。

注：后续课程将开始具体学习GPU架构细节和编程方法。

GPU计算：第20讲：线程束内同步 🧵

在本节课中，我们将学习GPU编程中的一项高级技术：线程束内同步。我们将探讨如何利用同一线程束内线程的特殊关系，通过共享数据和投票机制，实现比传统的块级同步更高效的协作。

课程回顾与概述

上一节我们介绍了各种并行模式及其优化。本节中，我们将开始探讨GPU程序员用于优化应用程序的一些高级特性，特别是线程束内的同步。

线程束是SM中的调度单位，通常包含32个线程。同一线程束内的线程遵循SIMD模型执行，这意味着它们几乎同时执行相同的指令。这种紧密的执行关系使得线程束内线程间的同步比使用 __syncthreads() 进行整个线程块的同步要高效得多。

CUDA提供了两类主要的线程束内同步内置函数：

数据洗牌指令：允许线程直接共享寄存器值，无需经过共享内存。
投票指令：允许线程就某个条件进行协作和判断。

数据洗牌指令 🔀

数据洗牌指令使同一线程束内的线程能够直接共享数据。传统上，线程块内共享数据需要通过共享内存和 __syncthreads()。如果只需在线程束内共享数据，使用洗牌指令可以避免这些开销，从而更快。

以下是几种主要的洗牌指令变体：

__shfl_sync(mask, var, srcLane): 从指定通道（线程）srcLane 复制数据。
__shfl_up_sync(mask, var, delta): 从相对ID较低（低 delta）的线程复制数据。
__shfl_down_sync(mask, var, delta): 从相对ID较高（高 delta）的线程复制数据。
__shfl_xor_sync(mask, var, laneMask): 通过自身通道ID与 laneMask 进行按位异或的结果来确定源通道。

这些函数名中的 _sync 表示它们会同步参与操作的线程，这在支持线程独立执行路径的现代GPU架构中很重要。参数 mask 指定了参与此次洗牌操作的线程。

在归约模式中应用洗牌指令

我们曾在归约和扫描模式中看到，一个线程需要读取另一个线程在前一迭代中产生的值。在传统的基于共享内存的实现中，每一步之后都需要 __syncthreads()。

当迭代进行到只剩一个线程束时，线程读取的数据都来自同一线程束内的其他线程。此时，我们可以用洗牌指令替代共享内存和全局同步。

以下是修改归约内核以使用洗牌指令的关键步骤：

在共享内存中进行部分归约：像往常一样进行归约树操作，但当步长 stride 减小到等于线程束大小时停止。

for (int stride = blockDim.x / 2; stride >= WARP_SIZE; stride >>= 1) {
    if (threadIdx.x < stride) {
        sdata[threadIdx.x] += sdata[threadIdx.x + stride];
    }
    __syncthreads();
}

将数据加载到寄存器：将共享内存中的部分结果加载到线程的私有寄存器中。
```
if (threadIdx.x < WARP_SIZE) {
    int sum = sdata[threadIdx.x] + sdata[threadIdx.x + WARP_SIZE];
}
```
使用洗牌指令完成归约：在寄存器中，使用 __shfl_down_sync 指令继续归约过程。
```
for (int stride = WARP_SIZE / 2; stride > 0; stride >>= 1) {
    sum += __shfl_down_sync(0xffffffff, sum, stride);
}
```
这里掩码 0xffffffff 表示线程束中的所有线程都参与操作。虽然某些线程可能读取到无效数据，但我们最终只使用线程0的结果。
存储最终结果：由线程0将寄存器 sum 中的最终结果写入全局内存。

通过这种优化，我们避免了在线程束内阶段使用共享内存和昂贵的全局同步，从而显著提升了性能。

线程束投票指令 🗳️

线程束投票指令使线程能够协作判断哪些线程满足特定条件。这在需要基于线程局部条件做出集体决策时非常有用。

主要的投票指令包括：

__all_sync(mask, predicate): 检查所有参与线程的 predicate 是否均为非零（真）。
__any_sync(mask, predicate): 检查是否有任意参与线程的 predicate 为非零（真）。
__ballot_sync(mask, predicate): 返回一个32位掩码，其中每个位对应一个线程，如果该线程的 predicate 非零则置1。
__activemask(): 返回一个掩码，指示在当前执行点哪些线程是活跃的（即处于同一执行路径）。这在控制流分叉时尤其有用。

应用示例：优化队列插入

考虑一个内核，其中多个线程检查条件并将满足条件的值原子性地添加到全局队列。这会导致所有线程争用同一个全局原子计数器。

我们可以优化此过程，让同一线程束内的线程协作，仅由其中一个线程代表整个线程束执行一次原子操作。

以下是实现步骤：

确定活跃线程：使用 __activemask() 获取当前执行分支中活跃线程的掩码。
```
unsigned int active = __activemask();
```
选举领导线程：选择第一个活跃线程作为领导线程。使用 __ffs()（查找第一个置位位）内在函数。
```
unsigned int leader = __ffs(active) - 1; // 索引改为从0开始
```
计算需插入的元素数量：使用 __popc()（人口计数）内在函数计算活跃掩码中置位位的数量。
```
unsigned int num_active = __popc(active);
```
领导线程分配空间：仅由领导线程执行一次原子加操作，为整个线程束预留队列空间。
```
unsigned int j;
if ((threadIdx.x % WARP_SIZE) == leader) {
    j = atomicAdd(&q_size, num_active);
}
```
广播起始位置：领导线程将其获得的起始索引 j 广播给线程束内的所有活跃线程。
```
j = __shfl_sync(active, j, leader);
```

计算各线程的偏移量：每个活跃线程需要计算自己在所分配空间内的偏移量。这等于在该线程之前活跃的线程数量。

unsigned int thread_idx_in_warp = threadIdx.x % WARP_SIZE;
unsigned int prev_threads_mask = (1 << thread_idx_in_warp) - 1;
unsigned int prev_active_threads = active & prev_threads_mask;
unsigned int offset = __popc(prev_active_threads);

写入队列：每个线程使用计算出的位置 j + offset 将其值写入队列。
```
q[j + offset] = value;
```

通过这种优化，我们将每个线程束的多次原子操作减少为一次，并避免了线程间的争用。虽然在这个简化的示例中可能因开销而未显优势，但在更复杂的、原子操作开销更大的上下文中，这种协作方式能显著提升性能。

总结 📝

本节课我们一起学习了GPU编程中的线程束内同步技术。

我们首先回顾了线程束作为基本执行单元的特性，并引出了利用其SIMD执行模式进行高效同步的概念。接着，我们深入探讨了两类核心的线程束内协作机制：数据洗牌指令和投票指令。

通过将洗牌指令应用于归约模式，我们看到了如何避免在线程束内阶段使用共享内存和全局同步，从而提升性能。随后，我们利用投票指令及相关内在函数（__activemask, __ffs, __popc）优化了一个队列插入场景，展示了如何让线程束内线程协作以减少昂贵的原子操作。

这些高级特性为优化GPU程序提供了强大的工具，特别是在需要线程紧密协作或减少全局通信开销的场景中。理解并恰当应用这些技术，是成为高效GPU程序员的重要一步。

GPU计算：第21讲：固定内存与流

在本节课中，我们将要学习两个紧密相关的主题：固定内存和CUDA流。我们将探讨如何通过使用固定内存来优化主机与设备之间的数据传输，以及如何利用CUDA流来实现数据传输与内核执行的并行化，从而进一步提升GPU程序的整体性能。

上一讲我们介绍了线程束内的同步，学习了如何利用同一线程束内线程的特殊关系进行快速同步。本节中，我们来看看如何优化主机与设备间的通信。

什么是固定内存？

首先，我们来谈谈固定内存。什么是固定内存？

回想一下我们之前编写的向量加法代码。我们分配了主机内存，将源向量复制到设备，调用加法内核，然后将结果向量复制回主机。当我们分析这段代码的性能时，会发现GPU的总执行时间主要由数据拷贝时间主导，而内核计算时间其实很短。

那么，我们能否优化这个拷贝过程，减少每次在主机和设备之间复制数据所花费的时间呢？

为了理解如何优化拷贝，我们首先需要了解调用 cudaMemcpy 时发生了什么。这个过程使用了直接内存访问技术。

直接内存访问：DMA允许硬件单元（如DMA引擎）在不涉及CPU的情况下直接访问内存。当从主机内存复制数据到GPU内存时，CPU并不参与逐字节的读写，而是由DMA引擎负责完成传输，这样CPU就可以去处理其他任务。

然而，DMA引擎使用物理地址来访问内存，而CPU使用虚拟地址。操作系统可以动态地将虚拟页面换出到磁盘，或重新映射到不同的物理页面。如果DMA引擎正在读取的物理页面被操作系统换出或移动，DMA引擎无法察觉，从而导致它读取到错误的数据。

为了避免数据损坏，任何被DMA访问的内存页面都必须被“锁定”或“固定”在物理内存中，防止操作系统移动或换出它们。这就是固定内存。

默认情况下，使用 malloc 分配的内存不是固定的。因此，当调用 cudaMemcpy 时，CUDA运行时会执行以下步骤：

将数据从用户分配的缓冲区复制到一个内部的、固定的临时缓冲区。
然后，DMA引擎从这个固定缓冲区将数据复制到设备内存。

可以看到，这实际上进行了两次拷贝，带来了额外的开销。

如何使用固定内存进行优化？

为了优化拷贝，我们可以避免这额外的拷贝。方法就是直接将我们的主机数组分配在固定内存中。这样，DMA引擎可以直接从我们的缓冲区读取数据。

CUDA提供了专门的API来分配和释放固定内存：

使用 cudaMallocHost 分配固定内存。
使用 cudaFreeHost 释放固定内存。

其语法与 cudaMalloc 和 cudaFree 类似。

重要提示：不应将所有内存都分配为固定内存。固定内存会减少操作系统可用于交换的虚拟页面数量，从而可能降低系统整体性能。固定内存应仅用于需要频繁在主机和设备之间传输的数据。

通过将主机数组分配在固定内存中，我们观察到向量加法示例的GPU执行时间显著减少，因为拷贝时间大幅缩短。

什么是CUDA流？

接下来，我们讨论CUDA流。这涉及到系统架构层面的并行性，是一种任务并行或流水线并行。

在一个典型系统中，我们有CPU、主机内存、GPU和设备内存。系统实际上能够同时执行以下操作：

在GPU上执行内核计算。
将数据从主机内存复制到设备内存。
将数据从设备内存复制回主机内存。

这是三种不同的硬件资源。然而，我们之前的代码是顺序执行的：先拷贝所有数据到设备，然后执行内核，最后将所有数据拷贝回主机。在每个阶段，其他资源都处于空闲状态。

为了更充分地利用硬件，我们可以采用流水线技术。例如，在向量加法中，我们可以将输入数组分成多个段。然后：

将第一段数据（A1, B1）拷贝到设备。
当第一段数据在设备上就绪后，启动内核计算 C1 = A1 + B1。同时，开始将第二段数据（A2, B2）拷贝到设备。
当第一段计算完成，开始将结果C1拷贝回主机。同时，对已就绪的第二段数据（A2, B2）执行内核计算，并开始拷贝第三段数据（A3, B3）到设备。

如此循环，我们就能重叠不同段的拷贝和计算操作，实现并行。

如何使用流实现流水线？

为了实现这种流水线，我们需要：

异步内存拷贝：允许主机发起拷贝操作后立即继续执行，而不等待拷贝完成。
CUDA流：流是一系列按顺序执行的操作（如内存拷贝、内核启动）的队列。不同流中的操作可以并行执行，而同一流中的操作则按顺序执行。

默认情况下，所有操作都放入“默认流”（stream 0）中，因此它们是串行的。

以下是使用流和异步操作的关键API：

cudaStreamCreate：创建一个新的流。
cudaMemcpyAsync：执行异步内存拷贝，可以指定目标流。
内核启动：在内核配置中，第四个可选参数就是流标识符。

以下是实现流水线向量加法的核心步骤代码示例：

// 1. 创建多个流
int num_streams = 32;
cudaStream_t streams[num_streams];
for (int i = 0; i < num_streams; ++i) {
    cudaStreamCreate(&streams[i]);
}

// 2. 计算分段大小
size_t segment_size = (n + num_streams - 1) / num_streams; // 向上取整

// 3. 为每个分段提交异步操作到对应的流
for (int s = 0; s < num_streams; ++s) {
    size_t start = s * segment_size;
    size_t end = min(start + segment_size, n);
    size_t seg_len = end - start;

    // 异步拷贝主机到设备 (指定流)
    cudaMemcpyAsync(&d_x[start], &h_x[start], seg_len * sizeof(float), cudaMemcpyHostToDevice, streams[s]);
    cudaMemcpyAsync(&d_y[start], &h_y[start], seg_len * sizeof(float), cudaMemcpyHostToDevice, streams[s]);

    // 在内核启动中指定相同的流
    vector_add_kernel<<<grid_dim, block_dim, 0, streams[s]>>>(&d_x[start], &d_y[start], &d_z[start], seg_len);

    // 异步拷贝设备到主机 (指定流)
    cudaMemcpyAsync(&h_z[start], &d_z[start], seg_len * sizeof(float), cudaMemcpyDeviceToHost, streams[s]);
}

// 4. 等待所有流中的操作完成
cudaDeviceSynchronize();

// 5. 销毁流
for (int i = 0; i < num_streams; ++i) {
    cudaStreamDestroy(streams[i]);
}

通过性能分析工具（如 NVIDIA Nsight Systems）查看时间线，可以清晰地看到不同流中的内存拷贝（H2D, D2H）和内核执行（Compute）操作实现了重叠。最耗时的操作（通常是H2D拷贝）的“尾部”时间掩盖了其他操作的延迟，从而提升了整体吞吐量。

总结

本节课中我们一起学习了优化GPU程序通信性能的两个关键技术。

首先，我们了解了固定内存。通过使用 cudaMallocHost 直接将数据分配在主机端的固定内存中，可以避免CUDA运行时在默认拷贝路径中的额外缓冲拷贝，从而显著减少主机与设备间的数据传输时间。

其次，我们深入探讨了CUDA流。通过创建多个流，并将计算任务分段，利用 cudaMemcpyAsync 和指定流的内核启动，我们可以实现数据传输与内核执行的流水线并行。这允许设备在计算一个数据段的同时，接收下一个数据段并传回上一个结果段，充分重叠了不同硬件资源（PCIe带宽、GPU计算单元）的操作，隐藏了操作延迟，最终提升了程序的整体执行效率。

掌握固定内存和流的使用，对于构建高性能的、数据密集型的CUDA应用程序至关重要。

GPU计算｜CMPS 297S396AA：22：动态并行性

在本节课中，我们将要学习CUDA编程中的一个高级特性：动态并行性。它允许GPU上的线程在运行时启动新的内核网格，这对于处理具有嵌套并行性或工作负载大小未知的应用程序非常有用。

上一节我们介绍了固定内存和CUDA流，它们用于优化主机与设备间的数据传输。本节中我们来看看动态并行性，它允许GPU线程自身启动新的内核，从而在GPU上实现更灵活的并行计算模式。

什么是动态并行性？

动态并行性指的是在GPU上执行的线程能够启动新的网格在GPU上运行的能力。如果一个网格中的线程正在GPU上运行，那么其中任何一个线程都可以去启动另一个线程网格。

为什么需要这样做？一个主要原因是，当我们的计算需要启动多个内核时，可以直接从GPU启动，而无需再与CPU同步。更重要的是，它适用于具有嵌套并行性的应用程序。

嵌套并行性

嵌套并行性是指应用程序中存在多个层次的并行性。例如，我们有一批可以并行执行的工作单元，而每个工作单元内部又包含更多可以并行执行的工作。当内层并行工作的数量在编译时未知时，动态并行性就特别有用。如果内层工作量未知，我们无法预先启动固定数量的线程。此时，让每个线程在运行时发现自己内部有多少并行性，然后自己启动更多线程，就变得非常有用。

以下是具有嵌套并行性模式的两种主要应用类型：

不规则嵌套工作：不同线程内部具有不同数量的嵌套并行工作。例如，在图算法中，每个顶点（对应一个线程）可能有不同数量的邻居需要并行访问。另一个例子是贝塞尔曲线绘制，每条曲线（对应一个线程）根据其曲率可能需要不同数量的点来绘制。
递归深度未知：所有线程可能具有相同数量的嵌套工作，但线程是否启动更多工作取决于运行时信息。例如，在四叉树/八叉树空间划分中，一个线程负责处理一个区域框，根据框内点的活动情况决定是否进一步细分（递归）。另一个例子是分治算法，如快速排序，根据分区属性决定是否继续划分。

CUDA SDK中包含了贝塞尔曲线、四叉树和快速排序的动态并行性示例。本节课我们将以广度优先搜索图算法为例进行讲解。

在BFS中应用动态并行性

作为回顾，我们之前实现的BFS使用“边界”法。在每次迭代中，前一边界中的每个顶点（对应一个线程）会顺序访问其所有邻居，并决定是否将邻居加入当前边界。

关键点在于这个顺序访问邻居的循环。我们可以使用动态并行性来并行化这个循环，即为每个邻居分配一个线程来处理。

以下是之前顺序处理邻居的核心循环代码：

for (unsigned int i = 0; i < num_neighbors; ++i) {
    // 加载邻居顶点索引
    // 尝试原子性地访问该邻居
    // 如果访问成功，将其加入当前边界
}

使用动态并行性后，父线程将启动一个子内核网格，网格中的每个线程负责处理一个邻居。子内核的代码大致如下：

__global__ void bfs_child_kernel(...) {
    unsigned int i = threadIdx.x + blockIdx.x * blockDim.x;
    if (i < num_neighbors) {
        // 加载第i个邻居的顶点索引
        // 尝试原子性地访问该邻居
        // 如果访问成功，将其加入当前边界
    }
}

父线程中启动子内核的代码：

// 计算网格配置
unsigned int threads_per_block = 1024;
unsigned int num_blocks = (num_neighbors + threads_per_block - 1) / threads_per_block;
// 启动子内核
bfs_child_kernel<<<num_blocks, threads_per_block>>>(...);

性能挑战与优化

直接为每个边界顶点线程启动一个子网格会导致性能严重下降（可能慢10倍以上）。主要原因包括：

待处理启动数限制：GPU运行时默认最多缓冲2048个动态启动的网格。超过此限制可能导致错误或挂起。可以使用 cudaDeviceSetLimit(cudaLimitDevRuntimePendingLaunchCount, ...) 提高限制，但这并非根本解决方案。
流与序列化：默认情况下，同一线程块内的所有线程共享同一个默认流。它们发起的动态内核会在这个流中序列化执行，即使GPU有空闲资源。可以通过编译器标志 --default-stream per-thread 改为每个线程拥有独立的默认流，从而允许更多并发。
大量小网格开销：许多顶点只有少量邻居，为它们启动网格（可能只有几个线程）的开销远大于并行收益，并且会因网格太小而无法充分利用SM。

针对以上问题，我们可以进行以下优化：

优化1：设置启动阈值：仅当邻居数量超过某个阈值（例如1024）时，才使用动态并行性启动子网格。否则，仍在父线程中顺序处理邻居。这显著减少了动态启动的数量和开销。
优化2：聚合启动（概念性）：让一个线程（例如一个Warp或一个块中的线程）收集多个线程的工作，然后代表它们发起一个更大的联合网格。这更复杂，通常需要线程间通信和协作。

应用启动阈值优化后，BFS性能从远差于静态版本转变为优于静态版本，展示了动态并行性在特定场景下的威力。

其他主题

卸载驱动代码：可以将主机端驱动计算循环（例如BFS中遍历各层的主循环）也放到GPU上一个线程中执行。这样主机只需启动一次内核即可完成整个计算，从而被释放出来执行其他任务。但这通常不是为了提升性能，而是为了释放CPU。
内存可见性：父线程在启动子网格前对全局内存的写入，对子网格中的线程是可见的。子网格对全局内存的写入，只有在父线程调用 cudaDeviceSynchronize() 同步后，才对父线程可见。局部内存和共享内存的指针不应传递给子网格，因为子网格可能在不同的SM上执行。
嵌套深度：硬件对动态启动的嵌套深度有限制，当前硬件通常为24层。超过此深度将无法启动新的网格。

本节课中我们一起学习了CUDA动态并行性的概念、应用场景（特别是嵌套并行性）、以及如何在广度优先搜索算法中实现并优化它。我们还了解了待处理启动限制、流的影响、以及通过设置阈值来避免过多小网格开销的优化方法。最后，我们简要探讨了卸载驱动代码、内存可见性和嵌套深度等其他相关主题。动态并行性是一个强大的工具，但需要谨慎使用并配合适当的优化策略。

GPU计算：23：课程总结与高级主题概览 🎓

在本节课中，我们将学习GPU计算课程的最后一个章节。由于课程时间有限，我们无法深入探讨所有主题。因此，本节课将作为一个“大杂烩”，简要介绍一些我们尚未覆盖但值得了解的高级功能和概念，以便你在未来的项目或阅读代码时能够识别它们。我们将快速浏览这些主题，并以一个总结结束本课程。

以下是本节课将要简要介绍的主题列表。

多GPU编程 🖥️🖥️

到目前为止，我们处理的模型是一个CPU连接一个GPU。但实际上，一个CPU可以连接多个GPU。在这种情况下，CPU需要协调各个GPU的工作。

在CUDA中，你可以使用 cudaGetDeviceCount 这样的API调用来获取节点上的GPU数量。然后，你可以使用 cudaSetDevice 来选择其中一个GPU。例如，如果你想在GPU 0上启动内核，可以设置设备为0，然后启动内核。

通常，在使用多个GPU和单个CPU时，我们也会使用多个CPU线程。每个CPU线程负责操作一个GPU，例如调用内核或执行数据拷贝。对于多线程，我们通常使用像OpenMP这样的库。

此外，我们还可以通过CPU在不同GPU之间进行同步，或者稍后我们会看到，也可以使用GPU之间的共享内存。

现在，这是一种同步方式。多GPU编程的一种形式是多个GPU连接到一个CPU上。

但我们也可以在网络中的多个不同节点上拥有多个GPU。这意味着我们有多个通过网络连接的CPU，每个CPU又可以连接多个GPU。在这种情况下，我需要协调的不仅仅是连接到一个CPU上的多个GPU，还包括网络中不同CPU上的多个GPU。

在这种情况下，我们通常使用像MPI（消息传递接口）这样的库。MPI允许我们编程通过网络连接的多个CPU。在MPI中，我们有“秩”的概念。一个MPI秩类似于CPU上的不同线程，但一个秩会在不同的CPU上运行。我们也可以在同一CPU上有多个秩，每个秩对应一个在该CPU上运行的进程，并且每个秩可以与一个GPU交互。

当我们进行网络上的多GPU编程时，一个重要的考虑因素是通信开销。当我们在一个CPU上编程多个GPU时，如果我想在它们之间共享数据，可以通过CPU内存进行，或者根据系统配置直接跨GPU进行。但通常，我可以使用CPU将数据从一个GPU复制到另一个GPU，这相对较快。

然而，如果我想跨网络共享数据，例如从一个网络节点上的GPU到另一个网络节点上的GPU，我不能直接复制。我需要将数据从GPU内存复制到CPU内存，然后通过网络传输到另一个CPU内存，再复制到另一个GPU。这个过程非常耗时。因此，在多GPU并行化时，一个重要的考虑是如何将网络通信与GPU上的计算重叠起来。

一个例子是，如果我们进行某种模板计算，并将计算空间划分到多个GPU上，每个GPU处理一部分空间。那么，我们需要在不同GPU之间共享“光环”元素。通常的做法是，先计算边界元素，然后在计算内部元素的同时，通过网络进行边界元素的通信。这就是如何将GPU之间的通信与GPU上内核执行重叠的一个例子。

以上简要介绍了如何编程多个GPU以及在此过程中面临的一些挑战。

互连技术 🔗

另一个我们未深入探讨的主题是互连技术。当一个CPU连接到多个GPU时，它们是如何连接的呢？通常通过一种互连技术连接，常见的是PCIe。PCIe是一种在许多系统中用于连接CPU和GPU的通用互连技术。

PCIe是一种非常通用的互连技术，我们用它来连接CPU和GPU，也可以用它连接CPU和其他设备，例如FPGA。由于PCIe通用性强，它在许多系统中被广泛使用。但也有其他类型的互连技术，例如NVLink。NVLink实际上提供了更快的互连速度，我们可以用它来连接多个GPU，甚至可以用来连接GPU和支持NVLink互连的某些CPU。PCIe得到广泛支持，大多数CPU和GPU都支持PCIe。然而，像NVLink这样的技术，只有部分CPU和特定代次的NVIDIA GPU（通常是较新的代次）支持。

因此，PCIe是通用的选择，而NVLink更专用于NVIDIA GPU及其较新代次，由部分CPU支持，能提供更高的带宽。但NVLink不能用于所有场景。如果我需要一个适用于任何CPU和GPU的通用方案，那么PCIe通常是人们使用的通用互连技术。

统一虚拟寻址 🗺️

下一个主题是统一虚拟寻址。到目前为止，我们一直假设CPU可以访问其主内存，GPU可以访问其自己的全局内存，并且我们将它们视为各自独立的地址空间。CPU主内存是一个独立的地址空间，GPU的全局内存是另一个独立的地址空间。

但实际上，后来的GPU系统或NVIDIA GPU系统引入了统一虚拟地址空间的概念。我们仍然有两个独立的物理内存：CPU的主内存和GPU的全局内存（或设备内存），它们有各自独立的物理地址空间。但我们有一个统一的虚拟地址空间，主机和设备使用相同的虚拟地址空间。

这意味着，指向CPU主内存支持地址范围内的内存地址，与GPU支持的地址范围是互斥的。同一个内存地址不能同时在两个地方都有效；它要么在CPU上有效，要么在GPU上有效。

当然，即使在没有统一虚拟内存的情况下，也不建议在CPU和GPU上使用相同的内存地址。但有了统一虚拟地址空间，地址范围是互斥的，这就不再是问题了。

那么，在CPU和GPU之间拥有统一的虚拟地址空间（即指针地址范围互斥）有什么优势呢？优势在于，当我知道CPU使用与GPU分离的虚拟地址范围时，我就可以从指针值本身判断出该指针是指向CPU内存还是GPU内存。数据的位置可以从指针值得知。

这样做的优势是，当我们执行 cudaMemcpy 时，我们不再需要指定是从主机到设备还是从设备到主机的拷贝。在统一虚拟地址空间下，我们可以简单地使用 cudaMemcpyDefault，运行时系统可以从指针值判断出该指针是指向设备还是主机。

这就是在CPU和GPU之间拥有统一虚拟地址空间（CPU主内存使用虚拟地址空间中的某个地址范围，GPU设备内存使用虚拟地址空间中另一个独立的地址范围）的优势。

零拷贝内存 📥

统一虚拟寻址的另一个优势，我稍后会展示。但让我先介绍另一个主题：零拷贝内存。

什么是零拷贝内存？零拷贝内存使得在GPU上运行的线程能够直接访问主机内存。

这意味着什么？在此之前，我们一直说GPU上的线程不能访问CPU主内存中的数据；我必须先将数据从CPU复制到GPU。但实际上，我们可以使用一种叫做零拷贝内存的技术。

“零拷贝”指的是我们不需要进行任何显式的 cudaMemcpy 操作。相反，我们可以让在GPU上运行的线程直接访问CPU内存中的数据。当我们这样做时，当GPU上运行的线程尝试访问指向CPU数据的某个内存地址时，运行时系统会按需将该数据从CPU内存复制到GPU内存，并提供给该线程。

因此，我们不需要显式地进行 cudaMemcpy 操作，我们可以直接分配数据（当然有特殊的方法），然后让GPU线程在访问时，运行时系统按需从CPU复制数据到GPU。

这里的一个重要点是，如果GPU线程要访问CPU内存中的数据，并且这些数据将按需复制，那么这意味着我们将使用DMA来复制数据。那么，这告诉我，我想要从GPU线程访问的、位于CPU内存中的数据必须位于哪里？必须是“固定”内存中。没错，为了使零拷贝内存工作，我希望GPU线程能够访问的、位于CPU内存中的数据应该分配在固定内存中。

我们使用 cudaHostAlloc 这个API调用来实现。cudaHostAlloc 命令在分配内存时，会在CPU上给我们一个固定的内存范围。这个内存不仅是固定的，而且如果我们尝试从GPU访问这些数据，运行时系统会自动为我们从CPU复制数据到GPU。

在没有统一虚拟内存的情况下，我们需要做的是：我们想要在CPU内存中访问的数据，也需要映射到GPU的地址空间。因此，我们可以使用一个叫做 cudaHostGetDevicePointer 的命令。它接受一个通过 cudaHostAlloc 分配的、指向CPU上固定内存的指针，并获取一个GPU上可以使用的对应指针，以便GPU访问这些数据。

然而，如果我们的系统支持统一虚拟寻址，那么上述步骤就是不必要的。因为如果我有统一虚拟寻址，就意味着CPU和GPU上的指针是互斥的。我从 cudaHostAlloc 得到的指针，可以直接交给GPU。当GPU解引用这个指针时，GPU可以从指针值知道这个内存实际上在CPU上，而不是在GPU上，因此它知道需要去CPU按需获取数据。但是，如果我不支持统一虚拟寻址，那么GPU就无法从指针值判断这是指向CPU内存还是GPU内存。这就是为什么我们需要使用像 cudaHostGetDevicePointer 这样的函数来将指向CPU内存的指针转换为GPU可用的指针，以便GPU知道需要去CPU获取该内存。

事件 ⏱️

下一个主题是事件。到目前为止，在我们编写的所有代码中，为了收集计时信息，我们一直在每次拷贝和内核调用后进行同步。我这样做只是为了保持代码简单，因为我不想从一开始就引入事件，否则我就必须从一开始就引入流，而我不想让你们从一开始就担心流的问题。

但在实践中，这种方法会干扰执行并导致阻塞。我真的不想每次进行拷贝或内核调用时都进行同步，特别是如果我使用多个流，并且这些内核本应彼此并行执行，或者拷贝本应彼此并行执行。

我们可以在不进行类似 cudaDeviceSynchronize 同步的情况下收集计时信息，方法是使用事件。事件可以通过简单地添加到流中来收集时间戳，而无需同步，以便我们在该流中的某个操作开始或完成时收集时间戳。

我们可以这样使用事件：例如，我们可以写 cudaEventCreate 来创建一个事件；我们可以写 cudaEventRecord 来记录一个事件，并提供我们想要记录到的流。如果我们将一个事件添加到一个流中，事件将捕获它在流中轮到的时间。然后，我们可以做类似 cudaEventSynchronize 的操作，等待某个事件发生。

如果我在一个流中放了多个操作，然后放了一个事件，接着又在事件后面放了多个操作。但后来我只想等待事件完成，而不关心事件后面的操作是否完成，我可以使用 cudaEventSynchronize，它等待流中特定事件完成，而无需等待该事件之后的操作完成。这比 cudaDeviceSynchronize 更好，因为后者基本上等待所有已提交的操作完成。

我们还可以获取事件之间的经过时间，使用 cudaEventElapsedTime 来获取两个事件之间的时间。例如，如果我想知道执行某个 cudaMemcpy 花了多长时间，但我不想在那次拷贝之后进行同步（它可能是一个异步拷贝），我只想在不进行前后同步的情况下计时，我可以做的是：在 cudaMemcpy 之前记录一个事件，在之后记录一个事件，然后获取这两个事件之间的经过时间。这告诉我那次拷贝花了多长时间，而无需我进行同步并干扰我试图实现的异步行为。

当然，我可以用 cudaEventDestroy 来销毁之前创建的事件。

这些信息中的很多也可以通过性能分析器收集。性能分析器也会告诉我们内核运行了多长时间，或者 cudaMemcpy 运行了多长时间。但如果出于某种原因，我也想在代码本身中收集这些信息（也许我想用它来做某种决策），我可以使用事件。

张量核心 ⚡

下一个主题是张量核心。张量核心随着Volta V100架构被引入，它们本质上是SM内部可编程的矩阵乘加单元。每个SM都有张量核心，可以用来在一条指令中完成一次小的矩阵乘加操作。

张量核心之所以重要，尤其是在深度学习兴起时，是因为很多深度学习工作负载主要就是做矩阵乘法。我们已经看到GPU在执行矩阵乘法方面有多出色，但实际上，通过使用张量核心，我们可以做得更好。张量核心是专用硬件，可以在一条指令中完成小的矩阵乘加操作，而不需要执行多条加载、加法和乘法指令来实现该操作。

因此，当我们使用矩阵乘法时，张量核心能给我们带来很大的加速，但我们必须使用特定的API来显式地使用它们。在Volta V100架构中，我们有640个张量核心，因为我们有80个SM，实际上每个SM有8个张量核心，每个张量核心能够执行一个4x4的矩阵乘法。它可以执行 D = A * B + C 操作，其中D、A、B、C都是4x4矩阵。

库 📚

下一个我们没机会讨论的主题是库。CUDA包含许多库，或者与CUDA一起流行的库，它们提供了常见的并行原语。我们在课程中实际使用的许多原语都在库中提供了。

例如，有一个叫做Thrust的库，它包含归约、扫描、筛选和排序等原语。我们在课程中手动实现了许多这些操作，这是一个重要的练习，因为我们想教你们这些不同原语的工作原理、它们在什么情况下性能好、什么情况下不好，以及优化它们背后的原理。但在实践中，如果你正在编写一个应用程序，你可能不会手动自己实现归约、扫描或排序操作作为应用程序的一部分；你可以直接调用像Thrust或CUB（Thrust使用的另一个库）这样的库。它们会为你执行归约、扫描、排序或筛选操作，你不需要自己手动实现。很多时候，我们在课程中实现的一些原语，你不需要担心自己实现；在一个成熟的应用程序中，你只需调用这些库来为你完成。但课程的目的不是教你们库，而是教你们关于如何推理并行化代码的概念，这就是为什么我们自己去实现它们。

另一个是cuBLAS，用于基本的密集线性代数，如向量加法、矩阵向量乘法、矩阵乘法。我们在课堂上实现了矩阵乘法，目的是教你们在并行化这类操作时面临的一些决策和挑战。但在实践中，如果你的应用程序中使用矩阵乘法或矩阵向量乘法，你可能不会从头实现；你可能会使用像cuBLAS这样的库，它已经有很多针对这些基本密集线性代数操作的高度优化实现。

还有NVBLAS，这是一个多GPU库，可以在多个GPU上执行这些线性代数操作，它构建在cuBLAS之上。

另一个是cuSPARSE，用于稀疏线性代数。在课程中，我们研究了如何执行SpMV，并研究了不同的稀疏矩阵存储格式。实际上，我们有一个库可以为我们执行各种稀疏线性代数例程。在实践中，如果你的应用程序需要执行稀疏矩阵操作，你可以直接使用这些内置函数来创建和操作这些矩阵，而不是从头开始自己实现。但同样，我们教你们这些是为了让你们了解在优化这类工作负载时需要经历的过程，这些原理也适用于你可能想要并行化的其他类型的工作负载。

另一个是cuSOLVER，用于分解和求解例程，如果你在进行矩阵分解。cuSOLVER实际上构建在cuBLAS和cuSPARSE之上。

还有cuDNN，在深度神经网络中很流行，它实现了许多对深度神经网络有用的例程。事实上，你们中可能有些人熟悉的许多框架（如果你之前从事过深度学习），比如Caffe、MXNet、TensorFlow、PyTorch，这些框架很多都给你提供了在GPU上运行的选项。这些框架能够获取你的模型并在GPU上运行，其底层原理就是调用cuDNN。因此，许多流行的、使用GPU来加速你的神经网络的深度学习框架，它们只是间接地使用cuDNN来加速你的代码。cuDNN本身是一个非常流行的库，被这些框架间接使用，而cuDNN本身也使用像cuBLAS这样的库，因为DNN本身会执行各种密集线性代数和其他操作。

另一个是NVGRAPH，用于图处理，它提供了图的不同表示和可以在图上执行的操作。cuFFT用于快速傅里叶变换，如果你对涉及FFT的任何操作感兴趣。还有NPP，用于各种信号、图像和视频处理。

所有这些库都非常有用且被广泛使用。每当你编写GPU应用程序时，我们在本课程中涉及的许多并行模式实际上都在这些库中实现了。很多时候，你不需要从头实现这些并行模式；你只需使用这些库来为你执行这些操作。构建一个应用程序通常涉及组合这些不同库中的许多函数。但同样，我们经历所有这些并行模式的原因是为了让你们了解在优化GPU代码时面临的挑战以及可以应用的优化原则，你们学到的这些原理也适用于如果你想实现一些不一定内置在任何这些库中的自定义应用程序。

其他编程接口 🔌

下一个主题是关于其他编程接口。在课程中，我们一直专注于CUDA。然而，CUDA并不是编程GPU的唯一接口。还有其他编程接口。

其中之一是OpenCL。OpenCL类似于CUDA，但它是开源的，并且在不同非NVIDIA GPU之间更具可移植性。还有其他硬件供应商，除了NVIDIA GPU还有其他GPU（我将在下一张幻灯片讨论）。但如果你对编写能够同时在NVIDIA GPU和非NVIDIA GPU上运行的代码感兴趣，OpenCL提供了一个更可移植的编程接口，而CUDA则非常特定于NVIDIA GPU。实际上，OpenCL代码不仅可移植到不同的GPU，还可移植到不同的设备。实际上，有些编译器可以将OpenCL编译到CPU，也有高级综合工具可以将OpenCL编译到FPGA。因此，如果你只对在NVIDIA GPU上运行感兴趣，也许CUDA比OpenCL更简单一些。但如果你希望你的代码能够跨不同种类的GPU移植，并且也希望移植到GPU和其他硬件（如FPGA和CPU），那么OpenCL是一个流行的编程接口。

另一个接口是OpenACC。OpenACC是一个基于指令的GPU编程接口。这里“基于指令”指的是类似于OpenMP的东西，你只需获取带有许多循环的常规顺序代码，你可以以某种方式注释这些循环，类似于使用OpenMP进行并行化注释的方式。通过简单地使用OpenACC以正确的方式注释你的循环，OpenACC编译器会自动获取你的顺序代码并将其并行化以在GPU上运行（当然，假设你注释的循环是可并行化的等等）。基于指令的编程在并行编程中非常流行，它在CPU上很流行，人们一直使用OpenMP来在CPU上并行化工作，而不是使用像Pthreads这样的低级线程库。类似地，如果你不想担心分配内存、调用内核、将工作划分到线程块等所有事情，你可以使用OpenACC来注释顺序代码中的循环，以便自动将该代码移植到GPU上运行。

另一个库是C++ AMP。C++ AMP是一个库，你可以用它来在C++内部实现GPU程序。在你的C++代码中，你可以使用Lambda函数和某些特殊语法等，以便直接在C++中实现在GPU上运行的代码，而无需所有这些CUDA扩展。

其他GPU硬件 🧩

下一个主题是其他GPU硬件。在课程中，我们一直专注于CUDA，也一直专注于独立的NVIDIA GPU。“NVIDIA GPU”指的是由NVIDIA制造的GPU，NVIDIA是一家流行的GPU制造公司。“独立”指的是GPU与CPU不在同一芯片上集成；你有一个与CPU分离的独立芯片的GPU。

但实际上，除了NVIDIA GPU还有其他GPU硬件。AMD制造GPU，你可能熟悉Radeon GPU，这些是由AMD制造的，它们在游戏和高性能计算中也很流行。NVIDIA GPU在游戏和高性能计算中很流行，AMD GPU在游戏和高性能计算中也很流行。

Arm设计GPU，你可能听说过Mali GPU。Arm在移动设备中非常流行，因此移动设备中使用的一些GPU遵循Arm架构。

此外，还有Intel GPU。Intel实际上也制造自己的GPU，尽管Intel以CPU闻名。很多这些GPU实际上集成在你的机器中。因此，如果你有一个Intel处理器，你的系统中可能也有一个Intel GPU，用于控制屏幕上的部分图形。

除了其他GPU硬件，正如我所说，我们一直专注于独立GPU，但实际上也有集成GPU，你可以将CPU和GPU放在同一芯片上。在这种情况下，CPU和GPU将共享相同的物理内存，这样做的好处是我们没有任何通过PCIe的数据传输；实际上，CPU和GPU都可以访问相同的物理内存，因为它们在同一芯片上。

关于CPU比较的说明 ⚖️

最后我想提一下的是关于本课程中所有与CPU比较的免责声明。到目前为止，在课程中，我们一直在将GPU与CPU进行比较，但我们使用的是单线程、非向量化的CPU代码。这样做夸大了我们看到的GPU加速效果。为了本课程的目的，也许这没关系。但一般来说，这种做法是不被鼓励的。例如，如果你要发表一篇论文，报告CPU和GPU之间的加速比，如果你想在GPU和CPU之间进行公平的比较，更公平的做法是将我们的GPU实现与使用向量化CPU代码的并行CPU实现进行比较。因为就像GPU可以并行运行许多线程块或线程束一样，CPU也可以并行运行许多线程，而且对于许多CPU，每个CPU核心都有一个向量单元，我们可以在那里向量化CPU代码，而不仅仅是在核心上运行单个线程。

我们在课程中没有这样做，原因是本课程的目标不是专注于CPU并行化，而是专注于GPU并行化。所以我不想让你们同时为CPU上的并行计算而烦恼。这可能是你们在本课程之前或之后在另一门课程中会涉及的内容。然而，我想说明这一点，因为我希望你们对我们一直做的CPU比较持保留态度。因为事实上，就像我们一直在GPU上进行并行化和向量化一样，我们也可以在CPU上做同样的事情。因此，更公平的CPU比较应该是与使用向量化的并行CPU代码进行比较。

总结 📝

本节课中，我们一起学习了GPU计算课程的最后一个章节。我们快速浏览了多个高级主题，包括多GPU编程、互连技术、统一虚拟寻址、零拷贝内存、事件、张量核心、常用CUDA库、其他GPU编程接口以及其他GPU硬件。我们还讨论了在进行CPU与GPU性能比较时应注意的公平性问题。

虽然我们未能深入探讨每个主题的细节，但希望本次概览能为你提供一个知识框架，让你知道在未来的学习或项目中，当遇到这些概念时可以去哪里寻找更多信息。GPU计算是一个广阔且快速发展的领域，本课程为你打下了坚实的基础。感谢你的参与，祝你未来在并行计算的道路上一切顺利！

GPU计算：02：数据并行编程

概述

在本节课中，我们将要学习数据并行编程的基本概念，并了解如何在GPU上实现它。我们将以向量加法为例，介绍CUDA编程模型的基础知识，包括内存管理、内核启动和线程组织。

回顾：GPU计算简介

上一节我们简要介绍了GPU计算的历史背景和基本原理。我们提到，尽管摩尔定律预测晶体管数量会持续增长，但自2005年起，处理器频率的提升遇到了瓶颈，这主要是由于功耗墙的限制。因此，为了获得更高的性能，行业转向了增加处理器核心数量的方向。

CPU和GPU代表了两种不同的设计哲学：

面向延迟的设计：旨在最小化完成单个任务所需的时间。CPU是这种设计的典型代表，它拥有少量但功能强大的算术逻辑单元、大容量缓存和复杂的控制逻辑。
面向吞吐量的设计：旨在最大化单位时间内完成的任务数量，即使这可能以增加单个任务的完成时间为代价。GPU是这种设计的典型代表，它拥有大量小型ALU、较小的缓存和简化的控制逻辑，并通过大量线程来隐藏延迟。

我们还介绍了CUDA平台的诞生，它使得通用GPU编程成为可能，并极大地推动了GPU在高性能计算、科学计算和机器学习等领域的应用。

数据并行编程与CUDA模型

本节中，我们来看看数据并行编程，并将在CUDA编程模型的背景下进行介绍。CUDA是当前最流行的GPU编程接口之一。

并行性主要有两种类型：

任务并行：指对相同或不同数据执行的不同操作可以并行进行。例如，文本编辑器在您编辑文本的同时，后台运行拼写检查。
数据并行：指对大量不同数据执行相同的操作可以并行进行。例如，计算屏幕上每个像素的颜色值。

应用程序中的任务并行性通常有限，因为不同任务的数量不会太多。而数据并行性则潜力巨大，因为只需对更大的数据集运行相同的程序即可释放更多并行性。这使得数据并行性非常适合在拥有大量并行计算单元的GPU上运行。

示例：向量加法

为了展示数据并行性及其在GPU上的实现，我们从一个简单的计算开始：向量加法。这可以看作是数据并行编程的“Hello World”。

向量加法即两个输入向量 X 和 Y 相加，得到输出向量 Z。其数学表示为：
Z[i] = X[i] + Y[i]，对于所有 i。

在CPU上顺序执行的代码如下：

for (int i = 0; i < n; ++i) {
    Z[i] = X[i] + Y[i];
}

CPU-GPU系统架构与内存管理

在讨论GPU并行实现之前，我们需要了解包含GPU的典型系统架构。

在CUDA术语中：

主机：指CPU及其主机内存。
设备：指GPU及其设备内存（也称为全局内存）。

CPU和GPU拥有独立的内存空间，通常不能直接访问对方的内存（除非使用高级特性如统一虚拟内存）。因此，要将计算任务卸载到GPU，需要遵循一个典型的操作序列：

以下是执行GPU计算的标准步骤：

在GPU上分配内存。
将数据从主机内存复制到设备内存。
在GPU上执行计算。
将结果从设备内存复制回主机内存。
释放GPU内存。

分配与释放GPU内存

使用 cudaMalloc 函数在GPU上分配内存。它返回一个指向设备内存的指针。

cudaError_t cudaMalloc(void** devPtr, size_t size);

对应的释放函数是 cudaFree。

cudaError_t cudaFree(void* devPtr);

主机与设备间数据传输

使用 cudaMemcpy 函数在主机和设备间复制数据。

cudaError_t cudaMemcpy(void* dst, const void* src, size_t count, cudaMemcpyKind kind);

其中 kind 指定复制方向，常用的是 cudaMemcpyHostToDevice 和 cudaMemcpyDeviceToHost。

在GPU上启动并行计算：内核与线程

现在数据已在GPU上，我们需要在GPU上并行执行向量加法。思路是：为向量中的每个元素分配一个GPU线程，每个线程负责计算一个输出元素。

在CUDA中，我们通过启动一个内核来在GPU上创建大量线程。这些线程被组织成一个两层结构：

网格：所有线程的集合。
线程块：网格中的线程被分组为多个线程块。同一个块内的线程可以相互协作，而不同块间的线程则相对独立。

启动内核时，需要指定网格大小（块的数量）和块大小（每个块中的线程数）。内核调用使用特殊的尖括号语法 <<< >>>。

// 假设每个块有512个线程
int threadsPerBlock = 512;
// 计算需要的块数（向上取整，以确保有足够线程覆盖所有元素）
int numBlocks = (n + threadsPerBlock - 1) / threadsPerBlock;
// 启动内核
vectorAddKernel<<<numBlocks, threadsPerBlock>>>(x_d, y_d, z_d, n);

编写GPU内核函数

内核函数是每个GPU线程将要执行的代码。它用 __global__ 关键字声明，表明这是一个由CPU调用、在GPU上执行的函数。

在内核函数内部，每个线程需要确定自己负责处理哪个数据元素。CUDA提供了内置变量来帮助线程识别自己：

blockIdx.x：线程所在块的索引。
threadIdx.x：线程在其块内的索引。
blockDim.x：每个块中的线程数量。

通过组合这些变量，线程可以计算出自己的全局索引：
int i = blockIdx.x * blockDim.x + threadIdx.x;

然后，线程使用这个全局索引 i 来访问数组中对应的元素：

__global__ void vectorAddKernel(float* x, float* y, float* z, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) { // 边界检查，防止访问越界
        z[i] = x[i] + y[i];
    }
}

注意 if (i < n) 这个边界条件。由于我们通过向上取整来计算块的数量，启动的线程总数可能略多于向量元素个数 n。这个检查确保多余的线程不会执行无效的内存访问。

编译与执行

CUDA程序使用 nvcc 编译器进行编译。nvcc 会自动将代码分离为主机代码（CPU）和设备代码（GPU），并分别进行编译和链接。

nvcc -o vector_add vector_add.cu

函数执行空间说明符

CUDA使用关键字来指定函数的执行位置：

__host__：默认，函数在CPU上执行。
__global__：内核函数，由CPU调用，在GPU上执行。
__device__：函数在GPU上执行，只能由其他 __device__ 或 __global__ 函数调用。
__host__ __device__：函数同时为CPU和GPU编译，可用于代码复用。

异步执行与同步

默认情况下，内核调用是异步的：CPU在启动内核后立即继续执行后续代码，而不会等待GPU内核完成。这允许CPU和GPU并行工作。

如果需要等待内核完成（例如，为了准确计时或确保数据可用），可以调用 cudaDeviceSynchronize() 函数。

错误检查

大多数CUDA API函数返回 cudaError_t 类型的错误码。良好的编程习惯是检查这些错误。

cudaError_t err = cudaMalloc(...);
if (err != cudaSuccess) {
    // 处理错误，例如打印错误信息
    printf("CUDA error: %s\n", cudaGetErrorString(err));
}

性能考量

在最初的完整示例中（包含内存分配、复制、计算、复制回和释放），GPU的加速比可能看起来并不惊人。这是因为数据迁移（在PCIe总线上复制）的开销很大。

然而，在典型的GPU加速应用中，数据通常只在计算开始前复制到GPU，在计算结束后复制回CPU。中间会执行大量长时间运行的内核。如果只计算内核执行时间，GPU相对于CPU的顺序代码通常会显示出巨大的性能提升。

总结

本节课中我们一起学习了数据并行编程的基础知识。我们以向量加法为例，完整走了一遍CUDA编程流程：

理解了CPU（主机）和GPU（设备）的分离式内存架构。
学会了使用 cudaMalloc、cudaMemcpy 和 cudaFree 管理设备内存。
掌握了如何通过 <<< >>> 语法配置并启动一个内核，在GPU上创建大量线程。
学会了编写 __global__ 内核函数，并使用 blockIdx、threadIdx 和 blockDim 等内置变量为每个线程分配唯一的工作。
了解了编译、异步执行、同步和错误检查等基本概念。

你现在已经掌握了编写一个简单CUDA程序的核心技能。本课程后续将深入探讨各种并行模式、优化策略以及GPU的高级特性，以帮助你充分发挥GPU硬件的强大计算能力。

GPU计算：03：多维网格与数据

在本节课中，我们将学习如何使用CUDA创建多维网格来处理多维数据。我们将通过三个具体的例子来演示：将RGB彩色图像转换为灰度图像、图像模糊处理以及矩阵乘法。这些例子将帮助我们理解如何有效地利用GPU的并行计算能力处理二维和三维数据。

回顾：数据并行编程

上一节课我们介绍了数据并行编程，并使用CUDA编程模型进行了演示。我们以向量加法为例，说明了任务并行与数据并行的区别。

任务并行：并行执行不同类型的任务，这些任务可能操作相同或不同的数据。通常能释放的并行度有限。
数据并行：对数据的不同部分执行相同的操作。随着数据量增大，可获得的并行潜力更高，无需为不同任务编写新代码。

在向量加法的例子中，我们学习了如何将计算卸载到GPU上。这涉及以下步骤：

在GPU上分配内存。
将数据从CPU复制到GPU。
在GPU上启动内核执行计算。
将结果从GPU复制回CPU。
释放GPU内存。

我们使用了cudaMalloc、cudaMemcpy和cudaFree等CUDA接口。内核启动时，我们配置了网格（grid）和线程块（block）。每个线程可以通过内置变量（如blockIdx、threadIdx、blockDim、gridDim）来区分自己，并通过公式 i = blockIdx.x * blockDim.x + threadIdx.x 计算全局索引。

CUDA代码使用NVIDIA的nvcc编译器进行编译，它会自动分离主机（CPU）代码和设备（GPU）代码。

多维网格与数据布局

本节中，我们将探讨如何创建多维网格以及如何存储和访问多维数据。许多应用中的数据都是多维的（如图像、矩阵），使用多维网格可以简化这类数据的处理。

示例1：RGB转灰度图像

这个计算的目标是将一个彩色图像（每个像素由红、绿、蓝三个分量组成）转换为灰度图像（每个像素只有一个表示强度的值）。

一个直观的数据并行化方法是为每个像素分配一个线程。图像是一个二维数组（宽度 x 高度），因此使用二维线程网格来处理会更加方便。

CUDA支持创建多维网格（一维、二维或三维）。这简化了多维数据的处理逻辑。我们可以将网格配置为二维的线程块数组，每个线程块本身也是二维的线程数组。

以下是配置多维网格的代码示例：

// 定义每个线程块的维度（例如 32x32 线程）
dim3 numThreadsPerBlock(32, 32); // Z维度默认为1

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/20daff6152991c4e587c4c7e35605976_27.png)

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/20daff6152991c4e587c4c7e35605976_29.png)

// 计算网格的维度（线程块的数量）
// 假设图像宽度为 `width`，高度为 `height`
int numBlocksX = (width + numThreadsPerBlock.x - 1) / numThreadsPerBlock.x;
int numBlocksY = (height + numThreadsPerBlock.y - 1) / numThreadsPerBlock.y;
dim3 numBlocks(numBlocksX, numBlocksY);

![](https://github.com/OpenDocCN/dsai-notes-zh/raw/master/docs/blt-gpu-comp/img/20daff6152991c4e587c4c7e35605976_31.png)

// 启动内核
rgbToGrayKernel<<<numBlocks, numThreadsPerBlock>>>(...);

在内核函数中，每个线程需要确定自己负责处理图像中的哪个像素。对于二维网格，内置变量（blockIdx, threadIdx, blockDim, gridDim）都有 .x, .y, .z 分量。

__global__ void rgbToGrayKernel(...) {
    // 计算线程对应的输出像素的行和列
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    // ...
}

数据的内存布局：行主序

在C语言中，动态分配的多维数组在内存中实际上是以一维数组的形式连续存储的。常见的存储约定是行主序，即先存储第一行的所有元素，接着是第二行，依此类推。

因此，给定行（row）和列（col），在线性数组中的索引计算公式为：
index = row * width + col

在我们的内核中，转换逻辑如下：

if (row < height && col < width) {
    int index = row * width + col;
    // 加权平均计算灰度值（示例权重）
    gray[index] = 0.3f * red[index] + 0.59f * green[index] + 0.11f * blue[index];
}

边界检查至关重要。因为图像尺寸可能不是线程块尺寸的整数倍，启动的线程数可能多于像素数。必须确保线程只在有效的像素范围内进行操作。

示例2：图像模糊处理

图像模糊的一种简单方法是，每个输出像素的值是其周围输入像素（例如一个3x3或5x5的区域）的平均值。

并行化策略是：为每个输出像素分配一个线程。每个线程读取其对应输入像素及其周围像素，计算平均值，并写入输出像素。

内核实现步骤与RGB转灰度类似，但包含嵌套循环来计算局部平均值：

计算输出像素的行列坐标（outRow, outCol）。
进行边界检查，确保输出坐标有效。
循环遍历以该像素为中心的周围区域（例如，从 outRow - blurSize 到 outRow + blurSize）。
在循环内部，必须对每个输入的访问进行边界检查，因为边缘像素的周围区域可能超出图像范围。
累加有效输入像素的值，最后除以总的有效像素数（或固定区域大小）得到平均值。

关键点：每次内存访问都应有对应的边界检查守卫。这是并行编程中常见错误来源，需要特别注意。

示例3：矩阵乘法

矩阵乘法 C = A x B 中，输出矩阵 C 的每个元素是矩阵 A 的一行与矩阵 B 的一列的点积。

一个高效的并行化方法是：为输出矩阵 C 的每个元素分配一个线程。这样，每个线程独立计算一个点积。

内核实现如下：

__global__ void matrixMulKernel(float* A, float* B, float* C, int N) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < N && col < N) {
        float sum = 0.0f;
        for (int i = 0; i < N; ++i) {
            // A中第row行，第i列的元素
            // B中第i行，第col列的元素
            sum += A[row * N + i] * B[i * N + col];
        }
        C[row * N + col] = sum;
    }
}

注意：

这里我们只并行化了输出矩阵的两个维度（行和列）。点积计算本身的循环（i）在每个线程内部串行执行。这是因为该循环存在跨迭代的依赖（累加和），并行化更复杂。
同样需要处理边界条件，特别是当矩阵维度不是线程块维度的整数倍时。
对于更一般的矩阵乘法（A: M x K, B: K x N, C: M x N），索引计算和边界检查需要相应调整。

总结

本节课我们一起学习了如何使用CUDA处理多维数据。主要内容包括：

创建多维网格：使用 dim3 类型配置二维或三维的线程块和网格维度，从而更自然地映射到像图像、矩阵这样的多维数据。
多维数据的内存布局：理解了C语言中多维数组通常以“行主序”方式在内存中线性存储，并掌握了通过 index = row * width + col 公式计算线性索引的方法。
边界条件处理：这是并行编程中的关键环节。必须确保每个线程的内存访问（无论是读取输入还是写入输出）都在数组的有效范围内，通常通过比较计算出的索引与数组维度来实现。
通过实例加深理解：我们通过RGB转灰度、图像模糊和矩阵乘法三个具体例子，实践了从问题分析、并行策略制定到内核代码实现和边界检查的完整流程。

处理多维数据是GPU计算中的常见任务。掌握多维网格的配置和数据索引技巧，是编写高效、正确CUDA程序的基础。在后续课程中，我们将学习更多优化技术来提升这类计算的性能。

GPU计算：04：GPU架构

概述

在本节课中，我们将暂时告别代码编写，深入探讨GPU的硬件架构。我们将了解GPU如何组织其计算单元，以及CUDA编程模型中的线程、块和网格是如何映射到这些硬件上执行的。理解这些概念对于编写高效的GPU程序至关重要。

从多维数据到硬件架构

上一节我们介绍了如何处理多维数据，并编写了三个数据并行的例子。本节中，我们来看看这些线程和块是如何在GPU的物理硬件上调度和执行的。

GPU的组成结构

一个GPU由多个流式多处理器组成，我们简称为SM。每个SM内部包含多个核心，这些核心是执行算术运算的单元。同一个SM内的核心共享控制逻辑和一部分内存。所有SM都可以访问同一个全局内存，也就是我们在第一讲中从CPU复制数据到GPU时使用的内存。

以本课程使用的Volta V100 GPU为例，它拥有80个SM，每个SM有64个核心，总计5120个处理核心。

线程块在SM上的调度

那么，我们编写的网格和线程块是如何在GPU上运行的呢？线程是以块为单位分配给SM的。同一个线程块内的所有线程必须被分配到同一个SM上执行。然而，一个SM可以同时执行多个线程块。

线程和线程块在执行时需要资源，例如寄存器和线程特定的控制数据。由于SM的硬件资源是有限的，因此它能同时支持的线程和线程块数量也是有限的。这解释了为什么线程块中的线程数存在上限（例如1024）。

如果一个网格启动的线程块数量超过了GPU能同时执行的数量，多余的线程块会排队等待，直到有SM空闲出来。

线程块内的协作与同步

线程在同一个块内被分配到同一个SM，这使得它们能够以不同块之间线程无法实现的方式进行协作。

以下是线程块内协作的主要方式：

屏障同步：线程可以调用 __syncthreads() 函数。这个屏障会确保块内的所有线程都到达代码中的这一点后，才允许任何线程继续执行。这有助于线程间协调工作。
共享内存：每个SM上都有一块共享内存，同一线程块内的线程可以访问这块共享内存，而不同块之间的线程则不能。这为线程间快速交换数据提供了可能。

这种协作能力是将网格划分为块的动机之一。当同一块内的线程在同一个SM上时，支持它们之间的同步和通信会更加高效。

透明可扩展性

线程块之间不进行同步或协作，这使得它们彼此独立。这种独立性带来了一个关键优势：透明可扩展性。

这意味着同一段代码可以在拥有不同数量SM的GPU设备上运行。在SM较少的设备上，线程块可能会被顺序执行；在SM较多的设备上，更多的线程块则可以并行执行。代码无需修改即可适应不同硬件规模的并行能力。

程序员不应编写试图在不同块之间进行同步的代码。因为块可能以任何顺序执行，甚至可能不会同时驻留在SM上，强行同步可能导致死锁。硬件不支持在块之间进行高效的上下文切换，因为这需要保存和恢复大量线程的上下文（如寄存器），开销巨大。

Warp：SM上的调度单元

当一个线程块被调度到SM上后，它会被进一步划分为更小的单元，称为 Warp。Warp是SM上进行调度的基本单位。

Warp的大小是设备特定的，但迄今为止一直是32个线程。例如，一个包含1024个线程的线程块会被分成32个Warp，每个Warp包含32个线程。

Warp的特殊之处在于，同一个Warp内的所有线程以单指令多数据的模式执行。这意味着硬件会取一条指令，然后让Warp内的所有线程执行这条相同的指令，但每个线程操作的是自己的数据。

SIMD的优势在于，它可以将取指、译码等控制逻辑的开销分摊到多个执行单元（核心）上，从而用更少的控制硬件支持更多的计算核心。

控制流分歧

SIMD模式的一个主要缺点是处理控制流分歧。当Warp内的线程遇到条件分支（如if-else语句）并需要执行不同的路径时，问题就出现了。

由于所有线程必须执行相同的指令，硬件会依次处理每条唯一路径：先执行“then”路径，此时条件为“假”的线程被禁用（不执行但占用资源）；然后执行“else”路径，此时条件为“真”的线程被禁用。最后，所有线程在分支结束后重新汇合。

这导致了硬件利用率下降，称为SIMD效率降低。在循环中，如果不同线程的循环次数不同，情况会更糟，因为所有线程都必须等待迭代次数最多的那个线程完成。

程序员可以通过组织数据或线程ID来尽量减少同一个Warp内的控制流分歧，例如对数据进行排序，让具有相似执行路径的线程聚集在同一个Warp里。

延迟隐藏与多线程

当Warp执行遇到长延迟操作（如未命中的内存访问或复杂的多周期算术运算）时，它需要等待。为了避免流水线停顿，GPU采用了多线程技术。

SM上会同时驻留许多Warp（远多于其核心数）。当一个Warp因长延迟操作而停滞时，调度器会迅速将其换出，并换入另一个就绪的Warp到核心上执行。通过这种方式，GPU的计算核心可以始终保持忙碌，从而隐藏操作延迟，实现高吞吐量。

例如，Volta V100的每个SM最多可支持2048个线程（即64个Warp），但其只有64个核心。这种高“线程-核心比”（32:1）确保了总有足够多的就绪Warp来填充因延迟而产生的空闲周期，这与CPU（通常为2:1）追求低延迟的设计哲学不同。

占用率

SM的占用率是指其实际活跃的Warp（或线程）数量与硬件支持的最大数量之比。通常，我们希望最大化占用率，因为更高的占用率意味着有更多可用的Warp来隐藏延迟。

然而，占用率可能受到多种因素的限制：

线程块大小：如果线程块大小设置不当，可能无法充分利用SM的线程容量。例如，若块大小太小，可能会受限于SM支持的最大块数量；若块大小不能被最大线程数整除，会导致线程资源浪费。
寄存器使用量：如果每个线程使用了大量寄存器，可能会在达到最大线程数之前就用完SM上的寄存器资源。
共享内存使用量：过度使用共享内存也会限制SM上可同时驻留的线程块数量。

因此，在选择线程块大小时，需要综合考虑这些约束，以优化占用率和最终性能。CUDA提供了 cudaGetDeviceProperties API，供程序员查询设备的这些具体限制。

总结

本节课我们一起学习了GPU的核心架构。我们了解了GPU由多个SM组成，SM内部包含众多核心。CUDA中的线程块被调度到SM上执行，并进一步划分为Warp，Warp以SIMD模式运行。我们探讨了控制流分歧带来的挑战、通过多线程和Warp调度来隐藏延迟的策略，以及占用率的概念及其对性能的影响。理解这些硬件原理是进行GPU性能优化的基础。

GPU计算：05：内存与分块优化

在本节课中，我们将学习GPU的内存架构、CUDA编程模型中的内存管理，以及一个重要的性能优化技术——分块（Tiling）。我们将通过矩阵乘法的例子，详细探讨如何利用共享内存来减少全局内存访问，从而提升程序性能。

概述：性能指标与瓶颈

上一节我们介绍了GPU的架构和线程调度。本节中，我们来看看GPU的内存系统。处理器设计者通常通过两个关键指标来告知用户其性能极限：

峰值浮点运算速率 (FLOPS)：处理器每秒能执行的浮点运算次数。这反映了计算核心的最大理论计算能力。
峰值内存带宽：内存每秒能向处理器核心提供的数据字节数。这反映了内存系统的最大理论数据传输能力。

例如，V100 GPU的峰值FLOPS约为14 TeraFLOPS，峰值内存带宽约为900 GB/s。这些是理论峰值，实际程序很难达到。它们为我们评估代码性能提供了参考基准。

根据程序对这两个资源的依赖程度，我们可以将其分为两类：

计算密集型：程序性能受限于处理器的计算能力。此时计算核心被充分利用，而内存相对空闲。
内存密集型：程序性能受限于内存带宽。此时计算核心经常因等待数据而空闲。

一个重要的概念是期望的计算与全局内存访问比，其公式为：
期望比率 = 峰值FLOPS / 峰值内存带宽
对于V100，这个比率约为15.6次浮点运算/字节。这意味着，为了充分利用计算核心，平均每从内存加载1字节数据，就需要执行约15.6次浮点运算。如果程序的实际比率远低于此，则很可能受内存带宽限制。

案例分析：向量加法与矩阵乘法

让我们分析两个内核的计算与内存访问比。

向量加法：每个线程执行 z[i] = x[i] + y[i]。

加载：2个单精度浮点数 = 8字节。
计算：1次浮点加法。
比率：1次运算 / 8字节 = 0.125 运算/字节。
这个比率极低，因此向量加法是典型的内存密集型应用，增加计算核心对其加速有限。

基础矩阵乘法：每个线程计算一个输出元素，通过循环累加。

每次循环迭代：加载2个浮点数（8字节），执行1次乘法和1次加法（2次运算）。
比率：2次运算 / 8字节 = 0.25 运算/字节。
这个比率仍然很低，但矩阵乘法本身具有很高的数据复用潜力。对于N×N的矩阵乘法，总运算量为2N³，而理论上最少只需加载8N²字节的输入数据。因此，其潜在比率为 (2N³) / (8N²) = 0.25N 运算/字节。当N较大时，这个值可以很高。我们当前的实现没有利用这种复用，导致了大量冗余的全局内存访问。

GPU内存架构与CUDA内存模型

在深入优化之前，我们需要了解GPU的内存层次结构。

GPU的流多处理器（SM）包含以下关键内存组件：

寄存器：速度最快，线程私有。访问延迟约为1个周期。
共享内存：块内线程共享，由程序员显式管理。访问延迟约数十周期。
L1缓存：硬件自动管理，用于缓存全局内存数据。
常量缓存：用于存储只读的常量数据。
L2缓存：所有SM共享的末级缓存。
全局内存：容量最大、速度最慢（访问延迟约数百周期），所有线程可访问。

在CUDA编程模型中，程序员可以通过以下关键字指定变量的存储位置：

__device__：变量位于全局内存。所有网格中的线程共享同一副本，生命周期为整个应用程序。
__constant__：变量位于常量内存。所有网格中的线程共享同一只读副本。
__shared__：变量位于共享内存。每个线程块拥有该变量的独立副本，块内线程共享，生命周期随块结束而结束。
局部变量：默认位于寄存器。每个线程拥有私有副本。

共享内存是程序员可以主动管理的快速内存，是进行分块优化的关键。

分块优化原理

矩阵乘法中存在大量数据复用机会。例如，计算输出矩阵中同一行的不同元素时，会重复使用输入矩阵A的同一行数据；计算同一列的不同元素时，会重复使用输入矩阵B的同一列数据。

在基础的并行实现中，每个线程独立地从全局内存加载所需的所有A行和B列元素，导致了大量重复加载。分块优化的核心思想是让同一个线程块内的线程协作，将输入数据的一块（Tile）加载到快速的共享内存中，然后所有线程从共享内存中重复访问这些数据。

具体步骤如下：

将输入和输出矩阵划分为大小相等的块（例如32×32）。
每个线程块负责计算输出矩阵中的一个块。
对于计算该输出块所需的输入数据，也对应地分成多个“片”。
线程块内的所有线程协作，将当前需要的一个输入A片和一个输入B片从全局内存加载到共享内存。
调用 __syncthreads() 确保所有线程完成加载。
每个线程从共享内存中读取数据，进行部分乘积累加计算。
调用 __syncthreads() 确保所有线程完成对当前共享内存片的使用。
重复步骤4-7，直到处理完所有相关的输入片。
将最终结果写回全局内存。

通过这种方式，每个数据元素从全局内存只加载一次（由某个线程负责），然后被同一个线程块内的多个线程从共享内存中多次访问，显著减少了全局内存访问次数。

分块矩阵乘法代码实现

以下是利用共享内存实现分块矩阵乘法的CUDA内核代码示例：

__global__ void matrixMultiplyShared(float* A, float* B, float* C, int n) {
    // 定义块大小（Tile维度）
    const int TILE_DIM = 32;

    // 声明共享内存数组，用于存储A和B的一个块
    __shared__ float As[TILE_DIM][TILE_DIM];
    __shared__ float Bs[TILE_DIM][TILE_DIM];

    // 计算当前线程负责的输出元素坐标
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    float sum = 0.0f;

    // 循环遍历所有输入块
    for (int t = 0; t < n / TILE_DIM; ++t) {
        // 协作加载A的一个块到共享内存
        // 每个线程加载一个元素：A中对应行，第t块中的对应列
        As[threadIdx.y][threadIdx.x] = A[row * n + (t * TILE_DIM + threadIdx.x)];
        // 协作加载B的一个块到共享内存
        // 每个线程加载一个元素：B中第t块中的对应行，对应列
        Bs[threadIdx.y][threadIdx.x] = B[(t * TILE_DIM + threadIdx.y) * n + col];

        // 等待块内所有线程完成共享内存加载
        __syncthreads();

        // 使用共享内存中的数据计算部分乘积和
        for (int i = 0; i < TILE_DIM; ++i) {
            sum += As[threadIdx.y][i] * Bs[i][threadIdx.x];
        }

        // 等待块内所有线程完成当前块的计算，然后再加载下一个块
        __syncthreads();
    }

    // 将最终结果写回全局内存
    if (row < n && col < n) {
        C[row * n + col] = sum;
    }
}

关键点说明：

__syncthreads() 是必要的屏障同步，确保加载或计算步骤的线程间一致性。
边界条件处理：在写入最终结果 C 时，需要检查线程坐标是否在矩阵有效范围内。但在加载共享内存时，参与加载的线程可能超出输出矩阵范围，这是允许且有时是必要的，以确保整个数据块被完整加载。
块大小（TILE_DIM）的选择需要权衡：更大的块能提高数据复用，但也会消耗更多共享内存，可能降低SM的线程块驻留数量（占用率）。

分块在CPU上的应用

分块优化并非GPU专属，它同样适用于CPU。在CPU上，分块主要利用的是硬件管理的多级缓存（L1、L2、L3）来提升数据局部性。其代码结构与GPU版本有很好的对应关系：外层循环遍历输出块和输入块，内层循环遍历块内的元素进行计算。通过分块，CPU版本的矩阵乘法也能获得显著的性能提升。

注意事项与总结

在使用共享内存和分块优化时，需要注意以下几点：

占用率影响：共享内存是SM上的有限资源。使用过多共享内存可能会限制SM上同时活跃的线程块数量，从而影响占用率和延迟隐藏能力。
动态共享内存：除了示例中的静态声明，CUDA也支持动态分配共享内存，这在编译时无法确定共享内存大小时非常有用。
复杂边界处理：当矩阵维度不是块大小的整数倍时，或处理非方阵乘法时，边界条件的处理会变得更加复杂，需要仔细编码。

本节课中我们一起学习了GPU的内存层次结构、CUDA内存模型，并深入探讨了利用共享内存进行分块优化的原理与实践。通过将数据从慢速的全局内存搬运到快速的共享内存中复用，我们能够显著降低内存密集型应用（如矩阵乘法）的全局内存访问压力，从而更有效地利用GPU强大的计算能力。分块是GPU编程中一项基础而强大的优化技术。

GPU计算：06：性能考量

在本节课中，我们将学习GPU计算中的更多性能考量。我们将回顾之前讨论过的内存架构和CUDA编程模型，并深入探讨两种关键的优化技术：内存合并与线程粗化。理解这些概念对于编写高效的GPU程序至关重要。

回顾：内存优化与矩阵乘法

上一节我们介绍了GPU架构中的内存层次结构和CUDA编程模型，并学习了如何使用共享内存分块技术来优化内存访问，并将其应用于矩阵乘法示例。

我们提到，处理器设计者使用各种指标来描述其性能，例如峰值浮点运算速率和峰值内存带宽。前者指处理器每秒能执行的浮点运算次数，后者指处理器内存每秒能提供的字节数。

应用程序根据其首先达到峰值浮点运算速率还是峰值内存带宽，可以分为计算密集型或内存密集型。在当前的处理器（如Volta V100）上，峰值浮点运算速率通常远高于内存带宽。这意味着，为了平衡两者，每加载一个字节的数据，大约需要执行15.6次运算（或每加载一个浮点值需要约60次运算）。如果程序的计算量低于这个比例，则更可能受限于内存带宽；反之，则更可能受限于计算能力。

在矩阵乘法内核的初始版本中，我们每加载8字节数据（A和B矩阵各4字节）执行2次浮点运算（一次乘法和一次加法）。这给出了一个计算与全局内存访问比为 0.25 ops/byte，远低于V100所需的 15.6 ops/byte。

为了改善这个比率，我们引入了共享内存。在GPU架构中，除了全局内存和硬件管理的缓存，每个流多处理器（SM）上还有一块更小但更快的共享内存。在CUDA模型中，同一线程块内的线程可以协作，将数据从全局内存加载到共享内存，然后多次从共享内存访问，从而最小化全局内存访问。

在矩阵乘法中，我们观察到处理同一输出块（Tile）的线程共享相同的输入行（来自A）和输入列（来自B）。通过分块，我们可以让线程块内的每个线程负责加载A和B输入块中的一个元素到共享内存。然后，所有线程从共享内存（而非全局内存）访问数据来进行部分点积计算。这种方法显著提升了性能。

在实现时，我们需要仔细处理边界条件。每个内存访问（加载A块、加载B块、存储输出C）都需要有其对应的边界检查。此外，共享内存和寄存器的使用会影响占用率，因为每个SM上的共享内存是有限的。

DRAM架构与内存合并

现在，让我们深入了解DRAM（动态随机存取存储器）的架构，这将引出我们今天要讨论的第一个优化：内存合并。

DRAM单元可以逻辑上视为一个存储电荷的电容器和一个三态器件。读取时，若值为1，电容器放电可被检测到；若值为0，则无放电。

多个DRAM单元连接在一根列线上，构成一个DRAM阵列。一个DRAM存储体由许多这样的垂直列线组成。在任何时刻，我们激活一行（通过行解码器），该行所有单元的电容器放电，其值通过感应放大器读取并暂存到列锁存器中，然后写回（因为读取会破坏数据）。这一整行数据被称为一个DRAM突发传输。

地址被分为行地址和列地址。行地址用于选择激活哪一行，列地址则通过多路复用器从整个突发传输中选择我们请求的特定部分。

访问DRAM阵列（包括放电、感应、写回）是整个过程最慢的部分。而从已读入列锁存器的突发传输中选择数据（多路复用）则快得多。

关键点在于：访问同一突发传输内的数据比访问不同突发传输的数据要快得多。

这引出了内存合并的概念。当同一个Warp中的线程访问连续的内存位置时，这些位置很可能位于同一个DRAM突发传输内。这些访问可以被合并，由一个DRAM事务服务。然而，如果同一个Warp中的线程访问分散在不同突发传输中的位置，则无法合并，需要多个内存事务，导致性能下降，这有时被称为内存发散。

在我们之前的例子中，向量加法和矩阵乘法的内存访问模式都天然地实现了内存合并。在向量加法中，同一个Warp中的线程访问数组X的连续元素（如 x[0], x[1], ..., x[31]）。在矩阵乘法中，由于Warp由连续的 threadIdx.x 线程组成，它们访问的B矩阵元素也是连续的，而访问的A矩阵元素则是相同的，因此访问也是合并的。

多存储体与占用率

通常，DRAM被组织成多个存储体。这有几个原因：制造单一巨大阵列不可行，并且多存储体允许我们在访问DRAM时实现并行性。

当拥有多个DRAM存储体时，可以在一个存储体读取突发传输的同时，让其他存储体进行下一次读取操作。这样，通过让多个存储体同时工作，可以隐藏内存访问延迟。

为了有效利用多个DRAM存储体，我们需要同时提供大量的内存访问请求。这需要许多线程同时运行并请求数据，以保持DRAM存储体繁忙。

这提醒我们高占用率的重要性。高占用率不仅有助于隐藏核心流水线的延迟，还能确保有足够多的内存访问请求来隐藏访问DRAM的延迟。

线程粗化

到目前为止，我们采用的并行化策略是让线程尽可能细粒度，即将线程分配给最小的并行单元。例如，向量加法中每个线程处理一个元素，矩阵乘法中每个线程处理一个输出元素。

细粒度线程的优势在于，它为硬件提供了尽可能多的线程，以充分利用所有计算资源，并实现透明可扩展性：硬件可以根据资源情况并行执行或低开销地串行化这些线程。

然而，细粒度线程的劣势在于，如果线程之间存在冗余操作（例如，多个线程块加载相同的数据），并且这些线程块最终被硬件串行执行，那么为并行化付出的冗余代价就白费了。

线程粗化优化正是为了解决这个问题。它让一个线程（或线程块）负责多个输出元素（或多个输出块），从而在串行执行时重用已加载的数据或已进行的计算，减少冗余。

让我们以矩阵乘法为例进行线程粗化。原本，每个线程块负责输出矩阵的一个分块。水平相邻的两个线程块会加载相同的A矩阵分块，但不同的B矩阵分块。如果这两个线程块在硬件上并行执行，这种冗余加载是值得的。但如果它们被串行执行，我们就可以让一个线程块顺序处理这两个水平相邻的输出块，从而只加载一次A分块并重复使用。

以下是实现线程粗化的关键代码修改步骤：

定义粗化因子：例如 coarsening_factor = 4，表示每个线程负责4个输出元素。
调整线程块网格：在x维度上，线程块的数量需要除以粗化因子。
计算起始列：每个线程负责的第一个元素的列索引需要根据粗化因子调整：col_start = blockIdx.x * blockDim.x * coarsening_factor + threadIdx.x。
使用数组存储部分和：每个线程现在需要维护一个大小为 coarsening_factor 的数组来存储多个输出元素的部分和。
嵌套循环：外层循环遍历A的输入分块。对于每个A分块，内层循环（粗化循环）遍历多个B分块（数量等于粗化因子）。在内层循环中，根据循环索引计算当前B分块对应的列，加载B分块，然后计算并累加到对应的部分和数组中。
结果写回：最后，通过一个循环将数组中的多个结果写回全局内存。

通过应用线程粗化，我们减少了冗余的A矩阵加载。性能提升取决于硬件和选择的粗化因子。因子太小可能优化不足，因子太大则可能导致资源利用不足或序列化过度，反而损害性能。

线程粗化的优势在于减少了为并行化付出的代价（如冗余内存加载、计算或同步）。劣势在于可能降低资源利用率（如果粗化过度），破坏了透明可扩展性（需要为不同设备调整因子），并且每个线程需要更多资源（如寄存器）。

优化清单与瓶颈分析

至此，我们学习了一系列常见的GPU优化技术：

调整资源以最大化占用率
最小化控制流发散
内存合并（确保统一的内存访问模式）
共享内存分块（捕获数据局部性）
线程粗化（减轻并行化开销）
私有化（处理输出竞争，将在后续课程介绍）

在应用这些优化时，需要注意它们之间可能存在冲突。例如：

最大化占用率可能与缓存抖动控制冲突。
使用大量共享内存进行分块可能会限制占用率。
线程粗化减少了冗余工作，但可能增加每个线程的资源使用，从而限制占用率。

因此，找到最佳的平衡点至关重要。关键在于识别应用程序的瓶颈。瓶颈可能是内存带宽、计算能力、占用率、控制流发散或原子操作竞争等。在优化之前，应先诊断瓶颈所在，然后选择针对性的优化策略。优化错误的瓶颈可能无法提升性能，甚至适得其反。

总结

本节课中，我们一起深入探讨了GPU性能优化的两个重要方面。
首先，我们通过理解DRAM的突发传输机制，学习了内存合并的原理及其对性能的关键影响。
其次，我们引入了线程粗化技术，通过让单个线程处理更多工作来减少并行化带来的冗余开销，并在矩阵乘法示例中实践了其实现方法。
最后，我们回顾了完整的优化清单，并强调了根据应用程序具体瓶颈来选择和平衡优化策略的重要性。掌握这些知识将帮助你编写出更高效的GPU程序。

GPU计算：07：性能剖析

在本节课中，我们将学习如何剖析CUDA应用程序，以了解其性能表现并识别瓶颈。我们将使用NVIDIA提供的性能剖析工具，通过一个向量加法的实例，学习如何收集和分析性能数据，从而指导我们进行有效的优化。

概述

上一节我们讨论了DRAM内存的组织结构、内存合并与分散访问、以及线程粗化等优化技术。我们了解到，不同的优化策略（如提高占用率、使用共享内存分块）之间可能存在权衡。因此，在应用优化之前，准确识别应用程序的性能瓶颈至关重要。本节我们将学习如何使用性能剖析工具来“了解你的瓶颈”。

性能剖析工具简介

CUDA平台提供了一个名为nvprof的性能剖析器。它可以测量应用程序中各种操作（如内核执行、内存拷贝）所花费的时间，并收集更详细的硬件性能指标。

基本时间线剖析

最简单的使用方式是直接运行nvprof并指定你的可执行文件。它会输出一个表格，列出所有CUDA API调用和GPU活动及其耗时。

nvprof ./vector_add

输出结果会显示：

GPU活动：如内核执行、设备内存拷贝。
API调用：如cudaMalloc，cudaMemcpy等CPU端调用。
每项活动的耗时、调用次数以及占总执行时间的百分比。

这种方式可以快速获得应用程序的时间分布概览。

导出数据以进行可视化分析

为了进行更深入的分析，我们可以将剖析数据导出到文件，并使用NVIDIA Visual Profiler (nvvp) 进行可视化。

首先，导出时间线数据：

nvprof -o timeline.prof ./vector_add

接着，导出详细的性能指标（如占用率、缓存命中率、内存带宽）。这需要使用硬件计数器，因此剖析器会多次运行内核以收集不同指标，导致总执行时间变长。

nvprof --metrics all -o metrics.prof ./vector_add

现在，我们可以启动可视化性能剖析器并导入生成的文件：

nvvp &

在nvvp中，选择“导入”功能，然后指定timeline.prof和metrics.prof文件。

分析向量加法内核

让我们以向量加法内核为例，使用nvvp进行分析。

时间线视图

时间线视图直观地展示了应用程序中各种活动的顺序和持续时间。我们可以看到：

cudaMemcpy（主机到设备）用于传输输入数组。
vector_add内核的执行。
cudaMemcpy（设备到主机）用于取回结果。
这个视图有助于理解整体执行流程和识别大的时间开销所在。

内核性能分析

点击分析单个内核，我们可以获得其配置和资源使用情况：

网格与块大小：例如，(N/512, 1, 1) 网格，(512, 1, 1) 块。
寄存器使用量：每个线程使用的寄存器数量。
共享内存使用量：每个线程块使用的共享内存量。
占用率：实际占用率与理论最大占用率。

对于我们的向量加法，由于块大小足够大（512），且未使用共享内存，寄存器使用量低，因此理论占用率可达100%，实际占用率也接近这个值（例如88.1%）。

识别性能瓶颈

执行“内核分析”功能，工具会显示计算单元和内存系统的利用率。

对于向量加法，分析结果显示：

内存系统利用率很高（>80%）。
计算单元利用率很低（约10%）。
结论提示：“内核性能受内存带宽限制”。

这验证了我们之前的分析：向量加法是内存密集型操作，计算与内存访问的比率很低，因此瓶颈在于全局内存带宽。

以下是更详细的分析维度：

内存带宽分析
此视图展示了内存层次结构中各部分的利用率。对于向量加法：

设备内存（全局内存） 利用率接近峰值。
L2缓存、统一缓存、共享内存 利用率很低。
图表清晰地显示，几乎所有请求都流向全局内存。

计算分析
此视图展示了核心内部各功能单元的使用情况。

负载/存储单元、单精度浮点单元、控制流单元的利用率都较低。
指令计数显示，执行了大量的加载/存储指令和整数指令（用于地址计算），相对较少的浮点运算指令。

延迟分析
此视图专注于占用率，这是隐藏延迟的关键。

对于初始配置（块大小512），占用率不是限制因素。
工具会显示，如果更改块大小或共享内存使用量，将如何影响理论占用率。

一个占用率受限的例子

如果我们故意将线程块大小设置为一个很小的值（例如32），重新编译并剖析：

nvprof -o timeline_occ.prof ./vector_add_small_block
nvprof --metrics all -o metrics_occ.prof ./vector_add_small_block

在nvvp中分析新数据：

占用率分析会显示，理论最大占用率降至50%（因为SM上最大线程数2048 / 每块线程数32 = 最大块数64，但SM的块数量限制可能更低，例如32，导致最大线程数为1024）。
延迟分析视图会明确提示“GPU利用率可能受块大小限制”，并用红色高亮块数量限制。
内核执行时间会显著增加，因为低占用率无法有效隐藏内存访问延迟。

这个例子说明了如何利用剖析工具来诊断由资源限制（此处为块大小）导致的低占用率问题。

总结

本节课我们一起学习了CUDA性能剖析的基本方法。我们了解到：

nvprof命令行工具可以快速获取应用程序的时间分布和性能指标。
NVIDIA Visual Profiler (nvvp) 提供了强大的可视化界面，用于深入分析性能瓶颈。
通过剖析，我们可以确定内核是受计算限制还是受内存带宽限制。
我们可以分析内存层次结构的利用率、计算单元利用率以及占用率，从而找到优化的方向。
性能剖析是理解优化权衡和做出正确优化决策的关键步骤。在应用共享内存分块或线程粗化等优化之前，务必先了解当前的瓶颈所在。

记住，优化的目标是找到在竞争性优化策略之间的最佳平衡点，而性能剖析正是照亮这条道路的灯塔。

GPU计算：08：卷积与常量内存

在本节课中，我们将开始学习并行模式，第一个要讨论的模式是卷积。我们还将利用卷积这个例子，来介绍GPU架构和CUDA编程模型中的一个新特性：常量内存。

课程回顾

在进入新的主题之前，我们先简要回顾一下之前学过的内容。

到目前为止，我们已经学习了GPU计算的基础知识。课程初期，我们探讨了由于频率停滞和功耗墙的限制，单线程性能发展停滞，这推动了并行计算的发展。我们比较了CPU（面向延迟的设计）和GPU（面向吞吐量的设计）的差异。

CPU（面向延迟）：专注于使单个任务尽可能快，使用强大的ALU、大容量缓存、复杂的控制流（如分支预测、数据转发）和少量多线程来隐藏延迟。
GPU（面向吞吐量）：拥有大量较简单的ALU，通过大规模并行来获得高吞吐量。缓存较小，控制逻辑更简单，从而能将更多芯片面积用于计算单元。通过海量线程来隐藏操作的高延迟。

我们了解了典型系统中GPU如何与CPU协同工作：CPU拥有主内存，GPU拥有设备内存。通常的流程是：在GPU上分配内存，将数据从CPU主内存复制到GPU内存，在GPU上执行内核（访问设备内存），然后将结果复制回CPU，最后释放GPU内存。虽然也有像统一内存这样的方式，但这是使用GPU的典型模式之一。

我们从简单的向量加法内核开始，学习了数据并行性，并看到它非常适合GPU。我们了解了网格如何组织成线程块和线程，以及如何计算线程索引，使每个线程对不同数据执行相同操作。

接着，我们学习了网格可以组织成多维线程数组（如2D、3D网格），以及多维数据通常按行主序存储。为了访问动态分配的多维数组，我们需要将线程获得的二维索引转换为一维索引。我们通过多个例子实践了这一点，包括RGB转灰度、图像模糊（这是卷积的一个特例）以及矩阵乘法。矩阵乘法的例子在后续几讲中持续被用作优化案例。

在学习了CUDA编程模型的基础后，我们转向讨论GPU架构。我们了解到GPU被组织成多个流式多处理器（SM），每个SM包含多个共享内存和控制单元的核，所有SM都能访问相同的全局内存。

网格在SM上调度时，被划分为线程块。调度以块为单位，这意味着同一个块内的所有线程会在同一个SM上同时运行。当然，一个SM上可以同时运行多个线程块。线程块在SM上被进一步划分为线程束（Warp），调度以线程束为单位。线程束的一个特点是，同一个线程束内的线程遵循SIMD模型执行，即所有线程在不同数据上执行相同的指令。

线程束大小自GPU引入以来一直是32，但这在未来可能会改变。由于线程束以SIMD方式执行，这引出了控制流发散的问题。如果同一个线程束内的线程需要走不同的控制路径（例如if语句），那么所有线程会先执行then部分（不满足条件的线程不活动），然后一起执行else部分（满足条件的线程不活动）。这些不活动的线程占用了计算资源但没有执行工作，导致硬件利用率低下。我们希望尽量减少这种情况，并将在未来的并行模式中看到相关例子。

我们还讨论了延迟隐藏。通过在SM上调度比核心数量更多的线程束和线程，当一个线程束遇到长延迟操作时，可以将其换出并换入另一个线程束执行。通过这种方式，我们可以隐藏延迟。为了隐藏延迟，我们希望有尽可能多的线程调度在SM上，这就引入了占用率的概念。每个SM的最大线程数取决于设备（例如V100是2048），但其他因素也可能限制占用率，例如每个SM能容纳的块数、每个线程使用的寄存器数量或SM的共享内存大小。我们需要关注资源使用情况，以确保最大化SM上的线程数。

接下来，我们研究了GPU架构中的内存。同一SM上的线程可以访问共享内存。我们可以利用共享内存来存放计划重用的数据。例如，在分块矩阵乘法中，我们不是让每个线程加载整行和整列，而是让同一个块内的线程协作加载A和B的一个数据块到共享内存中，然后进行计算。这减少了全局内存访问次数，提高了计算与全局内存访问的比率。

我们还了解了DRAM的组织方式。访问DRAM阵列本身很慢，但一旦开始DRAM突发传输，访问突发数据的一部分是很快的。因此，最好让线程在发出内存请求时访问同一个DRAM突发中的数据。如果访问不同的DRAM突发，就需要反复访问DRAM阵列。如果有多个DRAM阵列，我们可以在从一个DRAM阵列读取时，从另一个DRAM阵列的突发中服务数据，从而隐藏延迟。这也是最大化GPU占用率的另一个动机，因为高占用率不仅能帮助我们隐藏流水线延迟，还能提供大量内存访问请求，有助于隐藏内存访问延迟。

基于以上，我们总结了一个常见优化清单：

调整资源以最大化占用率。
最小化控制流发散。
采用合并访问的内存访问模式。
使用共享内存分块以捕获数据重用。
使用线程粗化以减轻并行化开销。
私有化（我们尚未涉及，将在后续讨论）。

选择应用哪种优化取决于应用程序的性能瓶颈是什么。优化通常是用一种资源换取另一种资源，因此需要确保用充裕的资源去换取瓶颈资源。这意味着我们需要选择合适的优化。为了了解瓶颈，CUDA提供了性能分析工具，可以帮助我们分析应用程序，找出哪些资源被充分利用，哪些没有，从而评估瓶颈所在。

以上是对课程第一部分——GPU计算基础知识的快速回顾。接下来，我们将进入课程的第二部分：并行模式。我们将花相当多的时间讨论这些模式，并希望随着每个模式的介绍，引入新的架构特性或优化技术。因此，学习并行模式的目标不仅是了解模式本身，也是利用这些模式来学习新的优化技术、方法或硬件特性。今天，我们将从卷积开始。

什么是卷积？

卷积是一种运算，它有一个输入。卷积运算的输出是：输出中的每个元素都是对应输入元素及其邻域元素的加权和。

我们之前看到的图像模糊就是卷积的一个特例，其中所有权重都相同。在图像模糊中，我们遍历输入元素并计算平均值。然而，在更一般的卷积形式中，我们可以为每个输入元素设置不同的权重。

这些权重由卷积掩码（Mask）决定。为了避免与CUDA内核函数混淆，这里我们称其为“掩码”而不是“核”。这个掩码由一组权重构成。这些权重被应用于输入，通过计算输入的加权平均（根据掩码中的权重）来得到输出。

卷积在信号处理、图像处理、视频处理等领域非常有用。它通常用于将信号（可能是一维信号、像素或图像等二维数据）转换为更理想的值。例如：

高斯模糊卷积：比我们之前看到的简单模糊更复杂，距离中心越远的像素权重越低。
图像锐化卷积。
边缘检测卷积。

卷积操作所实现的变换取决于掩码中的权重。本节课我们将以二维卷积为例，因为它易于可视化且足够复杂。当然，也存在一维和三维卷积。

如何并行化卷积？

基于我们目前所学，你认为并行化卷积的好方法是什么？一个简单的方法是：为每个输出元素分配一个线程。这个线程负责遍历对应的输入邻域元素和掩码，计算加权和。这是一种非常简单的卷积并行化方法，当然还有其他方法，但今天我们将采用这种方法。

这种方法很直接：每个线程负责一个输出元素，只需一个循环来遍历输入和掩码。

这里需要注意的一点是，掩码通常很小（例如5x5），并且对所有线程都是相同的。所有线程将访问输入的不同部分，但都会访问相同的掩码。因此，在这种情况下，由于掩码很小、所有线程都访问它、并且掩码在内核执行期间不会改变，我们可以将掩码存储在一个称为常量内存的特殊位置。

常量内存简介

之前介绍CUDA编程模型时，我们提到线程可以访问寄存器（同一块内的线程各有自己的寄存器），同一块内的线程可以访问共享内存，网格中的所有线程可以访问全局内存。实际上，网格中的所有线程还可以访问常量内存。

我们希望将掩码放入常量内存，因为访问常量内存比访问全局内存更快。稍后会解释原因。现在，让我们看看如何在CUDA编程模型中使用常量内存。

声明常量内存数组

要在常量内存中声明数组，可以使用 __constant__ 限定符。

__constant__ float mask_c[MASK_DIM][MASK_DIM];

这里，mask_c 是我们给常量内存数组起的名字（使用 _c 后缀是一种约定，表示常量内存）。MASK_DIM 是掩码的维度（例如5）。

从主机复制数据到常量内存

常量内存意味着在GPU执行时不能写入，但可以从CPU复制数据到其中。在主机代码中，我们使用一个特殊的CUDA复制函数 cudaMemcpyToSymbol。

cudaMemcpyToSymbol(mask_c, mask, MASK_DIM * MASK_DIM * sizeof(float));

mask_c：目标指针（常量内存中的数组）。
mask：源指针（主机上的掩码数组）。
第三个参数：要复制的数据大小（掩码元素个数 × 每个元素大小）。

常量内存的特点与优势

常量内存有一个限制：只能分配最多64KB。这一点很重要，否则我们可能会想把整个输入矩阵（如果它们在内核执行期间不变）都放进常量内存。但常量内存的权衡是：速度快，但容量小。

那么，为什么知道数据是常量就能让它更快？为什么容量有限？

缓存效率：为常量数据构建高效缓存更容易。因为不需要支持写入操作，所以无需跟踪数据变化（脏位）和写回机制，实现成本更低。
无需缓存一致性：在并行处理器中，如果多个线程可以写入同一数据，就需要缓存一致性协议来确保不同缓存间的数据一致性。对于常量数据，则不需要支持缓存一致性。
容量小，命中率高：由于总容量小，缓存可以容纳所有常量数据，从而最小化缓存失效，保持低缺失率。

从编程模型看，所有线程访问相同的常量内存。实际上，每个SM内部都有一个常量缓存，它与L1缓存或共享内存是不同的。

为什么不使用共享内存？ 当然可以将掩码放入共享内存，但这需要程序员手动管理（加载到共享内存）。而常量缓存是由硬件管理的，对程序员透明。

编写卷积内核

现在我们已经知道如何分配和使用常量内存，接下来编写卷积内核。它将与图像模糊内核非常相似。

基本思路是：每个线程负责计算一个输出元素。线程首先计算自己负责的输出元素的行列索引，然后检查边界条件（确保线程在输出图像范围内）。接着，线程初始化一个累加器 sum，然后循环遍历掩码的所有元素。对于每个掩码元素，计算对应的输入元素的行列索引，检查该输入索引是否在有效范围内，如果在，则从全局内存加载输入值，与常量内存中的掩码权重相乘，并累加到 sum 中。循环结束后，将 sum 写入输出数组的对应位置。

边界检查是必要的，因为对于边缘的输出元素，其对应的部分输入邻域可能超出图像边界。虽然边界检查会导致一些控制流发散（边缘线程与内部线程路径不同），但由于边缘线程占比较小，这种发散通常是可接受的，且优化收益不大。真正需要优化的是那些导致所有线程都发散的情况。

卷积性能与数据重用

卷积在GPU上性能很好。原因包括：

高度并行化：每个输出元素一个线程，并行度很高。
数据重用：相邻的线程（处理相邻的输出元素）会使用大量重叠的输入数据。例如，处理输出元素(i, j)的线程和使用元素(i, j+1)的线程，它们的输入邻域有很大重叠。

这种数据重用提示我们可以进行优化。主要的优化手段是共享内存分块。

分块卷积优化思路

我们可以让一个线程块协作加载一个输入数据块到共享内存中。每个线程加载输入块中的一个元素。然后，线程块内的线程可以基于共享内存中的数据块和常量内存中的掩码进行计算。

这里的分块与矩阵乘法的分块略有不同。在卷积中，输入块的大小大于输出块的大小。具体来说，如果输出块大小为 T x T，掩码大小为 M x M，那么需要的输入块大小至少为 (T + M - 1) x (T + M - 1)。

因此，线程块的大小（线程数）需要与输入块的大小匹配。在加载阶段，所有线程都参与将输入块加载到共享内存。但在计算阶段，只有一部分线程（对应输出块内的位置）会进行计算并将结果写回全局内存。换句话说，我们超额配置了线程，其中一部分线程仅负责加载数据而不进行计算。

另一种方法是让线程数等于输出块大小，然后让这些线程通过循环来加载更大的输入块，但代码会更复杂，性能提升可能有限。

计算与全局内存访问比率分析

我们通过计算与全局内存访问的比率来量化优化效果。

非分块版本：

每个线程的全局内存加载次数：M² 次（每次加载4字节浮点数）。
每个线程的浮点运算次数：2 * M² 次（每次乘加算2次操作）。
计算与内存访问比率：(2 * M² ops) / (M² * 4 bytes) = 0.5 ops/byte。

分块版本（以线程块为单位分析）：

设输入块大小为 T x T（元素个数）。
输出块大小为 (T - M + 1) x (T - M + 1)。
每个线程块的浮点运算次数：(T - M + 1)² * 2 * M²。
每个线程块的全局内存加载量：T² * 4 bytes（加载整个输入块）。
计算与内存访问比率：[(T - M + 1)² * 2 * M² ops] / [T² * 4 bytes] ≈ 0.5 * M² * (1 - (M-1)/T)² ops/byte。

当 M=5, T=32 时，比率约为 9.57 ops/byte，相比非分块版本有近19倍的提升。因此，分块能显著提高计算与内存访问的比率。

分块卷积的边界处理

在实现分块卷积时，需要小心处理边界。当加载输入块时，有些元素可能位于图像边界之外（称为“幽灵元素”）。有两种处理方法：

在从共享内存读取时进行边界检查。
更简单的方法：在将数据加载到共享内存时，对于越界的输入位置，直接存入0。这样，在后续基于共享内存的计算中，就无需再进行边界检查了。

总结

本节课中，我们一起学习了以下内容：

卷积的基本概念：输出元素是输入邻域的加权和，在图像处理等领域有广泛应用。
卷积的并行化：最简单的策略是为每个输出元素分配一个线程。
常量内存：我们介绍了CUDA中的常量内存，它适用于数据量小、恒定不变且被所有线程频繁访问的情况。我们学习了如何声明常量内存变量 (__constant__)，以及如何使用 cudaMemcpyToSymbol 从主机复制数据。常量内存通过硬件管理的缓存实现快速访问，其优势源于无需写回和缓存一致性协议，并且容量小、命中率高。
卷积内核实现：我们基于常量内存和简单的逐线程方法，勾勒出了卷积内核的代码框架。
数据重用与分块优化：我们分析了卷积中存在的数据重用模式，并提出了使用共享内存进行分块优化的思路。这种优化可以显著减少全局内存访问，提高计算与内存访问的比率。我们还讨论了分块时输入块与输出块大小不匹配的问题，以及边界条件（幽灵元素）的简化处理方法。

通过卷积这个并行模式，我们不仅学习了一种重要的算法模式，也深入了解了常量内存这一GPU硬件特性及其适用场景。在接下来的课程中，我们将继续探索其他并行模式和优化技术。

GPU计算｜CMPS 297S396AA - GPU Computing - Spring 2021：09：模板计算 🧮

在本节课中，我们将要学习一种新的并行模式——模板计算。我们将探讨其基本概念、与卷积运算的异同，并学习如何利用GPU的共享内存和寄存器进行优化。

概述

上一节课我们介绍了卷积运算，它引入了常量内存的使用。在卷积中，每个输出元素都是对应输入元素及其周围元素的线性组合，这个加权和由一个卷积核（或称掩码）定义。我们通过为每个输出元素分配一个线程来并行化计算，并利用共享内存来重用输入数据。

本节中，我们来看看模板计算。这是一种在结构化网格上进行的计算，其中网格点的值基于其邻居点的一个子集来计算。我们将从基础实现开始，逐步引入分块、线程协作和寄存器平铺等优化技术。

什么是模板计算？

模板计算指的是一类在结构化网格上进行的计算模式。网格中某个点的值是基于该点及其邻居点的一个子集计算得出的。

例如，在二维网格中，一个五点模板意味着某个网格点的下一个值是该点本身及其上下左右四个邻居点的组合。我们也可以在三维空间中进行模板计算，例如七点模板，其中点的值基于其在X、Y、Z三个维度上的邻居。

通常，这些网格点被方便地存储为多维数组。因此，我们可以将模板计算视为：输出是一个多维数组，其中每个输出网格点都基于对应的输入网格点及其在各个维度上的邻居。

模板计算与卷积的异同

模板计算实际上是卷积的一个特例。两者非常相似，但存在一些关键区别：

访问模式：卷积通常访问输入元素周围的整个块（例如3x3、5x5的卷积核），而模板计算通常只访问紧邻的邻居（例如上下左右、前后）。
优化机会：由于模板计算只访问特定的邻居子集，这为一些特定的优化提供了可能，而这些优化在通用的卷积中可能无法实现。
维度：本节课我们将重点实现三维模板计算，以练习使用三维网格和线程块，但概念同样适用于二维。

基础并行化策略

与卷积类似，最直接和细粒度的并行化方法是为每个输出元素分配一个线程。

为了简化边界条件的处理，我们将做一个贯穿课程的假设：只计算内部输出值，不计算边界值。这样做可以确保我们访问的输入值始终在边界内，从而简化代码。这个假设通常是合理的，因为在许多实际应用中，边界元素可能由其他GPU或进程负责计算。

对于三维模板计算，最自然的并行化方式是使用三维线程网格，为每个三维输出数组中的元素分配一个线程。

以下是实现的基础步骤：

配置内核：创建三维线程块和三维网格块。

dim3 blockDim(8, 8, 8); // 每个块有512个线程
dim3 gridDim((n + out_tile_dim - 1) / out_tile_dim,
             (n + out_tile_dim - 1) / out_tile_dim,
             (n + out_tile_dim - 1) / out_tile_dim);
stencil_kernel<<<gridDim, blockDim>>>(d_in, d_out, n);

计算线程索引：在内核中，每个线程计算其负责的输出元素在三维网格中的坐标（I, J, K）。

int i = blockIdx.z * blockDim.z + threadIdx.z; // Z维度
int j = blockIdx.y * blockDim.y + threadIdx.y; // Y维度
int k = blockIdx.x * blockDim.x + threadIdx.x; // X维度

边界检查：确保只处理内部元素（非边界）。

if (i >= 1 && i < n-1 && j >= 1 && j < n-1 && k >= 1 && k < n-1) {
    // 进行计算
}

执行计算：计算输出元素的线性索引，并根据模板公式计算其值。对于一个简单的七点模板（当前点 + 六个直接邻居），计算可能如下：

int idx = i * n * n + j * n + k; // 三维到一维的索引映射
d_out[idx] = c0 * d_in[idx] +
             c1 * (d_in[idx - 1] + d_in[idx + 1] +  // X方向邻居
                   d_in[idx - n] + d_in[idx + n] +  // Y方向邻居
                   d_in[idx - n*n] + d_in[idx + n*n]); // Z方向邻居

其中 c0 和 c1 是取决于具体应用的权重常数。

这种基础实现能获得一定的加速比，但仍有优化空间。

利用共享内存进行分块优化

观察发现，处理相邻输出元素的线程会加载许多相同的输入元素，即存在数据复用。与卷积类似，我们可以通过分块技术利用共享内存来优化。

核心思想是：一个线程块负责处理输出数据的一个“块”（Tile）。该块所需的输入数据会被集体加载到该线程块的共享内存中。随后，线程从快速的共享内存中读取数据进行计算，减少对全局内存的访问。

然而，这里存在一个与卷积相同的问题：输入块和输出块的维度不同。输出块的大小是 T x T x T，那么其所需的输入块大小则是 (T+2) x (T+2) x (T+2)（因为每个输出点需要其周围一层的邻居）。

解决方案是：

启动足够多的线程来加载整个输入块（(T+2)^3 个线程）。
在计算阶段，只让其中负责内部输出点的线程（T^3 个）保持活跃并进行计算。
在从全局内存加载数据到共享内存时，需要进行边界检查，防止加载不存在的边界外数据。

以下是关键代码步骤：

定义块尺寸并配置内核。

#define TILE_IN_DIM 10 // 输入块维度
#define TILE_OUT_DIM (TILE_IN_DIM - 2) // 输出块维度
// 线程块维度对应输入块维度
dim3 blockDim(TILE_IN_DIM, TILE_IN_DIM, TILE_IN_DIM);
// 网格维度基于输出块维度计算
dim3 gridDim((n + TILE_OUT_DIM - 1) / TILE_OUT_DIM, ...);

在内核中声明共享内存并加载数据。

__shared__ float tile_s[TILE_IN_DIM][TILE_IN_DIM][TILE_IN_DIM];
// 计算该线程应加载的全局内存坐标（考虑偏移）
int load_i = blockIdx.z * TILE_OUT_DIM + threadIdx.z - 1;
int load_j = blockIdx.y * TILE_OUT_DIM + threadIdx.y - 1;
int load_k = blockIdx.x * TILE_OUT_DIM + threadIdx.x - 1;
// 边界检查后加载到共享内存
if (load_i >= 0 && load_i < n && load_j >= 0 && load_j < n && load_k >= 0 && load_k < n) {
    tile_s[threadIdx.z][threadIdx.y][threadIdx.x] = d_in[load_i*n*n + load_j*n + load_k];
}
__syncthreads(); // 确保所有数据加载完毕

进行计算，仅活跃的内部线程从共享内存读取。

// 检查当前线程是否是负责内部输出点的活跃线程
if (threadIdx.x >= 1 && threadIdx.x < blockDim.x-1 &&
    threadIdx.y >= 1 && threadIdx.y < blockDim.y-1 &&
    threadIdx.z >= 1 && threadIdx.z < blockDim.z-1) {
    // 使用共享内存tile_s进行计算
    float result = c0 * tile_s[threadIdx.z][threadIdx.y][threadIdx.x] +
                   c1 * (tile_s[threadIdx.z][threadIdx.y][threadIdx.x-1] + ... ); // 访问邻居
    // 计算全局索引并写回结果
    int out_i = blockIdx.z * TILE_OUT_DIM + (threadIdx.z - 1);
    ... // 计算 out_j, out_k
    d_out[out_i*n*n + out_j*n + out_k] = result;
}

令人意外的是，这种初始的分块实现可能不会带来性能提升，甚至可能更慢。这是因为分块带来的同步、共享内存读写开销可能超过了数据复用减少的全局内存访问收益。性能提升与否取决于计算与内存访问的比率。

通过增大分块尺寸和线程协作进行优化

分析表明，分块优化的收益与输出块的尺寸 T 有关。T 越大，内部元素（有数据复用）与边界元素（无复用或少复用）的比例就越高，计算/内存访问比率就越好，潜在性能收益越大。

然而，直接增大 T 会遇到两个硬件限制：

线程数限制：线程块大小（(T+2)^3）不能超过硬件允许的最大值（通常是1024）。
共享内存限制：存储整个输入块所需的共享内存随 T^3 增长，可能耗尽资源或降低占用率（Occupancy）。

解决方案是结合线程协作（Thread Coarsening）。我们不再让一个线程只处理一个输出点，而是让一个线程处理多个输出点（例如，在Z维度上连续处理多个平面）。这样，我们可以在不增加线程数量的情况下，有效增大逻辑上的输出块尺寸。

具体策略是：

分配足够的线程来处理输出块中的一个二维平面（例如32x32个线程）。
这些线程协作处理整个三维输出块。在每次迭代中，它们处理一个Z深度的平面，并仅将当前计算所需的输入平面（前一个、当前、后一个）保存在内存中。
随着迭代进行，线程集体在Z维度上移动，处理下一个平面。

引入寄存器平铺以优化存储

进一步观察发现，在三维模板计算中，对于当前正在处理的平面（当前切片），其数据被该平面内的多个线程共享（例如，一个点被其左右邻居线程访问）。因此，当前切片适合放在共享内存中。

但是，对于“前一个”和“后一个”切片，每个点基本上只被一个对应的线程使用（用于计算其正上方和正下方的邻居贡献）。因此，这些数据不需要在线程间共享。

基于此，我们可以进行一项关键优化：寄存器平铺（Register Tiling）。

仅将需要共享的当前切片存储在共享内存中。
将不需要共享的前一个切片和后一个切片存储在寄存器中。每个线程将自己需要的那个点保存在其私有寄存器里。
在迭代过程中，数据在寄存器和共享内存之间移动：下一切片的数据先从全局内存加载到寄存器，当它变成“当前切片”时，再被存入共享内存；当“当前切片”变成“前一切片”时，其数据可以从共享内存移回寄存器（或通过指针交换逻辑避免实际拷贝）。

这种技术充分利用了GPU上通常比共享内存容量更大的寄存器文件，减少了对宝贵共享内存的占用，从而可能提高占用率和整体性能。

实际上，我们在矩阵乘法优化中已经无形使用了寄存器平铺：每个线程将累加结果（C矩阵的一个元素）存储在自己的寄存器中，整个C矩阵的块（Tile）实际上分布式存储在所有线程的寄存器集合里。

总结

本节课我们一起学习了模板计算这一重要的并行模式。

我们首先介绍了模板计算的基本概念，并将其与之前学过的卷积进行了对比。接着，我们实现了基础的三维模板计算内核。

然后，我们探讨了如何利用共享内存进行分块优化，并分析了其性能瓶颈。为了突破线程数和共享内存的限制，我们引入了线程协作的概念，以逻辑上增大分块尺寸。

最后，我们学习了高级优化技术——寄存器平铺，通过分析数据在模板计算中的共享特性，将不需要共享的数据存储在寄存器中，从而更高效地利用GPU的存储层次结构。

这些优化策略（共享内存分块、线程协作、寄存器平铺）是高性能GPU编程的核心技术，能够显著提升诸如模板计算这类具有数据复用模式的应用性能。

posted @ 2026-03-26 12:21 布客飞龙III 阅读(63) 评论(0) 收藏举报

刷新页面返回顶部

龙哥盟

贝鲁特-GPU-编程笔记-全-

贝鲁特 GPU 编程笔记（全）

GPU计算：10：归约（Reduction） 🧮

概述

回顾：模板计算（Stencil）

什么是归约？ 🤔

串行归约

并行归约：归约树 🌳

GPU上的归约挑战与策略

基础归约树实现（低效版本）

优化归约树实现（高效版本）

进一步优化：使用共享内存

高级优化：线程粗化（Thread Coarsening） ⚙️

边界条件处理

总结

GPU计算：第11讲：扫描（Kogge-Stone）🔍

概述

回顾：归约（Reduction）

扫描（Scan）操作介绍

分段扫描（Segmented Scan）

Kogge-Stone 并行扫描算法

优化：双缓冲（Double Buffering）

实现排除式扫描

工作效率（Work Efficiency）分析

总结

GPU计算：第12讲：扫描算法（Brent-Kung）🔬

概述

回顾：扫描模式与Kogge-Stone算法

Brent-Kung并行扫描算法

1. 归约阶段（Reduction Step）

2. 后归约阶段（Post-Reduction Step）

算法比较与分析

Brent-Kung算法实现要点

共享内存的使用

无需双缓冲

控制流分歧与线程重索引

从包含性扫描到排他性扫描

线程粗化优化

总结

GPU计算：13：直方图与原子操作

回顾：扫描操作

什么是直方图？

并行直方图的挑战

原子操作：解决数据竞争

优化一：私有化

优化二：线程粗化

总结

GPU计算：第14讲：合并模式 🧩

概述

顺序合并算法

并行合并的挑战与思路

寻找协同秩

并行合并内核实现

性能优化：使用共享内存

总结

GPU计算｜CMPS 297S396AA：15：排序算法

概述

基数排序 🧮

单比特基数排序演示

并行化基数排序的关键步骤

优化：改善合并访问

基数选择与线程聚合

归并排序 🔀

并行化归并排序

总结

GPU计算：16：稀疏矩阵计算（COO与CSR格式）

什么是稀疏矩阵？

稀疏矩阵存储格式概览

坐标格式

COO格式的原理

基于COO格式的SPMV实现

COO格式的优缺点分析

压缩稀疏行格式

CSR格式的原理

基于CSR格式的SPMV实现

CSR格式的优缺点分析

性能对比与总结

GPU计算：第17讲：稀疏矩阵计算（ELL与JDS格式）🚀

稀疏矩阵计算回顾