GPU GEMM 优化完全拆解（二）：寄存器分块、向量化访存与 Warp Tiling

动机：从 kernel 3 的病根出发#

kernel 3 的指令混合里，LDS（共享内存加载）占了大头，一次 FMA 配两次 LDS。要压低这个比例，只有一条路：让每个 LDS 的结果被多个 FMA 使用。换句话说，每个线程从”算 1 个输出元素”变成”算 8 个输出元素”，共享内存访问总量不变（还是每个元素读一遍），但每个结果摊到的 LDS 少了 8 倍。

设计：1D 分块，每线程一条”竖条”#

kernel 4 的配置：BM=BN=64，BK=8，TM=8。每个 block 负责 C 的一个 64×64 子块，256×256 个线程中的 512 个组成一个 block（blockDim = (BM×BN)/TM = 512），每个线程负责同一列上连续的 8 个元素（TM 沿 M 方向）。共享内存需求：$BM \times BK + BN \times BK = 64 \times 8 + 64 \times 8 = 1024$ 个 float，正好 4KB。

下图是它的访存与计算结构，作者把两个线程在内存循环中访问的值用橙色和红色标了出来：

kernel 4 的 1D Blocktiling 示意图（来源：siboehm 工作日志）：block 负责 64×64 的 C 子块；共享内存里缓存 A 的 64×8 子块与 B 的 8×64 子块；每个线程沿列方向算 8 个结果。

核心改动在内层循环：把点积循环（dotIdx）放到外面，先把 B 的一个元素缓存进寄存器变量 tmpB，内层 8 次乘加全部复用这一个值：

1
// 每个线程在寄存器文件里缓存 TM=8 个输出累加器
2
float threadResults[TM] = {0.0};
3

4
// 外层循环：沿 K 方向推进 block 子块
5
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
6
  // 加载共享内存（与 kernel 3 相同：每线程搬 1 个 A 元素 + 1 个 B 元素）
7
  As[innerRowA * BK + innerColA] = A[innerRowA * K + innerColA];
8
  Bs[innerRowB * BN + innerColB] = B[innerRowB * N + innerColB];
9
  __syncthreads();
10

11
  // 推进 block 子块
12
  A += BK;
13
  B += BK * N;
14

15
  // 计算：dotIdx 放外层，让 B 的元素只读一次、被 TM 次乘加复用
16
  for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
17
    float tmpB = Bs[dotIdx * BN + threadCol];      // B 元素缓存进寄存器
18
    for (uint resIdx = 0; resIdx < TM; ++resIdx) {
19
      threadResults[resIdx] +=
20
          As[(threadRow * TM + resIdx) * BK + dotIdx] * tmpB;
21
    }
22
  }
23
  __syncthreads();
24
}
25

26
// 写回：8 个结果各写一个 C 元素
27
for (uint resIdx = 0; resIdx < TM; ++resIdx) {
28
  C[(threadRow * TM + resIdx) * N + threadCol] =
29
      alpha * threadResults[resIdx] + beta * C[(threadRow * TM + resIdx) * N + threadCol];
30
}

为什么把 tmpB 单独提出来？内层循环里，8 次乘加共享同一个 B 元素、各自搭配一个不同的 A 元素。如果不缓存，Bs[...] 就要被读 8 次；缓存后只读 1 次，剩下的 7 次共享内存读取被省掉了。这就是寄存器分块的本质：把”读进寄存器”和”用寄存器计算”分开，一次读、多次算。

收益的精确账本：每结果的访存次数#

上一篇 kernel 3 中每个线程只算 1 个结果，K 循环迭代 $K/32$ 次（BK=32），每轮迭代做 2 次全局内存加载、$32 \times 2 = 64$ 次共享内存加载。摊到每个结果上：

kernel 4 中每个线程算 8 个结果，K 循环迭代 $K/8$ 次（BK=8），每轮迭代 2 次全局加载、$8 \times (1 + 8) = 72$ 次共享内存加载（8 步 dotIdx，每步 1 次读 B + 8 次读 A）。摊到每个结果上：

