matmul 复杂度多少？tiling 怎么做？
CPU 和 GPU 算 matmul 有什么区别？roofline 怎么用？
先弄清楚：每个概念是什么
- FLOPs：为什么是 2MNK
- C 有 \(M\times N\) 个元素，每个是长度 \(K\) 的 dot product，\(K\) 次乘 + \(K\) 次加。例：\(M=N=K=4096\) → \(2\cdot 4096^3\approx 1.37\times10^{11}\) FLOPs。
- Arithmetic intensity
- 每从 HBM 搬 1 byte 做几次 FLOP（基础见 GPU 与 NumPy）。matmul 数据是 \(n^2\) 级、计算是 \(n^3\) 级，本来很高。但 naive kernel 每次乘加都去 HBM 现读两个数，FP32 只有 0.25 FLOP/B。
- Tiling
- 把 C 切成 \(T\times T\) 小块，一个 thread block 负责一块。沿 K 每次把 A、B 各一个 tile 搬进 shared memory，里面每个数被用 \(T\) 次。例：FP32、\(T=32\)：intensity 从 0.25 涨到 8 FLOP/B。
- Roofline 与 ridge point
- 性能上限 = min(峰值算力, intensity × bandwidth)。先是一条斜线（memory-bound），再是一条平线（compute-bound），拐点叫 ridge point，intensity 在它左边时加算力没用。
- Tensor Core
- 专做小矩阵块乘加的单元，FP16 / BF16 输入、FP32 累加。例：A100 的 BF16 Tensor Core 312 TFLOPS，普通 FP32 只有 19.5 TFLOPS。elementwise 运算用不上它。
看图：roofline 上谁卡在哪，tiling 怎么搬数据
30 秒回答
Multiplying M×K by K×N takes MNK multiply-adds, so 2MNK FLOPs, or O(n³). A naive kernel rereads operands from memory for every multiply-add, so it's memory-bound. With tiling, each thread block owns a tile of C, walks along K, stages tiles of A and B in shared memory and reuses each value T times, so intensity grows with T. Then come register tiling, vectorized loads, double buffering and Tensor Cores. CPUs do the same with cache blocking and SIMD; GPUs use thousands of threads and hide latency by switching warps. Roofline: attainable FLOPs = min(peak, intensity × bandwidth).
推导 / 关键公式
tiling 为什么有用
一个 \(T\times T\) 的 C tile 沿 K 每走一步：搬 A、B 各 \(T^2\) 个数，做 \(2T^3\) FLOPs。FP32（4 bytes）：
\[ I_{\text{tile}} = \frac{2T^3}{4\cdot 2T^2} = \frac{T}{4} \]\(T\) 受 shared memory、register 大小限制。tile 太大则块数太少，填不满所有 SM。
roofline
\[ P = \min\big(P_{\text{peak}},\ I\times BW\big) \] \[ I^{*} = P_{\text{peak}} / BW \]整个 GEMM 理想情况下 A、B、C 各过一次 HBM（BF16，2 bytes）：
\[ I = \frac{2MNK}{2\,(MK + KN + MN)} \]手算例子：A100 上谁是 memory-bound
A100 80GB SXM：BF16 dense 312 TFLOPS，HBM 2,039 GB/s → ridge point \(I^{*}\approx 153\) FLOP/B。
| 操作（BF16） | intensity | 上限 | 瓶颈 |
|---|---|---|---|
| \(4096^3\) GEMM | ≈ 1365 | 312 TFLOPS | compute |
| GEMV（decode） | ≈ 1 | ≈ 2.0 TFLOPS（0.65%） | memory |
| elementwise add | 1/6 | ≈ 0.34 TFLOPS | memory → fusion |
CPU vs GPU
| CPU | GPU | |
|---|---|---|
| 并行 | 十来个强 core + SIMD（AVX） | 上千 thread，32 个一组成 warp |
| 分块放哪 | A 块放 L2，B 的窄条放 L1，C 的小块放 register | block tile 放 shared memory，thread tile 放 register |
| 藏 latency | 大 cache + prefetch | warp 等 memory 时切另一个 |
手撕代码（CUDA，shared memory tiling）
#define T 32
__global__ void matmul(const float* A, const float* B, float* C,
int M, int N, int K) {
__shared__ float As[T][T], Bs[T][T];
int tx = threadIdx.x, ty = threadIdx.y;
int r = blockIdx.y * T + ty, c = blockIdx.x * T + tx;
float acc = 0.f;
for (int k0 = 0; k0 < K; k0 += T) { // 沿 K 一块一块走
As[ty][tx] = (r < M && k0 + tx < K) ? A[r * K + k0 + tx] : 0.f;
Bs[ty][tx] = (c < N && k0 + ty < K) ? B[(k0 + ty) * N + c] : 0.f;
__syncthreads(); // tile 搬完再算
for (int k = 0; k < T; ++k) // 每个数在 SRAM 里用 T 次
acc += As[ty][k] * Bs[k][tx];
__syncthreads(); // 算完再搬下一块
}
if (r < M && c < N) C[r * N + c] = acc;
}读 A、B 时相邻 thread 读相邻地址（见 coalescing）。
必答点
- \(2MNK\) FLOPs，方阵 \(O(n^3)\)
- naive 是 memory-bound；tiling 让每个数复用 \(T\) 次，intensity \(\propto T\)
- 优化顺序：shared memory tiling → register tiling → vectorized load → double buffering → Tensor Core
- CPU：cache blocking + SIMD + 多线程；GPU：海量 thread、切 warp 藏 latency、Tensor Core
- roofline：A100 ridge ≈ 153 FLOP/B；大 GEMM compute-bound，GEMV、LayerNorm、elementwise memory-bound，靠 fusion 少过 HBM（如 FlashAttention）
追问
shared memory tiling 之后还能怎么优化？
register tiling 让每个 thread 算一小块 C（如 8×8）。float4 一次读 16 bytes，double buffering 让下一块的搬运和这一块的计算 overlap。Tensor Core（WMMA / MMA）的 FP16 维度取 8 的倍数（A100 上最好 64）。
有比 \(O(n^3)\) 更快的算法吗？
Strassen 是 \(O(n^{\log_2 7})\approx O(n^{2.81})\)，代价是很多额外的矩阵加法和临时矩阵。cuBLAS 的 GEMM 仍按 \(2MNK\) 计 FLOPs。
怎么判断一个 kernel 卡在哪？
算 achieved FLOP/s 和 achieved bandwidth 各占峰值多少，接近上限的那个就是瓶颈。Nsight Compute 能画出 kernel 在 roofline 上的位置。
来源NVIDIA Matrix Multiplication Background §2、§2.1–2.3 · NVIDIA GPU Performance Background §2、§4 · A100 datasheet（p.1）· Williams et al., CACM 2009（p.67）· Goto & van de Geijn 2008 §4、§6.1 · Nsight Compute Profiling Guide §2.9