前置讲 P2 · 并行编程与 GPU 编程模型 · Parallel Programming and the GPU Programming Model

Efficient ML Inference at Scale · 预备知识 2。个人学习讲义与配套实验,按主题整理并持续修订。

课程源码 ZIP ↓

个人整理的学习笔记,非官方讲义;课程来源与整理说明见课程介绍。配套示例未逐一运行验证。

本篇目录

定位:为没有并行编程背景的读者补充的前置讲,原课没有对应的课次。它支撑第 8 讲(硬件引擎)、第 10 讲(内核开发)、第 11 讲(内核优化),也是原课作业 4(写内核)的基础。建议先读 P1。 依据:NVIDIA CUDA C++ Programming Guide(编程模型、SIMT 架构、性能指南,以及各计算能力的规格表;本讲用 12.9 版);J. Austin et al., How to Scale Your Model 第 12 章(GPU);M. Harris 的两份经典材料 Optimizing Parallel Reduction in CUDA(NVIDIA, 2007)和 An Efficient Matrix Transpose in CUDA C/C++(NVIDIA 技术博客);V. Volkov, Better Performance at Lower Occupancy(GTC 2010)。 配套代码 Code:code/prereq2/(numpy;共享模块 simt_lab.py 是一台”用 numpy 模拟的 GPU”,+ 5 个脚本,最慢的约 3 秒) 预计阅读时间:3 小时(含运行代码)


0. 本讲导航 · How to read this lecture

0.1 学习目标 Learning goals

  1. 说清楚 GPU 为什么是”吞吐导向”的设计,以及它和 CPU 在晶体管、功耗、延迟处理上的取舍(承接 P1 §7)。
  2. 读懂一个最简单的 CUDA 内核:kernel、thread、block、grid,blockIdx * blockDim + threadIdx 这一行在算什么,为什么要边界检查,grid-stride 循环是什么。
  3. 解释硬件怎样执行内核:SM、warp、SIMT、分支发散、warp 调度与零开销切换。
  4. 判断一段访存是否合并(coalesced),会算一个 warp 的访问碰到几个 32 字节的 sector;用共享内存把不合并的访问变成合并的(矩阵转置)。
  5. 解释共享内存的 bank 冲突、为什么需要 __syncthreads()、漏写会怎样;读懂并行归约的几个版本各自消掉了什么开销。
  6. 用占用率计算器判断一个内核被什么资源卡住;用 Little 定律解释”占用率 vs 指令级并行”;知道波次量化。
  7. 用 Amdahl 定律和”内核启动开销”解释为什么小模型、小批的 LLM decode 会被 CPU 端拖慢,以及算子融合、CUDA Graphs、megakernel 的动机。

0.2 本讲地图 Map

§1  为什么要并行编程:数据并行;Amdahl 与 Gustafson
§2  从 CPU 到 GPU:吞吐导向的设计
§3  CUDA 编程模型:kernel、thread、block、grid;向量加法;二维网格;grid-stride;我们的模拟器      ← 代码 01
§4  硬件怎样执行:SM;warp 与 SIMT;分支发散;warp 调度与延迟隐藏
§5  全局内存与合并访存:存储空间一览;sector;访问模式;矩阵转置                                  ← 代码 02
§6  共享内存、同步与归约:bank 冲突;__syncthreads;并行归约的五个版本;warp shuffle              ← 代码 03
§7  占用率、延迟隐藏与波次                                                                        ← 代码 04
§8  Amdahl 定律与内核启动开销                                                                     ← 代码 05
§9  CUDA 之外:Triton、TPU 与其他编程模型
§10–14  小结、常见误区、练习、延伸阅读、和正课的联系

0.3 我们没有 GPU,怎么学?

simt_lab.py 是一台”纸上的 GPU”。你像写 CUDA 内核一样写一个 Python 函数,模拟器按网格和线程块的形状”启动”它,并且数出真实 GPU 上决定性能的那些事件:每个 warp 的访存碰到几个 sector、共享内存有没有 bank 冲突、有没有分支发散、漏写的同步造成了哪些数据竞争。

它和真实 GPU 有一个重要的区别:模拟器把一个线程块的所有线程放在 numpy 数组里,一条语句一次算完全部线程(所有线程严格同步地执行每一条语句)。所以它算出的结果总是确定的;为了弥补这一点,它另外记录”谁在哪个同步阶段读写过共享内存”,在真实 GPU 上会出问题的地方报告竞争(§6.2)。这种”整块向量化”的写法,其实很接近 Triton 和 TPU 的编程模型(§9)。

所有 P2 的数字都来自模拟(真实计算 + 事件计数),除了演示 05 的最后一部分是墙钟测量。模拟器的 GPU 参数是 H100(计算能力 9.0)的公开规格,几个标明的假设值(访存延迟、启动开销)见代码注释。


1. 为什么要并行编程 · Why parallel programming

1.1 数据并行

P1 §6 说过:单核频率不再上涨,性能只能靠并行。并行有很多种,但对 AI 来说最重要的是数据并行(data parallelism):对大量数据做同样的操作,每个数据的计算彼此独立。

  • 向量加法 ci=ai+bic_i = a_i + b_i:一百万个 ii,每个都可以同时算。
  • 矩阵乘法 C=ABC = AB:CC 的每个元素都是一次独立的点积。
  • LLM 的一步:批里每个 token 走同一个网络;同一个矩阵乘法里,每个输出元素独立。

数据并行的程序写起来有一个统一的模式:写一段”处理一个数据”的代码,然后让成千上万个执行者各处理一个。GPU 编程模型就是围绕这个模式设计的。

1.2 Amdahl 与 Gustafson

但总有一部分不能并行。P1 §6.4 的 Amdahl 定律:如果比例为 ss 的部分只能串行,

SAmdahl(n)=1s+(1−s)/n≤1s.S_\text{Amdahl}(n) = \frac{1}{s + (1-s)/n} \le \frac{1}{s} .

在 GPU 的尺度上(H100 有 16,896 个 FP32 lane),它非常残酷(代码 05 演示 1):

串行比例 ssn=132n = 132(SM 数)n=16,896n = 16{,}896(H100 lane 数)n→∞n \to \infty
1%5799100
0.1%1179441,000
0.01%1306,28210,000

