Lecture 10 · 内核开发：从朴素矩阵乘法到分块、流水与 Tensor Core · Kernel Development

对应 Maps to：Stanford CS349E / EE339 Efficient ML Inference at Scale（Autumn 2026），原课第 11 次课 · Tue Oct 27 · Kernel Development（Fredrik Kjolstad）；原课作业 4（Kernel Writing，10/19 – 10/26）。 依据说明：这次课没有指定阅读材料，课程网站也没有 slides。本讲以公开的经典材料为主线：S. Boehm 的博客 How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance: a Worklog（2022，下文简称 Boehm）的 kernel 1–10、NVIDIA CUDA C++ 编程指南、PTX ISA 文档（mma 的片段布局）、Triton 官方教程（矩阵乘法、融合 softmax）、GotoBLAS / BLIS 的论文。内容的取舍和顺序是我定的，不代表课堂上实际讲了什么。 配套代码 Code：code/lecture10/（numpy + 系统的 C 编译器；一个共享模块 + 6 个脚本，最长的约 11 秒） 先修：前置讲 P2（CUDA 编程模型、合并访存、共享内存与 bank、占用率）、第 8 讲（SM、Tensor Core、喂数据）、第 9 讲（分块与”访问次数”）。 预计阅读时间：3.5 小时（含运行代码）

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

0.1 学习目标 Learning goals

学完这一讲，你应该能够：

1. 说出什么时候值得自己写内核、一个 CUDA 内核从源码到在 SM 上运行经过哪些步骤。
2. 从最朴素的矩阵乘法内核出发，依次做合并访存、共享内存分块、一维 / 二维寄存器分块、向量化访存，并对每一步说清：它解除了哪个瓶颈、又撞上了哪个新瓶颈。
3. 数出一个内核每做 1000 次乘加要发多少条指令、碰多少个 cache line、共享内存要多少轮（bank 冲突），把这些数代入”多层屋顶”模型估计性能，并解释模型和实测的差距。
4. 解释寄存器分块的大小怎样由”共享内存带宽 ÷ 运算吞吐”决定，以及它和占用率的折中；为什么 A6000 和 A100 上的最优分块不同。
5. 读懂 Tensor Core 的 mma.m16n8k16 片段布局，解释 ldmatrix 的 bank 冲突和 XOR swizzle，算出喂饱 Tensor Core 需要的 warp 分块和流水线级数。
6. 用 Triton 风格的块级编程写同一个矩阵乘法，理解启动顺序对 L2 的影响和融合的收益；会给一个内核写测试（边缘形状、随 K 调整的容差）。

0.2 本讲地图 Map

§1  为什么要自己写内核
§2  一个内核的一生：源码 → PTX → SASS → 启动 → profiler
§3  目标：4096³ 的矩阵乘法；我们怎样在没有 GPU 的情况下"写内核"            ← 代码 01 演示 1
§4  kernel 1–2：朴素版与合并访存                                          ← 代码 01
§5  kernel 3：共享内存分块（同步、竞争、占用率、边界）                       ← 代码 02
§6  kernel 4–6：寄存器分块、向量化、bank 冲突；整条阶梯                     ← 代码 03
§7  延迟：Little 定律与双缓冲
§8  Tensor Core：mma 片段、ldmatrix 与 swizzle、喂数据的账、多级流水线       ← 代码 04
§9  同样的阶梯在 CPU 上：微内核与打包（真实测量）                           ← 代码 05
§10 Triton：块级编程、启动顺序与 L2、融合                                  ← 代码 06 演示 1–3
§11 测试、调试与 profiling                                                ← 代码 06 演示 4
§12 现状（2025–2026）
§13–17  小结、常见误区、练习（含答案）、延伸阅读、下一讲预告

0.3 本讲的代码

没有 GPU，我们分两步”写内核”：

• 正确性 + 事件计数：kernel_lab.py 在前置讲 P2 的 numpy GPU 模拟器（simt_lab）上加了向量访存（float4）、按相位计算的 bank 冲突、cache line 计数。六个矩阵乘法内核按 CUDA 的写法逐句翻译成 Python，在模拟器上真的执行、和 numpy 比较结果，同时数出每类指令和访存事件。
• 性能估计：把计数放大到 4096³，代入一个”多层屋顶”模型（FP32 运算、指令发射、L1、共享内存、显存），和 Boehm 在 RTX A6000 上的实测对照。

另外，05_cpu_microkernel.py 在你的 CPU 上真的测同样的阶梯（需要 C 编译器）；06_triton_style.py 用一个几十行的 mini-Triton 解释器演示块级编程。

1. 为什么要自己写内核 · Why write kernels at all

内核（kernel）是在 GPU 上运行的一个函数：CPU 发出一次”启动”，成千上万个线程一起执行它。PyTorch 里的每个算子（matmul、softmax、layer_norm……）背后都是一个或几个内核。

大多数时候你不需要自己写：矩阵乘法有 cuBLAS / cuBLASLt、卷积有 cuDNN、常见的融合有 CUTLASS 的模板和编译器（第 12–13 讲）。那为什么推理系统的工程师还在大量写内核？

1. 融合：两个算子之间的中间结果要写回显存再读出来。对受限于带宽的算子（逐元素运算、归一化、softmax），这一来一回就是全部的时间。把它们融合进一个内核，中间结果留在片上。FlashAttention（第 11 讲）是最著名的例子：注意力的 $N \times N$ 分数矩阵从来不写回显存。
2. 新算子：新的注意力变体（GQA、MLA、滑动窗口）、MoE 的路由和分组矩阵乘法、分页的 KV cache（第 7 讲）、投机解码的树形注意力（第 6 讲）——库里没有，或者有但不是为你的形状调的。
3. 新数据格式：INT4 / FP4 权重、块缩放（第 3 讲）。最快的写法是把”反量化”融合进矩阵乘法内核，读进来的是 4 位、乘的时候才展开。
4. 推理特有的形状：decode 时批很小，矩阵乘法退化成”细长”的矩阵 × 向量（第 1 讲），通用库的分块不一定适合。
5. 硬件的新特性：FP8、TMA、warp 专门化（第 8 讲）刚出来时，只有手写内核能用上。

写内核的工具也越来越多：CUDA C++（加上 CUTLASS / CuTe 模板库）、Triton（Python 写的块级内核，PyTorch 编译器的默认后端）、ThunderKittens（斯坦福 Hazy Research 的瓦片原语库）、JAX 的 Pallas、AWS 的 NKI（第 9 讲）、TileLang、Mojo 等。它们的抽象层次不同，但要解决的问题是同一个：让数据在离计算单元更近的地方被多用几次，并且让搬运和计算重叠。本讲用矩阵乘法把这件事从头走一遍。

2. 一个内核的一生 · The life of a kernel

你写的 .cu 文件
   │  nvcc：把主机代码交给 C++ 编译器；设备代码编译成
   ▼
PTX（虚拟指令集，类似"GPU 的汇编"，跨代兼容）
   │  ptxas：针对具体的架构（sm_80、sm_90……）分配寄存器、调度指令
   ▼
SASS（真正的机器码，cuobjdump -sass 可以看）  ──打包成 fatbinary，可以同时带 PTX 供将来的 GPU 即时编译
   │  运行时：kernel<<<grid, block, smem, stream>>>(参数)
   ▼
驱动把启动命令放进 stream 的队列 → GPU 的调度器按资源（寄存器、共享内存、线程数）把线程块分给各个 SM
   │  每个 SM 上的 warp 调度器每周期挑一个"准备好的" warp 发射一条指令（前置讲 P2 §4）
   ▼