对比 kernel 3：GMEM 访问减半，SMEM 访问降到原来的 9/16。K=4092 时，每个结果的平均共享内存访问从 8184 次降到 4603 次左右——指令流里的 LDS 密度大幅下降。实测性能 8474.7 GFLOPS（cuBLAS 的 36.5%），是 kernel 3 的 2.8 倍（作者文中按取整数字说成”约 2.2 倍”）。

Nsight Compute 的 warp stall 图印证了病根正在消退：

kernel 4 的 warp stall 分析（来源：siboehm 工作日志）。注意纵轴与 kernel 3 的图不同（作者特别提醒”轴变了”）：等待共享内存管线（MIO Throttle）的周期占比显著下降，说明 LDS 不再是唯一瓶颈。

插曲：循环顺序真的重要吗——编译器的力量#

一个值得展开的细节：作者一开始以为内层循环的顺序（dotIdx 在外还是 resIdx 在外）会影响性能，于是写了两版：

1
// 版本 A：dotIdx 在外（本篇上面的写法），显式缓存 tmpB
2
// 版本 B：resIdx 在外，每次乘加都直接读共享内存
3
for (uint resIdx = 0; resIdx < TM; ++resIdx) {
4
  for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
5
    threadResults[resIdx] +=
6
        As[(threadRow * TM + resIdx) * BK + dotIdx] * Bs[dotIdx * BN + threadCol];
7
  }
8
}

版本 B 直观上要差很多：内层两个循环共 $8 \times 8 \times 2 = 128$ 次共享内存访问，而版本 A 只有 $8 \times 9 = 72$ 次。但实测两者性能相同。看 SASS 汇编才明白原因：两个循环的循环次数都是编译期常量，编译器把它们完全展开，然后做公共子表达式消除（CSE）——重复读取 Bs 同一地址的指令被合并，最终生成的汇编和版本 A 一模一样。

更妙的是展开后的 SASS 里出现了一个伏笔：

1
LDS R26, [R35.X4+0x800]   // 32b 加载：从 As 读一个元素
2
LDS.128 R8, [R2]          // 128b 加载：从 Bs 一次读 4 个元素（float4 向量化）
3
LDS.128 R12, [R2+0x20]
4
LDS.128 R20, [R2+0x60]
5
LDS R24, [R35.X4+0x900]   // As 依然是 32b 加载
6
LDS.128 R16, [R2+0x40]
7
LDS.128 R4, [R2+0x80]
8
LDS R38, [R35.X4+0xd00]

编译器自动把 Bs 的连续读取合并成了 128 位的 LDS.128（一次读 4 个 float），但 As 的读取还是 32 位——因为 As 的布局是行优先（As[row*BK + col]），线程按列访问，地址不连续。如果能把 As 也转置成”列优先”，As 的读取同样可以向量化。 这个伏笔会在第六步兑现。

顺带一提：grid 映射的 30% 差距#

kernel 4 源码的注释里有一条很值钱的实测：把 cRow/cCol 与 blockIdx.x/blockIdx.y 的对应关系对调，大矩阵性能会掉约 30%。原因在 L2 缓存：当前配置（gridDim(N/BN, M/BM)，blockIdx.x 沿 N 方向）保证连续编号的 block 顺序访问 B 的连续列，同时共享 A 的同一行——空间局部性好，L2 命中率高；反过来则是共享 A 的列、乱序访问 B，L2 频繁失效。这个”把 block 到 C 子块的映射做成顺序友好”的思想，就是后面 kernel 10 里提到的 thread swizzling（线程交织）的雏形，虽然那次的实测结果是负的，后面会讲。

为什么 1D 还不够：同一个病，只是轻了#

kernel 4 从 12.8% 跳到 36.5%，但离 78% 还有距离。它的访存模式里，每个结果仍要摊到 $K/32$ 次全局访问和 $9K/8$ 次共享内存访问——只是比 kernel 3 少，不是没有。按第一篇的方法论：继续提高算术强度（每字节数据干的活），直到不再受内存管线限制。

方形 vs 条带：输入复用的几何学#