另一种看法是 Gustafson 定律:处理器多了,人们会去解更大的问题——并行部分变大、串行部分的时间不变,于是 SGustafson(n)=s+(1−s)nS_\text{Gustafson}(n) = s + (1-s) n,几乎线性(s=1%s = 1\%、n=16,896n = 16{,}896 时约 16,727)。

两者不矛盾,回答的是不同问题:Amdahl 问”同一个任务能快多少”(延迟),Gustafson 问”同样的时间能多做多少”(吞吐)。LLM 推理两者都要:单个请求的 decode 延迟受 Amdahl 限制(§8),批处理的吞吐更像 Gustafson。


2. 从 CPU 到 GPU:吞吐导向的设计 · Throughput-oriented design

P1 §7 列出了 CPU 和 GPU 对同一组物理限制的两种回答。这里再展开 GPU 这一边,因为后面所有的编程规则都从这里来:

  1. 晶体管花在运算单元上,而不是控制逻辑上。一个 H100 的 SM 有 128 个 FP32 单元、4 个 Tensor Core、4 个 warp 调度器(CUDA 编程指南),整块芯片 132 个 SM。没有复杂的分支预测和乱序执行。
  2. 用海量线程隐藏延迟,而不是用 cache 和预测避免延迟。P1 §4.4:一块 H100 要跑满 HBM 带宽,需要约 1.7 MB 的读请求同时在途。CPU 靠一个核里的乱序执行和预取凑出几十个在途请求;GPU 靠让成千上万个线程同时发请求——一个线程在等数据,硬件立刻切换去执行另一个线程。
  3. 一条指令驱动很多个运算单元(SIMT,§4.2),把取指、译码的开销摊薄(第 1 讲 §2.5)。
  4. 又宽又慢:频率约 2 GHz,靠单元数量取胜,更省电(P1 §6.5)。

代价是:GPU 只擅长数据并行、规整的计算。单个线程在 GPU 上比在 CPU 上慢得多;分支多、不规整的代码会浪费大量运算单元。程序员的工作,就是把问题组织成 GPU 喜欢的样子。


3. CUDA 编程模型 · The CUDA programming model

代码:01_cuda_programming_model.py

3.1 kernel、thread、block、grid

CUDA 是 NVIDIA 的 GPU 编程模型(C++ 的扩展)。几个概念:

  • host 与 device:CPU 和它的内存叫 host,GPU 和它的显存叫 device。程序在 host 上运行,把数据拷到 device,再让 GPU 执行计算。
  • kernel(内核):在 GPU 上执行的函数,用 __global__ 标记。你只写”一个线程要做什么”。
  • thread(线程):执行 kernel 的最小单位。启动一个 kernel 时会创建成千上万个线程,每个线程执行同一份代码,但用自己的编号去处理不同的数据。
  • block(线程块):线程被分成一组一组,每组最多 1024 个线程。同一个块的线程一定在同一个 SM 上运行,可以通过共享内存交换数据、用 __syncthreads() 同步。不同块之间不能这样合作。
  • grid(网格):一次 kernel 启动的所有块。块之间相互独立,可以按任意顺序、在任意 SM 上执行——这让同一个程序能在 SM 数量不同的 GPU 上运行。

每个线程通过几个内置变量知道”我是谁”:threadIdx(块内编号)、blockIdx(块编号)、blockDim(块的形状)、gridDim(网格的形状),各有 x、y、z 三个分量。

3.2 第一个内核:向量加法

__global__ void vector_add(const float* a, const float* b, float* c, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;   // 我是全局第几个线程?
    if (i < n) {                                     // 边界检查
        c[i] = a[i] + b[i];
    }
}
// host 端启动:⌈n / 256⌉ 个块,每块 256 个线程
int threads = 256;
int blocks  = (n + threads - 1) / threads;
vector_add<<<blocks, threads>>>(d_a, d_b, d_c, n);

逐行看:

  • blockIdx.x * blockDim.x + threadIdx.x:第 blockIdx.x 个块之前有 blockIdx.x × blockDim.x 个线程,再加上块内编号,就是全局编号。这一行几乎出现在每个 CUDA 内核里。
  • 块数用向上取整的除法算:⌈n/256⌉\lceil n / 256 \rceil。所以最后一个块通常有多余的线程。
  • if (i < n):多余的线程什么都不做。漏掉它,多余的线程就会读写数组外面的内存(§3.3)。
  • <<<blocks, threads>>>:CUDA 的启动语法,指定网格和块的形状。

在我们的模拟器里,同一个内核写成:

def vector_add(blk, a, b, c, n):
    i = blk.bx * blk.bdx + blk.tx          # 一个数组:块内每个线程一个值
    m = i < n                              # 边界检查变成一个布尔数组(mask)
    x = blk.load(a, i, mask=m)             # 只有 mask 为 True 的线程真的读
    y = blk.load(b, i, mask=m)
    blk.store(c, i, x + y, mask=m)

launch(vector_add, grid=(cdiv(n, 256),), block=(256,), args=(a, b, c, n))

if 变成了 mask:这正是 GPU 硬件执行分支的方式(§4.3)。

块大小怎么选? 代码 01 演示 1(n=10,000n = 10{,}000):

每块线程数块数warp 总数SIMT 效率合并效率
3231331399.8%100%
10010040078.1%86%
2564032099.8%100%
10241032099.8%100%

块大小 100 有两个代价:每块 4 个 warp,最后一个 warp 只有 4 个线程在干活(SIMT 效率 78%,§4.2);而且第 kk 个块从第 100k100k 个元素(第 400k400k 字节)开始,大多数块的起点没有对齐到 32 字节,访存会多碰 sector(合并效率 86%,§5)。块大小取 32 的倍数——CUDA 编程指南也这样建议。块大小在 128–1024 之间时,这个简单内核的计数几乎一样;块大小主要影响的是占用率(§7)。

3.3 忘了边界检查

把 if (i < n) 去掉、n=1000n = 1000、每块 256 个线程:最后一个块的第 232–255 号线程会读 a[1000]…a[1023]。模拟器直接报错(代码 01 演示 2);真实 GPU 上,越界读可能悄悄读到别的数据(结果出错但不报错),也可能触发非法地址错误让整个程序崩溃;越界写更危险。这是 CUDA 新手最常见的 bug,排查工具是 NVIDIA 的 compute-sanitizer。