所有块执行完 → 内核结束 → 同一个 stream 上的下一个内核才能开始

几个和性能直接相关的事实：

• 寄存器是编译时决定的。ptxas 决定每个线程用多少寄存器，这决定了一个 SM 能同时放几个块（占用率，§5.3）。寄存器不够时编译器会把变量溢出（spill）到”本地内存”（其实在显存里，经过 L1 缓存），非常慢。
• 启动有开销：每次几微秒（前置讲 P2 §8）。一个内核只做几微秒的活时，启动开销就是大头——这是第 11 讲 megakernel 的动机。
• 看 profiler，不要猜。NVIDIA Nsight Compute（ncu）会告诉你：计算单元和各级存储的利用率占峰值的百分比（“Speed of Light”）、每级存储搬了多少字节、warp 大部分时间在等什么（“stall reasons”：等全局内存叫 long scoreboard，等共享内存 / 特殊函数叫 short scoreboard 或 MIO throttle，等 __syncthreads 叫 barrier……）。Boehm 每一步都先看 profiler 再改代码，本讲引用的瓶颈判断就来自他的 profiler 截图。

我们没有 GPU，所以用模拟器代替”运行”、用计数器代替 profiler。

3. 目标：4096³ 的矩阵乘法 · The target

代码：01_naive_and_coalescing.py 演示 1

$C = A \cdot B$，$A$、$B$、$C$ 都是 4096×4096 的 FP32 矩阵，行主序存放。

强度上限远高于临界强度，所以这个矩阵乘法应该受限于计算（第 1 讲）：理想的内核把每个数从显存只读一次，然后在片上反复使用。下面所有的优化都是在逼近这个理想。

为什么选 A6000 和 FP32？因为 Boehm 的博客在这个设置下把每一步都测了，有一张完整的”阶梯”可以对照（表里的 GFLOP/s 是他的实测，矩阵大小 4092，下文按 4096 计算）：

（FP32 不走 Tensor Core。§8 再讲 Tensor Core 的内核——那里的原则一样，只是每一项数字都放大了一个数量级。）

我们的性能模型（kernel_lab.PerfModel）把一个 SM 每周期能处理的量列出来，乘以 SM 数和频率：

时间 = 最慢的那个资源（各资源完美重叠，roofline 的思想）。这是一个上限模型：它不模拟延迟和依赖，只数了访存、FFMA 和同步指令（没数地址计算和循环控制）。“L1 每周期处理一个 line”是根据 Nsight Compute 文档里 wavefront 的概念做的简化假设，不是 NVIDIA 的公开规格。用它的方式是：看每一步之后哪个资源最紧张，再和实测对照——实测接近模型，说明模型抓住了瓶颈；实测远低于模型，说明还有模型没看见的东西。

4. kernel 1–2：朴素版与合并访存 · Naive and coalesced

代码：01_naive_and_coalescing.py 演示 2–4

4.1 朴素版：一个线程算一个输出

最直接的写法：每个线程算 $C$ 的一个元素——$A$ 的一行和 $B$ 的一列的点积。

__global__ void sgemm_naive(int M, int N, int K, const float *A, const float *B, float *C) {
  const int x = blockIdx.x * blockDim.x + threadIdx.x;   // 行
  const int y = blockIdx.y * blockDim.y + threadIdx.y;   // 列
  if (x < M && y < N) {                                  // 边界检查（§5.4）
    float tmp = 0.0f;
    for (int i = 0; i < K; ++i) tmp += A[x * K + i] * B[i * N + y];
    C[x * N + y] = tmp;
  }
}
// 启动：dim3 grid(ceil_div(M, 32), ceil_div(N, 32)); dim3 block(32, 32);

代码是对的，但只有 309 GFLOP/s——cuBLAS 的 1.3%。问题出在一个 warp 的 32 个线程同时访问了什么地址（前置讲 P2 §5）。块是 32×32，threadIdx.x 变化最快，所以一个 warp 是 threadIdx.x = 0..31、threadIdx.y 相同的 32 个线程：它们的 x（行号）连续、y（列号）相同。

• 读 A[x*K + i]：32 个线程读 32 个不同的行，地址相差 $K \times 4 = 16$ KB，碰到 32 个不同的 cache line；
• 读 B[i*N + y]：32 个线程读同一个地址（广播），1 个 line。

模拟器数出来：朴素版每个 warp 每步碰 33 个 128 字节的 line、33 个 32 字节的 sector，每次乘加要经过 L1 的字节是 33 字节。

4.2 合并访存：只改一行映射

// 块改成一维 1024 个线程
const int row = blockIdx.x * 32 + threadIdx.x / 32;
const int col = blockIdx.y * 32 + threadIdx.x % 32;

现在一个 warp 是同一行、连续 32 列：读 A[row*K + i] 是同一个地址（广播，1 个 line），读 B[i*N + col] 是连续 32 个 float = 128 字节（1 个 line）。每步 2 个 line、5 个 sector，每次乘加经过 L1 的字节从 33 降到 5。两个内核做的乘加、读的数完全一样，只是”哪个线程算哪个元素”不同。

4.3 模型 vs 实测

• 朴素版：模型说瓶颈是 L1——每 32 次乘加要处理 33 个 line，L1 每周期只能处理一个。按这个算出约 290 GFLOP/s，和实测的 309 很接近。Boehm 的 profiler 显示朴素版的显存吞吐只有约 15 GB/s：慢的不是显存（显存那一列只要 11.5 ms），而是 L1 里一个个 line 地排队。
• 合并版：L1 的压力降了约 16 倍，模型的上限跳到 4.8 TFLOPS，实测只有 2.0（模型的 41%）。剩下的差距来自模型没算的东西：每次乘加仍要发两条 load 指令，每条要等几十个周期才回来，加上地址计算。这一步之后瓶颈已不在 L1；要再快，就得少发 load——把数据先搬到片上，让很多线程共用（§5）。

直觉：GPU 的内存系统以 32 字节的 sector、128 字节的 line 为单位工作。一个 warp 的 32 次访问如果落在一个 line 里，就是一次事务；落在 32 个 line 里，就是 32 次。合并访存是 GPU 编程的第一条规矩，而它只取决于”线程号怎样映射到地址”。

5. kernel 3：共享内存分块 · Shared-memory tiling

代码：02_shared_memory_tiling.py

5.1 想法：先搬进来，再用 32 遍

合并版里，同一个块的 1024 个线程会反复从全局内存读同样的数：同一行的 32 个线程都要 $A$ 的那一行，同一列的 32 个线程都要 $B$ 的那一列。这些重复的读虽然大多在 L1 / L2 命中，但每次都要发一条 load 指令、占一个 L1 的周期。

共享内存（shared memory）是每个 SM 上一块由程序员管理的 SRAM（A6000 每 SM 最多 100 KB），同一个块的线程都能读写，延迟和带宽接近 L1。把 $A$ 的 32×32 块和 $B$ 的 32×32 块先搬进共享内存，块里的线程就从共享内存读：

