solutions/c4_gemm.cu 逐行读C4 讲解：GEMM 三连，复用度一级级往上提
naive → shared memory tiling → register tiling，三个 kernel 做的是同一件事：每从慢的存储读一个数，就多做几次乘加。逐段读参考解：每个版本的数据复用怎么算、两次 __syncthreads 各自防什么、线程怎么分到 4×4 的输出块、外积为什么省读取，以及小矩阵为什么反而变慢。
对照 drillsc4_gemm.cu 读。第 ④ 节的 register tiling 是这题的核心，建议拿纸跟着“线程 (ty=2, tx=5)”的例子把下标算一遍。
① 先看全貌：一个数能被用几次
C[M,N] = A[M,K] @ B[K,N]，C[i][j] = Σₖ A[i][k] · B[k][j]。一共 M·N·K 次乘加（FMA），而 A、B 只有 MK + KN 个数。n = 4096 时，每个输入数在数学上被用了 4096 次。问题是：每用一次，要从哪里把它拿过来？
| 版本 | 每次 FMA 平均从显存（global）读几个数 | 从 shared memory 读几个数 | 复用发生在哪 |
|---|---|---|---|
| v0 naive | 2 | 0 | 没有（全靠 L1/L2 缓存碰运气） |
| v1 smem，T = 32 | 2 / 32 = 0.0625 | 2 | 32×32 的 tile 里，每个数被 32 个线程用 |
| v2 regtile，64×64×16，4×4 | 2048 / 65536 ≈ 0.031 | 8 / 16 = 0.5 | 再加一层：寄存器里每个数用 4 次 |
（v2 的 global 一栏：每轮 k0 读 (BM + BN)·BK = 2048 个数，做 BM·BN·BK = 65536 次 FMA。shared 一栏：每个 k 读 TM + TN = 8 个数，做 TM·TN = 16 次 FMA。）
naive：每写一句话就跑一趟图书馆借两本书，查完还回去。smem：小组 32 个人一次把一摞资料复印到会议桌上（shared memory），大家围着桌子查，查完一摞再去复印下一摞。regtile：每个人再把自己要反复用的几页抄到手边的便签上（寄存器），一张便签上的数字能连用四次，连会议桌都不用伸手去够。
② v0 naive：每个线程一个输出
drillsmatmul_naive_kernel（15–23 行）
const int row = blockIdx.y * blockDim.y + threadIdx.y;
const int col = blockIdx.x * blockDim.x + threadIdx.x; // threadIdx.x → 列
if (row >= M || col >= N) return;
float acc = 0.0f;
for (int k = 0; k < K; ++k) acc += A[(int64_t)row * K + k] * B[(int64_t)k * N + col];
C[(int64_t)row * N + col] = acc;
为什么让 threadIdx.x 走列：一个 warp 是同一个 row、32 个连续 col。
- 读
B[k][col]：32 个连续 float，合并。 - 读
A[row][k]：32 个线程读同一个地址，硬件广播一次就够。 - 写
C[row][col]：连续，合并。
如果反过来让 threadIdx.x 走行，读 A 和写 C 都会变成跨步。问题出在复用上：A 的第 row 行被同一行的 N 个线程各读一遍，B 的每一列被 M 个线程各读一遍，全靠缓存兜底。
③ v1 shared memory tiling：两次 barrier
drillsmatmul_smem_kernel（28–53 行）
block = 32×32 个线程，负责 C 的一个 32×32 tile。沿 K 方向每次取 A 的一个 32×32 tile 和 B 的一个 32×32 tile：
逐行：
- drills40–41 行：每个线程搬 A 的
(row, k0 + tx)和 B 的(k0 + ty, col)。越界的位置填 0：0 乘任何数都是 0，内层循环就不需要任何边界判断。 - drills44–49 行：内层循环
acc += As[ty][k] * Bs[k][tx]。同一个 warp（相同 ty，tx = 0..31）读As[ty][k]是同一个地址（广播），读Bs[k][0..31]是 32 个连续 float，落在 32 个不同的 bank，都没有 bank conflict。 - 每个从显存读进来的数，在 shared memory 里被 32 个线程各用一次：global 读取减少到 1/32。
block 用了 32×32 = 1024 个线程，正好是上限（drills149 行）。
④ v2 register tiling：每个线程算 4×4
drillsmatmul_regtile_kernel（56–119 行）
参数（drills57–58 行）：block 负责 C 的 BM×BN = 64×64，K 方向每次 BK = 16，每个线程负责 TM×TN = 4×4 = 16 个输出，所以只要 (64/4)·(64/4) = 256 个线程。
线程到输出块的映射
const int tx = tid % (BN / TN); // tid % 16：第几个 4 列宽的竖条
const int ty = tid / (BN / TN); // tid / 16：第几个 4 行高的横条
64×64 的 block tile 被切成 16×16 个 4×4 小块，每个线程一块。例如 tid = 37：ty = 2，tx = 5，负责 block tile 里的第 8–11 行、第 20–23 列：
每一格是一个线程（格子里是 tid），对应 C 的一个 4×4 小块。第 ty 行第 tx 列的线程负责 C 的行 4ty…4ty+3、列 4tx…4tx+3（都是相对 block 左上角）。
协作加载：e 循环和转置存放的 As
for (int e = tid; e < BM * BK; e += THREADS) { // A 的 64×16 = 1024 个数，每线程 4 个
const int m = e / BK, k = e % BK; // k 是快维度
As[k][m] = (gm < M && gk < K) ? A[gm * K + gk] : 0.0f; // 注意：转置存
}
- 为什么 e 循环里 k 是快维度：A 是行优先，同一行里 k 连续。相邻线程（相邻 e）读相邻的 k，读显存是合并的。B 的循环里 n 是快维度，同理。
- 为什么存成
As[k][m]（转置）：内层循环需要“固定 k，取连续 4 个 m”（As[k][ty*4 + i]）。转置之后这 4 个数在 shared memory 里是挨着的，读起来快，也方便以后改成一次float4读取。 - 代价：转置写入时相邻线程 m 相同、k 不同，地址相差 64 个 float，正好落在同一个 bank 上（16 路冲突）。题目页的“测完之后”建议用 padding 或
float4加载去优化它。
外积：16 次 FMA 只读 8 个数
for (int k = 0; k < BK; ++k) {
for (int i = 0; i < TM; ++i) a_frag[i] = As[k][ty * TM + i]; // A 的一列里的 4 个数
for (int j = 0; j < TN; ++j) b_frag[j] = Bs[k][tx * TN + j]; // B 的一行里的 4 个数
for (int i = 0; i < TM; ++i)
for (int j = 0; j < TN; ++j) acc[i][j] += a_frag[i] * b_frag[j]; // 4×4 外积
}
从 shared memory 读 4 + 4 = 8 个数进寄存器，做 4 × 4 = 16 次 FMA，每个数被用 4 次。16 个格子分别累加进 acc[i][j]，K 方向全部走完后 acc 就是这 16 个 C 元素的最终值。
acc[TM][TN]、a_frag、b_frag 都是固定大小的局部数组，配合 #pragma unroll 全部展开后下标是常量，编译器就能把它们全部放进寄存器（16 + 4 + 4 = 24 个 float）。如果下标不是编译期常量，数组会被放到 local memory（实际在显存里），性能会断崖式下跌。
外积就是填一张乘法表：上面一排写 4 个数，左边一列写 4 个数，表里 16 个格子各是行头乘列头。你只需要记住 8 个数，就能填满 16 个格子。表越大（TM、TN 越大），每记一个数能填的格子越多，但手里要同时记的格子（寄存器）也越多，记不住就得写到纸上（寄存器溢出到 local memory）。
写回和 __launch_bounds__
- drills109–118 行：每个线程把自己的 16 个结果写回 C，逐个检查
gm < M、gn < N。 - drills60 行
__launch_bounds__(THREADS)：告诉编译器“这个 kernel 每个 block 最多 256 个线程”。编译器据此决定每个线程最多能用多少个寄存器，保证 256 个线程放得下。不写的话编译器按最坏情况（1024 个线程）保守估计。 - host 端（drills162–163 行）：block 是一维的 256 个线程，二维映射是 kernel 里自己用
tid % 16、tid / 16算的。
⑤ 为什么 n = 512 时 regtile 反而最慢
题目页的实测（RTX 5090，TFLOP/s）：n = 512 时 naive 5.2、smem 5.2、regtile 2.5；n = 4096 时 regtile 33.0，是 smem 的 3.6 倍。
n = 512 时 regtile 的 grid 是 (512/64)² = 64 个 block，而 5090 有 170 个 SM：106 个 SM 完全闲着。smem 版本是 (512/32)² = 256 个 block，至少每个 SM 都有活干。
公司有 170 辆出租车（SM）。大订单（大 tile）一单能拉很多人，效率高，但今天只来了 64 个大订单，106 辆车停在车库里。把订单拆小一点（小 tile），每单效率低一些，但所有车都跑起来了，总运量反而更大。
⑥ 常见坑与自测
- 只写了第一个
__syncthreads()：大矩阵上偶尔出错，而且很难复现。 - 边界不填 0：K 不是 32 的倍数时，最后一个 tile 里读进了越界的垃圾值，加进了结果。
- acc 用了运行时下标：例如
acc[i][j]的循环没法展开，数组被放到 local memory。用ncu看local load/store是否为 0。 threadIdx.x走行：读写都不合并，比 naive 还慢。
自测 1：v1 里一个 block 每轮 k0 从显存读多少个数？做多少次 FMA？每个数被用了几次？
每轮读 A、B 各 32 × 32，共 2048 个数；做 32 × 32 × 32 = 32768 次 FMA。每次 FMA 用 2 个操作数，共 65536 次“使用”，除以 2048 次读取，每个读进来的数被用了 32 次。换个说法，每次 FMA 平均读 2048 / 32768 = 0.0625 = 2/32 个数，和第 ① 节表里一致。
自测 2：TM = TN = 8 会怎样？
每个线程 64 个 acc，加上 16 个 fragment，寄存器至少 80 个 float。每次 k 读 16 个数做 64 次 FMA，shared memory 读取再减半。但 block 里只需要 (64/8)² = 64 个线程，或者把 BM、BN 加大到 128。寄存器占用高会降低 occupancy，太高还会溢出。这是 GEMM 调参的核心权衡，Triton 的 autotune（T4）就是在搜这些参数。
源码与资料索引
本仓库：drillssolutions/c4_gemm.cu · drillsmine/c4_gemm.cu · drillstests/test_c4_gemm.py · drillscontent/c4.md
对照实现
- cuda-samples 的 shared memory 版（和 v1 同构）：cudasamplesMatrixMulCUDA，两次 barrier 在 cudasamples109 行 和 cudasamples123 行
- CuTe 版本的同一个算法：cutlassgemm_device，循环体是 cutlasscopy → __syncthreads → gemm，G2 会讲 CuTe 的布局代数
- Tensor Core：cudasampleswmma::mma_sync 一次算 16×16×16；FlashInfer 的 PTX 封装 flashinfermma.cuh
文档与文章
- Simon Boehm · How to Optimize a CUDA Matmul Kernel：从 naive 一路到接近 cuBLAS，每一步都有 profile，和这题的 v0 → v2 正好对应
- Programming Guide · launch_bounds · Shared Memory · Bank Conflicts · WMMA