3.4 二维网格

二维问题(矩阵)用二维的块和网格更自然:

int r = blockIdx.y * blockDim.y + threadIdx.y;   // 行
int c = blockIdx.x * blockDim.x + threadIdx.x;   // 列
if (r < M && c < N) X[r * N + c] += bias[r];

一个关键细节:线程按 x 最快变化的顺序编号(tid=tx+ty⋅blockDim.x\text{tid} = t_x + t_y \cdot \text{blockDim.x}),再每 32 个打包成一个 warp。16×16 的块里,一个 warp 是两行、每行 16 个线程。所以 x 方向应该对应内存里连续的那一维(按行存储时是列)。代码 01 演示 3 把 x 方向错配给行:结果仍然正确,但合并效率从 95% 掉到 35%,访存 sector 数从 2,750 涨到 7,455。

3.5 grid-stride 循环

另一种常见写法让网格大小和数据大小脱钩:

for (int i = blockIdx.x * blockDim.x + threadIdx.x; i < n; i += gridDim.x * blockDim.x)
    y[i] = alpha * x[i];

每个线程从自己的全局编号出发,每次跳过”全网格的线程总数”。代码 01 演示 4:4 个块 × 128 线程、16 × 256、391 × 256 三种网格处理 10 万个元素,结果都对,warp 级 load 指令总数都是 3,125——工作量由数据决定,网格只决定谁来做。实践中常把网格设成”SM 数的几倍”,让每个 SM 都有活干(§7.3)。


4. 硬件怎样执行 · How the hardware executes a kernel

4.1 SM

GPU 由几十到上百个 SM(streaming multiprocessor,流式多处理器)组成,H100 有 132 个。按 CUDA 编程指南,计算能力 9.0(H100)的每个 SM 有:128 个 FP32 单元、64 个 FP64 单元、64 个 INT32 单元、4 个 Tensor Core、16 个特殊函数单元(SFU,算 exp、sin 等)、4 个 warp 调度器;一个 65,536 个 32 位寄存器的寄存器堆(256 KB);一块 256 KB 的 SRAM,可以按需在 L1 cache 和共享内存之间划分(共享内存最多 228 KB)。

kernel 启动后,硬件把块分配到 SM 上:一个 SM 可以同时驻留多个块(最多 32 个块、64 个 warp、2048 个线程,还要看寄存器和共享内存够不够,§7.1)。块一旦分配到某个 SM,就一直在那里执行到结束。

4.2 warp 与 SIMT

SM 不是一个线程一个线程地执行,而是32 个线程一组,叫一个 warp。一个块的线程按编号每 32 个分成一个 warp(第 0–31 号是 warp 0,以此类推)。warp 调度器每个周期从它负责的 warp 里挑一个”准备好了”的,发射一条指令,这条指令由 warp 里的 32 个线程同时执行。

NVIDIA 把这叫 SIMT(single instruction, multiple threads)。它和 P1 §6 的 SIMD 很像——一条指令、多个数据——区别在于编程模型:SIMD 要程序员(或编译器)显式地把数据装进向量寄存器;SIMT 让你写”一个线程”的代码,硬件自动把 32 个线程打包执行。每个线程有自己的寄存器、自己的程序计数器(从 Volta 开始),看起来像独立的线程,执行起来像 SIMD。

SIMT 效率:一条 warp 指令执行时,32 个线程里真正在干活的比例。线程数不是 32 的倍数时,最后一个 warp 不满(§3.2 的块大小 100);分支发散时,一部分线程闲着(§4.3)。

4.3 分支发散

if (threadIdx.x % 2 == 0) { A(); } else { B(); }

一个 warp 只有一条指令流。如果 warp 里有的线程要走 A、有的要走 B,硬件只能先执行 A 路径(走 B 的线程闲着),再执行 B 路径(走 A 的线程闲着)。两条路的时间加在一起——这叫分支发散(branch divergence)。模拟器的 blk.branch(cond) 返回走 if 的那部分线程的 mask,同时记录发生发散的 warp。

避免发散的办法:让分支条件以 warp 为单位一致。CUDA 编程指南的原话大意是:如果分支条件只依赖 threadIdx / warpSize,就不会有 warp 发散。§6.3 的归约会看到,仅仅改变”哪些线程干活”,就把发散消掉了 90%。

4.4 warp 调度与延迟隐藏

一条访存指令发出去,要几百个周期才回来。这期间这个 warp 不能执行依赖这个数据的指令。GPU 的办法很简单:换一个 warp。

关键在于切换是零开销的:每个 warp 的寄存器一直留在寄存器堆里(这就是为什么 SM 的寄存器堆有 256 KB 那么大),切换 warp 不需要保存和恢复任何状态,调度器下一个周期直接从另一个 warp 发射指令。CPU 的线程切换要保存寄存器、刷新流水线,开销是微秒级的。

所以 GPU 隐藏延迟的能力取决于:每个 SM 上有多少个 warp 可以轮流执行(占用率,§7.1),以及每个 warp 自己能同时发出多少个独立的请求(指令级并行,§7.2)。这又回到了 P1 §4 的 Little 定律。


5. 全局内存与合并访存 · Global memory and coalescing

代码:02_memory_coalescing.py

5.1 GPU 的存储空间

存储空间位置谁能访问容量(H100)延迟量级
寄存器SM 内单个线程私有每 SM 256 KB(每线程最多 255 个)几乎为 0
共享内存SM 内(和 L1 共用 256 KB)同一个块的线程每 SM 最多 228 KB约 20 ns
L1 cacheSM 内硬件自动管理上一项剩下的部分约 20 ns
L2 cache芯片上,所有 SM 共享硬件自动管理50 MB一两百 ns
全局内存(HBM)芯片旁所有线程、host 也能拷入拷出80 GB几百 ns

(延迟是 Gimlet 博客引用 Jarmusch et al. 的测量量级,只作参考。)

程序员显式管理的是寄存器(局部变量)、共享内存(__shared__)和全局内存(cudaMalloc 分配的数组)。L1、L2 是硬件自动管理的 cache。