template <int BS>   // BS = 32
__global__ void sgemm_smem(int M, int N, int K, const float *A, const float *B, float *C) {
  __shared__ float As[BS * BS], Bs[BS * BS];
  const int tr = threadIdx.x / BS, tc = threadIdx.x % BS;
  A += blockIdx.x * BS * K;                      // 这个块要用的 A 的第一行
  B += blockIdx.y * BS;                          // 这个块要用的 B 的第一列
  C += blockIdx.x * BS * N + blockIdx.y * BS;
  float tmp = 0.0f;
  for (int bk = 0; bk < K; bk += BS) {
    As[tr * BS + tc] = A[tr * K + tc];           // 每个线程搬一个 A、一个 B（合并访存）
    Bs[tr * BS + tc] = B[tr * N + tc];
    __syncthreads();                             // ① 等大家都搬完
    A += BS;  B += BS * N;
    for (int d = 0; d < BS; ++d)
      tmp += As[tr * BS + d] * Bs[d * BS + tc];
    __syncthreads();                             // ② 等大家都用完，下一轮才能覆盖
  }
  C[tr * N + tc] = tmp;
}

这正是第 9 讲的分块：外层循环按 32 的块走 $K$，块在”片上缓冲”（共享内存）里被复用 32 次。

5.2 计数：瓶颈挪到了共享内存

全局访问少了 32 倍。可每次乘加仍要从共享内存读两个数（As[tr*BS + d] 一个 warp 读同一个地址——广播；Bs[d*BS + tc] 是连续 32 个字——没有 bank 冲突），共享内存每周期只能送 128 字节。模型说瓶颈挪到了共享内存，上限 4.7 TFLOPS；实测 3.0（模型的 64%）。Boehm 的 profiler 里最多的停顿是 “Stall MIO Throttle”——共享内存指令的队列满了，和模型一致。

块大小有用吗？ 代码 02 演示 2 把 BS 换成 8、16、32：L1 的 line 数 ∝ 1/BS（31.25 → 7.81 → 1.95 每 1000 次乘加），可共享内存的轮数三种都是 62.5。只要”一个线程只算一个输出”，每次乘加就要从共享内存拿两个操作数，这个天花板和块大小无关。要突破它，必须让寄存器也复用数据（§6）。

5.3 同步、竞争与占用率

两次 __syncthreads() 缺一不可（代码 02 演示 3）：

• 去掉第一次：有的 warp 还没把数写进共享内存，别的 warp 就开始读——写后读（RAW）竞争，模拟器报告 256 次；
• 去掉第二次：有的 warp 还在读这一轮的块，跑得快的 warp 已经开始写下一轮——读后写（WAR）竞争，报告 4 次。

有意思的是，模拟器里两种错误写法算出来的数都是对的：模拟器让所有线程严格同步地一条条执行。真实 GPU 上 warp 之间不同步，这样的程序可能这次对、下次错——测试通过不代表没有竞争。所以模拟器不看结果，而是记录”谁在哪个同步阶段读写过哪个位置”；真实 GPU 上对应的工具是 compute-sanitizer --tool racecheck。

占用率（occupancy）= 一个 SM 上同时驻留的 warp 数 / 上限（前置讲 P2 §7）。A6000 每 SM 最多 48 个 warp、16 个块、64K 个寄存器、100 KB 共享内存，四项里最紧的那项决定能放几个块：

1024 线程的块只放得下一个：67%，和 Boehm 的 profiler 一致，他也说占用率不是 kernel 3 的主要问题。BS = 16 能到 100%，可模型说它并不更快：瓶颈是共享内存的吞吐，不是延迟。占用率只在”在途的请求不够多、延迟藏不住”时才是问题（§7）。

5.4 边界

朴素版和合并版都有 if (x < M && y < N)：网格按 $\lceil M/32 \rceil \times \lceil N/32 \rceil$ 个块启动，总线程数比元素多，多出来的线程必须什么都不做。代码 01 演示 4 用 $100 \times 70 \times 50$ 测试：有边界检查时结果正确，边缘块里一部分线程被屏蔽（SIMT 效率 73%–78%）；去掉检查，模拟器立刻报告”线程 4（块 3,0）越界访问下标 5000，数组大小 5000”。真实 GPU 上，越界读可能悄悄读到别的数组、算出错误结果，也可能非法访问后崩溃——这是 LLM 写内核时最常见的错误之一（第 14 讲）。

k3 之后的内核和 Boehm 的原版一样假设形状是块的整数倍。生产代码的做法：(1) 在搬运时加边界检查、越界读成 0（Triton 的 mask，§10）；(2) 把输入补齐（padding）到块的倍数；(3) 为”主体”和”边缘”各写一个内核。

6. kernel 4–6：寄存器分块与向量化 · Register tiling and vectorization

代码：03_register_tiling.py

6.1 一维分块：一个线程算一列 8 个输出

想法：让一个线程算多个输出。从共享内存读进寄存器的一个数，被这个线程用好几次。kernel 4 的块算 64×64 的 $C$，512 个线程，每个线程算一列里连续 $TM = 8$ 个元素：

float res[TM] = {0.0f};                               // 8 个累加器，在寄存器里
for (int bk = 0; bk < K; bk += BK) {                  // BK = 8
  As[iAr * BK + iAc] = A[iAr * K + iAc];              // 64×8 的 A 块、8×64 的 B 块，每个线程各搬一个数
  Bs[iBr * BN + iBc] = B[iBr * N + iBc];
  __syncthreads();
  A += BK;  B += BK * N;
  for (int d = 0; d < BK; ++d) {
    float tb = Bs[d * BN + tcol];                     // 读 1 个 B……
    for (int r = 0; r < TM; ++r)
      res[r] += As[(trow * TM + r) * BK + d] * tb;    // ……和 8 个 A 各乘一次
  }
  __syncthreads();
}

每次乘加的共享内存读从 2 次降到 $(1 + 8)/8 = 1.125$ 次（Boehm 的说法是”9K/8 次共享内存访问 / 每个结果”）。LDS 指令从每 1000 次乘加 62.5 条降到 35.2 条，模型的上限 8.4 TFLOPS，实测 8.5——几乎完全一致。

6.2 二维分块：外积

再进一步：每个线程算一个 $TM \times TN = 8 \times 8$ 的输出块。每个 $k$ 步从共享内存读 8 个 $A$（regM）、8 个 $B$（regN），做 64 次乘加——一个外积：

float res[TM * TN] = {0.0f}, regM[TM], regN[TN];      // 64 个累加器 + 16 个操作数
...
for (int d = 0; d < BK; ++d) {
  for (int i = 0; i < TM; ++i) regM[i] = As[(trow * TM + i) * BK + d];
  for (int j = 0; j < TN; ++j) regN[j] = Bs[d * BN + tcol * TN + j];
  for (int i = 0; i < TM; ++i)
    for (int j = 0; j < TN; ++j) res[i * TN + j] += regM[i] * regN[j];
}

每次乘加的共享内存读 = $(TM + TN)/(TM \cdot TN) = 16/64 = 0.25$ 次。块变成 128×128、256 个线程，从全局内存读的数也少了一半（每 1000 次乘加的全局 load 从 0.98 条降到 0.49 条：每个搬进来的数参与 128 次乘加）。

6.3 寄存器的代价：占用率

64 个累加器 + 16 个操作数 + 地址，每个线程要 100 多个寄存器。按每线程 104 个估计，256 线程的块每 SM 只放得下 2 个，占用率 33%。这就是寄存器分块的折中：块越大，共享内存的压力越小，但每线程寄存器越多、占用率越低、延迟越难藏。块再大（16×8），寄存器 176 个，占用率只剩 17%。

分块多大才够？（代码 03 演示 4）每次乘加要从共享内存读 $4(TM + TN)/(TM \cdot TN)$ 字节，硬件每次乘加能供应的是”共享内存 128 字节/周期 ÷ 每周期的乘加数”：

