SGEMM 实测与三本账：八个指标怎么算出来的

一、三句话↑

二、cuBLAS 这个分母怎么定的↑

K0 = 2.401 ms / 57242.4 GFLOP/s / 占峰值 69.32%，抽样校验通过。 两项设置让它只跑纯 FP32：源码里设 cublasSetMathMode(CUBLAS_PEDANTIC_MATH)，环境变量设 NVIDIA_TF32_OVERRIDE=0。 如果 TF32 没禁掉，tensor core 会让结果超过 FP32 峰值。所以 69.32% 只可能是纯 FP32 的成绩。

// cuBLAS 是列主序，而我们的矩阵是行主序。不想转置数据，就转置运算：
//   行主序的 C = A·B   等价于   列主序的 Cᵀ = Bᵀ·Aᵀ
cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
            N, M, K, &alpha,
            B, N,          // ← B 在前
            A, K, &beta,   // ← A 在后
            C, N);

恒等式和那行调用的关系，要从一个事实讲起： 同一段字节，按行主序读和按列主序读，得到的矩阵互为转置。 所以这里的"转置"不是一个操作。

所以恒等式可以直接当参数表读：

             Cᵀ    =    Bᵀ    ·    Aᵀ
             ↓          ↓          ↓
        写进 C 缓冲区  传 B 指针   传 A 指针

// 右边从左往右念一遍 —— Bᵀ 在前、Aᵀ 在后 —— 就是参数的先后顺序。
// 而 Bᵀ 和 Aᵀ 根本不用构造：把 B、A 的指针原样递过去，
// cuBLAS 用列主序一读，看到的就已经是 Bᵀ 和 Aᵀ 了。

维度和 ld 也自动对上。 cublasSgemm(h, …, m, n, k, …) 算的是 (m×k) · (k×n)。这里的结果是 N×M 的 Cᵀ， 所以 m=N, n=M，中间那一维仍是 k=K。 参数顺序 N, M, K 就是这样来的。

参数            值      为什么
─────────────────────────────────────────────────────────
m               N       结果 Cᵀ 是 N×M，它的行数是 N
n               M       Cᵀ 的列数是 M
k               K       中间那一维没动
第一个矩阵       B       cuBLAS 读成 Bᵀ（N×K），正好是 m×k ✓
它的 ld         N       Bᵀ 在列主序下有 N 行 → ld = N
第二个矩阵       A       cuBLAS 读成 Aᵀ（K×M），正好是 k×n ✓
它的 ld         K       Aᵀ 有 K 行 → ld = K
输出            C       cuBLAS 写进去一个 N×M 的列主序矩阵
它的 ld         N       同上

// 最后一步：cuBLAS 往 C 缓冲区里写的是列主序的 Cᵀ（N×M）；
// 我们回过头按行主序去读同一段字节，看到的就是 M×N 的 C。闭环。

一、六条规矩↑

二、两个口径，必须一起换↑

FP32 峰值   128 SM × 128 core × 2 × 2.52 GHz  = 82,575 GFLOP/s   // 设备自报的 clockRate ← 主表用这个
            128 SM × 128 core × 2 × 2.64 GHz = 86,508 GFLOP/s   // 本批 1 Hz 采样的中位数（稳态 2625–2655）

显存带宽    10501 MHz × 2（DDR 双沿）× 384 bit ÷ 8 = 1008 GB/s

主表的"占 FP32 峰值"按 2.52 GHz 那一行算。 这个时钟来自 cudaDeviceProp，换一台 4090 也读到同一个数。

标称值的意思是"保证不低于"。实际频率由 GPU Boost 按功耗墙（450 W）和温度动态调整，本批稳定在 2640 MHz。 上限是 3105，但稠密 FP32 GEMM 功耗太高，达不到。 要按实测口径引用，把整列乘 0.955（K6 从 49.57% 变成 47.32%）。 换口径时，分子和分母要一起换。 这个时钟每台机器不同，不能跨批次沿用。

"占 cuBLAS"那一列不受口径影响： 分子和分母都是 GFLOP/s，峰值约掉了。 所以"K6 = cuBLAS 的 71.5%" 在哪套口径下都成立，这也是本报告把它当主结论的原因。

一、一个 warp 读 shared memory 有两种姿势↑

一条 load 是整个 warp 的指令，它的代价由 32 个地址的分布决定。 32 个线程同时执行这条指令，各自把自己的 tx 代进下标，算出一个地址。 所以"一条 load"实际上是 32 个地址。

二、主例 K3：从代码一步步数出全部八个↑

K3 的内层循环只有两行，八个指标全从这两行数出来：

for (int k = 0; k < TILE; ++k)
    acc += Mds[ty][k] * Nds[k][tx];
//         ↑ 下标没有 tx       ↑ 下标含 tx 且系数为 1
//           = 广播               = 32 路

