02 · CUDA from a Systems Engineer’s Eyes · 报告 01
同一块 RTX 4090、同一个可执行文件、同一次开机： 一个 FP32 矩阵乘从 naive 改到了 cuBLAS 的 71.5%。 本报告只讲代码：每一级改了哪一行、为什么这样改、原理图是什么样。 每一级只引入一个新概念， 所以每一步的收益归谁很清楚。 最重要的发现： 收益最大的两步（8.209× 和 3.359×），算术强度完全没变。
本报告只讲代码怎样一步步改成现在的样子，以及每一步的原理。 数字的来源、测量方法和八个指标的算法，见 第 02 份报告。
教科书只给一本账：算术强度 → roofline。 实测中，这本账在七个 kernel 里失效五次， 达成率最高到 2112.7%。 达成率 2112.7% 说明这个 kernel 不在带宽墙上： 它读的字节大部分命中了 4090 的 75.5 MB L2，没有走到显存。
| 账 | 量什么 | 什么时候它才是墙 | 这条链里管用几次 |
|---|---|---|---|
| 字节账 | 算术强度 → roofline 上限 | 显存带宽真的饱和时 | 七次里失效五次 |
| 请求账 | 一条指令触发几个内存事务 | 访存不合并时 | 只在 K1→K2 成立一次 |
| 条数账 | 指令条数 × 分类单价 | 前两者都不是墙时 | K3 之后全靠它 |
| 级 | 引入的概念 | 动了哪本账 | 算术强度 | load/FMA | 比上一级快 | 占 cuBLAS |
|---|---|---|---|---|---|---|
| K1 naive | 一线程一输出，索引映射写反 | 基线 | 0.250 | 2.000 | 基准 | 1.13% |
| K2 合并 | 把 x/y 的映射换回来，其余一字不改 | 请求账 | 0.250 不变 | 2.000 不变 | 8.21× | 9.30% |
| K3 shared memory 分块 | 片上暂存，一块一块搬 | 字节账 | 7.969 ↑32× | 2.062 | 1.16× | 10.80% |
| K3c 书上的粗化 | Fig. 6.13：一线程 4 个输出，复用 M tile | 字节账 | 12.720 ↑1.6× | 2.039 | 1.08× | 11.69% |
| K4 1D 寄存器分块 | 一线程 8 行，N 元素进寄存器用 8 次 | 条数账 | 7.969 和 K3 相同 | 1.188 ↓ | 3.10× | 36.29% |
| K5 2D 寄存器分块 | 8 行 × 4 列 = 32 个输出，两个方向一起复用 | 条数账+字节账 | 12.720 | 0.414 ↓ | 1.66× | 60.28% |
| K6 float4 | 128 位访存，4 条并 1 条 | 条数账 | 12.720 不变 | 0.104 ↓ | 1.19× | 71.50% |
要细看的两处
① K2 的算术强度和 load/FMA 都和 K1 相同，却快 8.209 倍。 两本账都说"没变"，时间却变了。第三本账（请求账）就是从这里来的。
② K4 的算术强度和 K3 相同（都是 7.969），却比 K3 快 3.359 倍，是整条链上第二大的一步（仅次于 K2）。K4 是从 K3 改出来的，所以这里和 K3 比；总表「比上一级快」一列的 3.10× 是和上一行 K3c 比。 算术强度真正涨了的 K3c（7.969 → 12.720）， 只快了 8%。 教科书"提高算术强度就能更快"的因果链，在这里断得很干净。
// ---------------------------------------------------------------------------
__global__ void sgemm_k1(int M, int N, int K, float alpha,
const float *A, const float *B, float beta, float *C) {
// 这两行是 CUDA 里最常见的一句话："我是谁"。
// blockIdx 是我这个 block 在 grid 里的编号，blockDim 是 block 的宽度，
// threadIdx 是我在 block 内的编号 —— 三者拼出全局唯一的坐标。
const int x = blockIdx.x * blockDim.x + threadIdx.x; // 行（M 方向）
const int y = blockIdx.y * blockDim.y + threadIdx.y; // 列（N 方向）
if (x < M && y < N) {
float acc = 0.0f;
for (int i = 0; i < K; ++i)
acc += A[x * K + i] * B[i * N + y]; // A 走一行，B 跨一列
C[x * N + y] = alpha * acc + beta * C[x * N + y];
}
}
算法和结果都正确，唯一的问题是 x（行）和 y（列）的分工。
同一个 warp 里，threadIdx.x 连号，threadIdx.y 相同。
K1 把 threadIdx.x 给了行，所以 32 个线程落在 A 的 32 个不同行上。
它们的地址间隔 K 个 float，落在 32 个不同的 128 字节段里。
row*4+col），不是元素的值。
一条 LDG 被拆成 32 个请求
请求账：每个 FMA 33 个内存请求
A 的那条 load → 32 个请求（32 个线程各落一段）
B 的那条 load → 1 个请求（32 个线程落同一段，合并）
----
33
// 而字节账完全看不到这件事：它只数"逻辑上要读多少字节"，
// 不管这些字节是一次取回来的，还是分 32 次、每次浪费 124 字节取回来的。
实测 211.911 ms，只有 cuBLAS 的 1.13%。 这是整条链的起点。
// ---------------------------------------------------------------------------
__global__ void sgemm_k2(int M, int N, int K, float alpha,
const float *A, const float *B, float beta, float *C) {
const int y = blockIdx.x * blockDim.x + threadIdx.x; // 列（N 方向）
const int x = blockIdx.y * blockDim.y + threadIdx.y; // 行（M 方向）
if (x < M && y < N) {
float acc = 0.0f;
for (int i = 0; i < K; ++i)
acc += A[x * K + i] * B[i * N + y];
C[x * N + y] = alpha * acc + beta * C[x * N + y];
}
}
和 K1 逐字对比，只差两行
K1: const int x = blockIdx.x * blockDim.x + threadIdx.x; // 行 const int y = blockIdx.y * blockDim.y + threadIdx.y; // 列 K2: const int y = blockIdx.x * blockDim.x + threadIdx.x; // 列 const int x = blockIdx.y * blockDim.y + threadIdx.y; // 行
threadIdx.x 赋给了不同的变量，
一个 warp 的 32 个地址就从"撒在近 500 KB 里"变成"挤在 128 B 里"。
代码没变，形状变了。边界测试意外抓到的东西
N=4096（2 的幂） N=4097（打破对齐） K1 211.911 ms 63.675 ms ← 快了 3.33 倍 K2 25.815 ms 26.615 ms ← 只慢 3% K1/K2 8.209× 2.392× // 4097 上两者都不受跨度冲突影响，剩下的 2.392× 才是"不合并"本身的代价。 // 于是 8.209 ÷ 2.392 = 3.43 就是 2 的幂跨度的额外惩罚。 // 回乘检验：2.392 × 3.43 = 8.209，和实测逐位相同。
N=4096 时，A 的行跨度正好是 16384 字节，是 2 的幂。 这在"不合并"之上又叠了一层冲突。N=4097 打破了对齐，这层冲突就消失了。
但我们没能指认出具体是哪个硬件结构在冲突（L2 slice 哈希？DRAM row buffer？）。 把地址逐位算一遍，两种 N 在"行号"这一级的模式完全相同， 所以"32 个地址撞同一个 L2 set"这个解释站不住。 要判断到具体结构，需要 ncu 的 DRAM 指标，而这一批数据没有（见 报告 02 · §7）。
对照组：同样从 4096 改到 4097，其余 kernel 只慢 2.8%–7.1%， cuBLAS 慢 15.7%。所以不是"4097 让所有 kernel 都变快"，而是只有 K1 被单独惩罚了。
所以"合并访存是性价比最高的一步"这句话要加个注： 那 8.209 倍里只有 2.392 倍来自合并本身， 另外 3.43 倍来自顺带治好的跨度冲突。 书里和多数教程都没提这一点。
// ---------------------------------------------------------------------------
__global__ void sgemm_k3(int M, int N, int K, float alpha,
const float *A, const float *B, float beta, float *C) {
__shared__ float Mds[TILE][TILE];
__shared__ float Nds[TILE][TILE];
const int tx = threadIdx.x, ty = threadIdx.y;
const int Col = blockIdx.x * TILE + tx;
const int Row = blockIdx.y * TILE + ty;
float acc = 0.0f;
for (int ph = 0; ph < (K + TILE - 1) / TILE; ++ph) {
const int aCol = ph * TILE + tx;
const int bRow = ph * TILE + ty;
Mds[ty][tx] = (Row < M && aCol < K) ? A[Row * K + aCol] : 0.0f;
Nds[ty][tx] = (bRow < K && Col < N) ? B[bRow * N + Col] : 0.0f;
__syncthreads();
#pragma unroll
for (int k = 0; k < TILE; ++k)
acc += Mds[ty][k] * Nds[k][tx];
__syncthreads();
}
if (Row < M && Col < N)
C[Row * N + Col] = alpha * acc + beta * C[Row * N + Col];
}
下图标的"行 96–127 × 列 64–95"是从 blockIdx 算出来的，不是随便挑的。
先把它落到一个具体的 block 上，图才好懂。
// K3 沿用 K2 的映射（见 §3）：blockIdx.x → 列，blockIdx.y → 行 const int Col = blockIdx.x * TILE + tx; const int Row = blockIdx.y * TILE + ty; // 取 block(2,3)，也就是 blockIdx.x = 2、blockIdx.y = 3： Row = 3 * 32 + ty = 96 + ty // ty = 0..31 → 行 96 … 127 Col = 2 * 32 + tx = 64 + tx // tx = 0..31 → 列 64 … 95 // 所以这个 block 负责 C 的 行 96–127 × 列 64–95 这一块 32×32 的输出。 // 注意 (2,3) 里的 2 是 x（列）、3 是 y（行）—— 和"行在前列在后"的习惯正好相反， // 这是 K1 那个索引 bug 的同一个来源（§2·二）。
再看这个 block 里的线程 (tx=5, ty=7)，
它负责 C[103][69]
（Row = 96+7 = 103，Col = 64+5 = 69）。
下表是它在 128 个阶段里搬什么、算什么：
| 这个线程搬什么 | 然后用什么算 | |
|---|---|---|
| ph = 0 | Mds[7][5] ← A[103][5] Nds[7][5] ← B[7][69] |
用 Mds[7][0..31] × Nds[0..31][5]= A[103][0..31] · B[0..31][69] |
| ph = 1 | Mds[7][5] ← A[103][37] Nds[7][5] ← B[39][69] |
用 Mds[7][0..31] × Nds[0..31][5]= A[103][32..63] · B[32..63][69] |
| … | … | … |
| ph = 127 | Mds[7][5] ← A[103][4069] Nds[7][5] ← B[4071][69] |
= A[103][4064..4095] · B[4064..4095][69] |
算术强度从 0.250 涨到 7.969，涨了 32 倍。 原因：在 32×32 的 tile 里，每个搬进来的元素被 32 个线程用了 32 次。
但实测只快 1.162×
字节账涨了 32 倍，时间只省了 14%。为什么？
K2 的 roofline 达成率 = 2112.7% // 达成率超过 100% 二十倍，说明 K2 根本没在带宽墙上 —— // 它那些"逻辑上要读的字节"绝大部分命中 L2，没走到 DRAM。 // 一个不受带宽限制的 kernel，你把它的 DRAM 算术强度提高 32 倍，换不来时间。
这 14% 可以用条数账算出来（三个单价的标定见报告 02 · §3）：
K2 = 2.000 条 LDG × 1.94 周期 = 3.88 周期
K3 = 1 条 32 路 LDS × 2.686
+ 1 条广播 LDS × 0.534
+ 0.062 条 LDG × 1.94 = 3.34 周期
预测 3.88 ÷ 3.34 = 1.162× 实测 1.162× ✓
在这台卡上，K3 的 shared memory 分块不是靠省带宽赢的。 它赢在把 2 条 global memory 的 load（每条 1.94 周期）换成 1 条 32 路和 1 条广播的 shared memory load （2.686 + 0.534）。
// ---------------------------------------------------------------------------
__global__ void sgemm_k3c(int M, int N, int K, float alpha,
const float *A, const float *B, float beta, float *C) {
__shared__ float Mds[TILE][TILE];
__shared__ float Nds[TILE][TILE];
const int tx = threadIdx.x, ty = threadIdx.y;
const int Row = blockIdx.y * TILE + ty;
const int colStart = blockIdx.x * TILE * COARSE + tx; // Fig. 6.13 第 13 行
float p0 = 0.0f, p1 = 0.0f, p2 = 0.0f, p3 = 0.0f; // 代替 Pvalue[4]
for (int ph = 0; ph < (K + TILE - 1) / TILE; ++ph) {
const int aCol = ph * TILE + tx;
const int bRow = ph * TILE + ty;
// ---- M tile：每阶段只搬一次 —— 粗化省下的就是这里 ----
Mds[ty][tx] = (Row < M && aCol < K) ? A[Row * K + aCol] : 0.0f;
for (int c = 0; c < COARSE; ++c) { // 粗化循环
const int col = colStart + c * TILE;
// ---- N tile：随列变，每个 c 各搬一次 ----
Nds[ty][tx] = (bRow < K && col < N) ? B[bRow * N + col] : 0.0f;
__syncthreads(); // ① 写后读
float acc = 0.0f;
#pragma unroll
for (int k = 0; k < TILE; ++k)
acc += Mds[ty][k] * Nds[k][tx]; // 内层和 K3 一字不差
__syncthreads(); // ② 读后写
// 寄存器轮转：p0 加上本轮结果后排到队尾
const float t = p0 + acc;
p0 = p1; p1 = p2; p2 = p3; p3 = t;
}
}
if (Row < M) {
const int c0 = colStart, c1 = c0 + TILE, c2 = c1 + TILE, c3 = c2 + TILE;
if (c0 < N) C[Row * N + c0] = alpha * p0 + beta * C[Row * N + c0];
if (c1 < N) C[Row * N + c1] = alpha * p1 + beta * C[Row * N + c1];
if (c2 < N) C[Row * N + c2] = alpha * p2 + beta * C[Row * N + c2];
if (c3 < N) C[Row * N + c3] = alpha * p3 + beta * C[Row * N + c3];
}
}
粗化后，一个 block 横跨 TILE × COARSE = 128 列。
接着看 §4·二 的例子：
K3 的 block(2,3) 只管行 96–127 × 列 64–95；
K3c 的 block(0,3) 管列 0–127 共四个 tile，
其中第三块就是 K3 的那个 block。
Nds。
一个线程每阶段搬的格数变多了（2 → 5），但它负责的输出从 1 个变成 4 个。 相除后，每个输出摊到的格数从 2 降到 1.25， 算术强度从 7.969 升到 12.720，涨到 1.6 倍。
又是同一个剧本
字节账涨到 1.6 倍，实测只快 1.082×。 因为条数账几乎没动（2.062 → 2.039）。 K3c 省的是"搬进 shared memory"那一侧；内层循环里每个 FMA 仍然要发 1 条 32 路 + 1 条广播，一条都没少。 K3 → K3c 是整条链上唯一一次"只动字节账"的干净实验，结果只快 8%。
顺带：书上 Fig. 6.13 的写法会掉进 local memory
书里用数组 float Pvalue[COARSE_FACTOR]，
上面的代码改用 p0…p3 四个标量加轮转。
数组下标 c 在运行期变化，
而寄存器不能动态寻址，编译器只好把 Pvalue 放进 local memory。
它物理上就在显存里。
于是内层每做一次
Pvalue[c] +=，就要多一次片外读和一次片外写
，本来只该碰两条片内 LDS 的一行，凭空多出两条走显存的线。
换成标量后，累加全在寄存器里完成，ptxas 报 0 bytes spill。
// ---------------------------------------------------------------------------
__global__ void sgemm_k4(int M, int N, int K, float alpha,
const float *A, const float *B, float beta, float *C) {
__shared__ float Mds[TILE][TILE];
__shared__ float Nds[TILE][TILE];
const int tx = threadIdx.x, ty = threadIdx.y;
const int Col = blockIdx.x * TILE + tx; // 列：tx 连续
const int i0 = ty * CF; // tile 内第 i0..i0+7 行归我
const int Row0 = blockIdx.y * TILE + i0; // 全局第 Row0..Row0+7 行
float a0 = 0.0f, a1 = 0.0f, a2 = 0.0f, a3 = 0.0f,
a4 = 0.0f, a5 = 0.0f, a6 = 0.0f, a7 = 0.0f;
for (int ph = 0; ph < (K + TILE - 1) / TILE; ++ph) {
const int aCol = ph * TILE + tx;
// ---- 协同加载：128 个线程 × 每人 16 格 = 两块 32×32 的 tile ----
for (int r = 0; r < CF; ++r) {
const int row = Row0 + r;
const int bRow = ph * TILE + i0 + r;
Mds[i0 + r][tx] = (row < M && aCol < K) ? A[row * K + aCol] : 0.0f;
Nds[i0 + r][tx] = (bRow < K && Col < N) ? B[bRow * N + Col ] : 0.0f;
}
__syncthreads(); // ① 写后读
// ---- 内层：1 个 N 进寄存器，用 8 次 ----
for (int k = 0; k < TILE; ++k) {
const float n = Nds[k][tx];
a0 += Mds[i0 ][k] * n;
a1 += Mds[i0 + 1][k] * n;
a2 += Mds[i0 + 2][k] * n;
a3 += Mds[i0 + 3][k] * n;
a4 += Mds[i0 + 4][k] * n;
a5 += Mds[i0 + 5][k] * n;
a6 += Mds[i0 + 6][k] * n;
a7 += Mds[i0 + 7][k] * n;
}
__syncthreads(); // ② 读后写
}
if (Col < N) {
if (Row0 < M) C[(Row0 ) * N + Col] = alpha * a0 + beta * C[(Row0 ) * N + Col];
if (Row0 + 1 < M) C[(Row0 + 1) * N + Col] = alpha * a1 + beta * C[(Row0 + 1) * N + Col];
if (Row0 + 2 < M) C[(Row0 + 2) * N + Col] = alpha * a2 + beta * C[(Row0 + 2) * N + Col];
if (Row0 + 3 < M) C[(Row0 + 3) * N + Col] = alpha * a3 + beta * C[(Row0 + 3) * N + Col];
if (Row0 + 4 < M) C[(Row0 + 4) * N + Col] = alpha * a4 + beta * C[(Row0 + 4) * N + Col];
if (Row0 + 5 < M) C[(Row0 + 5) * N + Col] = alpha * a5 + beta * C[(Row0 + 5) * N + Col];
if (Row0 + 6 < M) C[(Row0 + 6) * N + Col] = alpha * a6 + beta * C[(Row0 + 6) * N + Col];
if (Row0 + 7 < M) C[(Row0 + 7) * N + Col] = alpha * a7 + beta * C[(Row0 + 7) * N + Col];
}
}
继续用 §4·二 的 block(2,3)：行 96–127 × 列 64–95。
K3 里线程 (tx=5, ty=7) 只算 C[103][69] 一格；
K4 把同一列上的 8 格 C[96..103][69] 交给一个线程。
B[k][69] 就被读了 8 遍；
K4 交给 1 个线程，读一次进寄存器，供 8 次乘加。
每 8 个 FMA 的 shared load：8 条 32 路 + 8 条广播 → 1 条 32 路 + 8 条广播。
K3c 和 K4 是同一招的两个方向
横着排（同一行、不同列）共用 A；竖着排（同一列、不同行）共用 B。 K3c 选横向，把省下的 A 留在 shared memory 里； K4 选纵向，把省下的 B 留在 寄存器 里。 差别不在方向，而在省下的那份数据放在哪一层。 两者分别快 1.082× 和 3.359×。
一条 load 是整个 warp 的指令：32 个线程各把自己的 tx 代进下标，算出 32 个地址。
这 32 个地址的分布决定这条指令便宜还是贵。
K4 砍掉的那 7 条全是贵的一类。
判定法
看 shared memory 下标里 tx 前面乘了什么：
没有 tx → 广播，硬件取一次、复制给全 warp，
本批实测单价 0.534 周期；系数是 1 → 干净的 32 路，32 个 bank 各服务一个地址，
单价 2.686 周期；系数 ≥ 2 → 开始撞 bank，更贵。
和 §3·二 的合并判据是同一套，
只是对象不同：那边判 global memory 的事务数，这边判 shared memory 的 bank 占用。
K3 每 8 个 FMA：8 条广播 + 8 条 32 路 = 8×0.534 + 8×2.686 = 25.76 周期
K4 每 8 个 FMA：8 条广播 + 1 条 32 路 = 8×0.534 + 1×2.686 = 6.96 周期
↓ 3.70×
// 注意广播那一列：K4 和 K3 按比例算是同一个数（bcast 都是 1.000），一点没省。
// 那 8 条广播不是 K4 的额外开销 —— K3 是"每配 1 个 FMA 就再搭一条贵的 32 路"，
// K4 是"8 个 FMA 才搭一条"。省下来的全在 32 路那一侧。
tile 还是 32×32，但 block 只剩 32×4 = 128 个线程。
32 行按 ty 切成 4 条，每条 8 行。
线程 (ty, tx) 算自己那一条里第 tx 列的 8 个输出；
搬运时也在这一条里，沿第 tx 列给 Mds 和 Nds 各搬 8 格。
为什么是"每线程管一列上的 8 行"，而不是"一行上的 8 列"
因为 warp 按 threadIdx.x 连号切分。
block 为 32×4 时，一个 warp 是同一个 ty 的 32 个线程，
tx 从 0 到 31。所有 global memory 地址都以 tx 结尾，32 个连号，全部合并。
沿行方向粗化就得把 block 改成 4×32，
这时一个 warp 跨 8 个 ty，合并被劈成 8 段。
而且那样共用的是 A，A 本来就走广播，不值得省。
算术强度没变（7.969 → 7.969），dist 从 1.000 降到 0.125。
节拍从 3.34 降到 0.99，实测快 3.359×。
节拍比上面的 shared memory 账多出约 0.12 周期，是 global memory load 那一项（每个 FMA 0.0625 条 × 1.94 周期）：3.22 + 0.12 = 3.34，0.87 + 0.12 = 0.99。这是整条链上第二大的一步（仅次于 K2），字节账全程旁观。
// ---------------------------------------------------------------------------
__global__ void sgemm_k5(int M, int N, int K, float alpha,
const float *A, const float *B, float beta, float *C) {
__shared__ float Mds[TILE][TILE]; // 4 KB
__shared__ float Nds[TILE][TILE * TN]; // 16 KB
const int tx = threadIdx.x, ty = threadIdx.y;
const int i0 = ty * CF;
const int Row0 = blockIdx.y * TILE + i0;
const int colBase = blockIdx.x * TILE * TN;
float acc[CF][TN];
#pragma unroll
for (int r = 0; r < CF; ++r)
#pragma unroll
for (int j = 0; j < TN; ++j) acc[r][j] = 0.0f;
for (int ph = 0; ph < (K + TILE - 1) / TILE; ++ph) {
const int aCol = ph * TILE + tx;
// ---- 搬 A 的 tile：每线程 8 格（和 K4 一样）----
#pragma unroll
for (int r = 0; r < CF; ++r) {
const int row = Row0 + r;
Mds[i0 + r][tx] = (row < M && aCol < K) ? A[row * K + aCol] : 0.0f;
}
// ---- 搬 B 的 tile：每线程 8 行 × 4 段 = 32 格，段内仍以 tx 结尾 → 合并 ----
#pragma unroll
for (int r = 0; r < CF; ++r) {
const int bRow = ph * TILE + i0 + r;
#pragma unroll
for (int j = 0; j < TN; ++j) {
const int col = colBase + j * TILE + tx;
Nds[i0 + r][j * TILE + tx] =
(bRow < K && col < N) ? B[bRow * N + col] : 0.0f;
}
}
__syncthreads(); // ① 写后读
// ---- 内层：一次 8×4 的外积 ----
#pragma unroll
for (int k = 0; k < TILE; ++k) {
float m[CF], n[TN];
#pragma unroll
for (int r = 0; r < CF; ++r) m[r] = Mds[i0 + r][k]; // 8 条广播
#pragma unroll
for (int j = 0; j < TN; ++j) n[j] = Nds[k][j * TILE + tx]; // 4 条 32 路
#pragma unroll
for (int r = 0; r < CF; ++r)
#pragma unroll
for (int j = 0; j < TN; ++j)
acc[r][j] += m[r] * n[j]; // 32 个 FMA
}
__syncthreads(); // ② 读后写
}
#pragma unroll
for (int r = 0; r < CF; ++r) {
const int row = Row0 + r;
if (row >= M) continue;
#pragma unroll
for (int j = 0; j < TN; ++j) {
const int col = colBase + j * TILE + tx;
if (col < N) C[row * N + col] = alpha * acc[r][j] + beta * C[row * N + col];
}
}
}
K4 只把 n 读进寄存器，所以只能在一个方向上复用。
K5 两个方向一起复用：每步 k 先读 8 个 m 和 4 个 n，
再做一次 8×4 的外积。
12 条 load 换 32 个 FMA。
m 沿自己那一行横着用 4 次，每个 n 沿自己那一列竖着用 8 次。
这就是"二维"的含义。
32 个累加器全程留在寄存器里，128 个阶段跑完才写回 C 一次。
acc[7][2] 对应的就是 §4·二 的那一格 C[103][69]。
归一到每 32 个 FMA： K3 32 条广播 + 32 条 32 路 ≈ 103.0 周期 K4 32 条广播 + 4 条 32 路 ≈ 27.8 周期 K5 8 条广播 + 4 条 32 路 ≈ 15.0 周期 ← 两个方向一起省 // K4 砍的是 32 路那一类（32 → 4）；K5 接着砍广播那一类（32 → 8）。 // 到这一步两类 load 都已经很少，下一堵墙是 FP32 管线本身的 0.25 周期。
block 还是 32×4 = 128 个线程，没有改；变的是它负责的输出块大小。
输出块从 32×32 拉成 32×128，于是覆盖 C[103][69] 的 block 变成
block(0,3)，正是 §5·二 里 K3c 的那个 block。
tx"，32 个线程连号。
连着取，内层读 Nds 会变成跨步 4，
只用到 1/4 的 bank，撞成 4 拍。
右边是同一个 block 在 A、B、C 上各占的窗口：
阶段循环只推进 K 方向，M 由 blockIdx.y 覆盖，N 由 blockIdx.x 覆盖。
这一步是整条链上唯一一次两本账一起改善
B 的搬运量变成 4 倍，输出也变成 4 倍，两者打平；A 一格没多。
于是每个输出摊到的格数从 2.00 降到
(1 + TN) / TN = 1.25，算术强度从 7.969 升到 12.720。
K4 那一步只动条数账，
K3c 那一步只动字节账，只有 K5 两样都动了。
tile 横向拉长后，搬运分工也跟着变：A 那边一格没多（还是 8 格）， B 那边从 8 格变成 32 格。 每线程每阶段搬 40 格，是 K4（16 格）的 2.5 倍； 负责的输出则是 K4 的 4 倍。
__syncthreads() 也不能省。
每线程每阶段 K4 K5 为什么 搬 A 8 格 8 格 输出块的行数没变，还是 32 行 搬 B 8 格 32 格 输出块从 32 列拉成 128 列，B 的 tile 跟着大 4 倍 合计 16 格 40 格 搬运量多 2.5 倍 负责的输出 8 个 32 个 但输出多了 4 倍 摊到每个输出 2.00 格 1.25 格 ← 省下的就是 A 的那 3/4 算术强度 7.969 12.720 // 和 K3c 一模一样 —— 殊途同归
| kernel | 每线程寄存器 | 占用率 | GFLOP/s | 比上一级快 |
|---|---|---|---|---|
| K4 | 54 | 75.0% | 20770.6 | — |
| K5 | 72 | 33.3% | 34506.4 | 1.661× |
| K6 | 111 | 33.3% | 40928.8 | 1.186× |
下面逐行读代码，仍用同一个例子：block(0,3) 的线程
(ty=0, tx=5)，第 ph=1 阶段。
注释里的数字都属于这个线程：
tid = 5、行 96..103、列 20..23。
下面三张图也用这组数字。
01 __global__ void sgemm_k6(int M, int N, int K, float alpha,
02 const float *A, const float *B, float beta, float *C) {
03
04 __shared__ __align__(16) float Ads[TILE][TILE]; // 转置存：Ads[k][行]（K5 是 Mds[行][k]）
05 __shared__ __align__(16) float Bds[TILE][TILE * TN]; // 尺寸和 K5 的 Nds 一样，但列的排法变了
06 // __align__(16)：128 位读写必须 16 字节对齐
07 const int tx = threadIdx.x, ty = threadIdx.y; // block 仍是 32×4
08 const int tid = ty * TILE + tx; // 新增：0..127，加载改成按"第几个 float4"分工
09 const int i0 = ty * CF; // 我在 tile 里的行段起点：ty=0 → 0（同 K5）
10 const int Row0 = blockIdx.y * TILE + i0; // 我的 8 个输出行：96..103（同 K5）
11 const int colBase = blockIdx.x * TILE * TN; // 本 block 的 128 列起点：0（同 K5）
12 const int col0 = colBase + tx * TN; // 变了：我的 4 列是"连着的"20..23
13 // K5 是 5/37/69/101，每隔 32 列一个
14 float acc[CF][TN]; // 32 个累加器 —— 和 K5 一字不差
15 #pragma unroll
16 for (int r = 0; r < CF; ++r)
17 #pragma unroll
18 for (int j = 0; j < TN; ++j) acc[r][j] = 0.0f;
19
20 const int phases = (K + TILE - 1) / TILE;
21 for (int ph = 0; ph < phases; ++ph) { // 阶段循环：和 K5 逐字相同
22
23 // ---- 搬 A：整块 32×32 = 256 个 float4，128 个线程各 2 条 ----
24 #pragma unroll
25 for (int s = 0; s < (TILE * TILE / 4) / 128; ++s) { // = 2 轮
26 const int idx4 = tid + s * 128; // 我是第几个 float4：5、133
27 const int r = idx4 / (TILE / 4); // tile 内行号 = idx4/8 → 0、16
28 const int c = (idx4 % (TILE / 4)) * 4; // tile 内列号 = 20（必是 4 的倍数）
29 const int row = blockIdx.y * TILE + r; // 全局行 96、112 —— 不按我的输出行分
30 const int aCol = ph * TILE + c; // 全局列 52
31 float4 a4 = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
32 if (row < M && aCol + 3 < K) // 整条一起判：不能只读半条 float4
33 a4 = *reinterpret_cast<const float4 *>(&A[row * K + aCol]);
34 Ads[c + 0][r] = a4.x; Ads[c + 1][r] = a4.y; // 转置写：读进来是横着的 4 个，
35 Ads[c + 2][r] = a4.z; Ads[c + 3][r] = a4.w; // 写出去是竖着的 4 个 → 4 条 store
36 }
37
38 // ---- 搬 B：整块 32×128 = 1024 个 float4，128 个线程各 8 条 ----
39 #pragma unroll
40 for (int s = 0; s < (TILE * TILE * TN / 4) / 128; ++s) { // = 8 轮
41 const int idx4 = tid + s * 128; // 5、133、261 … 901
42 const int r = idx4 / (TILE * TN / 4); // 行号 = idx4/32 → 0、4 … 28
43 const int c = (idx4 % (TILE * TN / 4)) * 4; // 列号 = tx*4 = 20（恒定）
44 const int bRow = ph * TILE + r; // 全局行 32、36 … 60
45 const int col = colBase + c; // 全局列 20 = 我自己的 col0
46 float4 b4 = make_float4(0.0f, 0.0f, 0.0f, 0.0f);
47 if (bRow < K && col + 3 < N)
48 b4 = *reinterpret_cast<const float4 *>(&B[bRow * N + col]);
49 *reinterpret_cast<float4 *>(&Bds[r][c]) = b4; // 直着写：一条 128 位 store
50 }
51
52 __syncthreads(); // ① 写后读：和 K5 同一道屏障，位置、作用都没变
53
54 #pragma unroll
55 for (int k = 0; k < TILE; ++k) { // 内层 32 步 —— K6 的全部收益只发生在这里
56 const float4 m0 = *reinterpret_cast<const float4 *>(&Ads[k][i0 ]); // 广播：行 96..99
57 const float4 m1 = *reinterpret_cast<const float4 *>(&Ads[k][i0 + 4]); // 广播：行 100..103
58 const float4 n0 = *reinterpret_cast<const float4 *>(&Bds[k][tx * TN]); // 32 路：列 20..23
59 const float m[CF] = { m0.x, m0.y, m0.z, m0.w, // 拆成 8 个标量：纯粹给寄存器改名，
60 m1.x, m1.y, m1.z, m1.w }; // 编译后不产生任何指令
61 const float n[TN] = { n0.x, n0.y, n0.z, n0.w };
62 #pragma unroll
63 for (int r = 0; r < CF; ++r)
64 #pragma unroll
65 for (int j = 0; j < TN; ++j)
66 acc[r][j] += m[r] * n[j]; // 还是那 32 个 FMA —— 但这次只用了 3 条 load
67 }
68
69 __syncthreads(); // ② 读后写
70 }
71
72 #pragma unroll
73 for (int r = 0; r < CF; ++r) { // 写回：8 行，每行一条 float4
74 const int row = Row0 + r;
75 if (row >= M || col0 + TN > N) continue; // 整条要么全写要么不写 → 边界粒度变粗成 4
76 const float4 old = *reinterpret_cast<const float4 *>(&C[row * N + col0]); // beta 的读也向量化
77 float4 out;
78 out.x = alpha * acc[r][0] + beta * old.x;
79 out.y = alpha * acc[r][1] + beta * old.y; // 四个分量对应 acc 的 j=0..3
80 out.z = alpha * acc[r][2] + beta * old.z; // 顺序不能错：acc[r][j] ↔ C 的第 col0+j 列
81 out.w = alpha * acc[r][3] + beta * old.w;
82 *reinterpret_cast<float4 *>(&C[row * N + col0]) = out; // 一条 128 位写
83 } // 共 8 条写（K5 是 32 条），字节数一样多
84 }
一条 128 位 load 一次搬 4 个 float。 K5 每步 k 发 12 条 32 位读，K6 发 3 条 128 位读。 两者搬的都是 48 个字节，换来的都是 32 个 FMA。
K6 和 K5 的区别，看这七行就够了。
行的算法一字未改，K6 相对 K5 的改动只发生在列方向。
变的只有两处：多了一个 tid；col0 的算法变了。
tid 让搬运改成按"第几个 float4"分工；
col0 把输出列从"隔 32 一个"并成"连着的 4 个"。
后面所有向量化都建立在这一步上。
也因为这一步，§4·二 那一格 C[103][69] 换了主人。
这一改的方向，和 K5 当初的要求正好相反
K5 特意让 4 列"每隔 32 取一个"，为的是避开读 Nds 时的 bank 冲突。
K6 必须改回连着的 4 列，否则 n[0..3] 在 shared memory 里不连续，合不成一条 float4。
同一件事，两版的最优解相反：
K5 守的是 32 位访问的 bank 规则，K6 守的是 128 位访问的对齐规则。
这类「相反要求」的实际代价只能靠 ncu 量，推理定不了量。
向量化只有一个前提：要用的那几个数在内存里必须连着。
K5 的 Mds[行][k] 不满足这一条。
所以 K6 不能只"加个 float4"，要先改数据的摆放。
Mds[行][k] 存，它们每隔 32 个 float 才出现一次；
转置成 Ads[k][行] 后才连成一片，2 条 float4 读完。
代价是搬进 shared memory 时要横读竖写，这一步会撞 bank。
怎么撞的、有哪三条退路，见 §8·五；
它在反汇编里留下的痕迹（STS 从 40 条降到 16 条）见 §8·八。
K4、K5 的搬运分工都是"我搬我自己那几行 / 那几列"，K6 不是。
向量化要求连续 4 个 float 归同一个线程，并且地址 16 字节对齐。
所以 K6 改成按 tid 线性铺 float4。
一个线程搬的数据，和它自己要算的输出不再对应。
这种分法最容易数清楚：给 tile 里每一条 float4 编个号，128 个线程按号领。
float4，128 个线程各领 2 条。
三步映射每一步都是一一对应（编号 → 位置 → Ads 的格子）。
所以不会漏掉一格，也不会有两个线程撞在同一格上。
不用逐个去想每个线程搬哪几格：
只要总数对得上、映射是一一的，覆盖就自动成立。
Ads[c+i][r] 时 bank = r mod 32，而一个 warp 只有 4 个不同的 r，
所以 32 个线程挤在 4 个 bank 上。
B 侧没有这个问题：一个 warp 的 32 个 lane 各拿 4 列，
合起来正好是一整行 512 字节。
转置写入撞 bank：这是 K6 最可能出问题的地方，本报告没有量
根因是两个方向垂直：global memory 那边要"按行连续读"才能合并，shared memory 这边要"按列写"才能转置。 这和 §4 里 corner turning 那对矛盾是同一类。
// 三条候选退路，按省事程度排 A padding：__shared__ float Ads[TILE][TILE + 4]; // 必须加 4 不是 1 —— 加 1 行首就不再 16 字节对齐，内层那两条 float4 读会挂。 // 加 4 之后 bank = (4(c+i) + r) mod 32，冲突从 8 路降到约 4 路，没有完全消除。 B 换 A 的加载分工，让一个 warp 的 r 取 32 个不同值 // 但那样全局读就变成跨步 K，直接丢掉合并 —— 得不偿失。 C A 侧干脆不转置、不向量化，内层退回 8 条 32 位广播读，只向量化 B 和写回 // 广播本来就便宜（0.534 周期），A 侧向量化的收益最小 —— 这是最省事的退路。
不要先猜，先量：
ncu --metrics l1tex__data_bank_conflicts_pipe_lsu_shared.sum
K5、K6 各跑一次。如果 K6 明显更高，再按 A → C 的顺序试。
这台机器上 ncu 起不来，所以上面整段是推断，没有证据。
K6 的全部收益只发生在这 14 行里。
前面的改造（转置、连续 4 列、线性 tid）都只为让这三条读成立。
Ads 的第 k 行取 2 条（广播），Bds 的第 k 行取 1 条（32 路），共 48 字节。
第 59–61 行那两句"拆成标量"看着像在拷贝 12 个数，但不产生任何指令。
float4 的四个分量本来就在 4 个连号寄存器里，
这两句只是给它们起名字，方便写下面的双重循环。
前提仍是两层 #pragma unroll：
不展开，m[r] 的下标就不是常量，一样会掉进 local memory。
写回除了 alpha、beta 的缩放，没有别的计算，只是把寄存器里已经算完的 32 个数，按初始化时定好的坐标送回去。
所以这一段最容易出的错不是算错，而是坐标公式和前面三处
（搬 B、读 Bds、写回）不一致。不一致就会写到别人的格子上。
out 的四个分量 → 一条 128 位写。
顺序不能错：acc[r][j] 对应 C 的第 col0+j 列。
下半是同一个线程在两版里负责的格子：
K5 是四条细竖线，K6 是一条 4 列宽的竖带，合起来覆盖的 128 列完全一样。
注意边界粒度：K6 不能只写半条 float4，
判的是 col0 + 4 > N。这就是它要求 N % 4 == 0 的原因。
没算的一笔小账：beta 的那次读
C = alpha*acc + beta*C 里的 beta*C 要先把 C 读回来。
所以写回其实是"读 → 改 → 写"，每个输出是 4 B 读 + 4 B 写。
本报告的 bytesPerOut 只算了写的那 4 B。
漏掉的 4 B 摊到 644 B 上不到 1%；但以后如果跑 beta ≠ 0 的基准，要补上这一笔。
这也说明了"写回只发生一次"的价值：
128 个阶段里 C 只被碰这一次。如果每个阶段都往 C 里累加，这 8 B 就要乘 128 倍，
算术强度会被打回原形。
"float4 编成了 128 位"是一条纯靠推理、没有直接证据的断言。
编译器完全可以把 float4 拆成 4 条 32 位 load，
而 ptxas -v 只报寄存器和 spill，看不到指令形态。
所以用反汇编数一遍（cuobjdump -sass，静态条数）：
kernel FFMA LDS.128 LDS.32 LDG.128 LDG.32 STS STG.128 STG.32 TOTAL k5 1056 64 128 - 72 40 - 32 1856 k6 1056 96 - 18 - 16 8 - 1520 // ✓ 证实：K6 有 96 条 LDS.128 + 18 条 LDG.128，float4 确实编成了 128 位访存。 // ✗ 推翻：K5 竟然也有 64 条 LDS.128 —— 编译器自己把它三分之一的 shared load 向量化了。
第二条改写了 K5→K6 的解释
原来的说法是"标量 vs 向量化"，实际是"部分向量化 vs 完全向量化"。
访存指令合计 K5 = 192 + 72 = 264 条 K6 = 96 + 18 = 114 条 比值 2.32× // 而手推的条数账说 load/FMA 从 0.414 降到 0.104，比值 3.98× —— 高估了。 // 根源就是手推时假设 K5 全是标量 load，没料到编译器已经替它向量化了三分之一。 // 两者 FFMA 条数完全相同（1056），说明展开程度一样，静态条数可以直接比。
另有一个没预测到的旁证：K1 和 K2 的 SASS 逐项相同 （都是 30 FFMA / 59 LDG.32 / 216 条总计）。 两者只差一个索引映射，指令构成完全一样。 所以 K1→K2 那 8.209× 里， 指令条数的贡献严格为零，100% 来自访存行为。 这是"请求账必须独立于条数账"最干净的证据。
| 这一步 | 字节账 算术强度 |
请求账 请求/FMA |
条数账 load/FMA | 实测 |
|---|---|---|---|---|
| K1→K2 | 不变 0.25 | 33 → 2 | 不变 2.000 | 8.209× |
| K2→K3 | 0.25 → 7.97 | 不变 | 2.000 → 2.062 | 1.162× |
| K3→K3c | 7.97 → 12.72 | 不变 | 2.062 → 2.039 | 1.082× |
| K3→K4 | 不变 7.97 | 不变 | 2.062 → 1.188 | 3.359× |
| K4→K5 | 7.97 → 12.72 | 不变 | 1.188 → 0.414 | 1.661× |
| K5→K6 | 不变 12.72 | 不变 | 0.414 → 0.104 | 1.186× |
把六级压成六句
K1 索引映射写反，一条 LDG 被拆成 32 个请求。
K2 把映射换回来，快 8.209 倍：其中 2.392 倍来自合并，3.43 倍来自顺带治好的 2 的幂跨度冲突。
K3 用 shared memory 分块。算术强度涨 32 倍，但 K2 本来就不在带宽墙上，所以只快 16%；
真正的收益来自 LDG 换成了更便宜的 shared memory load。
K3c 是书上的粗化，只动字节账的干净实验，快 8%。
K4 让一个 N 元素进寄存器后用 8 次。算术强度没变，dist 降到 1/8，快 3.359 倍。
K5 两个方向一起复用。占用率砍到 33.3%，反而快 1.661 倍。
K6 用 float4。字节不变，条数 ÷4，快 1.186 倍。但 SASS 显示编译器本来就把 K5
向量化了三分之一，所以真实的条数比只有 2.32 倍，不是 4 倍。
最终：K1 → K6 共快 63.1 倍，做到 cuBLAS 的 71.5%。 siboehm 在 A6000 上做到 93%。差的 21 个点主要是 warptiling 和 double buffering。
| 要什么 | 去哪 |
|---|---|
| 八个 kernel 的完整源码 + 测试脚手架 | src/sgemm_all.cu（含 K0 cuBLAS 基准） |
| 测试脚本 | Scripts/ —— run_p0.sh 一条命令跑完全套 |
| 原始输出与解析后的 CSV | Results/Round1/ |
| 八个指标各自怎么算出来的 | 报告 02 · §3 |
| 形状扫描（prefill / decode） | 报告 02 · §5 |
这份报告的代码从哪来
K1 出自 PMPP 第 3 章 Fig. 3.11，K2 是第 4 章的修正， K3 是第 5 章的分块，K3c 是第 6 章 Fig. 6.13。 K4 / K5 / K6 书上没有，是照着 siboehm 那篇的思路自己写的。 所有 kernel 合并进同一个可执行文件，在同一次开机里跑。 因为"K6 是 cuBLAS 的百分之多少"这类比值，只有在同一个二进制下才成立。