A6000（GA102）每 SM 有 128 个 FP32 单元，要 8×8 才喂得饱；A100（GA100）每 SM 只有 64 个，4×4 就够了。这正好解释了 Boehm 的自动调优结果：A6000 上最好的是 $TM = TN = 8$（BM = BN = 128、BK = 16），A100 上最好的是 $TM = TN = 4$（BM = BN = 64、BK = 16）——在 A100 上用更大的分块没有好处，只会多占寄存器。最优点随硬件变，这就是要自动调优的原因（第 13 讲）。

6.4 bank 冲突

模拟器数出 k5 的 LDS 指令只有 7.8 条/1000 次乘加，可共享内存的轮数是 23.4——比指令数多 3 倍。原因是 bank 冲突：共享内存分 32 个 bank，每个 bank 每周期送一个 4 字节的字；一个 warp 的访问里，同一个 bank 上有 $k$ 个不同的字，就要 $k$ 轮（同一个字被多个线程读只算一次——广播）。

两种常用的解法：padding——每行多留一个字，让相邻行错开一个 bank（浪费一点空间）；改映射——让同一相位的线程读连续地址。Boehm 的 kernel 10（warp 分块）就是后者：把每个线程的 8 列拆成两组 4 列（tcol·4 和 64 + tcol·4），LDS.128 就没有冲突了。

向量访问的”相位”：一条 32 位的 LDS，整个 warp 一个相位；64 位分两个相位（每个 16 个线程）；128 位（float4）分四个相位（每个 8 个线程），冲突只在相位内部计。所以 LDS.128 即使没有冲突也要 4 轮（32 个线程 × 16 字节 = 512 字节，每轮 128 字节）。向量化省的是指令，不是共享内存的带宽。（这是 NVIDIA 论坛和 profiler 文档里常见的描述，NVIDIA 没有给出完整的官方规格；kernel_lab.smem_wavefronts 按这个模型计算。）

6.5 kernel 6：向量化访存

kernel 6 的分块和 kernel 5 一样，改了三处：

float4 t = reinterpret_cast<const float4 *>(&A[iAr * K + iAc * 4])[0];   // LDG.128：一条指令读 16 字节
As[(iAc * 4 + 0) * BM + iAr] = t.x;                                       // A 在共享内存里转置存放（As[k][m]）
As[(iAc * 4 + 1) * BM + iAr] = t.y;
As[(iAc * 4 + 2) * BM + iAr] = t.z;
As[(iAc * 4 + 3) * BM + iAr] = t.w;
reinterpret_cast<float4 *>(&Bs[iBr * BN + iBc * 4])[0] =
    reinterpret_cast<const float4 *>(&B[iBr * N + iBc * 4])[0];            // B 直接 128 位拷贝
// 内层：regM 的 8 个数现在是连续的（As[d][trow*8 .. +7]），和 regN 一样用两条 LDS.128 读

为什么转置 $A$？k5 的 regM[i] = As[(trow*TM + i)*BK + d] 的 8 个数相隔 BK 个字，不连续，不能合并成一条 LDS.128；转置之后变成 As[d*BM + trow*TM + i]，连续了。

6.6 整条阶梯

每 1000 次乘加的账本（模拟器在 128³ 上数，256³ 上验证——主循环里的计数完全相同）：

FFMA 一直是 31.25 条（一条指令 32 个线程各做一次乘加）；优化的是 FFMA 之外的那些指令和轮数。k6 的 warp 指令总数已经接近 FFMA 本身。

代入模型（A6000，4096³）：

三点值得细看：

1. 模型对瓶颈的判断和 Boehm 的 profiler 一致：k1 卡在不合并的访问，k3 卡在共享内存（MIO throttle），之后共享内存的压力逐步下降。k1、k4 的实测和模型几乎一样。
2. k5 的模型低于实测。一个吞吐上限模型不应该比实测还慢——说明某个假设错了。最可能的是：regN[j] = Bs[d*BN + tcol*TN + j] 的 8 个数是连续的、对齐也能被证明，nvcc 把它们合并成了两条 LDS.128（我们的模拟器是逐句翻译源码的）。按这个假设，模型是 18.8 TFLOPS，实测是它的 85%。我们看不到 SASS，这只是推断——在真机上用 cuobjdump -sass 一看便知。这是用模型的正确姿势：模型和实测矛盾时，去找哪个假设错了。
3. k6 之后剩下的差距（18.2 → cuBLAS 23.2）来自：Bs 的 2 路冲突、没有双缓冲（§7）、以及 warp 分块带来的更好的寄存器复用。Boehm 的 kernel 10 到了 cuBLAS 的 94%。

warp 分块（kernel 10）多加了一层：块 → warp → 线程。每个 warp 负责块里的一个子块（比如 64×32），warp 里的线程再各算一小块。好处是：一个 warp 读的共享内存地址更集中（容易避开 bank 冲突），而且这一层正好对应 Tensor Core 的”一个 warp 一条 mma”（§8）。

7. 延迟：Little 定律与双缓冲 · Latency, Little’s law and double buffering

上面的模型只看吞吐。真实的 GPU 上，一条全局 load 要几百个周期才回来，共享内存几十个周期。Little 定律：要让某个资源一直忙，在途的请求数 = 吞吐 × 延迟（前置讲 P1、第 8 讲 §8.2）。

GPU 用两种办法凑够在途的请求：

• 线程级并行（TLP）：很多 warp 轮流执行，一个 warp 在等，就换另一个（占用率）；
• 指令级并行（ILP）：一个线程里有很多互不依赖的指令——比如 k5 的 64 个累加器，64 条 FFMA 之间没有依赖，可以连着发。

寄存器分块降低了占用率，但增加了 ILP；所以 k5 只有 33% 的占用率仍然很快。代码 05 演示 1 在 CPU 上演示了同一件事：一个核上只有 1 条 FMA 依赖链时，每条都要等上一条的结果；链越多越快，约 16 条时到顶（§9）。

双缓冲（double buffering）解决的是另一种等待：k3–k6 每一轮都是”搬一块 → __syncthreads → 算 → __syncthreads”，搬的时候所有线程都在等。双缓冲在共享内存里放两份缓冲：算第 $t$ 块的同时，异步地搬第 $t+1$ 块。

// 伪代码：Ampere 的 cp.async（全局 → 共享内存的异步拷贝，不经过寄存器）
issue_async_copy(stage[0], tile 0);  commit();
for (int t = 0; t < n_tiles; ++t) {
  if (t + 1 < n_tiles) { issue_async_copy(stage[(t + 1) % 2], tile t + 1);  commit(); }
  wait_until_at_most_1_pending();  __syncthreads();      // 第 t 块到了（允许第 t+1 块还在路上）
  compute(stage[t % 2]);
  __syncthreads();                                       // 大家用完第 t 块，它的缓冲区才能被第 t+2 块覆盖
}

（真实的 API 是 PTX 的 cp.async / cp.async.commit_group / cp.async.wait_group，或 CUDA 的 cuda::memcpy_async 和 cuda::pipeline。）级数可以多于 2——这就是多级流水线，§8.4 计算 Tensor Core 内核要几级。第 8 讲 §8.4 的流水线模拟器给出了”级数 vs 计算单元忙碌比例”的曲线。

8. Tensor Core 内核 · Tensor Core kernels

代码：04_tensor_cores.py

8.1 一条 mma 指令：warp 的 32 个线程合作

Tensor Core 不是”更快的 FFMA”，而是一种warp 级的指令。Ampere 的 mma.sync.aligned.m16n8k16（FP16 输入、FP32 累加）：一个 warp 的 32 个线程合作完成 $D_{16 \times 8} = A_{16 \times 16} \cdot B_{16 \times 8} + C_{16 \times 8}$，2048 次乘加——一条 FFMA 只有 32 次。