第 1 步 · dist（32 路条数 ÷ FMA）

循环体里含 tx 的 shared 读：Nds[k][tx] —— 1 条
一次循环做几个 FMA：       1 个

dist = 1 条 ÷ 1 个 FMA = 1.000

除以 FMA 数，是为了把不同 kernel 放到同一把尺子上。 不同 kernel 一次循环做的 FMA 个数不同 （K3 做 1 个，K5 做 32 个）。除以之后，指标的含义统一为 "每做一次乘加，平均要发几条这类 load"。 这个商没有量纲，可以横向比较，这是"条数账"成立的前提。

第 2 步 · bcast（广播条数 ÷ FMA）

循环体里不含 tx 的 shared 读：Mds[ty][k] —— 1 条
bcast = 1 条 ÷ 1 个 FMA = 1.000

第 3 步 · ldg（全局 load 条数 ÷ FMA）

一个阶段里，每个线程从 global 搬 2 个元素进 shared（A 一个、B 一个）—— 2 条 LDG
同一个阶段里，内层循环跑 TILE = 32 步，每步 1 个 FMA        —— 32 个 FMA

ldg = 2 条 ÷ 32 个 FMA = 0.062

三步合起来：load/FMA = dist + bcast + ldg = 1.000 + 1.000 + 0.062 = 2.062。

第 4 步 · bytesPerOut（每个输出要过多少字节）

一个 block 算 32×32 = 1024 个输出，沿 K 走 N/32 = 128 个阶段
每个阶段搬进来：A 的 32×32 + B 的 32×32 = 2048 个 float = 8192 字节
128 个阶段合计：128 × 8192 = 1,048,576 字节
摊到 1024 个输出：1,048,576 ÷ 1024 = 1024 字节
再加上写回 C 自己那 4 字节：            1028 字节/输出

第 5 步 · 算术强度与 roofline 达成率

算术强度回答"这个 kernel 每搬一个字节，能换来多少次浮点运算"。 它把第 4 步的 bytesPerOut 换算成一个 和矩阵规模无关的比值，所以不同 kernel 可以直接比较。

算术强度 = 每个输出的浮点运算 ÷ 每个输出的字节
         = 2K ÷ bytesPerOut
         = 2×4096 ÷ 1028 = 7.969 FLOP/byte

// 分子的 2K：C 的每个元素要做 K 次乘加，一次乘加算 2 个浮点运算（1 乘 1 加）。
// 分母的 1028：第 4 步数出来的，每个输出摊到的显存字节。

roofline 上限 = 算术强度 × 峰值带宽
             = 7.969 × 1008 = 8033 GFLOP/s

// 意思是：如果这个 kernel 把显存带宽跑满，它最快只能到 8033。
// 1008 GB/s 是 4090 的理论峰值带宽，口径见 §2·二。

下面的 6184.0 是「实测吞吐」，是八步里唯一一个量出来的数； 其余七个都从代码数出来。 所以下面把它单独拆开。

实测 GFLOP/s = 2·M·N·K ÷ 耗时(秒) ÷ 1e9

// ① M·N·K 是这次乘法一共要做多少次乘加（FMA）
//    C 有 M×N 个元素，每个元素要沿 K 方向累加 K 次
   M·N·K = 4096 × 4096 × 4096 = 68,719,476,736 次 FMA

// ② 乘以 2：一次 FMA 在 FLOP 的记法里算 2 个运算（1 次乘 + 1 次加）
//    这是整个行业的统一口径 —— cuBLAS、Nsight、论文都这么算，
//    否则跟别人的数字没法比
   2 × 68,719,476,736 = 137,438,953,472 FLOP ≈ 137.44 GFLOP

// ③ 除以耗时。22.225 ms 就是主表里那个中位耗时 ——
//    和第 6 步算节拍用的是【同一个数】，同一次测量，没有第二次计时
   137,438,953,472 ÷ 0.022225 s = 6,183,979,908,751 FLOP/s

// ④ 换成 GFLOP/s（除以 1e9）
   6,183,979,908,751 ÷ 1e9 = 6184.0 GFLOP/s

达成率 = 实测 ÷ 上限 = 6184.0 ÷ 8033 = 77.0%

≈ 100%  贴着屋顶   带宽真跑满了 —— 这才是"受带宽限制"，想更快只能提高算术强度
≪ 100%  离屋顶远   带宽没跑满 → 瓶颈不是带宽，在别处（指令吞吐 / 延迟 / 占用率）
> 100%  超出屋顶   分母错了 —— 字节没全走显存，命中了 L2。模型在这里被证伪

K3 的 77.0% 没超过 100%，所以字节账在 K3 上没有被证伪，但这不等于字节账成立。 K2 的 2112.7% 才是直接证伪。

K1 257.4%   K2 2112.7%   K3 77.0%   K3c 52.2%   K4 258.6%   K5 269.1%   K6 319.2%