kernel 4 的 1D 分块里，8 个结果共享 1 个 B 元素、需要 8 个不同的 A 元素——9 个输入喂 8 次乘加。为什么不干脆让每个线程算一个 8×8 的方形？那样 8 个 A 元素和 8 个 B 元素可以任意两两配对，64 次乘加只需要 16 个输入：

为什么要算”方形”结果而不是”条带”结果（来源：siboehm 工作日志）：8 个 B 元素和 8 个 A 元素两两配对可以产生 64 个结果，而 1 个 B 元素配 8 个 A 元素只有 8 个结果。

用算术强度（每字节输入对应的 FMA 次数）表达更清楚：

方形分块的算术强度是条带的 4.5 倍。这是 GEMM 优化里最重要的几何直觉：输出是输入的二维函数，把输出分块成方形，输入复用的机会最大。下面这张图把”算更多结果 → 算术强度上升”的链条画了出来：

算更多结果如何提高算术强度（来源：siboehm 工作日志）：所有 kernel 的总 FLOPs 相同，区别只在全局内存访问总量——每线程算的结果越多，需要的输入搬运次数越少。

实现一：合作加载——每线程搬多个元素#

kernel 5 的配置：BM=BN=128，BK=8，TM=TN=8。每个 block 负责 128×128 的 C 子块，blockDim = (128×128)/(8×8) = 256 个线程，共享内存 $128 \times 8 + 8 \times 128 = 2048$ 个 float（8KB）。

每轮迭代要搬 2048 个 float 进共享内存，256 个线程每人要搬 8 个。加载循环用步长（stride）方式展开：每个线程沿列方向每隔 stride 行取一个元素，保证每个线程的访问在全局内存里是连续的（合并访问）：

kernel 5 的 GMEM 加载示意（来源：siboehm 工作日志）：线程们合作填满共享内存，每个线程用 stride 循环搬运多个元素，同时保持全局内存访问的连续性。

1
// 每个线程要搬运的元素行数步长
2
const uint strideA = numThreadsBlocktile / BK;   // 256 / 8 = 32
3
const uint strideB = numThreadsBlocktile / BN;   // 256 / 128 = 2
4

5
for (uint loadOffset = 0; loadOffset < BM; loadOffset += strideA) {
6
  As[(innerRowA + loadOffset) * BK + innerColA] =
7
      A[(innerRowA + loadOffset) * K + innerColA];
8
}
9
for (uint loadOffset = 0; loadOffset < BK; loadOffset += strideB) {
10
  Bs[(innerRowB + loadOffset) * BN + innerColB] =
11
      B[(innerRowB + loadOffset) * N + innerColB];
12
}
13
__syncthreads();

注意这里加载 A 时线程的映射：innerRowA = threadIdx.x / BK，innerColA = threadIdx.x % BK——线程沿 K 方向排布，同一 warp 的 32 个线程访问的是 A 的同一行、连续 K 列（合并访问）。B 的加载同理按 BN 排布。

实现二：计算循环——寄存器缓存 + 外层积#

共享内存填好后，每个线程要完成 8×8 的乘加。计算循环变成三层：外层 dotIdx（沿 K 步进），中间层把需要的 A、B 元素读进两个寄存器数组 regM[8]、regN[8]，内层做 8×8 的外层积（outer product）累加进 threadResults[64]：

kernel 5 的 2D Blocktiling 结构（来源：siboehm 工作日志）：上方是沿 K 方向推进的外层循环（A 向右、B 向下移动子块），下方是三个内层循环——dotIdx 步进、TM 维、TN 维。

1
float threadResults[TM * TN] = {0.0};   // 64 个输出累加器，全在寄存器里
2
float regM[TM] = {0.0};                 // A 的寄存器缓存
3
float regN[TN] = {0.0};                 // B 的寄存器缓存
4