每个线程只拿着输入和输出的一小片（片段，fragment），哪个线程拿哪些元素是 PTX ISA 规定死的。记 groupID = lane >> 2、t = lane % 4：

规律：同一行相邻的两个 FP16 放在一个 32 位寄存器里（所以列号成对出现）；4 个相邻 lane 按列分，8 组按行分。$C$ 的布局和 $A$ 的左半边相同——上一条 mma 的输出可以直接当下一条的输入，FlashAttention 里 $P = \mathrm{softmax}(S)$ 接着乘 $V$ 就用到这一点（第 11 讲）。代码 04 演示 1 按这个布局把随机矩阵拆成片段、模拟 mma、再拼回来，100 条 mma 和 FP64 的结果最大相差 $2.3 \times 10^{-6}$（FP16 × FP16 的乘积在 FP32 里是精确的，误差只来自累加）。

CUDA 里最简单的写法是 WMMA API（片段的布局对程序员不透明）：

#include <mma.h>
using namespace nvcuda;
wmma::fragment<wmma::matrix_a, 16, 16, 16, half, wmma::row_major> a;
wmma::fragment<wmma::matrix_b, 16, 16, 16, half, wmma::row_major> b;
wmma::fragment<wmma::accumulator, 16, 16, 16, float> acc;
wmma::fill_fragment(acc, 0.0f);
for (int k = 0; k < K; k += 16) {
  wmma::load_matrix_sync(a, A + row * K + k, K);     // 一个 warp 一起读一个 16×16 的块
  wmma::load_matrix_sync(b, B + k * N + col, N);
  wmma::mma_sync(acc, a, b, acc);                    // 一个 warp 一起算
}
wmma::store_matrix_sync(C + row * N + col, acc, N, wmma::mem_row_major);

高性能的内核（CUTLASS、FlashAttention）直接用 PTX 的 mma 和 ldmatrix，自己控制片段的布局和寄存器的复用。

8.2 ldmatrix 与 swizzle

数据先从全局内存搬进共享内存，再用 ldmatrix.x4 读成片段：一个 warp 读 4 个 8×8 的 FP16 小矩阵，每个 lane 提供一行的地址（16 字节）。每个 8×8 小矩阵是一个相位：8 行 × 16 字节 = 128 字节，正好铺满 32 个 bank——前提是这 8 行落在不同的 bank 上。

如果共享内存里的数据块按行主序存放、每行 $BK = 64$ 个 FP16（128 字节），同一列的 8 行在共享内存里相隔 128 字节，全都从 bank 0 开始：8 路冲突。代码 04 演示 2 在一个块级的 Tensor Core 矩阵乘法（128³、每块 64×64、4 个 warp、FP16 输入 / FP32 累加）里数出来：

XOR swizzle：把地址看成一串 128 字节的 line，每个 line 8 个 16 字节的 chunk；物理 chunk 号 = 逻辑 chunk 号 XOR（line 号 mod 8）。同一列的第 0–7 行于是落在 bank 0、4、8、……、28 开始的 8 个 chunk 上，冲突消失，而且不浪费一个字节（padding 也能解冲突，但会让整块拷贝变复杂）。swizzle 只改变物理地址，写和读都经过同一个映射，结果完全一样（表里四种情况的误差都是 $5.3 \times 10^{-6}$）。CUTLASS 的 CuTe 把它写成 Swizzle<3,4,3>；Triton 编译器自动选择共享内存的 swizzle 布局；Hopper 的 TMA 硬件直接支持几种 swizzle 模式。

8.3 喂数据的账

Tensor Core 把每 SM 每周期的乘加从 64–128 次提高到 512–2048 次（第 8 讲 §7），共享内存的带宽没变（128 字节/周期）：

一个 warp 做 $k = 16$ 的一步、算 $W_M \times W_N$ 的输出块，复用 = $W_M W_N / (2(W_M + W_N))$：

H100 上一个 warp 至少要算 64×64，每个线程 128 个累加器——寄存器就快用完了。Hopper 的回答是 WGMMA：4 个 warp（一个 warpgroup）合作一条指令，$B$ 直接从共享内存读、不经过寄存器，块可以大到 64×256，累加器分摊到 128 个线程上（第 8 讲 §7.2）。

再往外一层，块的大小决定显存流量。4096³ FP16 在 H100 上（波次模型）：64×64 的块读显存要约 340 µs，比按 989 TFLOPS 计算的约 140 µs 还长——整个内核会受限于带宽；128×128 降到约 100 µs，128×256 约 60 µs。所以 Hopper 上的 GEMM 和 FlashAttention-3 都用很大的块。

8.4 多级流水线要几级

一级流水线 = 一个 $B_M \times B_N \times B_K$ 的输入块。按 Little 定律，除了正在算的那一级，还要有足够多的级”在路上”：

块越大，每级算得越久、需要的级数越少，但每级占的共享内存越多。实际的 Hopper GEMM 一般 3–5 级。每 SM 需要的输入带宽约 90–120 GB/s，132 个 SM 合起来远超显存的 3.35 TB/s——这些数据大部分来自 L2（同一个 $A$ 块被同一行的多个块共用）。

Hopper 的写法（第 8 讲 §8.3–§8.5）：TMA（一个线程发一个”张量描述符”，硬件搬整块、顺便做 swizzle）+ mbarrier（按到达的字节数计数，取代 __syncthreads）+ warp 专门化（一个 producer warp 只发 TMA，consumer warpgroup 只做 WGMMA）。CUTLASS 3 和 FlashAttention-3 就是这样写的。

9. 同样的阶梯在 CPU 上 · The same ladder on a CPU

代码：05_cpu_microkernel.py（真实测量；下面是作者机器 Apple M4 Pro、Apple clang 21 上一次运行的结果，你的数字会不同）

GPU 内核的每一步在 CPU 的高性能矩阵乘法（GotoBLAS、BLIS、OpenBLAS）里都有对应。这一节是真的测的。

Little 定律在一个核上：只有 1 条 FMA 依赖链时约 9 GFLOP/s；链越多越快，约 16 条时到顶（约 121 GFLOP/s，单核 FP32 向量的实测峰值）。和 GPU 上”足够多的 warp 或足够多的独立累加器”是同一件事。

阶梯（单核，FP32）：

• 寄存器微内核：每次算 $C$ 的 $M_R \times N_R$ 小块，累加器留在寄存器里走完整个 $k$——和 k5 的”每线程 8×8 个累加器”一模一样。比 ikj 快 3 倍，到了实测峰值的约 84%。
• 打包（packing）：n = 1024 时只多约 11%；n = 2048 时不打包的微内核掉到 51 GFLOP/s，打包后仍是 111。不打包时微内核读的 $B$ 是一条 16 列宽、$n$ 行高的竖条，相邻两行隔 $4n$ 字节，矩阵一大，这些散开的 cache line 和页就装不下了；打包把它们拷成一段连续内存。这就是 GPU 上”先搬进共享内存、按计算的顺序重排”（k6 转置 $A$）的 CPU 版本。
• 寄存器分块的大小：4×16 最快（101），8×8 次之（84）；8×16 和 4×32 要 32 个累加器寄存器，加上 $B$ 和 $a$ 超过了 NEON 的 32 个向量寄存器，编译器溢出到内存，掉到 42–51。GPU 上同样的约束表现为”占用率下降”，CPU 上是直接溢出。
• 厂商库：numpy 调用的 Accelerate 单线程就有约 1.7 TFLOPS，是我们向量微内核的 15 倍。据公开的分析，它在 M 系列芯片上使用专门的矩阵协处理器（AMX / SME）——和 Tensor Core 对 CUDA core 的关系一样：同样的算法，换成专用的矩阵单元（第 8 讲）。