// 七个里没有一个贴近 100% —— 最近的就是 K3 的 77.0%。
// 也就是说：这条链上没有一个 kernel 真的受带宽限制。
// 五个被直接证伪；K3/K3c 虽然没被证伪，但离屋顶太远，同样说明带宽不是它们的墙。

K3 真正的瓶颈是指令条数。 第 8 步用条数账算出 3.340，和实测 3.34 一致； §4 再把七个 kernel 并排对比。

第 6 步 · 节拍（把实测耗时换算成周期）

节拍 = 32 × SM数 × 时钟 × 耗时(秒) ÷ (M·N·K)
     = 32 × 128 × 2.52e9 × 22.225e-3 ÷ 4096³
     = 3.338 周期/warp圈        （程序打印 3.34）

第 7 步 · 每条 load 多少周期

每条 load 周期 = 节拍 ÷ load/FMA = 3.34 ÷ 2.062 = 1.62

这是诊断量，不是预测量。 它把"这个 kernel 慢"换成"它每发一条 load 要付几个周期"。 如果这个数在 kernel 之间差很多，说明 load 的类型构成变了。 例如 K1 是 15.91：它每个 FMA 的两条 load 一共触发 33 个请求（A 32 个、B 1 个）。

第 8 步 · 模型预测

预测节拍 = 2.69 × dist + 0.53 × bcast + 1.94 × ldg
        = 2.69 × 1.000 + 0.53 × 1.000 + 1.94 × 0.062
        = 3.340        实测 3.34   ✓

K3 是标定点之一，所以这里必然吻合。 真正的检验是用同样三个单价外推 K3c 和 K5，见 §3 · 六。

三、K1 / K2：没有 shared memory 的特例，以及为什么需要请求账↑

下面按 K1 → K2 → K3c → K4 → K5 → K6 的顺序补齐其余六个 kernel，每个都走同样的八步。 K3 是主例，上面已经逐步算过。 只有 K1/K2 多一步"请求账"，因为从 K3 起访存全部合并，这本账归零，不再变化。

K1 · 九步

const int x = blockIdx.x * blockDim.x + threadIdx.x;   // 行 ← 从 .x 来
const int y = blockIdx.y * blockDim.y + threadIdx.y;   // 列 ← 从 .y 来
float acc = 0.0f;
for (int i = 0; i < K; ++i)
    acc += A[x * K + i] * B[i * N + y];                // ← 1 个 FMA + 2 条 LDG（全局！）
C[x * N + y] = alpha * acc + beta * C[x * N + y];

第 1 步  这个线程负责几个输出？
  一个 acc → 1 个输出                              // 和 K3 相同

第 2 步  内层每轮：几个 FMA、几条 shared load、各属哪类？
  acc += A[x * K + i] * B[i * N + y];
        └─ 1 个 FMA ─┘  └ 2 条 LDG ┘                    // 是 LDG（全局），不是 LDS
  一条 shared load 都没有  →  dist  = 0 ÷ 1 = 0.000
                            bcast = 0 ÷ 1 = 0.000

第 3 步  全局 load 怎么摊到每个 FMA 上？
  内层每轮 2 条全局 load、1 个 FMA，没有阶段可摊
  →  ldg = 2 条 ÷ 1 个 FMA = 2.000
  // K3 是 0.062，差 32 倍 —— 那 32 倍正是"搬一次进 shared、被 32 个线程各用一次"换来的

第 4 步  bytesPerOut：每个输出摊多少字节？
  每个输出要读 A 的一整行 + B 的一整列 = 2K = 8192 个 float
  8192 × 4 字节 = 32768，再加写回 C 的 4 字节
  →  bytesPerOut = 32772 字节                  // 代码：2.0 * K * 4.0 + 4.0

第 5 步  AI = 2K ÷ bytesPerOut
  8192 ÷ 32772 = 0.250 FLOP/byte           // 七个 kernel 里最低
  roofline 上限 = 0.250 × 1008 GB/s = 252 GFLOP/s
  达成率 = 648.6 ÷ 252 = 257.4%       // 程序打印的就是这个数

  // 648.6 是这九步里唯一一个【量】出来的数，其余全是从代码数出来的：
  //   实测 GFLOP/s = 2·M·N·K ÷ 耗时(秒) ÷ 1e9
  //   分子 2 × 4096³ = 137,438,953,472 FLOP（一次 FMA 记 2 个运算）
  //   分母就是主表那个中位耗时 211.911 ms = 0.211911 s
  //   137,438,953,472 ÷ 0.211911 ÷ 1e9 = 648.6 GFLOP/s      —— 逐项拆解见 §3·二

第 6 步  节拍（把实测耗时换算成周期）
  32 × 128 × 2.52e9 × 0.211911 ÷ 4096³ = 31.83 周期/warp圈