5
for (uint bkIdx = 0; bkIdx < K; bkIdx += BK) {
6
  // ... 合作加载共享内存（见上），__syncthreads() ...
7
  A += BK;      // A 子块右移 BK 列
8
  B += BK * N;  // B 子块下移 BK 行
9

10
  for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
11
    // 把需要的共享内存元素读进寄存器缓存
12
    for (uint i = 0; i < TM; ++i) {
13
      regM[i] = As[(threadRow * TM + i) * BK + dotIdx];
14
    }
15
    for (uint i = 0; i < TN; ++i) {
16
      regN[i] = Bs[dotIdx * BN + threadCol * TN + i];
17
    }
18
    // 外层积：regM 的每个元素与 regN 的每个元素相乘累加
19
    for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
20
      for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
21
        threadResults[resIdxM * TN + resIdxN] +=
22
            regM[resIdxM] * regN[resIdxN];
23
      }
24
    }
25
  }
26
  __syncthreads();
27
}

为什么要在中间层显式加一层”读进寄存器”？因为内层 64 次乘加要反复使用同样的 16 个值——每次乘加都去读共享内存的话，64 次乘加配 128 次 LDS，比例又回去了。先取到寄存器，64 次乘加只需要 16 次共享内存访问（8 次读 A + 8 次读 B）。下图展示了 dotIdx 循环推进时，哪些共享内存元素被装进寄存器：

寄存器缓存的装载过程（来源：siboehm 工作日志）：图中的尺寸为便于绘制做了缩减，真实 kernel 中 BK=TM=TN=8。每一步 dotIdx，线程从 As 装载一列（8 个元素）到 regM、从 Bs 装载一行（8 个元素）到 regN，然后做 64 次乘加。

账本与性能#

kernel 5 的访存账本（每线程每轮迭代：全局加载 $2048/256 = 8$ 次；共享内存加载 $8 \times (8 + 8) = 128$ 次；每线程算 64 个结果）：

对比 kernel 4（$K/32$ 和 $9K/8$）：每结果的全局访问再砍一半，共享内存访问降到 1/18。实测 15971.7 GFLOPS（68.7%），又是一次接近 2 倍的提升。注意一个细节：这个版本的共享内存读取 As[(threadRow*TM+i)*BK + dotIdx] 是按列访问（threadRow 连续 → 行号连续 → 地址间隔 BK×4B=32B → 全部落在同一个 bank 上），正是第一篇埋下的 bank conflict 伏笔，kernel 6 之后会专门处理。

优化一：加载时转置 As，让共享内存读取也能向量化#

第四步的 SASS 里，Bs 的读取被编译器合并成了 128 位的 LDS.128，As 却不行——因为 As 按行优先存储，而计算时线程按列访问。解法很优雅：在从全局内存往共享内存搬运的时候就把 As 转置。存储时把 A 的”一行 8 个元素”拆开放成”一列”，让 As 变成列优先布局。这样计算循环里 As[dotIdx * BM + threadRow * TM + i] 的下标 i 连续 → 地址连续 → 编译器可以生成 LDS.128：

As 转置的效果（来源：siboehm 工作日志）：同样的三个内层循环，As 在共享内存中的排布从行优先变成列优先后，线程按连续地址读取，编译器自动生成 128 位向量化加载。

1
// 加载 A 时顺手转置：tmp 是全局内存里连续 4 个元素（一次 128 位读取）
2
float4 tmp = reinterpret_cast<float4 *>(&A[innerRowA * K + innerColA * 4])[0];
3
As[(innerColA * 4 + 0) * BM + innerRowA] = tmp.x;   // 拆开放进"转置后"的位置
4
As[(innerColA * 4 + 1) * BM + innerRowA] = tmp.y;
5
As[(innerColA * 4 + 2) * BM + innerRowA] = tmp.z;
6
As[(innerColA * 4 + 3) * BM + innerRowA] = tmp.w;

这一步单独带来约 500 GFLOPS（相对 kernel 5 的 16 TFLOPs 约 3%）。收益不大，但它让 As、Bs 的读取统一成了 128 位向量，为后面的自动调优清掉了障碍（向量化访问要求分块尺寸满足整除性约束，见第七步）。

优化二：float4 向量化全局内存访问#

全局内存的搬运同样向量化：加载和写回都用 float4（4 个连续 float = 128 位）替代单个 float：

1
// B 的搬运：一次 128 位加载 + 一次 128 位存储
2
reinterpret_cast<float4 *>(&Bs[innerRowB * BN + innerColB * 4])[0] =
3
    reinterpret_cast<float4 *>(&B[innerRowB * N + innerColB * 4])[0];