5.2 合并访存:以 sector 为单位

GPU 的全局内存以 32 字节的 sector 为单位搬数据。一个 warp 执行一条 load 指令时,32 个线程各给出一个地址,硬件把这些地址落在的 sector 合并起来一起搬。所以一条 warp load 的代价取决于它碰到了几个 sector——这叫合并访存(memory coalescing)。

代码 01 / 02 的模拟器正是这样计数的。一个 warp 读 32 个数的各种模式(代码 02 演示 1):

访问模式每条指令的 sector 数实际搬的字节合并效率
连续:线程 ii 读第 ii 个 float324128 B100%
错位 1 个元素5160 B80%
跨 2 个元素8256 B50%
跨 8 个元素(32 B)及以上321024 B12.5%
广播:所有线程读同一个数132 B—
随机下标≈321024 B12.5%
连续读 float648256 B100%
连续读 16 字节(向量化读取)16512 B100%

**规则:相邻的线程访问相邻的地址。**跨度一大,每个线程独占一个 sector,搬 1024 字节只用 128 字节——这和 P1 §2.5 “跨度达到一整行,cache line 只用了一点点”是同一件事,只是单位从 cache line 换成了 sector。更宽的类型和向量化读取(每个线程一次读 16 字节)在连续访问时仍然 100% 有效,而且一条指令搬的字节更多,对隐藏延迟有帮助(§7.2)。

5.3 案例:矩阵转置

转置 B=A⊤B = A^\top 是”合并访存”最经典的教学案例(Harris 的 NVIDIA 博客)。AA 和 BB 都按行存储。块是 32×8 个线程,每块负责 32×32 的一块,每个线程处理 4 个元素。

朴素转置:B[x][y] = A[y][x]。读 AA 时相邻线程读同一行的相邻元素(合并);写 BB 时相邻线程写同一列的相邻元素,地址相隔一整行(不合并)。

共享内存分块:

__shared__ float tile[32][32];
// 1. 按行把 A 的一块读进 tile(合并)
for (int j = 0; j < 32; j += 8) tile[threadIdx.y + j][threadIdx.x] = A[(y + j) * N + x];
__syncthreads();                          // 2. 等整块都写进来
// 3. 换成转置后的块位置,按列从 tile 读、按行写到 B(合并)
x = blockIdx.y * 32 + threadIdx.x;  y = blockIdx.x * 32 + threadIdx.y;
for (int j = 0; j < 32; j += 8) B[(y + j) * N + x] = tile[threadIdx.x][threadIdx.y + j];

“转置”这件不规整的事,挪到了片上的共享内存里做;全局内存的读写都是连续的。但第 3 步按列读 tile 时,32 个线程读 tile[0][c]、tile[1][c]……,地址相隔 32 个字——全部落在同一个 bank 上(§6.1),32 路冲突。解决办法:tile[32][33],每行多填一个元素。

代码 02 演示 2(512×512):

内核写 sector合并效率共享内存每条指令几轮竞争估计 HBM 时间(相对 copy)
copy(基线)32,768100%–01.00×
朴素转置262,14422%–04.50×
共享内存分块32,768100%16.501.00×
分块 + 填充一列32,768100%1.001.00×
分块 + 填充,但忘了 sync32,768100%1.01,0241.00×

对照 Harris 在真实 GPU 上的实测(1024×1024,有效带宽 GB/s):

内核Tesla M2050Tesla K20c
copy105.2136.0
朴素转置18.855.3
共享内存分块51.397.6
分块 + 填充(无 bank 冲突)99.5144.3

顺序完全一致:朴素转置只有 copy 的 18%–41%;分块后提到 49%–72%(还被 bank 冲突拖着);消除 bank 冲突后达到最快 copy 的约 95%。

两点要注意:

  1. 模拟器估计朴素转置的 HBM 时间是 copy 的 4.5 倍;Harris 实测在 M2050 上慢 5.6 倍、在 K20c 上只慢 2.5 倍。真实 GPU 的 L2 cache 会把相邻 warp 写同一个 sector 的部分合并起来,模拟器没有建模这一点。计数器是用来比较写法、找出问题的,不是用来预测时间的。
  2. 忘了 sync 的版本,模拟器算出的结果是对的,但报告了 1,024 次竞争——§6.2 解释。

6. 共享内存、同步与归约 · Shared memory, synchronization, and reduction

代码:03_shared_memory_sync_reduction.py

6.1 bank 与 bank 冲突

共享内存是 SM 上一块程序员自己管理的 SRAM(可以看成一个”手动 cache”)。为了让一个 warp 的 32 个线程能同时访问,它被分成 32 个 bank,每个 bank 每个周期能服务一个 4 字节的字;第 kk 个字在第 k mod 32k \bmod 32 个 bank。

  • 32 个线程访问 32 个不同的 bank:一轮完成。
  • 多个线程访问同一个字:广播,也是一轮。
  • mm 个线程访问同一个 bank 里的不同的字:要分 mm 轮——mm 路 bank 冲突。

代码 03 演示 1:

访问模式(下标)需要的轮数
tid(连续)1
2·tid / 4·tid / 8·tid / 16·tid2 / 4 / 8 / 16
32·tid(按列读 32×32 的 tile)32
33·tid(按列读 32×33 的 tile)1
3·tid(奇数跨度)1
所有线程读同一个字1

跨度是 2k2^k 时冲突 2k2^k 路;跨度是奇数时没有冲突(奇数和 32 互质)。这和 P1 §2.5 的 cache 冲突缺失是同一个数学问题:地址按 2 的幂次方取模分组,2 的幂次方的跨度让所有访问落到同一组;解法也一样——填充一个元素错开。

6.2 同步:__syncthreads()

一个块里的线程不是严格同步执行的:8 个 warp 各自推进,谁快谁慢不确定。所以当一个线程要读另一个线程写进共享内存的数据时,必须保证对方已经写完——这就是 __syncthreads() 的作用:块里所有线程都到达这里之后,才一起继续;之前的共享内存读写对之后的所有线程可见。

NVIDIA 博客里的经典小例子:把一个块里的数组倒过来。