第 7 步  每条 load 周期
  load/FMA = 0.000 + 0.000 + 2.000 = 2.000
  31.83 ÷ 2.000 = 15.91 周期

第 8 步  模型预测
  2.69 × 0 + 0.53 × 0 + 1.94 × 2.000 = 3.880
  实测 31.83，误差 -87.8%  ← 模型在这里彻底失效
  // 换哪组单价都救不回来：差的不是系数，是这个 kernel 根本不在条数账的适用区间里

第 9 步  请求账 —— 只有 K1 / K2 走这一步
  warp 按 threadIdx.x 每 32 个切一刀 → warp 内 x 在变、y 不变。而 K1 里 x 是行、y 是列：
  A[x*K+i]   32 个线程读 32 个不同行，地址隔 K×4 字节      32 个请求
  B[i*N+y]   warp 内 y 相同，读同一地址                     1 个请求
  →  reqPerFMA = 33 个/FMA
  31.83 ÷ 33 ≈ 0.96 周期/请求
  // 换成"请求"当单位，数量级立刻就回来了。

K2 · 九步 —— 前四步逐字相同，从第 5 步的达成率开始分家

const int y = blockIdx.x * blockDim.x + threadIdx.x;   // 列 ← 改成从 .x 来
const int x = blockIdx.y * blockDim.y + threadIdx.y;   // 行 ← 改成从 .y 来
...
    acc += A[x * K + i] * B[i * N + y];                // ← 和 K1 逐字相同

第 1 步  几个输出？   1 个                                        // 同 K1
第 2 步  dist / bcast  没有 shared → 0.000 / 0.000       // 同 K1
第 3 步  ldg          2 条 ÷ 1 个 FMA = 2.000                 // 同 K1
第 4 步  bytesPerOut  2K × 4 + 4 = 32772 字节                 // 同 K1
第 5 步  AI           8192 ÷ 32772 = 0.250，上限 252 GFLOP/s   // 同 K1
             达成率 = 5324.0 ÷ 252 = 2112.7%          ← 第一个不同的数
  // 前四步是人从代码里数出来的，代码没变所以数也没变；
  // 达成率带上了实测 GFLOP/s，所以从这里开始分家。

  // 5324.0 和 K1 的 648.6 是同一个公式代不同的耗时：
  //   实测 GFLOP/s = 2·M·N·K ÷ 耗时(秒) ÷ 1e9
  //   分子一模一样：2 × 4096³ = 137,438,953,472 FLOP —— 两版要算的乘加一次不多一次不少
  //   分母才是变的：K1 211.911 ms → K2 25.815 ms
  //   137,438,953,472 ÷ 0.025815 ÷ 1e9 = 5324.0 GFLOP/s      —— 逐项拆解见 §3·二
  //   所以 5324.0 ÷ 648.6 = 8.21，正好就是耗时之比 —— 吞吐涨的全部来自"跑得快"，
  //   而不是"算得多"。这也是为什么分母 252 不变，达成率却从 257.4% 跳到 2112.7%。

第 6 步  节拍
  32 × 128 × 2.52e9 × 0.025815 ÷ 4096³ = 3.88
  // K1 是 31.83 —— 同样的 2.000 条 load/FMA，差 8.21 倍

第 7 步  每条 load 周期
  3.88 ÷ 2.000 = 1.94 周期            // K1 是 15.91

第 8 步  模型预测
  2.69 × 0 + 0.53 × 0 + 1.94 × 2.000 = 3.880   // 和 K1 同一个数
  实测 3.88，误差 +0.0%   ← 模型对 K2 是准的（它是标定点）

第 9 步  请求账 —— 差别全在这里
  warp 内 tx 在变、ty 不变。K2 把 y（列）接到 tx、x（行）接到 ty：
  A[x*K+i]   x 由 ty 定 → warp 内行相同，32 个线程读同一地址   1 个请求（广播）
  B[i*N+y]   y 由 tx 定 → 列连号，32 个地址首尾相接            1 个请求（合并）
  →  reqPerFMA = 2 个/FMA                     // K1 是 33
  3.88 ÷ 2 ≈ 1.94 周期/请求
  // 和 K1 的 0.96 同一量级，说明"按请求算钱"是对的。

              第 1–5 步（全是人算）        第 6 步   第 7 步   第 8 步      第 9 步
                                          节拍   每条load  模型/误差    请求账
K1    1 输出 · dist 0 · bcast 0 ·        31.83    15.91   3.880 / -87.8%   33
K2    ldg 2.000 · 32772 B · AI 0.250     3.88    1.94    3.880 / +0.0%    2
      ↑ 一模一样，五个量一个都没变

// 指令条数都是 2 条，但硬件真正要处理的事务数差了 16.5 倍。

四、K3c / K4：同叫"粗化"，动的却是两本不同的账↑