10. Triton：块级编程 · Block-level programming with Triton

代码：06_triton_style.py 演示 1–3

10.1 写”一个块做什么”

CUDA 里你写一个线程做什么，再自己安排线程怎样合作（共享内存、同步、寄存器分块、swizzle、流水线）。Triton（Tillet、Kung、Cox，MAPL 2019；OpenAI 维护）换了一个抽象层次：你写一个块（Triton 叫一个”程序实例”）做什么，操作的是整块的张量：

@triton.jit
def matmul_kernel(a_ptr, b_ptr, c_ptr, M, N, K, stride_am, stride_ak, stride_bk, stride_bn, stride_cm, stride_cn,
                  BLOCK_M: tl.constexpr, BLOCK_N: tl.constexpr, BLOCK_K: tl.constexpr, GROUP_M: tl.constexpr):
    pid = tl.program_id(0)
    num_pid_m, num_pid_n = tl.cdiv(M, BLOCK_M), tl.cdiv(N, BLOCK_N)
    group_id = pid // (GROUP_M * num_pid_n)                        # 分组顺序（§10.2）
    first_m = group_id * GROUP_M
    group_size_m = min(num_pid_m - first_m, GROUP_M)
    pid_m = first_m + (pid % (GROUP_M * num_pid_n)) % group_size_m
    pid_n = (pid % (GROUP_M * num_pid_n)) // group_size_m

    offs_m = pid_m * BLOCK_M + tl.arange(0, BLOCK_M)               # 这个块负责的行号（一个向量）
    offs_n = pid_n * BLOCK_N + tl.arange(0, BLOCK_N)
    offs_k = tl.arange(0, BLOCK_K)
    a_ptrs = a_ptr + offs_m[:, None] * stride_am + offs_k[None, :] * stride_ak    # 一整块的指针
    b_ptrs = b_ptr + offs_k[:, None] * stride_bk + offs_n[None, :] * stride_bn
    acc = tl.zeros((BLOCK_M, BLOCK_N), dtype=tl.float32)
    for k in range(0, tl.cdiv(K, BLOCK_K)):
        a = tl.load(a_ptrs, mask=(offs_m[:, None] < M) & (offs_k[None, :] < K - k * BLOCK_K), other=0.0)
        b = tl.load(b_ptrs, mask=(offs_k[:, None] < K - k * BLOCK_K) & (offs_n[None, :] < N), other=0.0)
        acc = tl.dot(a, b, acc)                                    # 块乘法：编译器会用 Tensor Core
        a_ptrs += BLOCK_K * stride_ak
        b_ptrs += BLOCK_K * stride_bk
    c_ptrs = c_ptr + offs_m[:, None] * stride_cm + offs_n[None, :] * stride_cn
    tl.store(c_ptrs, acc.to(tl.float16), mask=(offs_m[:, None] < M) & (offs_n[None, :] < N))

（结构和官方教程 Matrix Multiplication 一样；教程里 $M$、$N$ 方向用取模代替 mask，这里统一写成 mask。）

没有 threadIdx、没有 __syncthreads、没有共享内存。这些交给编译器：它决定块内的线程怎样分工、哪些数据经过共享内存、用什么 swizzle、tl.dot 用哪种 Tensor Core 指令、流水几级（num_stages）、每个程序用几个 warp（num_warps）。程序员保留的是分块（BLOCK_M/N/K）和启动顺序——恰好是第 9 讲说的”最重要的那部分”。Triton 论文的主张是这样写出的矩阵乘法和卷积能接近厂商库；PyTorch 2 的编译器（第 12 讲）默认生成 Triton 内核。

代码 06 用一个几十行的 mini-Triton 解释器（kernel_lab.TL）逐个执行程序实例：$300 \times 200 \times 150$、$1 \times 1000 \times 64$、$129 \times 65 \times 33$ 这些不规整的形状都正确——mask 让越界的位置不读（当作 0）、不写。块越大，每读一个字节做的运算越多：512³ 时从 16×16 块的 4 FLOP/字节到 128×256 的约 43，总读取 = $M K \cdot N/B_N + K N \cdot M/B_M$。

10.2 启动顺序与 L2

每个程序实例读 $A$ 的 $B_M$ 行、$B$ 的 $B_N$ 列。同时在跑的那一波程序实例如果碰到的是同一些行和列，读一次就能在 L2 里被大家共用。Triton 教程用一个 9×9 的例子说明：

代码 06 演示 2 用本讲的波次模型估计 4096³ FP16 的显存流量：128×128 的块、GROUP_M = 8 时降到行主序的 69%（A100）–77%（H100）；128×256 的块只有 GROUP_M = 16 时降到 87%–89%；GROUP_M = 4 反而比行主序多（一波的块跨过两个组时，碰到的块行和块列都变多）。这个优化不碰内核的计算部分，只改”程序号 → 块坐标”的映射（CUDA 里叫 threadblock swizzle）。是否有用取决于问题大小和 L2：Boehm 在 kernel 10 上试过，L2 命中率已经约 80%，没有收益。

10.3 融合：softmax

Triton 教程的第二个例子是融合的 softmax。逐个算子地算（PyTorch eager 的写法）：求每行最大值、减、exp、求和、除，5 个内核，每个都把 $M \times N$ 的中间结果写回显存、下一个再读回来。代码 06 演示 3 的计数正好是教程给的公式：

softmax 每个元素只有几次运算，完全受限于带宽，时间正比于流量——教程据此估计融合有约 4 倍的理论加速。融合的前提是一行放得进一个程序实例的片上存储；一行太长时，就要用”在线 softmax”分块地算——那正是 FlashAttention（第 11 讲）。

11. 测试、调试与 profiling · Testing, debugging and profiling

代码：06_triton_style.py 演示 4

一个只在边缘出错的 bug：把矩阵乘法内核里 $K$ 方向的 mask 去掉。在 $64^3$、$128 \times 64 \times 96$ 这些”维度都是 32 的倍数”的形状上，它全部通过——最后一个 $K$ 块正好不越界；只有 $K$ 不是 32 的倍数时才出错（读到了下一行的数，误差 13–66）。所以测试集一定要有边缘形状：1、质数、比块大小多一点或少一点。

容差要随 $K$ 放大：$K$ 个乘积相加，累加的舍入误差随 $K$ 增长。用固定的 $10^{-6}$ 去比，大 $K$ 的正确内核也会”失败”；用太松的容差，又会放过真正的 bug。KernelBench 用 atol = rtol = $10^{-2}$（第 14 讲会讨论这个选择的后果）。

累加精度：FP16 输入、64×64 的输出，$K$ 从 64 到 4096：

FP16 累加的误差是 FP32 的约 8000 倍，差不多就是两种格式机器精度之比 $2^{13}$（前置讲 P3 §4.5）。Tensor Core 的 GEMM 默认 FP32 累加就是这个原因。

测试清单：

1. 和高精度的参考实现比较（FP64 的 numpy / PyTorch）；
2. 随机形状 + 边缘形状；不同的数据分布（全 0、很大的数、有离群值）；
3. 随 $K$ 调整的容差；
4. 多次运行：竞争导致的错误不一定每次出现（§5.3）；
5. compute-sanitizer：memcheck 查越界，racecheck 查共享内存竞争，synccheck 查同步的误用。