对应的 SASS 从 32 位的 LDG.E/STG.E 变成 128 位的 LDG.E.128/STG.E.128——指令数减到 1/4。作者在这里困惑过一个问题：为什么不能直接写 4 个单元素赋值让编译器自己合并？

1
Bs[innerRowB * BN + innerColB * 4 + 0] = B[innerRowB * N + innerColB * 4 + 0];
2
Bs[innerRowB * BN + innerColB * 4 + 1] = B[innerRowB * N + innerColB * 4 + 1];
3
Bs[innerRowB * BN + innerColB * 4 + 2] = B[innerRowB * N + innerColB * 4 + 2];
4
Bs[innerRowB * BN + innerColB * 4 + 3] = B[innerRowB * N + innerColB * 4 + 3];

答案在硬件对齐上：LDG.E.128 要求地址 16 字节对齐。kernel 的参数是指针（float *B），编译器无法在编译期验证传入的指针满足对齐要求，所以不敢生成向量化指令。reinterpret_cast<float4 *> 的实质作用不是类型转换，而是向编译器承诺”这个地址一定 16 字节对齐”——一个编译期约束，代价是要求调用方真的保证对齐（生产代码里通常配合 cudaMalloc 的 256 字节对齐保证）。对比共享内存：它是运行时管理的片上存储，编译器不需要用户承诺就能自动向量化共享内存访问。

kernel 6 总计 18237.3 GFLOPS（78.4%）。Nsight Compute 还能看到几个”已知但未解决”的问题：共享内存访问仍有 bank conflict（cuBLAS 是避免的）、占用率高于必要水平、没有实现双缓冲（double buffering）。作者选择了先摘更容易的果子——自动调优。

5 个模板参数与合法配置的约束#

到 kernel 6 为止，代码里累积了 5 个模板参数：BM、BN、BK（控制全局内存 → 共享内存的分块）和 TM、TN（控制共享内存 → 寄存器的分块）。它们不是独立的——合法性受一堆整除性约束限制，比如：

• 共享内存加载要向量化：BM × BK 必须是 4 × NUM_THREADS 的整数倍（每个线程每轮搬运 4 个 float）；
• 全局内存加载要避免”最后一行只加载一半”的量化问题：NUM_THREADS × 4 必须被 BK 整除；
• 寄存器分块不能留下死角：BM 必须是 16 × TM 的整数倍。

作者把这些约束写成 static_assert，编译期就把非法组合挡掉（编译器在模板实例化时检查）。约束滤掉大部分组合后，仍然剩约 400 个合法配置。于是写了个 bash 脚本逐个编译、跑基准，找出最快的那组——这就是自动调优（autotuning）。

结果：最优参数随 GPU 而变#

kernel 6 的初始参数是 BM=BN=128、BK=TM=TN=8。自动调优在 A6000 上找到的最优解是 BM=BN=128、BK=16、TM=TN=8——BK 从 8 翻倍到 16，性能从 19 TFLOPs 提到 20 TFLOPs（表格口径 18237.3 → 19721.0 GFLOPS，84.8%）。

更有意思的是作者顺手在另一块卡（A100 SMX4 40GB）上跑了同样一组参数：A6000 的最优配置在 A100 上只有 12.0 TFLOPs，而 A100 自己的最优配置（BM=BN=64、BK=16、TM=TN=4，块更小、每线程结果更少）能到 12.6 TFLOPs——差了 6%。A100 的 FP32 峰值（19.5 TFLOPs）本来就低于 A6000（38.7 TFLOPs），所以绝对数字不可比，但**“最优参数随 GPU 型号改变”这件事是确定的**。作者坦言无法解释为什么这些特定参数在各自卡上最优——寄存器压力、占用率、L2 行为、bank 分布交织在一起，已经超出直觉能把握的范围。

这解释了高级库的工程形态：cuBLAS 内部存的是一张 {GPU 型号, 矩阵尺寸, 数据类型} → 最优 kernel 的查找表；编译器级的 Triton 也内置自动调优例程。在参数空间巨大且无法理论推导最优解的情况下，搜索比思考更可靠。 作者引用了 nelhage 的话来表达这种”不满足”：我们应当保持”计算机是可以被理解的”信念，但接受目前还解释不了的事实。