两个 kernel 都让一个线程多算几个输出，但省下来的东西完全不同。 判断方法是看多出来的输出摊薄了哪一项。

K3c · 每线程 4 个输出（同一行上，跨 4 个 tile）

float p0 = 0.0f, p1 = 0.0f, p2 = 0.0f, p3 = 0.0f;   // ← 4 个累加器 = 4 个输出
for (int ph = 0; ph < K / TILE; ++ph) {             // ← 128 个阶段
    Mds[ty][tx] = A[Row * K + aCol];                // ← M tile 每阶段只搬一次（1 格）
    for (int c = 0; c < COARSE; ++c) {              // ← 粗化循环，4 轮
        Nds[ty][tx] = B[bRow * N + col];            // ← N tile 随列变，4 格
        __syncthreads();
        for (int k = 0; k < TILE; ++k)
            acc += Mds[ty][k] * Nds[k][tx];         // ← 内层和 K3 一字不差！
        ...
    }
}

第 1 步  几个输出？
  p0…p3 四个累加器 → 4 个输出（同一行上，跨 4 个 tile）

第 2 步  内层每轮：几个 FMA、几条 LDS？
  acc += Mds[ty][k] * Nds[k][tx];       // 和 K3 一字不差
  1 个 FMA、2 条 LDS  →  dist = 1 ÷ 1 = 1.000   bcast = 1 ÷ 1 = 1.000
  ← 条数账一分没省，这是 K3c 的全部问题所在

第 3 步  全局 load 摊到每个 FMA？
  每阶段每线程搬 1 格 M + 4 格 N = 5 格
  每阶段做 4 个 c × 32 轮 = 128 个 FMA
  →  ldg = 5 ÷ 128 = 0.039                     // K3 是 0.062

第 4 步  bytesPerOut
  每阶段 5 格 ÷ 4 个输出 = 1.25 格/输出                 // K3 是 2 格
  128 阶段 × 1.25 × 4 字节 + 4 = 644 字节
  // 末尾那个 +4 不用再除以 4 —— 输出数和写回格数是 1:1 的：
  //   每线程 4 个输出 → 写回 4 格 = 16 字节 → 摊到每个输出仍是 4 字节。

第 5 步  AI
  8192 ÷ 644 = 12.720 FLOP/byte          // K3 是 7.969，提高 1.6 倍
  上限 = 12.720 × 1008 = 12823 GFLOP/s
  达成率 = 6690.6 ÷ 12823 = 52.2%     // AI 越高，这个百分比反而越低
  // 6690.6 从哪来 —— 程序里就一行：
  //     double gflops = 2.0 * M * N * K / (ms / 1000.0) / 1e9;
  //     = 2 × 4096³ ÷ 0.020542 s ÷ 1e9，即第 6 步那个耗时的另一种写法

第 6 步  节拍
  32 × 128 × 2.52e9 × 0.020542 ÷ 4096³ = 3.09   // K3 是 3.34

第 7 步  每条 load 周期
  load/FMA = 1.000 + 1.000 + 0.039 = 2.039
  3.09 ÷ 2.039 ≈ 1.51 周期          // K3 是 1.62

第 8 步  模型预测
  2.69 × 1.000 + 0.53 × 1.000 + 1.94 × 0.039 = 3.295
  实测 3.09，误差 +6.6%   // 它没参与标定，这是一次外推

每阶段、每个输出摊到的格数 = (1 + COARSE) ÷ COARSE
                              ↑       ↑
                        M 搬 1 格  N 搬 COARSE 格

  COARSE = 4：(1 + 4) ÷ 4 = 1.25 格        // K3 相当于 COARSE=1：(1+1)÷1 = 2 格

// 分子那个 1 就是粗化的全部收益：同一行上的 COARSE 个输出共用同一块 M tile，
// 每阶段只搬一次。不粗化的话（4 个独立的 K3 block），这一格要被各搬一次，共 4 次。
// N 那一侧省不掉 —— COARSE 个输出的列号各不相同，每人都得要自己那一格。

再加下去还剩多少？让 COARSE 一路涨：
  COARSE → ∞ 时 (1 + COARSE) ÷ COARSE → 1     // M 那 1 格被摊成 0，只剩 N 的 1 格
  bytesPerOut → 128 阶段 × 1 格 × 4 字节 + 4 = 516 字节
  AI 上限     → 8192 ÷ 516 = 15.876          // ≈ TILE ÷ 2 = 16

  COARSE = 4 已经拿到 12.720 = 上限的 80.1%

K4 · 每线程 8 个输出（同一列上，连续 8 行）

float a0…a7 = 0.0f;                          // ← 8 个累加器 = 8 个输出
for (int k = 0; k < TILE; ++k) {
    const float n = Nds[k][tx];              // ← 贵的那条（32 路）只读一次
    a0 += Mds[i0    ][k] * n;                // ↓ 8 个 FMA 共用它
    a1 += Mds[i0 + 1][k] * n;                //   Mds 这 8 条全是广播（下标无 tx）
    ...
    a7 += Mds[i0 + 7][k] * n;
}

