CUDA C++ Best Practices Guide Notes 1

原文地址：https://docs.nvidia.com/cuda/cuda-c-best-practices-guide/

1. 序言

1.2. Access（评估）、Parallelize（并行化）、Optimize（优化）、Deploy（部署）

Access、Parallelize、Optimize、Deploy（APOD）是一个循环过程，即优化一个 CUDA 程序是一个持续迭代的过程。即初始投入最小化、测试与部署、循环迭代和持续改进，最终实现程序的最大性能。

1.2.1. Access

对于现有项目，优化的第一步就是评估，以定位哪一部分代码消耗了最多的执行时间。有了这些信息之后开发者就可以评估这些瓶颈是否适合并行化，以及如何通过 GPU 加速来优化。

优化的过程中需要理解用户的需求和限制，然后应用阿姆达尔定律（Amdahl’s Law，CSAPP 里面提到过）和古斯塔夫森定律（Gustafson’s Law）来评估并行化的潜力。这些定律可以帮助开发者确定哪些部分的代码可以并行化，以及并行化的潜力有多大。

古斯塔夫森定律评估程序通过并行计算获得的加速 $S$，公式如下：
$S = s + P \times N$
$= s + (1 - s) \times N$
$= N + (1 - N) \times s$
其中，$S$ 表示具有并行性的程序的理论加速；
$N$ 是处理器数量；
$s$ 和 $P$ 是分别程序串行部分执行时间和并行部分执行时间的占总执行时间的比例，在并行系统中，$s + P = 1$。
所以，$S$ 也可以使用 $P$ 来表示：
$S = (1 - P) + P \times N$
$= 1 + (N - 1) \times P$

1.2.2 Parallelize

识别出热点代码和确定了优化目标之后，开发者就需要并行化代码了。根据原始代码的情况，这可能简单到调用现有的 GPU 优化库（如 cuBLAS、cuFFT 或 Thrust），也可能简单到添加一些预处理器指令作为并行化编译的提示。

另一方面，一些程序的设计需要做大量的重构工作才能把它内在的并行性暴露出来。正如即使是 CPU 架构也需要暴露并行性来提升或仅仅维持顺序应用程序的性能一样，CUDA 系列并行编程语言（如 CUDA C++、CUDA Fortran 等）的目标是尽可能简化并行性的表达，同时支持在专为最大并行吞吐量设计的 CUDA GPU 上运行。

1.2.3. Optimize

在对程序的每一轮并行化完成之后，开发者就可以把工作重点移到优化上来提升性能了。由于有很多种优化方法，因此越理解应用的需求，APOD 的循环就越顺利。然而，正如 APOD 整体一样，程序优化也是一个迭代的过程。这意味着程序员不需要在获得良好的加速效果之前花费大量时间记忆所有的优化策略。相反，策略可以在学习过程中逐步应用。

优化可以在多个层次上进行，从重叠数据传输与计算到微调浮点运算序列。可能的性能分析工具在这一过程中非常宝贵，因为它们可以帮助开发者确定哪些优化策略是最有效的。

1.2.4. Deploy

在对程序的一部分或者更多的部分完成 GPU 加速之后，就可以将结果与最初的预期进行比较了。回想一下，初始的评估过程确定了优化热点代码后的潜在的加速上限。

在解决其它热点代码以进一步提高总加速效果之前，开发者应该考虑将部分并行化的实现部署到生产环境中。

1.3 推荐和最佳实践

1.4 评估你的程序

套话。

2. 异构计算

CUDA 编程涉及在两个不同的平台上同时运行代码：一个拥有一个或多个 CPU 的 Host 系统和一个拥有一个或多个支持 CUDA 的 NVIDIA GPU 设备。

2.1. *Host 和 Device 之间的不同

线程资源

Host 系统上的执行流水线可以支持有限数量的并发执行线程。比如一个有 32 核的 CPU 可以并发跑 64 个线程。作为比较，现代的 NVIDIA GPU 的每个 Multiprocessor 最多可以同时支持 2048 个活动线程。这就导致，在拥有 80 个 Multiprocessor 的 GPU 上，可以同时支持超过 160,000 个活动线程。

线程