__shared__ float s[256];
s[t] = d[t];
__syncthreads();            // 下一行读的是别的线程写的元素
d[t] = s[255 - t];

漏写 __syncthreads(),代码 03 演示 2 的模拟器给出:结果是对的(模拟器让所有线程严格同步),但报告了竞争:“warp 0 读共享内存 [255],它在本阶段刚被 warp 7 写过,中间没有 sync”。在真实 GPU 上,warp 0 可能在 warp 7 写入之前就去读,得到旧值。这类 bug 最可怕的地方是它可能 99 次碰巧正确、第 100 次出错。(NVIDIA 的 compute-sanitizer --tool racecheck 能检查。)

模拟器检查两种竞争:写后读(读到别的 warp 本阶段写的数据)和读后写(覆盖别的 warp 本阶段刚读过的数据——比如分块矩阵乘法的循环里,下一轮装载新 tile 之前漏了第二个 __syncthreads(),第 10 讲会遇到)。

还要记住:__syncthreads() 只同步一个块内的线程。块之间没有这种同步;要交换数据,通常只能结束这个 kernel、再启动下一个(或者用原子操作、协作组等更复杂的机制)。

6.3 并行归约:五个版本

求和(归约,reduction)看起来简单,却是理解 GPU 优化的最好例子。Mark Harris 2007 年的报告 Optimizing Parallel Reduction in CUDA 一步一步优化了这个问题,下面是其中几个版本(代码 03 演示 3,2182^{18} 个数、每块 256 线程,每块先算出一个部分和):

版本 1:交错寻址,分支发散。第 ss 步(s=1,2,4,…s = 1, 2, 4, \dots),编号是 2s2s 倍数的线程做 sdata[tid] += sdata[tid + s]。if (tid % (2*s) == 0) 让每个 warp 里一半、四分之三……的线程闲着,严重发散。

版本 2:让连续的线程干活。第 ss 步线程 tid 处理 index = 2*s*tid。干活的线程是连续的一段,整个 warp 要么全干、要么全不干,发散几乎消失;但 2*s*tid 的跨度是 2s2s,带来 bank 冲突。

版本 3:顺序寻址。第 ss 步(ss 从 128 减半到 1)前 ss 个线程做 sdata[tid] += sdata[tid + s]。干活的线程连续(无发散),访问的地址也连续(无 bank 冲突)。

版本 4:加载时先加一次。版本 3 的第一步就有一半线程闲着。让每个块负责 512 个元素、每个线程读两个数先加起来,块数减半。

版本 5:最后一个 warp 用 shuffle。归约到剩 64 个数后,第 0 个 warp 用 __shfl_down_sync 在寄存器之间直接交换数据:16、8、4、2、1 五步,不用共享内存、不用同步。

版本块数warp 指令相对版本 1SIMT 效率发散的 warpbank 每指令轮数sync 次数
1 交错寻址(发散)1024350,2081.0054%48,1281.009,216
2 连续线程(bank 冲突)1024206,8480.5991%5,1203.339,216
3 顺序寻址1024206,8480.5991%5,1201.009,216
4 加载时先加一次512111,6160.3292%2,5601.004,608
5 + warp shuffle51256,3200.1699%01.001,536

每一步消掉了一种浪费:发散 → bank 冲突 → 闲置的线程 → 最后几步的同步和共享内存往返。最终的指令数只有版本 1 的 16%。

这些优化为什么重要?归约的强度极低(每读 4 字节只做 1 次加法),理想情况下应该是 memory-bound:读完 1 MB 输入至少要约 320 ns。模拟器估计版本 1 光发射指令就要约 335 ns——它受限于指令,而不是带宽;后面的版本把指令开销压到读数据的时间以下,内核才真正受限于带宽。Harris 的报告也是以”接近峰值带宽”为优化目标的。

一个历史注脚:Harris 2007 年的原版在最后一个 warp 里不做同步,依赖”warp 内 32 个线程严格同步执行”。从 Volta 架构开始,warp 内的线程也可以独立推进,那个技巧不再安全;现代写法用 __syncwarp() 或者 shuffle(版本 5)。


7. 占用率、延迟隐藏与波次 · Occupancy, latency hiding, and waves

代码:04_occupancy_latency_waves.py

7.1 占用率

占用率(occupancy)= 一个 SM 上同时驻留的 warp 数 ÷ 最大 warp 数(H100 是 64)。一个块被放到 SM 上以后,它所有线程的寄存器、它的共享内存一直占着,直到块结束。四种资源会卡住能放的块数:

  1. 线程(warp):每 SM 最多 64 个 warp;
  2. 块:每 SM 最多 32 个块;
  3. 寄存器:每 SM 65,536 个,按 warp 以 256 个为单位分配;
  4. 共享内存:每 SM 228 KB,每块另有 1 KB 系统保留。

代码 04 演示 1 的占用率计算器:

配置每 SM 块数活跃 warp占用率瓶颈
256 线程/块,32 寄存器/线程864100%warp 数
256 线程,64 寄存器43250%寄存器
256 线程,128 寄存器21625%寄存器
256 线程,255 寄存器1812.5%寄存器
128 线程,48 KB 共享内存41625%共享内存
128 线程,100 KB 共享内存2812.5%共享内存
32 线程/块323250%块数
230 KB 共享内存00—无法启动

寄存器是最常见的瓶颈:复杂的内核(矩阵乘法、注意力)每线程常常要上百个寄存器来放累加器——第 1 讲 §3.9 说”输出块的累加器要一直待在片上”,就是待在这里。

7.2 延迟隐藏:占用率 vs 指令级并行

占用率越高越好吗?不一定。真正需要的是:在途的请求足够多(P1 §4 的 Little 定律)。一个 SM 分到的 HBM 带宽约 25 GB/s,访存延迟按 500 ns 算(假设值),需要约 12.7 KB 的请求同时在途。一个 warp 能贡献的在途字节 = 它同时发出的读请求数 × 每个请求的字节数。代码 04 演示 2:

活跃 warp(占用率)4 B/线程,1 个请求4 B/线程,2 个请求16 B/线程,1 个请求16 B/线程,2 个请求
16(25%)16%32%65%100%
32(50%)32%65%100%100%
64(100%)65%100%100%100%
  • 每个线程一次只读 4 字节、每个 warp 只有 1 个请求在途时,就算 100% 占用率也只能拿到约 65% 的带宽。
  • 让每个 warp 同时发出多个请求(指令级并行,循环展开就能做到),或者每个线程一次读 16 字节(向量化读取),只要 25% 的占用率就能喂饱带宽。

这就是 Volkov 在 GTC 2010 的演讲 Better Performance at Lower Occupancy 的观点:**占用率只是隐藏延迟的手段之一,每个线程多做几件独立的事(ILP)可以替代很多个线程(TLP)。**高性能的矩阵乘法和注意力内核通常就是”低占用率、每个线程干很多活、高 ILP”的风格。

7.3 波次量化

所有 SM 同时能跑 132 × (每 SM 块数) 个块,这叫一”波”(wave)。块更多时,要等前一波结束才能开始下一波。块数恰好比一波多一点时,最后一波几乎是空的(代码 04 演示 3,每 SM 放 1 个块):

块数波数整体 SM 利用率
1321100%
133250.4%
2642100%
265366.9%
10,0007699.7%

这在分块矩阵乘法里形成一个取舍:块越大,每个块的算术强度越高(第 1 讲 §3.9,约 tmtn/(tm+tn)t_m t_n / (t_m + t_n)),但块数越少:

输出 CC块大小块数SM 利用率块的强度
4096×4096256×25625697%128
2048×2048256×2566448%128
512×1024128×1283224%64
512×1024256×25686%128

大矩阵用大块没问题;小矩阵用大块,132 个 SM 只用了几个。这时宁可用小块,或者把 KK 维也切开(split-K)让更多 SM 参与。LLM decode 的矩阵乘法正是”很瘦”的情形(一个维度是批大小,可能只有几到几十),块数不够填满 GPU 是 decode 内核的常见问题——第 10–11 讲的 split-K、flash-decoding 就是为它设计的。


8. Amdahl 定律与内核启动开销 · Launch overhead

代码:05_amdahl_and_launch_overhead.py

8.1 每启动一次内核都要付钱

每次启动一个 kernel,CPU 端要准备参数、经过驱动把命令提交给 GPU,GPU 再把块分配到 SM 上——这段固定开销通常是几微秒。如果一个 kernel 本身要跑几毫秒,这点开销可以忽略;但 LLM decode 在 batch = 1 时,每一步要启动几百个小 kernel。

以一层 Transformer 为例,一种有代表性的划分(代码 05 演示 2,是假设,不对应某个具体引擎):

  • 未融合(每层 12 个):RMSNorm、QKV 投影、RoPE、写 KV cache、注意力、输出投影、残差相加、RMSNorm、gate/up 投影、SiLU × 乘、down 投影、残差相加。
  • 适度融合(每层 8 个):残差+RMSNorm、QKV 投影、RoPE+写 KV、注意力、输出投影、残差+RMSNorm、gate/up 投影+SiLU×乘、down 投影。

再加上最后的 RMSNorm、输出层、采样 3 个。GPU 真正的工作时间取第 1 讲的 roofline 下界(H100、BF16、ctx = 1024),每次启动开销扫描 2、4、8 µs(假设值):

模型GPU 工作下界内核划分内核数/步启动开销 4 µs 时每步(不重叠)开销占比
Llama-3.1-8B4.52 ms未融合3876.07 ms26%
Llama-3.1-8B4.52 ms适度融合2595.56 ms19%
Llama-3.2-1B0.75 ms未融合1951.53 ms51%
Llama-3.2-1B0.75 ms适度融合1311.27 ms41%
任意–megakernel1≈ GPU 时间≈0%

“不重叠”是上界;如果 CPU 提前把 kernel 排进队列、和 GPU 的执行重叠,下界是 max(GPU 时间, 内核数 × 启动开销)。真实系统在两者之间。即使完全重叠,1B 模型在启动开销 8 µs 时,CPU 发 kernel 的速度也跟不上 GPU(每步从 0.75 ms 变成 1–1.6 ms)。

模型越小、批越小,启动开销的占比越高。这就是 Amdahl 定律在推理里的样子:GPU 再快,CPU 端串行的部分也压不下去。三种对策:

  1. 算子融合(operator fusion):把几个相邻的小 kernel 合成一个,减少 kernel 个数,还能少读写中间结果(第 11、12 讲)。
  2. CUDA Graphs:把一整串 kernel 的启动”录”下来,之后一次性提交,大幅降低每个 kernel 的 CPU 端开销。主流推理引擎在 decode 阶段都用它。
  3. megakernel:把整个前向写成一个常驻 GPU 的 kernel,在内部自己调度各个算子,彻底消灭 kernel 边界。第 11 讲的阅读材料——Hazy Research 的 Designing a Low-Latency Megakernel for Llama-1B——针对的正是这里的 1B 模型。

8.2 在你的电脑上看到同样的现象

每次 numpy 调用也有固定开销(Python 解释器、参数检查……),是 kernel 启动开销的 CPU 版本。代码 05 演示 3 的实测(编写时的 M4 Pro,只看数量级):对 1 个元素做 np.add 每次约 0.16 µs,几乎全是固定开销;到 26 万个元素时固定开销只占约 1%。同样把 26 万个元素相加,切成 128 次小调用比一次大调用慢约 5 倍。


9. CUDA 之外:其他编程模型 · Beyond CUDA

CUDA 的”每个线程一份代码”(SIMT)不是唯一的写法。你会在后面的讲次里遇到:

  • Triton(OpenAI):用 Python 写 GPU 内核,但编程单位是一个块而不是一个线程——你对一整块数据写向量化的操作(像 numpy),编译器负责把它映射到线程、处理合并访存和共享内存。我们的模拟器其实就是这个风格。第 10、13 讲。
  • CUTLASS / CuTe:NVIDIA 的 C++ 模板库,用”布局代数”描述数据怎样切块、怎样映射到线程和 Tensor Core,是写高性能矩阵乘法和注意力的主力工具。
  • TPU 与 Pallas:TPU 没有”线程”的概念,程序员(或 XLA 编译器)面对的是大的向量和矩阵单元、以及手动管理的片上存储(VMEM);Pallas 是 JAX 里写 TPU/GPU 内核的接口,也是块级编程。第 8 讲。
  • AWS NKI:AWS Trainium/Inferentia 芯片的内核接口,同样是块级、显式管理片上存储(SBUF、PSUM)。第 9 讲的阅读材料。