第 1 步  几个输出？   a0…a7 → 8 个输出（同一列上，连续 8 行）

第 2 步  内层每轮：几个 FMA、几条 LDS、各属哪类？
  一轮 k 里：1 条 Nds[k][tx]（含 tx → 32 路）+ 8 条 Mds[..][k]（无 tx → 广播）
            8 个 FMA
  →  dist  = 1 条 ÷ 8 个 FMA = 0.125            ← 这才是 K4 的全部收益
     bcast = 8 条 ÷ 8 个 FMA = 1.000           // 和 K3 一样，没省

第 3 步  ldg
  每阶段每线程搬 8 格 M + 8 格 N = 16 格
  每阶段 32 轮 × 8 个 FMA = 256 个 FMA
  →  ldg = 16 ÷ 256 = 0.062                      // 和 K3 相同

第 4 步  bytesPerOut
  每阶段 16 格 ÷ 8 个输出 = 2 格/输出                    // 和 K3 一模一样
  128 阶段 × 2 × 4 字节 + 4 = 1028 字节

第 5 步  AI
  8192 ÷ 1028 = 7.969 FLOP/byte
  ← 和 K3 的 7.969 逐字相同！字节账一个字都没动
  上限 8033，达成率 = 20770.6 ÷ 8033 = 258.6%
  // 达成率从 K3 的 77.0% 冲到 258.6% —— 不是因为字节变了，是因为速度变了

  // 20770.6 怎么来的（唯一一个量出来的数）：
  //   实测 GFLOP/s = 2·M·N·K ÷ 耗时(秒) ÷ 1e9
  //   137,438,953,472 ÷ 0.006617 ÷ 1e9 = 20770.6        —— 逐项拆解见 §3·二
  //   分子和 K3 是同一个数，分母从 22.225 ms 缩到 6.617 ms：
  //   20770.6 ÷ 6184.0 = 3.36，正好等于耗时之比 —— 这一步"算得"一点没多，纯粹是跑得快。

第 6 步  节拍
  32 × 128 × 2.52e9 × 0.006617 ÷ 4096³ = 0.99   // K3 是 3.34

第 7 步  每条 load 周期
  load/FMA = 0.125 + 1.000 + 0.062 = 1.188
  0.99 ÷ 1.188 ≈ 0.84 周期

第 8 步  模型预测
  2.69 × 0.125 + 0.53 × 1.000 + 1.94 × 0.062 = 0.990
  实测 0.99，误差 -0.0%              // 它是标定点之一，必然吻合

五、K5 / K6：两个方向一起复用，再换成 128 位↑

K5 · 8 行 × 4 列 = 32 个输出

float m[8], n[4], acc[8][4];
for (int k = 0; k < TILE; ++k) {
    for (int r = 0; r < 8; ++r) m[r] = Ads[k][i0 + r];      // 8 条广播
    for (int c = 0; c < 4; ++c) n[c] = Bds[k][tx + c * 32]; // 4 条 32 路
    for (int r = 0; r < 8; ++r)
        for (int c = 0; c < 4; ++c)
            acc[r][c] += m[r] * n[c];                      // 32 个 FMA
}

第 1 步  几个输出？   acc[8][4] → 32 个输出（8 行 × 4 列）

第 2 步  内层每轮：几个 FMA、几条 LDS、各属哪类？
  一轮 k 里：4 条 Bds[k][tx+c*32]（含 tx → 32 路）
            8 条 Ads[k][i0+r]（无 tx → 广播）
            32 个 FMA
  →  dist  = 4 条 ÷ 32 个 FMA = 0.125          // 和 K4 相同
     bcast = 8 条 ÷ 32 个 FMA = 0.250         ← 这一项 K4 是 1.000，现在也摊薄了
  // 每个 m 用 4 次、每个 n 用 8 次 —— 两个方向一起复用，这是 K5 相对 K4 的全部改动

第 3 步  ldg
  每阶段每线程搬 8 格 A + 32 格 B = 40 格
  每阶段 32 轮 × 32 个 FMA = 1024 个 FMA
  →  ldg = 40 ÷ 1024 = 0.039                     // 和 K3c 相同

第 4 步  bytesPerOut
  每阶段 40 格 ÷ 32 个输出 = 1.25 格/输出                // 和 K3c 一样
  128 阶段 × 1.25 × 4 + 4 = 644 字节