warp 为什么值得显式化#

到目前为止，我们的代码里只有 block 和 thread 两层概念。但 GPU 硬件实际执行的粒度是 warp：32 个连续 threadId 的线程共享一条指令流。warp 在 CUDA 代码里”隐形”（软件没有直接对应物），但它的行为直接影响性能：

• 调度单位：SM 上的 warp 调度器（A6000 每个 SM 有 4 个）以 warp 为单位分发指令；
• bank conflict 的边界：冲突只发生在同一 warp 内部，不同 warp 互不影响；
• 寄存器缓存局部性：较新的 GPU 上有寄存器缓存（operand collector），同一 warp 内寄存器访问有局部性收益。

warp 与 warp 调度器（来源：siboehm 工作日志，作者想象画法）：A6000 的每个 SM 有 4 个 warp 调度器，每条指令按 warp 为单位发射。

每个线程可以用一个表达式算出自己属于哪个 warp：warpId = threadIdx.x % 32（warpSize 在所有 NVIDIA GPU 上都是 32）。

三级分块：block → warp → thread#

kernel 10 在原有的两层分块之间插入一层 warp tiling，把并行层次彻底显式化：

• Block tiling：不同 block 并行跑在不同 SM 上；
• Warp tiling：不同 warp 并行跑在不同 warp 调度器上，或并发排队在同一个调度器上；
• Thread tiling：同一线程内的少量指令级并行（ILP）。

kernel 10 的循环层次（来源：siboehm 工作日志）：新增的 warptiling 循环位于 blocktiling 与 threadtiling 之间。

A6000 上的配置：128 线程（4 个 warp），BM=BN=128，BK=16，WM=WN=64，WNITER=4，TM=8，TN=4。block 的 128×128 子块切成 2×2 个 warp 分块，每个 warp 负责 64×64；每个 warp 的 64×64 再沿 N 方向切成 WNITER=4 个 64×16 的 warp 子块，每个线程负责其中 8×4 的部分。各层尺寸的约束关系：

WMITER 是每个 warp 子块沿 M 方向切分的份数，这里恰好为 1（每个 warp 子块由整 warp 的 32 个线程一次算完）。代码结构如下（核心循环）：

1
// dotIdx 沿共享内存推进
2
for (uint dotIdx = 0; dotIdx < BK; ++dotIdx) {
3
  // 把本线程在 warp 子块里需要的 A 元素装进寄存器
4
  for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
5
    for (uint i = 0; i < TM; ++i) {
6
      regM[wSubRowIdx * TM + i] =
7
          As[dotIdx * BM + warpRow * WM + wSubRowIdx * WSUBM +
8
             threadRowInWarp * TM + i];
9
    }
10
  }
11
  // 把本线程需要的 B 元素装进寄存器（WNITER 个 warp 子块）
12
  for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
13
    for (uint i = 0; i < TN; ++i) {
14
      regN[wSubColIdx * TN + i] =
15
          Bs[dotIdx * BN + warpCol * WN + wSubColIdx * WSUBN +
16
             threadColInWarp * TN + i];
17
    }
18
  }
19
  // 对每个 warp 子块做外层积累加
20
  for (uint wSubRowIdx = 0; wSubRowIdx < WMITER; ++wSubRowIdx) {
21
    for (uint wSubColIdx = 0; wSubColIdx < WNITER; ++wSubColIdx) {
22
      for (uint resIdxM = 0; resIdxM < TM; ++resIdxM) {
23
        for (uint resIdxN = 0; resIdxN < TN; ++resIdxN) {
24
          threadResults[(wSubRowIdx * TM + resIdxM) * (WNITER * TN) +
25
                        (wSubColIdx * TN) + resIdxN] +=
26
              regM[wSubRowIdx * TM + resIdxM] *
27
              regN[wSubColIdx * TN + resIdxN];
28
        }
29
      }
30
    }
31
  }
32
}

作者把三级分块画成了下面这张大图（这是 siboehm 工作日志里信息量最大的一张图）：