CPU 上的线程一般是重量级的实体（相较于 GPU 线程）。操作系统必须把线程换入换出 CPU 的执行通道（Execution Channal）来提供多线程的支持。因此这种换入换出带来的上下文切换（Context Switches）是缓慢且昂贵的。作为比较，GPU 的上线程则是极度轻量的。在一个典型的 GPU 系统上，上千个线程排队等待工作（以每组 32 个线程的 Warp 为例）。如果 GPU 必须等待一个 Warp 的线程，它会立即开始执行另一个 Warp 的工作。由于所有的活动线程都分配了独立的寄存器，因此在执行 GPU 线程时，不需要进行寄存器或其它状态的交换。资源会一直分配给每个线程，直到其执行完成。简而言之，CPU 核心的设计目标是最小化少量线程的延迟，而 GPU 的设计目标是处理大量并发的轻量级线程，以最大化吞吐量（Throughput）。

Warp 这个词源自纺织业，指的是织布机上的经纱（Warp Threads）。在织布过程中，经纱是纵向排列的纱线，它们被固定在织布机上，而纬纱（Weft Threads）则横向穿过经纱，形成布料。经纱是织布的基础，而它们以统一的方式被拉动和操作，这与 GPU 中线程的行为非常相似，这就是为什么 Warp 这个词被用来描述 GPU 中的线程组。

RAM

Host 和 Device 各自拥有独立的物理内存。由于 Host 内存和 Device 内存是分离的，Host 内存中的数据必须偶尔在 Device 内存和 Host 内存之间进行传输。

这些就是涉及并行编程时 CPU 和 GPU 之间的主要硬件区别。考虑到这些差异而组成的应用程序可以将 Host 和 Device 视为一个有凝聚力的异构系统，其中每个处理单元都将被用来做它最擅长的工作：Host 的串行代码和 Device 的并行代码。

2.2. 在启用 CUDA 的设备上运行什么？哪些任务适合跑在 CUDA 设备上？

在决定这个问题时，需要考虑一下问题：

CUDA 设备是非常适合处理大量数据的并行计算任务的。这通常设计对大型数据集（如矩阵）的算术计算，其中可以同时对数千个元素执行相同的操作。并行运行大量线程的支持在上述的 GPU 轻量级线程模型中已经介绍过了。
为了使用 CUDA，数据必须从 Host 传输到 Device。这些传输性能方面代价高昂，应尽量减少，这个代价有几个后果：
1. *传输数据到 Device 和从 Device 传输数据的操作的复杂性应该被证明是合理的。传输数据给少量的线程使用的代码（相比于传输大量数据的代码）几乎没有任何优势。理想的情况是许多线程同时执行大量的工作。
  1. 举一个矩阵相加的例子。假设要在 Device 上执行两个 $N \times N$ 矩阵的加法，那么相加需要 $N^2$ 次计算操作，传输数据有 $3N^2$ 次移动。因此，计算和传输的比值是 $1:3$ 或者 $O(1)$。当这个比值较高时，可以更容易地取得性能收益。比如对于还是同样的矩阵，矩阵乘法需要 $N^3$ 次计算操作，因此，计算和传输的比值是 $N:3$，也就是 $O(N)$，在这种情况下，矩阵越大，收益越大。当然操作的类型是一个额外的因素，因为加法的复杂性和三角函数不同。总之，在确定究竟在 Host 还是 Device 上执行操作时，数据的传输开销是一个重要的考虑因素。
2. 如果数据要在 Device 上计算，那么数据应该就被尽可能长地保留在 Device 上。这个很好理解，尽量建设数据在 Host 和 Device 之间的传输。
为了取得最好的性能，Device 上相邻的线程访问的内存数据应该时连续的。一定的内存访问模式可以让硬件将多个数据项的读取或写入合并到一个操作中。无法因布局以实现操作合并的数据，或者没有足够的局部性来有效使用 L1 或者纹理缓存的数据，在 GPU 上的计算中往往只有较小的加速。
1. 但有一个完全例外的模式值得注意，那就是完全随机的内存访问模式。一般来说，任何架构都只能以较低的效率处理这些内存访问模式。然而，于基于缓存的架构（Cache Based Architecture，如 CPU）相比，延迟隐藏架构（Latency Hiding Architecture，如 GPU）往往能更好地应对完全随即的内存访问模式。
  1. 注：延迟隐藏架构是一种通过并行性和高效调度来隐藏内存访问延迟的设计。例如当一个 Warp 的线程在等待内存访问时，GPU 会立即切换到另一个 Warp 执行。

3. 应用程序性能分析

3.1. Profile

许多代码都是一小段代码的执行会占据实际工作的大部分，因此这部分代码需要被好好分析。

3.1.1. Creating the Profile

可以使用开源的工具，如gprof生成性能报告。注意，需要在使用gcc编译时加上-pg选项。

