HOW CUDA PROGRAMMING WORKS

STEPHEN JONES, GTC 2022

目录

• 为什么CUDA会是现在这个样子
• 使用GPU是为了性能
• NVIDIA Ampere GPU 架构
• FLOPS不是问题所在 - 带宽才是

• 深入探讨随机存取存储器（RAM）

  • DRAM工作原理
  • DRAM读取过程
  • 物理原理介入与HBM性能

• 数据访问模式的重要性

• CUDA的GPU执行层级

  • CUDA 线程块
  • SIMT模型：每个线程运行相同的程序
  • 线程块的组织与数据访问
  • GPU 上的 Warp 执行
  • 内存访问模式：并行与吞吐量

• 充分利用GPU资源

  • Ampere A100 GPU 关键规格与架构
  • CUDA的GPU执行层次结构

• 流式多处理器内部

  • SM“满载”的含义
  • A100 SM 资源概览
  • 线程块的结构与资源需求
  • GPU 如何将块放置在 SM 上

• 占用率

  • 占用率概念
  • 占用率量化
  • 占用率是程序调优最强大的工具

• 填补空白：优化资源利用

  • 资源需求 (蓝色网格
  • 资源需求 (绿色网格
  • 结合蓝色和绿色网格填充资源

• 并发性

  • 并发：同时做多件事
  • 并发：依赖性
  • 并发：多个依赖性任务
  • 并发：独立性
  • 并发：本质是超额订阅
  • 并发：超额订阅的重要性

• 总结与CUDA设计理念

为什么CUDA会是现在这个样子？ (Page 3-4)

本演讲探讨了CUDA编程为何会呈现出当前的形式。核心原因在于：

• 物理定律 (Physics)：GPU存在的根本原因是为了性能。使用CUDA在GPU上编写高性能代码，但性能始终受到物理定律的限制。因此，CUDA的设计方式是受物理定律驱动的。

使用GPU是为了性能 (Page 8-9)

• 使用GPU是为了追求性能。
• 这意味着要尽可能利用所有可用的GPU资源。

NVIDIA Ampere GPU 架构 (Page 10-12)

Ampere GPU 架构提供了以下资源：

参数	数值
SMs	108
总线程数	221,184
峰值FP32 TFLOP/s	19.5
峰值FP64 TFLOP/s (非张量)	9.7
峰值FP64 TFLOP/s (张量)	19.5
张量核心精度	FP64, TF32, BF16, FP16, I8, I4, B1
每SM共享内存	160 kB
L2 缓存大小	40960 kB
内存带宽	1555 GB/sec
GPU Boost 时钟	1410 MHz
NV/Link 互联	600 GB/sec

值得注意的是，Ampere GPU 的峰值 FP64 TFLOP/s（非张量）达到 9.7，这一性能甚至超越了2001年排名Top500第一的ASCI White超级计算机（其性能为7.9 teraflop/s）。

FLOPS不是问题所在 - 带宽才是 (Page 13-15)

尽管GPU拥有强大的FLOPs能力，但实际性能瓶颈往往在于带宽。

• Ampere A100 GPU 架构分析：
  • GPU中包含108个SM，运行在1410 MHz的提升时钟频率下。
  • 每个SM可以在每个时钟周期加载64字节数据。
  • 峰值内存请求速率 = 64B x 108个SMs * 1410 MHz = 9750 Gigabytes/sec。
  • 而HBM2内存带宽仅为1555 GBytes/sec。
  • 这意味着请求的带宽与提供的带宽之比为 9750 / 1555 ≈ 6.3倍。

• 基于内存速度的FP64 FLOP/s计算：
  • 如果仅考虑内存带宽作为限制，FP64 FLOP/s = 1555 GB/SEC / 8 BYTES = 194 GFLOP/S。
  • 这与1996年排名Top500第一的Hitachi SR2201超级计算机（性能为220 gigaflop/s）处于同一数量级。

这表明，即使GPU具备极高的浮点运算能力，若内存带宽不足，实际的有效计算能力仍会受到严重限制，甚至可能退化到与多年前的超级计算机相当的水平。

深入探讨随机存取存储器（RAM）

DRAM工作原理

幻灯片深入探讨了随机存取存储器（RAM）的工作原理，展示了1比特和0比特在DRAM单元中（由一个晶体管和一个电容器组成）的存储方式。当电容器充电时存储“1”比特，放电时存储“0”比特。

此幻灯片展示了DRAM存储器阵列的整体架构，包括行解码器（Row Decoder）、列解码器（Column Decoder）、感应放大器（Sense Amplifiers）和数据缓冲区（Data Buffers）。读取地址被清晰地分解为行地址和列地址部分，分别由行解码器和列解码器处理。

DRAM读取过程

DRAM读取过程的第一步是：

1. 激活行并将数据拉入感应放大器。这个过程会因为电容器放电而破坏该行中的数据。
   

DRAM读取过程的第二步是：
2. 从感应放大器中保存的页面（page）中读取指定列索引的数据。此操作不会破坏感应放大器中的数据。

DRAM读取过程的第三步是：
3. 可以从感应放大器中保存的同一个页面（page）重复读取不同列索引的数据。

第三步的进一步阐释：
3. 可以从感应放大器中保存的同一个页面重复读取不同列索引的数据。“突发”读取（"Burst" reads）可以一次加载多个列。

DRAM读取过程的第四步是：
4. 在获取新页面之前，必须将旧行数据写回，因为之前的读取操作已经破坏了原始数据。

物理原理介入与HBM性能

幻灯片解释了DRAM访问中的物理原理。
示例HBM值：
- 读取新列的时间：C_L = 16 周期
- 加载新页面的时间：T_RCD = 16 周期
- 写回数据的时间：T_RP = 16 周期
- 页面（行）大小 = 1KB

每次读取的成本 (C_L = 16 周期)。
切换页面的成本是其3倍 (T_RP + T_RCD + C_L = 48 周期)。
这是因为切换页面需要对电容器进行充电/放电，这涉及到物理RC时间常数。
电容器电压公式： $V_C = V_S(1 - e^{-t/RC})$

这意味着什么？我们预期对于地址合并（coalesced reads）与分散读取（scattered reads）会存在显著的性能差异。
在A100上，对于广泛间隔读取的内存带宽为：111 / 1418 = 峰值带宽的8%。这相当于峰值带宽的1/13！

图表：HBM内存吞吐量随地址发散（8字节读取，A100）
- 峰值带宽：1418 GB/sec
- 突发大小：64 Bytes
- HBM页面大小：1KB
- 在步长间隔为8192字节时，性能下降到111 GB/sec。

此幻灯片内容与上一页相同，但图中叠加了一张CM-5超级计算机的图片，CM-5是1993年洛斯阿拉莫斯国家实验室首次达到第一的超级计算机，性能为59.7 gigaflop/s。

数据访问模式的重要性

此幻灯片强调了数据访问模式的重要性。
伪代码展示了行主序数组遍历：

for(y=0; y<M; y++) {
    for(x=0; x<N; x++) {
        load(array[y][x]);
    }
}

右侧网格图展示了行主序数组布局，并标示了列读取延迟 C_L 的方向。

此幻灯片对比了行主序和列主序数组遍历对性能的影响。
列主序数组遍历的伪代码：

for(x=0; x<N; x++) {
    for(y=0; y<M; y++) {
        load(array[y][x]);
    }
}

指出“行读取延迟 T_RAS = T_RP + T_RCD + C_L”比列访问慢13倍，再次强调了数据访问模式对性能的巨大影响。

使用GPU的性能原因在于：充分利用所有GPU资源，这意味着需要管理内存访问模式。

但这与CUDA有什么关系？

CUDA的GPU执行层级

此幻灯片介绍了CUDA的GPU执行层级结构：
- 工作网格（Grid of work）
- 被划分为多个块（Divide into many blocks）
- 每个块包含多个线程（Many threads in each block）

CUDA 线程块 (THE CUDA THREAD BLOCK)

一个线程块（Thread block）包含固定数量的线程，这些线程保证在同一个SM（流多处理器）上同时运行。

SIMT模型：每个线程运行相同的程序 (EVERY THREAD RUNS EXACTLY THE SAME PROGRAM)

这就是“SIMT”模型（单指令多线程）。一个线程块拥有固定数量的线程。

CUDA 内核函数 euclidean_distance 示例：

__global__ void euclidean_distance(float2 *p1, float2 *p2, float *distance, int count) {
    // Calculate the index of the point my thread is working on
    int index = threadIdx.x + (blockIdx.x * blockDim.x);

    // Check if thread is in-range before reading data
    if (index < count) {
        // Compute the Euclidean distance between two points
        float2 dp = p2[index] - p1[index];
        float dist = sqrtf(dp.x * dp.x + dp.y * dp.y);
        // Write out the computed distance
        distance[index] = dist;
    }
}

其中，int index = threadIdx.x + (blockIdx.x * blockDim.x); 这一行代码是计算每个线程处理的数据索引的关键。

SIMT 的核心在于条件检查 if (index < count)，它确保了在所有线程执行相同程序时，只有在数据范围内的线程才进行实际计算，从而实现并行控制。

线程块的组织与数据访问 (Thread Block Organization and Data Access)

线程块被分解为 32 个线程的“warp”（线程束）。一个“warp”是 GPU 的矢量元素。
例如，一个线程块可以包含多个 warp，每个 warp 负责处理一定范围的线程 ID：
* warp 1: thread 0...31
* warp 2: thread 32...63
* warp 3: thread 64...95
* warp 4: thread 96...127
* warp 5: thread 128...159
* warp 6: thread 160...191

在代码层面，数据访问通常通过一个计算出的索引完成，例如 float2 dp = p2[index] - p1[index];，其中 index 是基于线程 ID 计算的。

float2 数据类型表示一个包含两个浮点数（float x, float y）的结构，其大小为 sizeof(float2) = 2x 4 bytes = 8 bytes。

因此，一个 warp（32 个线程）在进行数据加载时，会加载 32 x 8 bytes 的数据，总计 256 bytes。

GPU 上的 Warp 执行 (WARP EXECUTION ON THE GPU)

多个 warp 被调度到 GPU 的流多处理器 (SM) 上执行。
下图展示了 warps 如何映射到 A100 SM 块图上。

从一个 SM 发出的内存请求：例如，如果 4 个 warps 同时发出请求，那么总共的内存请求量为 4 warps x 256 bytes per warp = 1024 bytes。

这种并行内存请求与内存硬件设计紧密相关。内存页大小通常为 1024 bytes。这意味着来自一个 SM 的 4 个 warps 的内存请求（总计 1024 字节）可以一次性地加载一个完整的内存页，从而实现高效的内存访问。

内存访问模式：并行与吞吐量 (Memory Access Patterns: Parallelism and Throughput)

对于单个线程，其内存访问可能看起来像随机地址内存读取。

然而，由于 GPU 上是 128 个（例如 4 个 warp x 32 线程/warp）线程并行执行，这实际上变成了对整个内存页的相邻读取。这种“合并访问”（coalesced access）是 GPU 实现高内存吞吐量的关键。

上图展示了 HBM 内存吞吐量与地址发散（连续读取之间的步幅间隔）的关系（8 字节读取，A100）。
* 当步幅间隔（Stride Interval）较小（例如 8 字节，即连续访问）时，Achieved Bandwidth（实现带宽）最高，接近 2048 GB/sec。
* 随着步幅间隔的增加，实现带宽会急剧下降。当步幅间隔达到 64 字节或更大时，带宽会显著降低并趋于稳定在较低水平，例如 128 GB/sec 左右。
这表明，为了获得最佳性能，GPU 程序应尽量采用连续或小步幅的内存访问模式，以最大化 HBM 内存吞吐量。

充分利用GPU资源

为了获得最佳性能，需要充分利用所有可用的GPU资源。这意味着要有效地管理内存访问模式，并考虑“占用率”（occupancy）这一概念。

Ampere A100 GPU 关键规格与架构

Ampere A100 GPU 的详细规格已在“NVIDIA Ampere GPU 架构”章节中列出。其架构包括多个SM（例如，从SM0到SM107），每个SM包含寄存器（256K）和L1/共享内存（192K）。所有SM都连接到一个40MB的L2缓存，L2缓存再连接到80GB的HBM显存。

CUDA的GPU执行层次结构

CUDA的GPU执行层次结构定义了工作如何被分解和分配到GPU上：

1. 起始工作：首先，有一项需要处理的工作，例如一张图片。
   

2. 划分为等大小的块（Grid of Work）：这项工作被划分为一系列等大小的块，这构成了“工作网格”（Grid of Work）。
   

3. 独立处理每个块：每个块都将独立进行处理。CUDA不保证执行顺序，并且块之间不能直接交换数据。
   

4. 块被放置到SM上：每个块都会被放置到一个流处理器（SM）上。
   

5. 持续放置块直至每个SM“满载”：当一个块完成其工作并退出后，一个新的块会被放置到其位置，这个过程持续进行，直到整个网格的工作完成。
   
   
   

流式多处理器内部 (LOOKING INSIDE A STREAMING MULTIPROCESSOR)

SM“满载”的含义 (WHAT DOES IT MEAN FOR AN SM TO BE “FULL”?)

SM的“满载”状态指的是其内部资源被充分利用。一张幻灯片展示了多个流式多处理器（SM 0到SM 11）的概览，其中一些SM显示为部分或完全被占用（以绿色填充程度表示），而另一些则为空白或处于未激活状态。一个SM的内部结构由四个处理单元组成，每个单元包含L1指令缓存、L0数据缓存、调度器、寄存器文件（F0，256 x 32b）和Tensor Core。SM的“满载”状态与这些内部资源的利用率有关。

A100 SM 资源概览

以下是A100流式多处理器（SM）的关键资源列表：
* 每个SM的最大线程数：2048
* 每个SM的最大块数：32
* 每个SM的总寄存器数：65536
* 每个SM的总共享内存：160 kB
* 每个warp的线程数：32
* 并发活动的warp数：4
* 每个SM的FP32核心数：64
* 每个SM的FP64核心数：32
* 最大L1缓存大小：192 kB
* 每个SM的加载带宽：90 GB/秒
* GPU Boost Clock：1410 MHz

幻灯片也强调了“每个SM的最大线程数”、“每个SM的总寄存器数”、“每个SM的总共享内存”和“每个warp的线程数”等资源。

线程块的结构与资源需求 (ANATOMY OF A THREAD BLOCK)

GPU上的工作负载（如对图像进行处理）被组织成一个“工作网格”（Grid of work），该网格进一步划分为许多“块”（Blocks）。每个块又包含许多“线程”（Threads）。每个块中的线程数量是固定的，并且每个线程都运行完全相同的程序，这被称为SIMT（Single Instruction, Multiple Thread）模型。

一个计算欧几里得距离并原子性累加平均距离的CUDA C++代码片段展示了线程块的工作方式：

__shared__ float mean = 0.0f;
__device__ float mean_euclidian_distance(float2 *p1, float2 *p2) {
    // compute the euclidian distance between two points
    float2 dp = p2[threadIdx.x] - p1[threadIdx.x];
    float dist = sqrtf(dp.x * dp.x + dp.y * dp.y);

    // Accumulate the mean distance atomically and return distance
    atomicAdd(&mean, dist / blockDim.x);
    return dist;
}

一个网格中的所有块都运行相同程序，使用相同数量的线程，这会产生3个资源需求：
1. 块大小（Block size）：必须并发执行的线程数量。

2. 共享内存（Shared memory）：块中所有线程共用的内存。在上述代码中，__shared__ float mean = 0.0f;就是共享内存的示例。

3. 寄存器（Registers）：取决于程序复杂度的资源。寄存器是每个线程的资源，因此总预算是：(每个块的线程数 x 每个线程的寄存器数)。

GPU 如何将块放置在 SM 上 (HOW THE GPU PLACES BLOCKS ON AN SM)

GPU将线程块放置在流式多处理器（SM）上进行执行。一个块总是在单个SM上运行，并且具有固定数量的线程。

A100 SM的关键资源如下：
* 2048 线程
* 65536 寄存器
* 160 kB 共享内存

例如，一个线程块的资源需求可能为：
* 每个块的线程数：256
* 每个线程的寄存器数：64
* 每个块的寄存器数：(256 * 64) = 16384
* 每个块的共享内存：48 kB

GPU会根据这些资源需求将块分配给SM。例如，当放置第一个块（Block 0）时，它会占用SM上的相应线程、寄存器和共享内存资源。

接着，GPU会放置第二个块（Block 1），它也会消耗类似的资源。

然后是第三个块（Block 2）。

在示例中，如果每个块需要48 kB的共享内存，那么在SM的总160 kB共享内存中，只能放置3个这样的块（3 * 48 kB = 144 kB）。当尝试放置第四个块（Block 3）时，由于其48 kB的共享内存需求将使总消耗达到192 kB，超过了160 kB的SM共享内存容量，因此Block 3无法被放置。这表明共享内存可能成为一个瓶颈。

然而，如果线程块的共享内存需求更低，例如每个块仅需32 kB共享内存，那么同一个SM就可以放置更多的块。在这种情况下，160 kB的共享内存足以支持最多5个块（5 * 32 kB = 160 kB）。因此，减少每个块的共享内存需求可以提高SM的利用率。

SM 上的线程块放置概览 (Page 76-78)
上述详细步骤说明了GPU如何根据资源限制在SM上放置线程块。总结来说：
- GPU 在流式多处理器 (SM) 上放置线程块的方式。
- 一个块（Thread block）包含固定数量的线程，始终运行在一个 SM 上。
- A100 SM 关键资源：
- 2048 线程
- 65,536 寄存器
- 160 KB 共享内存

• 示例块资源需求：

  • 256 线程/块
  • 64 寄存器/线程
  • (256 * 64) = 16384 寄存器/块
  • 32 KB 共享内存/块

• 下图展示了线程块如何利用 SM 资源，图中显示了在线程资源上可以放置多达 8 个块 (Block 0-7)，在寄存器上可以放置 4 个块 (Block 0-3)，在共享内存上可以放置 5 个块 (Block 0-4)。
  
  

• 此图以三维网格形式可视化了线程块在 SM 上的放置，暗示了资源维度的约束或布局。
  

占用率 (Occupancy)

占用率概念 (Page 79)

占用率描述了 SM 上并行运行的活跃块的数量，受限于不同的资源。
对比两种情况：
- 共享内存受限情况（Shared memory limited case）：SM 主要因为共享内存不足而限制了可运行的块数量。
- 寄存器受限情况（Register limited case）：SM 主要因为寄存器不足而限制了可运行的块数量。
图示展示了两种情况下 A100 SM 关键资源（线程、寄存器、共享内存）的分配情况。

占用率量化 (Page 80)

根据资源限制，量化了两种情况下的占用率：
- 共享内存受限情况：占用率 3 块/SM。
- 寄存器受限情况：占用率 4 块/SM。

占用率是程序调优最强大的工具 (Page 81)

提高占用率可以显著提升程序性能。
从 3 块/SM 的占用率提升到 4 块/SM，可以使程序运行速度提升 33%。

填补空白：优化资源利用 (Filling in the Gaps)

资源需求 (蓝色网格) (Page 82)

介绍了 A100 SM 关键资源以及“蓝色网格”的资源需求。
资源需求 (蓝色网格)：
- 256 线程/块
- 64 寄存器/线程
- (256 * 64) = 16384 寄存器/块
- 48 KB 共享内存/块
图示展示了 SM 资源被部分“蓝色网格”块占据。

资源需求 (绿色网格) (Page 83)

在“蓝色网格”的基础上，引入了“绿色网格”的资源需求。
资源需求 (绿色网格)：
- 512 线程/块
- 32 寄存器/线程
- (512 * 32) = 16384 寄存器/块
- 0 KB 共享内存/块
此页仅显示了两种网格的资源需求，未展示填充情况。

结合蓝色和绿色网格填充资源 (Page 84)

演示了如何通过结合不同资源需求的块（蓝色和绿色）来更有效地“填补空白”，最大化 SM 资源利用率。
绿色网格的块不使用共享内存，因此可以与使用共享内存的蓝色网格块并发运行，以更好地利用 SM 资源。

并发性 (Concurrency)

并发：同时做多件事 (Page 85)

并发性是指能够同时执行多个任务。
图中显示 GPU 芯片可以同时处理图像（"Process Flower"）和进行内存拷贝（"Copy Memory"）。

并发：依赖性 (Page 86)

说明了操作之间的依赖关系。
处理一个花的图像通常是一个顺序过程：
1. 复制到 GPU (Copy to GPU)
2. 处理花 (Process Flower)
3. 从 GPU 复制 (Copy from GPU)
这些操作之间存在数据依赖性，必须按顺序执行。

并发：多个依赖性任务 (Page 87)

展示了两个独立的图像处理任务链，但它们是顺序排列的，暗示了在没有明确并发机制的情况下，即使任务独立也可能按顺序执行。
第一个任务链：复制到 GPU -> 处理花 -> 从 GPU 复制 (白色花)。
第二个任务链：复制到 GPU -> 处理花 -> 从 GPU 复制 (粉色花)。

并发：独立性 (Page 88)

引入了独立流的概念，允许独立的任务链并发执行。
通过将任务分成不同的“流”（Stream 1 和 Stream 2），可以使它们在没有数据依赖的情况下同时运行，从而提高效率。
Stream 1：处理白色花。
Stream 2：处理粉色花。

并发：本质是超额订阅 (Page 89)

并发的实现本质上是通过对 GPU 资源进行超额订阅（oversubscription）。
通过运行多个独立的流，GPU 可以同时调度来自不同任务的线程块，从而更充分地利用其 SM 资源，即使总资源需求超过了单个任务的限制，也能通过调度机制实现高效利用。
图中左侧显示了 Stream 1 和 Stream 2 独立运行，右侧显示了 A100 SM 关键资源被来自不同流（如 Flower Block 和 Copy Block）的块所占据，实现了资源的更密集利用。

并发：超额订阅的重要性 (Page 90)

再次强调了并发的重要性，它通过超额订阅来实现更高的 GPU 利用率。
图中展示了两个独立的任务流（Stream 1 和 Stream 2），它们可以并发执行，共同利用 GPU 资源。
此页与 Page 88 视觉上相似，但标题明确指出了超额订阅是并发的核心。

总结与CUDA设计理念

本节对CUDA编程的设计原理及其工作方式进行总结性探讨。

CUDA设计理念解析

幻灯片提出了一个核心问题：为什么CUDA会是现在这个样子？其设计深受物理定律驱动，旨在最大化性能。

NVIDIA Ampere A100 GPU架构的内存层次结构和数据流是关键。SM、L2缓存和HBM内存之间的协同工作，以及不同层级间的带宽优化，是实现高性能的基础。

内存访问的物理机制，特别是DRAM单元的充放电过程，直接影响了读写时序和性能。理解电容器电压公式 V_c = V_0(1-e^(-t/RC)) 有助于深入理解DRAM的物理限制。

内存访问模式对HBM内存吞吐量有显著影响。高效的内存访问模式，如连续地址访问，对于充分利用HBM的带宽至关重要。图表再次强调了HBM内存吞吐量随地址步长变化的显著性，连续访问可达峰值，而分散访问则会急剧下降。

CUDA编程模型中的线程组织方式，特别是线程块和warp的概念，是实现并行性的基础。32个线程组成的warp作为GPU的向量元素，是指令执行的基本单元。

数据局部性在GPU编程中至关重要。通过优化数据在内存中的布局和线程对数据的访问模式，可以有效利用缓存，提升处理效率。

A100 SM的关键资源分配直接影响并行执行能力。SM可以并行处理多个线程块，并需细致管理线程、寄存器和共享内存等资源，以确保高效率运行。

GPU资源的抽象和分配可以通过三维模型来表示，这有助于理解计算单元、内存区域和逻辑资源如何被线程块高效利用。

CUDA流（Streams）机制是实现并发执行不同操作的关键，它通过允许独立的任务链同时运行来提高GPU利用率。不同的流可以共享和超额订阅SM资源，从而最大化GPU的吞吐量。

最终，所有这些对硬件和编程模型的深入理解和优化，都服务于高性能计算需求，特别是在物理模拟等复杂领域。

本演讲回顾并强调了CUDA编程的工作原理及其设计哲学。其核心理念是基于物理定律，通过精妙的硬件架构和编程模型协同，实现卓越的并行计算性能。