它们的共同点都在这一讲里:把数据切成放得进片上存储的块,让访存连续、让计算单元满负荷、让同步尽量少。


10. 小结 · Summary

  1. 数据并行是 GPU 编程的核心模式:写”一个线程做什么”,让成千上万个线程各处理一个数据。Amdahl 定律限制了固定规模问题的加速,Gustafson 定律说明扩大规模后吞吐几乎线性。
  2. GPU 是吞吐导向的:晶体管花在运算单元上,用海量线程(而不是 cache 和预测)隐藏延迟,一条指令驱动 32 个线程,又宽又慢更省电。代价是只擅长规整的数据并行计算。
  3. CUDA 编程模型:kernel 由 grid 中的 block 中的 thread 执行;blockIdx * blockDim + threadIdx 算全局编号;向上取整的网格 + 边界检查;块大小取 32 的倍数;二维时 x 方向对应内存连续的那一维。
  4. 硬件执行:SM 以 warp(32 线程)为单位发射指令(SIMT);分支发散让两条路径都要执行;warp 切换零开销,靠足够多的 warp 和足够多的在途请求隐藏延迟。
  5. 合并访存:全局内存以 32 字节 sector 为单位搬运,相邻线程访问相邻地址才不浪费。不规整的访问(转置)先在共享内存里重排。
  6. 共享内存有 32 个 bank,2 的幂次方跨度造成冲突,填充一个元素解决;读别的线程写的数据前必须 __syncthreads(),漏写是不确定的 bug。归约的优化依次消掉发散、bank 冲突、闲置线程和同步开销。
  7. 占用率受寄存器、共享内存、线程数、块数限制;但隐藏延迟要的是”在途请求足够多”,指令级并行和向量化读取可以替代高占用率。块数和 SM 数对不上时有波次损失;小矩阵用大块会让大部分 SM 闲着。
  8. 启动开销:每个 kernel 几微秒的固定开销,在小模型、小批的 decode 里占比可达一半;融合、CUDA Graphs、megakernel 是对策。

11. 常见误区 · Common pitfalls

  1. “GPU 线程和 CPU 线程一样。” GPU 线程是 warp 里的一个 lane,没有自己独立的指令流水线;单个线程在 GPU 上很慢,快是因为同时有几万个。
  2. “块里的线程是同时执行每条语句的,所以不用同步。” 不同 warp 之间推进顺序不确定(从 Volta 起 warp 内也不保证);读别的线程写的共享内存数据前必须同步。
  3. “只要结果对,内核就没问题。” 竞争类 bug 可能大多数时候碰巧正确;越界访问也可能悄悄读到别的数据而不报错。
  4. “访存次数决定速度。” 决定速度的是搬了多少个 sector(合并与否),以及在途请求够不够多,而不是 load 指令的条数。
  5. “占用率越高越好。” 占用率只是隐藏延迟的一种手段;很多最快的内核占用率很低,靠的是每个线程的指令级并行和大块复用。
  6. “块越大,算术强度越高,所以越快。” 块太大会让块数少于 SM 数,大部分 SM 闲着(波次量化)。
  7. “共享内存总是比全局内存快。” 有 bank 冲突时,共享内存的速度会降到 1/32;而且用共享内存要付同步的代价。
  8. “GPU 越快,推理就越快。” 小模型、小批的 decode 可能被 CPU 端的启动开销拖住,GPU 有一半时间在等。

12. 练习 · Exercises

带 💻 的题可以用 simt_lab 验证。参考答案在每题后面。

练习 1(网格大小)。n=1,000,000n = 1{,}000{,}000,每块 256 线程,要多少个块?最后一个块有几个线程是多余的?如果每块 1000 个线程呢?

参考答案

⌈106/256⌉=3907\lceil 10^6 / 256 \rceil = 3907 块,3907×256−106=1923907 \times 256 - 10^6 = 192 个多余线程。每块 1000 个线程:正好 1000 块,没有多余线程;但每块 1000 个线程 = 31 个满 warp + 1 个只有 8 个线程的 warp,每块都浪费了一个 warp 的大部分(SIMT 效率约 97.7%),不如取 1024。

练习 2(二维下标)。一个 1000×1000 的矩阵(按行存储),用 32×8 的块。网格的形状是什么?线程 (bx=3, by=5, tx=7, ty=2) 负责哪个元素(行、列、一维下标)?

参考答案

网格:x 方向 ⌈1000/32⌉=32\lceil 1000/32 \rceil = 32 个块,y 方向 ⌈1000/8⌉=125\lceil 1000/8 \rceil = 125 个块。行 =5×8+2=42= 5 \times 8 + 2 = 42,列 =3×32+7=103= 3 \times 32 + 7 = 103,一维下标 =42×1000+103=42,103= 42 \times 1000 + 103 = 42{,}103。

练习 3(合并访存)💻。一个 warp 的线程 ii 读 float32 数组的第 (i×3)(i \times 3) 个元素。碰到几个 sector?合并效率多少?用 02_memory_coalescing.py 的 pattern_kernel 加一行验证。

参考答案

地址 0,12,24,…,3720, 12, 24, \dots, 372 字节,跨越 00–375375 字节,即 12 个 sector(⌈376/32⌉=12\lceil 376 / 32 \rceil = 12)。需要 128 字节,搬 384 字节,效率约 33%。

练习 4(bank 冲突)💻。共享内存里一个 float32 的 tile,一个 warp 的 32 个线程读同一列的 32 行(tile[r][c],r=0..31r = 0..31)。tile 每行 64 个元素时几路冲突?每行 65 个呢?每行 66 个呢?用 03 脚本的 bank_kernel 验证。

参考答案

每行 64 个:相邻行相隔 64 个字,64 mod 32=064 \bmod 32 = 0,32 个线程全在同一个 bank,32 路冲突。每行 65 个:相隔 65 个字(奇数),32 个线程落在 32 个不同的 bank,无冲突。每行 66 个:相隔 66 个字,66 mod 32=266 \bmod 32 = 2,只用到 16 个偶数号 bank,第 rr 行和第 r+16r + 16 行撞在一起,2 路冲突。所以填充量要让行距和 32 互质,填 1 个最常见。