profiling 的顺序（Nsight Compute）：先看 “Speed of Light”——计算和存储哪个更接近峰值？再看 Memory Workload Analysis——每级存储搬了多少，L1 / L2 命中率，共享内存的 bank 冲突（“shared store/load bank conflicts”）；再看 Warp State Statistics——warp 大部分时间在等什么；最后看 Source Counters——停顿集中在哪几行代码。Boehm 每一步都按这个顺序决定下一步改什么。

12. 现状（2025–2026）· State of the field

• 一条指令越来越大：Hopper 的 WGMMA 由一个 warpgroup 发出；Blackwell 的第五代 Tensor Core（tcgen05 指令）由单个线程发出，累加器放在专门的 Tensor Memory 里而不是寄存器，还能让两个 SM 合作一条 MMA。写内核的重心进一步从”怎样算”转向”怎样搬”。
• 抽象层次的分化：一端是 CUTLASS / CuTe（C++ 模板，2025 年起也有 Python 的 CuTe DSL），把布局（layout）和 swizzle 写成代数对象；一端是 Triton（以及 2025 年起暴露更底层控制的 Gluon）、ThunderKittens、TileLang 这类”瓦片级”语言；JAX 的 Pallas、AWS 的 NKI 面向各自的硬件。FlashAttention 的新版本、DeepSeek 开源的 FlashMLA / DeepGEMM 都是手写或半手写的内核。
• 自动调优和编译器：同一个内核的最优分块随硬件和形状变（§6.3），Triton 的 @triton.autotune、cuBLASLt 的启发式和 PyTorch 的 max-autotune 都在搜索这个空间（第 13 讲）。
• LLM 写内核：KernelBench 之类的基准和”智能体 + 编译器 + profiler”的循环（第 14–15 讲）。本讲 §5.4、§11 的那些 bug——越界、竞争、只在边缘形状出错——正是它们最常犯、评测最需要抓住的。

13. 小结 · Summary

1. 自己写内核的理由：融合、新算子、新数据格式、推理特有的形状、新硬件特性。所有优化做的都是同一件事——让数据在离计算单元更近的地方被多用几次，并让搬运和计算重叠。
2. 合并访存：只改”线程号 → 地址”的映射，一个 warp 碰的 cache line 从 33 个变成 2 个，快 6 倍。
3. 共享内存分块：全局访问少 32 倍，但每次乘加仍要从共享内存读 2 个数——瓶颈挪到共享内存，这个天花板和块大小无关。
4. 寄存器分块：每个线程算 $TM \times TN$ 个输出（外积），每次乘加的共享内存访问 = $(TM + TN)/(TM \cdot TN)$。要多大由”共享内存带宽 ÷ 运算吞吐”决定（A6000 要 8×8，A100 4×4 就够），代价是寄存器和占用率。
5. bank 冲突和向量化：冲突让轮数比指令数多几倍，用 padding、改映射或 swizzle 消除；LDS.128 省指令不省带宽。
6. 一个简单的多层屋顶模型对每一步瓶颈的判断和 profiler 一致；模型和实测的矛盾（k5）指向错误的假设。
7. 延迟靠 TLP（占用率）和 ILP（独立的累加器）凑够在途的请求；双缓冲 / 多级流水线让搬运和计算重叠。
8. Tensor Core：warp 级的 mma，片段布局由指令规定；ldmatrix 需要 swizzle 避免冲突；H100 要 64×64 以上的 warp 块，所以有了 WGMMA、TMA、warp 专门化和 3–5 级流水线。
9. CPU 上是同一个阶梯：寄存器微内核 + 打包 + 缓存分块；厂商库靠专用的矩阵单元再快一个数量级。
10. Triton 让你只写”一个块做什么”：分块和启动顺序留给程序员，其余交给编译器。测试要有边缘形状和随 $K$ 调整的容差。

14. 常见误区 · Common pitfalls

1. “合并访存就是让每个线程读连续的地址”——要连续的是同一个 warp 里相邻线程同时访问的地址；单个线程沿着自己的循环连续读，对合并没有帮助（k1 的每个线程都在连续读 $A$ 的一行，§4.1）。
2. “用了共享内存就快了”——如果每次运算仍要从共享内存取两个操作数，瓶颈只是从 L1 挪到了共享内存（§5.2）。
3. “占用率越高越好”——占用率是藏延迟的手段之一；寄存器分块用更低的占用率换来了更少的共享内存访问和更多的 ILP，反而更快（§6.3、§7）。
4. “float4 让共享内存的带宽翻 4 倍”——LDS.128 分四个相位，至少 4 轮；它省的是指令（§6.4）。
5. “模型算出来比实测慢，说明硬件比想象的好”——吞吐上限模型不该比实测慢；这说明模型的某个假设（比如编译器生成的指令）错了（§6.6）。
6. “测试都通过了，内核就是对的”——没有边缘形状的测试抓不住 mask 错误（§11）；模拟器 / 某次运行结果正确也不代表没有竞争（§5.3）。
7. “Tensor Core 只是更快的乘法器”——它是 warp 级的指令，数据必须按规定的片段布局分布在线程之间；喂不饱它的内核拿不到它的峰值（§8）。
8. “FP16 输入就可以 FP16 累加”——累加误差大几千倍（§11）。
9. “写好一个内核就能用在所有 GPU 上”——最优的分块随每 SM 的运算吞吐、寄存器数、共享内存大小而变（§6.3）；要么为每种硬件调，要么交给自动调优。

15. 练习 · Exercises

练习 1（合并访存） 朴素内核的块改成 dim3 block(32, 32)，但把 x、y 的计算对调（x 用 threadIdx.y、y 用 threadIdx.x）。一个 warp 每步碰几个 cache line？

答案：一个 warp 现在是同一行（x 相同）、连续 32 列（y 连续）：读 A[x*K + i] 是同一个地址（1 个 line），读 B[i*N + y] 是连续 128 字节（1 个 line），共 2 个——和 kernel 2 一样。kernel 2 只是用一维的块做了同样的事。

练习 2（共享内存分块的全局访问） kernel 3 的块大小 BS = 16（256 个线程）。每 1000 次乘加碰多少个 L1 line？

答案：一个 warp 是 tile 的 2 行 × 16 列：读 $A$ 的块碰 2 个 line（每行 64 字节）、读 $B$ 也是 2 个，每块 4 个；这个 warp 每块做 $32 \times 16 = 512$ 次乘加。$4 / 512 \times 1000 = 7.81$（代码 02 演示 2 的数）。

练习 3（寄存器分块的大小） 一块 GPU 每 SM 每周期 96 次 FP32 乘加，共享内存每周期 128 字节。$TM \times TN$ 至少多大（取 $TM = TN$）？

答案：每次乘加能分到 $128/96 = 1.33$ 字节；需要 $4 \cdot 2T / T^2 = 8/T \le 1.33$，$T \ge 6$。取 2 的幂就是 8×8（或者非方的 8×4：$4 \cdot 12/32 = 1.5$，不够）。

练习 4（占用率） 256 个线程的块，每线程 104 个寄存器，A6000（每 SM 64K 个寄存器、48 个 warp，寄存器按每 warp 256 个的粒度分配）。每 SM 能放几个块？占用率多少？

答案：每 warp $\lceil 104 \times 32 / 256 \rceil \times 256 = 3328$ 个寄存器，每块 8 个 warp = 26,624 个，$65536 / 26624 = 2.46$ → 2 个块、16 个 warp，占用率 $16/48 = 33\%$（代码 03 演示 3 表里 k5 / k6 的占用率）。