kernel 10 的 Warp Tiling 全景图（来源：siboehm 工作日志）：每个 warp 计算 WSUBN×WNITER 宽、WSUBM×WMITER 高的分块；每个线程计算 WNITER×WMITER 个 TM×TN 小块。

寄存器压力：用占用率换算术强度#

看配置就能猜到 kernel 10 的寄存器账单很重：仅输出累加器就要

个寄存器，加上 regM（8 个）与 regN（16 个）缓存，光数据就 152 个寄存器，再算上地址与循环变量，轻松超过 160 个——Ampere 单线程 255 个寄存器的上限被占掉一大半。每个 block 128 线程，SM 上能同时驻留的 warp 数必然很低。

但第一篇讲过的 Cusp Behavior 在这里给出了设计依据：计算密集的 kernel 不需要高占用率。指令流几乎全是 FMA 时，少量 warp 就能把 FMA 流水线喂满；寄存器和共享内存省下来的空间反而让 block 能开得更大。现代 GEMM 的普遍做法就是”牺牲占用率换每线程计算量”。

实测 21779.3 GFLOPS（cuBLAS 的 93.7%）。作者在文中写”自动调优后从 19.7 TFLOPs 提到 21.7 TFLOPs”并标注为 A100，但这个数字与基准表（A6000 口径 19721 → 21779）一致，而与 A100 上 cuBLAS 只有 14.7 TFLOPs 的实测矛盾（21.7 不可能超过 cuBLAS 的 147%），应当按 A6000 理解——原文的”on an A100”疑为笔误。另一组明确的 A100 数据来自下面的尺寸扫描图（作者注明”这张图在 A100 上生成，绝对 FLOPs 数值与 A6000 不同”）。

与 cuBLAS 的尺寸扫描：为什么小矩阵差这么多#

kernel 10 vs cuBLAS 的尺寸扫描（来源：siboehm 工作日志，A100）：在 2048 和 4096 尺寸上只差几个百分点；小矩阵差距巨大。

在 2048、4096 的大矩阵上，kernel 10 与 cuBLAS 只差几个百分点；但小矩阵（128-1024）差得非常多。原因不在 kernel 本身，而在 cuBLAS 的工程形态：cuBLAS 不是一个 GEMM 实现，而是几百个 GEMM 实现的集合（编译产物约 500MB）。作者用 cuobjdump --list-text 数出了 16 个不同的 SGEMM kernel，再用 Nsight Systems 追踪每次调用实际启动的 kernel：

注意 256 那一行：cuBLAS 调用了两个 kernel——一个分块 GEMM 加上一个归约 kernel。这是 Split-K：矩阵太小时，单个 block 的计算量不足以喂饱整个 GPU，于是把 K 维度切分给多个 block 并行算部分和，再用一个 reduce kernel 累加（代价是需要额外的中间结果存储空间）：

Split-K 示意（来源：siboehm 工作日志，作者自绘并注明”这是我的想象画法，不确定”）：每个 block 只算 C 子块的一部分（沿 K 切分），部分和需要额外存储，再由归约 kernel 合并。

作者的结论是：要写出对所有形状和尺寸都快的库，就必须为不同形态准备专门化的 kernel，运行时按尺寸派发——这正是 cuBLAS 16 个 kernel 存在的意义，也是我们这一系列学到的”一个 kernel 打天下”行不通的原因。

插曲二：thread swizzling 的负结果#

kernel 10 的作者还试了另一项优化：thread swizzling（线程交织）——利用 block 按 blockIdx 递增顺序启动的特性，重排 block 到 C 子块的映射，提高 L2 局部性（L2 是整块 GPU 共享的全局内存缓存）。结果：没有任何性能提升，作者最终移除了这段代码，理由是L2 命中率已经高达 80%，没有可优化的空间了。

又一个负结果，但和 kernel 7/8 的负结果一样有价值：它界定了优化的边界——当某项资源（这里是大到 80% 的 L2 命中率）已经接近饱和时，围绕它的优化注定是空转。先测量，再决定优化哪里，而不是把教科书上的优化挨个试一遍。