第 5 步  AI
  8192 ÷ 644 = 12.720，上限 12823 GFLOP/s
  达成率 = 34506.4 ÷ 12823 = 269.1%

  // 上限 = 算术强度 × 显存峰值带宽（口径见 §2·二）
  //   带宽 = 10501 MHz × 2（DDR 双沿）× 384 bit ÷ 8 = 1008.1 GB/s
  //   12.720 × 1008.1 = 12823 GFLOP/s
  //   （程序乘的是未截断的 8192÷644 和 1008.1；
  //     拿印出来的 12.720 × 1008 手算会差一个零头，那是四舍五入）

  // 34506.4 是实测吞吐，唯一一个量出来的数：
  //   实测 GFLOP/s = 2·M·N·K ÷ 耗时(秒) ÷ 1e9
  //   137,438,953,472 ÷ 0.003983 ÷ 1e9 = 34506.4        —— 逐项拆解见 §3·二
  //   分子和 K4 仍是同一个数，分母从 6.617 ms 缩到 3.983 ms：
  //   34506.4 ÷ 20770.6 = 1.66，正好等于耗时之比。

第 6 步  节拍
  32 × 128 × 2.52e9 × 0.003983 ÷ 4096³ = 0.60   // K4 是 0.99

第 7 步  每条 load 周期
  load/FMA = 0.125 + 0.250 + 0.039 = 0.414
  0.60 ÷ 0.414 ≈ 1.45 周期

第 8 步  模型预测
  2.69 × 0.125 + 0.53 × 0.250 + 1.94 × 0.039 = 0.545
  实测 0.60，误差 -9.2%   ← 第二个外推点，同样在 ±10% 内

K5 同时动了两本账： 字节账 7.969 → 12.720（和 K3c 同一个手法）， 条数账 1.188 → 0.414（把广播也摊薄了）。 占用率从 75.0% 掉到 33.3%， 速度反而快 1.66×：寄存器分块和占用率本来就互相对立。

K6 · 把 12 条 32 位读换成 3 条 128 位读

float4 m01 = ((float4*)&Ads[k][i0])[0];      // 2 条广播覆盖 8 个 m
float4 m23 = ((float4*)&Ads[k][i0])[1];
float4 nv  = ((float4*)&Bds[k][tx * 4])[0];  // 1 条 32 路覆盖 4 个 n

第 1 步  几个输出？   32 个                                  // 和 K5 相同

第 2 步  内层每轮
  1 条 128 位 32 路 + 2 条 128 位广播，仍是 32 个 FMA
  →  dist  = 1 ÷ 32 = 0.031       // K5 是 0.125
     bcast = 2 ÷ 32 = 0.062      // K5 是 0.250

第 3 步  ldg   每阶段 10 条 float4 ÷ 1024 个 FMA = 0.010

第 4 步  bytesPerOut
  10 条 × 4 个 float = 40 个 float —— 和 K5 的 40 格一模一样
  →  644 字节，和 K5 同一个数
  // 搬的字节一个没少，只是打包方式变了：40 个 float 从 40 条指令变成 10 条

第 5 步  AI   8192 ÷ 644 = 12.720   ← 和 K5 逐字相同
             达成率 = 40928.8 ÷ 12823 = 319.2%   // 七个里最高

  // 上限 12823 和 K5 也是同一个数 —— 它只由算术强度和带宽决定，和跑多快无关：
  //   12.720 × 1008.1 GB/s = 12823 GFLOP/s        // 带宽口径见 §2·二

  // 40928.8 是实测吞吐，唯一一个量出来的数：
  //   137,438,953,472 ÷ 0.003358 ÷ 1e9 = 40928.8        —— 逐项拆解见 §3·二
  //   分子和 K5 仍是同一个数，分母从 3.983 ms 缩到 3.358 ms：
  //   40928.8 ÷ 34506.4 = 1.19，正好等于耗时之比。

  // 所以这一节的三个"逐字相同"（AI、上限、分母）之外，只有耗时在动 ——
  // 达成率从 269.1% 涨到 319.2%，涨的全是分子。

第 6 步  节拍   32 × 128 × 2.52e9 × 0.003358 ÷ 4096³ = 0.50
  // 4090 的物理地板是 0.25 —— K6 是它的 2 倍

第 7 步  每条 load 周期
  load/FMA = 0.031 + 0.062 + 0.010 = 0.104
  // 每线程每阶段：内层 32 步 × 3 条 128 位读 + 全局 10 条 float4，共 32×32 = 1024 个 FMA
  //   dist  = 32 步 × 1 条（Bds[k][tx*4]）      ÷ 1024 = 32 ÷ 1024 = 0.031
  //   bcast = 32 步 × 2 条（Ads[k][i0]、[i0+4]） ÷ 1024 = 64 ÷ 1024 = 0.062
  //   ldg   = A 的 2 条 + B 的 8 条            ÷ 1024 = 10 ÷ 1024 = 0.010
  //   合计 106 ÷ 1024 = 0.1035 → 0.104（三个数各自四舍五入再相加会得 0.103，差的是零头）
  0.50 ÷ 0.104 ≈ 4.85 周期
  ← 突然跳到 4.85，因为"一条"的含义变了：现在一条搬 4 个 float