练习 5（bank 冲突） 一个 warp 的 32 个线程读共享内存 buf[lane * 2]（32 位的 float）。几路冲突？读 buf[lane * 33] 呢？

答案：lane * 2：字号 0、2、……、62，bank = 字号 mod 32，每个偶数 bank 有 2 个不同的字 → 2 路。lane * 33：bank = lane mod 32，各不相同 → 无冲突（这就是”每行补一个元素”的原理）。

练习 6（LDS.128 的轮数） 一个 warp 的 32 个线程各用 LDS.128 读 16 字节，地址是 lane * 4 个字（连续的 512 字节）。几轮？换成每个线程用 4 条 LDS.32 读同样的数呢？

答案：LDS.128 分 4 个相位，每个相位 8 个线程读连续 32 个字，无冲突，共 4 轮、1 条指令。4 条 LDS.32：第 $j$ 条读 lane*4 + j，32 个字步长 4 → 只用 8 个 bank、每个 bank 4 个字 → 4 路冲突，4 条 × 4 轮 = 16 轮。向量化在这里既省指令，也避免了标量写法的冲突。

练习 7（Tensor Core 的 warp 分块） A100 每 SM 每周期 1024 次 FP16 乘加。一个 warp 算 $W_M \times W_N$，$W_M = W_N = W$，$W$ 至少多大才不受共享内存带宽限制？

答案：需要 $W^2 / (4W) = W/4 \ge 1024/128 = 8$，$W \ge 32$。H100 要 $W/4 \ge 16$，$W \ge 64$（§8.3）。

练习 8（swizzle） 共享内存块每行 64 个 FP16（128 字节），不 swizzle。ldmatrix 读一个 8×8 小矩阵（第 0–7 行的第 8–15 列）几路冲突？用 XOR swizzle（物理 chunk = 逻辑 chunk XOR (行号 mod 8)）后，第 3 行第 8 列落在哪个字节？

答案：每行的第 8 列在行内偏移 16 字节，8 行都落在 bank 4–7 上 → 8 路冲突。第 3 行：line 3，逻辑 chunk 1，物理 chunk $1 \oplus 3 = 2$，字节 $3 \times 128 + 2 \times 16 = 416$（bank 8–11）。

练习 9（流水线级数） H100 上一个 128×256×64 的 FP16 块，每 SM 每周期 2048 次乘加、1.83 GHz。算一级要多久？如果数据的延迟是 1 µs，至少要几级？227 KB 的共享内存放得下几级？

答案：$128 \times 256 \times 64 / 2048 = 1024$ 周期 = 0.56 µs。$1 + \lceil 1 / 0.56 \rceil = 3$ 级。每级 $(128 + 256) \times 64 \times 2 = 48$ KiB，$227 / 48 = 4$ 级（代码 04 演示 4）。

练习 10（分组顺序） Triton 教程的 9×9 例子里，如果 GROUP_M = 9，前 9 个程序实例要读多少个块？

答案：GROUP_M = 9 时组内列主序：前 9 个程序实例是第 0 列的 9 个块，读 $A$ 的全部 9 行（81 块）和 $B$ 的 1 列（9 块）= 90——和行主序一样差。最好的是接近方形的分组。

练习 11（融合的收益） 逐元素的 GELU：$y = 0.5x(1 + \tanh(\sqrt{2/\pi}(x + 0.044715x^3)))$。PyTorch eager 逐个算子大约要 7 个内核（立方、乘、加、乘、tanh、加 1、乘），每个读一个或两个整张量、写一个。和融合成一个内核相比，流量大约多几倍？

答案：大致每个内核读 1–2 个、写 1 个张量，7 个内核共约 $7 \times 2.5 \approx 17$ 个张量大小的流量；融合后读 $x$ 写 $y$，2 个。约 8 倍。（前置讲 P4 §6 用编译器做过 GELU 的融合；第 12 讲的图编译器自动做这件事。）

练习 12（Little 定律在 CPU 上） 一个核每周期能发 4 条向量 FMA，每条 FMA 的延迟是 4 个周期。至少要几条独立的 FMA 链？如果一个微内核的累加器块是 2×16（每个向量 4 个 float），够吗？

答案：$4 \times 4 = 16$ 条。2×16 的块 = 8 个向量累加器 = 8 条链，只够一半——代码 05 演示 3 里 2×16 只有约 56 GFLOP/s，4×16（16 条链）约 100。

练习 13（FP16 累加） 为什么 FP16 累加的误差大约是 FP32 累加的 $2^{13}$ 倍？

答案：两种累加的误差都主要来自每次加法的舍入，正比于”当前部分和 × 单位舍入误差”。FP16 有 11 位有效数字（单位舍入 $2^{-11}$），FP32 有 24 位（$2^{-24}$），比值 $2^{13} = 8192$。代码 06 测出 6,600–9,100。

16. 延伸阅读 · Further reading

主线

• S. Boehm, How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance: a Worklog, 2022. https://siboehm.com/articles/22/CUDA-MMM ；代码：https://github.com/siboehm/SGEMM_CUDA
• NVIDIA, CUDA C++ Programming Guide（编程模型、共享内存、异步拷贝、WMMA）；PTX ISA（mma、ldmatrix、cp.async、wgmma 的片段布局）；Nsight Compute 文档（Speed of Light、warp stall reasons、wavefront）。
• Triton 官方教程：Matrix Multiplication、Fused Softmax。https://triton-lang.org/main/getting-started/tutorials/
• P. Tillet, H. T. Kung, D. Cox, Triton: An Intermediate Language and Compiler for Tiled Neural Network Computations, MAPL 2019.

更多

• NVIDIA CUTLASS 文档（Efficient GEMM in CUDA、CuTe 的 layout 与 swizzle）。https://github.com/NVIDIA/cutlass
• K. Goto, R. van de Geijn, Anatomy of High-Performance Matrix Multiplication, ACM TOMS 2008；F. Van Zee, R. van de Geijn, BLIS: A Framework for Rapidly Instantiating BLAS Functionality, ACM TOMS 2015（CPU 上的打包与微内核）。
• B. Spector et al., ThunderKittens: Simple, Fast, and Adorable AI Kernels, 2024（瓦片级原语，Hazy Research）。
• J. Shah et al., FlashAttention-3: Fast and Accurate Attention with Asynchrony and Low-precision, 2024（WGMMA、TMA、warp 专门化、ping-pong 的完整例子）。

本仓库

• 前置讲 P2：CUDA 编程模型、合并访存、bank 冲突、占用率；simt_lab 模拟器。
• 第 8 讲：Tensor Core 的演进、数据块强度、TMA 与流水线。
• 第 9 讲：分块、循环顺序与”访问次数”的模型。

17. 下一讲预告 · Next lecture

这一讲的矩阵乘法是”计算密集”的：难点是把 Tensor Core 喂饱。推理里另外两类内核的难点不一样。注意力的中间结果（$N \times N$ 的分数矩阵）比输入大得多，FlashAttention 用”分块 + 在线 softmax”让它从不离开片上存储——这是本讲 §10.3 的融合推到极致；decode 时批很小，FlashAttention 又占不满 GPU，flash-decoding 沿 KV 的长度再切一刀。另一类是整个模型的一步 decode：几百个小内核，每个只跑几微秒，启动开销和内核之间的空隙成了主要的时间——Hazy Research 的 megakernel 把整个 Llama-1B 的前向写成一个内核。第 11 讲讲这两件事。

→ 第 11 讲：内核优化技术