$ gcc -O2 -g -pg myprog.c
$ gprof ./a.out > profile.txt
Each sample counts as 0.01 seconds.
  %   cumulative   self              self     total
 time   seconds   seconds    calls  ms/call  ms/call  name
 33.34      0.02     0.02     7208     0.00     0.00  genTimeStep
 16.67      0.03     0.01      240     0.04     0.12  calcStats
 16.67      0.04     0.01        8     1.25     1.25  calcSummaryData
 16.67      0.05     0.01        7     1.43     1.43  write
 16.67      0.06     0.01                             mcount
  0.00      0.06     0.00      236     0.00     0.00  tzset
  0.00      0.06     0.00      192     0.00     0.00  tolower
  0.00      0.06     0.00       47     0.00     0.00  strlen
  0.00      0.06     0.00       45     0.00     0.00  strchr
  0.00      0.06     0.00        1     0.00    50.00  main
  0.00      0.06     0.00        1     0.00     0.00  memcpy
  0.00      0.06     0.00        1     0.00    10.11  print
  0.00      0.06     0.00        1     0.00     0.00  profil
  0.00      0.06     0.00        1     0.00    50.00  report

3.1.2. 分析热点

在上述例子中，我们可以看到函数genTimeStep占据了程序运行总时间的$\frac{1}{3}$。

3.1.3. *理解 Scaling

在 CUDA 上能取得的性能收益完全取决于并行化的程度。不能充分并行化的代码应该运行在 Host，除非这么做会导致过度的数据传输。

通过理解程序如何 scale（拓展?），可以设定期望并规划并行化策略。

3.1.3.1. 强拓展性和阿姆达尔定律

强扩展性是衡量在固定整体问题规模的情况下，随着处理器数量的增加，求解时间如何减少的指标。如果一个应用程序表现出线性强扩展性，则其加速比等于所使用的处理器数量。

强拓展性通常等同于阿姆达尔定律，该定律规定了并行化串行程序部分所能预期的最大加速。本质上，它指出程序的最大加速比为：

$$ S = \frac{1}{(1 - P) + \frac{P}{N}} $$

这里 $P$ 是可以并行化的代码部分占总串行执行时间的比例，而 $N$ 是运行代码并行部分的处理器数量。

$N$越大，$\frac{P}{N}$越小。可以把$N$看作一个很大的数，这个这个等式就可以简化为$S = \frac{1}{1 - P}$。现在，如果$\frac{3}{4}$的代码可以并行化，那么最大加速比就是$S = \frac{1}{1 - \frac{3}{4}} = 4$。

3.1.3.2. 弱拓展性和古斯塔夫森定律

弱扩展性是衡量在每个处理器的问题规模固定的情况下，随着处理器数量的增加，求解时间如何变化的指标；也就是说，随着处理器数量的增加，整体问题规模也会增加。

弱拓展性通常等于古斯塔夫森定律，该定律指出，在实践中，问题规模会随处理器数量的增加而拓展。因此，程序的最大加速比为：

$$ S = N + (1-P)(1-N) \ = 1 - P + N \times P $$

其中，$P$ 是可以并行化的代码部分占总串行执行时间的比例，$N$ 处理器数量。

看待古斯塔夫森定律的另一种方式是，当我们扩大系统时，不是问题规模保持固定，而是执行时间保持固定。注意，古斯塔夫森定律假设串行与并行代码执行时间的比值时固定的，这反映了设置和处理更大问题的额外成本。

公式表明，随着$N$的增加，加速比趋近于$N \times P$，即并行部分主导性能提升。

3.1.3.3. 应用强弱拓展性

4. 并行化应用程序

确认热点代码并且分析代码确定优化目标之后，就可以并行化代码了。根据原始代码的情况，这可能简单到调用现有的 GPU 优化库（如 cuBLAS、cuFFT 或 Thrust），也可能简单到添加一些预处理器指令作为并行化编译的提示。

5. 开始并行化

5.1. 并行库

最直接的方式是利用线程的并行库，CUDA 工具包有很多为 CUDA 设备优化好的并行库，比如cuBLAS、cuFFT等。

5.2. 并行化编译器

5.3. 并行化编程模型

6. 获取正确的答案

并行系统中可能会遇到传统编程从不会错遇到的问题，例如线程问题、浮点值计算方式导致的意外值等。

6.1. 验证

6.1.1. 参考比较

对任何现有程序的修改，正确性验证的一个关键方面就是建立某种机制，使得可以从代表性输入中获得的已知良好的参考输出与新结果进行比较。

人话，建立数据集，检测修改是否会导致回退。

6.1.2. 单元测试

还有一种方法是将代码本身结构化，使其在单元测试中易于验证。例如，可以将 CUDA kernel 函数编写成许多短小的__Device__函数，而不是一个庞大的__global__函数。