练习 5(同步)💻。在 transpose_tiled 里,如果把 blk.sync() 放到第一个循环里面(每装载 8 行就 sync 一次),结果、sync 次数、竞争数会怎样?这样做有什么坏处?

参考答案

结果仍然正确,竞争数为 0(最后一次 sync 之后才读 tile),sync 次数从每块 1 次变成 4 次。坏处是多了 3 次不必要的同步:每次同步都要等块里最慢的 warp,降低了 warp 之间互相填补空隙的机会。同步要放在”必须的地方”,不多也不少。

练习 6(归约的指令数)。为什么版本 2 和版本 3 的 warp 指令数完全一样(206,848),但估计的发射时间不同?

参考答案

两者”哪些线程在第几步干活”完全一样(都是前 ss 个线程),所以发射的 warp 指令数相同。区别在共享内存的地址:版本 2 访问 2*s*tid,有 bank 冲突,每条共享内存指令平均要 3.33 轮;版本 3 访问 tid 和 tid+s,没有冲突。多出来的轮数在估计模型里算作额外的发射时间(真实硬件上是共享内存流水线的重放)。

练习 7(占用率)💻。一个内核每块 128 线程、每线程 96 个寄存器、每块 32 KB 共享内存。用 occupancy() 算每 SM 块数、占用率和瓶颈。把共享内存降到 16 KB 会提高占用率吗?

参考答案

每 warp 寄存器 =⌈96×32/256⌉×256=3072= \lceil 96 \times 32 / 256 \rceil \times 256 = 3072,寄存器允许 ⌊65536/3072⌋=21\lfloor 65536 / 3072 \rfloor = 21 个 warp → 5 个块(每块 4 warp);共享内存每块 3333 KB,允许 ⌊233472/33792⌋=6\lfloor 233472 / 33792 \rfloor = 6 个块;所以每 SM 5 块、20 个 warp、占用率 31%,瓶颈是寄存器。共享内存降到 16 KB 不会提高占用率——寄存器仍然先卡住。优化前先找对瓶颈。

练习 8(Little 定律与向量化)。每个 SM 需要约 12.7 KB 在途。一个内核每 SM 只有 16 个活跃 warp,每个线程读 float32、每个 warp 一次只有 1 个请求。能拿到多少带宽?改成每个线程读 float4(16 字节)并展开循环让每个 warp 同时有 2 个请求呢?

参考答案

16×128=204816 \times 128 = 2048 字节,约 16%。改后 16×2×512=16,38416 \times 2 \times 512 = 16{,}384 字节 ≥ 12.7 KB,约 100%。不提高占用率,也能喂饱带宽。

练习 9(波次)💻。一个矩阵乘法 C[1024,1024]C[1024, 1024],块大小 128×128,每 SM 能放 2 个块。要几波?SM 利用率多少?换成 64×64 的块(每 SM 放 4 个)呢?

参考答案

128×128:64 个块,每波容量 264,1 波,利用率 64/264≈24%64 / 264 \approx 24\%。64×64:256 个块,每波容量 528,1 波,利用率约 48%。小矩阵要更小的块或 split-K 才能用满 GPU。

练习 10(启动开销)。一个 decode 步骤有 300 个 kernel,GPU 工作时间 2 ms,每次启动开销 5 µs。不重叠时每步多久?完全重叠时呢?用 CUDA Graphs 把每次启动开销降到 1 µs 后呢?

参考答案

不重叠:2+300×0.005=3.52 + 300 \times 0.005 = 3.5 ms;完全重叠:max⁡(2,1.5)=2\max(2, 1.5) = 2 ms。CUDA Graphs 后:不重叠 2.3 ms,完全重叠 2 ms。


13. 延伸阅读 · Further reading

  • NVIDIA, CUDA C++ Programming Guide:第 5 章(编程模型)、“Hardware Implementation”(SIMT 架构,也是第 8 讲的阅读材料)、“Performance Guidelines”(合并访存、bank、占用率)、“Compute Capabilities”(各代规格表)。
  • J. Austin et al., How to Scale Your Model, 第 12 章 How to Think About GPUs:从 roofline 角度看 SM、Tensor Core、存储层次和 NVLink。
  • M. Harris, Optimizing Parallel Reduction in CUDA(NVIDIA, 2007);An Efficient Matrix Transpose in CUDA C/C++、Using Shared Memory in CUDA C/C++(NVIDIA 技术博客)。
  • V. Volkov, Better Performance at Lower Occupancy(GTC 2010):用 ILP 代替占用率。
  • D. Kirk and W. Hwu, Programming Massively Parallel Processors(教材):系统地讲 CUDA 编程和并行模式(归约、扫描、矩阵乘法、卷积)。
  • 本仓库 CS349D 前置讲 P1:另一个角度的 GPU 组织、分块与复用、逐算子 roofline。

14. 和正课的联系 · Where this shows up

  • 第 1 讲:§3.9 的”分块与累加器”在本讲 §7.1 有了具体的约束(寄存器、共享内存、占用率);§4.2 的批处理在 GPU 上就是让矩阵乘法的块数和强度都足够大。
  • 第 8 讲(硬件引擎):SM、warp 调度器、Tensor Core、TMA(异步拷贝)、warp specialization 的硬件细节;TPU 的另一种设计。
  • 第 10 讲(内核开发):用本讲的全部概念从零写一个分块矩阵乘法:合并访存 → 共享内存分块 → 寄存器分块 → 双缓冲 → Tensor Core;原课作业 4。
  • 第 11 讲(内核优化):FlashAttention 是”分块 + 在线 softmax”,flash-decoding 是 split-K(§7.3),megakernel 是消灭启动开销(§8)。
  • 第 12–13 讲(编译器):算子融合(§8)、调度语言(怎样自动选择分块和线程映射)、Triton 的块级编程模型(§9)。
  • 第 14 讲(智能体生成内核):KernelBench 要检查的正是本讲的那些 bug:越界、竞争、结果不对但碰巧对。
← 返回课程目录