第 8 步  模型预测
  2.69 × 0.031 + 0.53 × 0.062 + 1.94 × 0.010 = 0.136
  实测 0.50，误差 -72.8%
  // 三个单价都是用 32 位 load 标定的，套到 128 位 load 上必然不准

六、七个 kernel 的全部指标，和单价标定↑

三个单价由 K2 / K3 / K4 三个实测点解三元一次方程组得到。

K2:  ldg?  × 2.000                                    = 3.88   →  LDG = 1.940
K3:  dist? × 1.000 + bcast? × 1.000 + 1.940 × 0.062   = 3.34
K4:  dist? × 0.125 + bcast? × 1.000 + 1.940 × 0.062   = 0.99
                     ↓ 两式相减，bcast 消掉
     dist? × 0.875 = 3.34 − 0.99 − 0
                                        →  32 路 = 2.686
                                        →  广播  = 0.534

// 09-21 那一批解出的是 2.61 / 0.53 / 1.89 —— 换了机器、换了驱动，
// 三个单价的相对关系稳得住，绝对值差 3% 左右。

一、按账本拆开看每一步↑

二、要改写的那句结论↑

一、尺寸：为什么主表必须用 4096↑

1024 测的不是显存带宽：它的 footprint 只有 12.6 MB，整个装得进 4090 的 75.5 MB L2。 到 4096（201 MB）才明显超出 L2。 主表取 4096 是折中：够大，能脱离 L2；又不会大到 K1 单次就要约 1.7 秒（8192 上的实测）。

二、形状：换成真实推理负载还剩多少↑

主表全是方阵，但生产里的 GEMM 几乎从来不是方阵。 prefill 是瘦长矩阵；decode 的 M 等于 batch size，常常只有 1。 N/K 取自 Qwen2.5-7B 的 FFN 维度（18944 × 3584）。

方阵 4096³     K1 → K6   63.4×
prefill        K1 → K6   49.2×
decode b32     K1 → K6   23.8×
decode b1      K1 → K6   3.4×   ← 只剩零头

cuBLAS   decode b1   M=1    0.302 ms   449.6 GFLOP/s
         decode b32  M=32   0.313 ms   13882.7 GFLOP/s

  活多了 32×（M 从 1 变成 32），耗时只多 1.04×
  →  实测吞吐比 = 13882.7 ÷ 449.6 = 30.9×        // = 32 ÷ 1.04

一、K1 的 211.911 ms 里，四分之三和"不合并"无关↑

            N=4096（2 的幂）    N=4097（打破对齐）
K1           211.911 ms          63.675 ms     ← 快了 3.33 倍
K1/K2          8.209×             2.392×

8.209 ÷ 2.392 = 3.43   ← 2 的幂跨度的额外惩罚
回乘检验：2.392 × 3.432 = 8.209（3.43 是 3.432 四舍五入后的值）

N=4096 时，A 的行跨度恰好是 16384 字节 = 2¹⁴；N=4097 时变成 16388。 两次运行只改了这一个变量。 所以"合并访存是性价比最高的一步"要加注： 那 8.21 倍里，只有 2.392 倍来自合并本身，另外 3.43 倍来自顺带治好的跨度冲突。 对照组：同样从 4096 改到 4097，K2–K5 只慢 2.8%–7.1%，cuBLAS 慢 15.7%。K6 用 float4，要求 N 是 4 的倍数，所以测的是 4100，慢 7.5%。

// 最顺手的解释是"warp 内 32 个地址撞同一个 L2 set"。把地址逐位算一遍：
              行号(addr÷128)        行内偏移      sector 效率
N=4096        base + 128t           恒为 0         12.5%
N=4097        base + 128t           0,4,8,…        12.5%
              ↑ 完全相同                      ↑ 完全相同

// 16388 = 2¹⁴ + 2²，多出来的只影响第 2–6 位 —— 那是 128 字节行的【内部】，
// 不参与任何 set / bank 的选择。所以那个解释解释不了 3.33 倍。

二、编译器已经替 K5 做了一半向量化↑

kernel     FFMA  LDS.128   LDS.32  LDG.128   LDG.32    TOTAL
k1           30        -        -        -       59      216
k2           30        -        -        -       59      216   ← 和 k1 逐项相同
k5         1056       64      128        -       72     1856   ← 预测说这里该是 0
k6         1056       96        -       18        -     1520

一、六个指标里，只有两个是真的拿不到↑

二、"合并效率"图的三条替代证据↑

最后那个指标要证明"K1 的 8.21 倍劣势来自访存不合并"。没有它，这件事照样能证明：

实测 nsys 在这种容器里能跑（CUDA 时间线可用）。 但它的 --gpu-metrics-devices 走同一条权限路径，同样报 ERR_NVGPUCTRPERM。