更进一步，如果大多数函数被定义为__Host__ __Device__，而不仅仅是__Device__，那么这些函数可以在 CPU 和 GPU 上同时进行测试，从而增强对函数正确性的信心。

如果__Host__ __Device__是在 Host 代码（普通的函数）中调用，那它就执行在 CPU 上，如果被在 Device 代码（kernel 函数，例如global）中调用，那它就执行在 GPU 上。

6.2. 调试

CUDA-GDB。

6.3. 数值准确性与精度

7. 优化 CUDA 程序

在每一轮应用程序并行化完成后，开发者可以转向优化实现以提高性能。

8. *性能指标

8.1. Timing

8.1.1. 使用 CPU 计时器

CUDA API 函数都是异步的，因此控制流会在 CUDA API 完成它们的工作之前返回到 Host 线程，因此在计时之前，需要使用cudaDeviceSynchronize()将 CPU 线程和 GPU 进行同步。cudaDeviceSynchronize()会阻塞 CPU 线程，直到所有在调用之前的 CUDA API 函数都完成。

当然也可以将 CPU 线程与 GPU 上的特定 Stream 或者事件进行同步，但这些同步函数并不适合用于计时默认流以外的流中的代码（？）。cudaStreamSynchronize() 会阻塞 CPU 线程，直到之前发布到给定流中的所有 CUDA 调用完成。cudaEventSynchronize() 会阻塞，直到 GPU 记录了特定流中的给定事件。

在接触到的项目中，使用的几乎都是cudaStreamSynchronize()，Stream 几乎是伴随着数据流的。

由于 Default Stream 在设备上表现出序列化的行为（Default Stream 中的操作只能在先前所有 Stream 的调用完成后开始；并且直到 Default Stream 中的操作完成之后，其它 Stream 中的后续操作才能开始，也就是独占的），因此可以将需要计时的代码放在 Default Stream 中，然后使用cudaStreamSynchronize()和cudaEventSynchronize()以实现精确计时。而其它 Stream 由于并行性，计时可能会收到其它 Stream 的影响。

为什么我自己写的 sample code，使用多个 custom Stream 然后统计其中一个 stream，得到的耗时，比使用 Default Stream 得到的耗时要短？

8.1.2. 使用 CUDA GPU 计时器

GPU 计时器的值以毫秒表示，精度约为 0.5 微秒。该值是在 GPU 时钟上测量的，精度和操作系统无关。

Sample code

8.2. 带宽

带宽（即数据传输速率）是性能最重要的限制因素之一。几乎所有代码的修改都应在考虑它们如何影响带宽的背景下进行。正如本指南的**内存优化**部分所述，带宽会受到数据存储内存的选择、数据的布局方式、访问顺序以及其他因素的显著影响。

为了准确测量性能，计算理论带宽和有效带宽是非常有用的。当后者远低于前者时，设计或实现细节可能会降低带宽，而提高带宽应成为后续优化工作的主要目标。

8.2.1. 理论带宽计算

理论带宽可以通过硬件规格表中的内存频率和内存总线宽度来计算。例如 NVIDIA Tesla V100 使用的内存为 HBM2（双倍数据速率），其频率为 877 MHz（$877\times 10^6$），内存总线宽度为 4096 位。使用这些数据可以计算得到 NVIDIA Tesla V100 的理论带宽为 898 GB/s：

$$ (0.877 \times 10^9 \times (4096 / 8) \times 2) \div 10^9 = 898 \text{ GB/s} $$

由于是 HBM2，因此结果乘 2，将结果除以 $10^9$，以将单位转化为 GB。

8.2.2. 有效带宽计算

有效带宽是通过测量实际应用程序的数据传输速率来计算的。使用这个公式：

$$ Effective\ Bandwidth = \frac{(B_r + B_w)\div 10^9}{time} $$

有效带宽的单位是 GB/s，$B_r$ 和 $B_w$ 分别是读取和写入的字节数，time 是数据传输的时间。

例如，计算一个$2048 \times 2048$的矩阵复制的有效带宽：

$$ Effective\ Bandwidth = \frac{(2048^2 \times 4 \times 2) \div 10^9}{time} $$

其中，$2048^2$ 是矩阵的大小，乘 4 是每个元素（float）的字节数，乘 2 是读取和写入，time 是数据传输的时间，除以 $10^9$ 是为了将单位转化为 GB。

8.2.3. Visual Profiler 报告的吞吐量

对于 Compute capability 为 2.0 或更高的设备，Visual Profiler 可用于收集多种不同的内存吞吐量测量值。