Content
GPU 访存与计算优化
> 本文为 LCPU Infra Seminars 第二部分课程笔记。回放：[2.1-Memory Abstraction & Hierarchy（主讲：林若瑜）](https://www.bilibili.com/video/BV1gqGA6ZEfg)，[2.2-FP32 GEMM Quick Walkthrough（主讲：周宇轩）](https://www.bilibili.com/video/BV1L8GA6YEAH)
>
> 注：AI 参与部分笔记整理工作
下文沿用[CUDA 的层级](https://aerisun.top/posts/51#cuda-的层级)：一个 warp 包含 32 个线程，lane 是线程在 warp 内的编号；block 也称 CTA，是线程协作和共享内存的边界；SM 是执行这些 warp 的硬件单元。
Roofline 与算术强度
计算吞吐与显存带宽的上界
Roofline 模型用计算吞吐和内存带宽估算性能上界：
$$
P \leq \min\left(P_{\text{peak}},\ BW_{\text{mem}}\times I\right)
$$
其中，$P$ 是 kernel 的浮点运算吞吐，$P_{\text{peak}}$ 是对应数据类型和计算管线的峰值吞吐，两者单位均为 FLOP/s；$BW_{\text{mem}}$ 是显存带宽，单位为 Byte/s；算术强度 $I$ 表示每搬运一个字节完成多少次浮点运算，单位为 FLOPs/Byte。
*斜线表示带宽给出的上界，水平线表示计算吞吐给出的上界。*
$$
I^*=\frac{P_{\text{peak}}}{BW_{\text{mem}}}
$$
同一个 kernel 在不同 GPU 上可能落在不同区域。
> 若实测点远低于上界，仅凭 Roofline 还不能定位原因。请求并发不足、访问分散、指令依赖和同步，都可能让硬件资源没有充分利用。
>
GEMM 的理想复用与实际访存
对 $C=A\times B$，设 $A\in\mathbb{R}^{M\times K}$、$B\in\mathbb{R}^{K\times N}$，则 $C\in\mathbb{R}^{M\times N}$。若每个输入只从 DRAM 读取一次，输出只写回一次：
$$
I_{\text{ideal}}\approx
\frac{2MNK}{4(MK+KN+MN)}
$$
当 $M=N=K=n$ 时，代入得到 $2n^3/(12n^2)=n/6$。这个理想值随矩阵规模增加，前提是实现能够充分复用输入。
“一个线程算一个输出”的朴素实现会让相邻线程重复读取同一行 $A$ 和同一列 $B$。若暂不计缓存复用，单个线程执行 $K$ 次 FMA，读取 $2K$ 个 float，再写回一个结果：
$$
I_{\text{naive}}\approx\frac{2K}{8K+4}\approx 0.25\ \text{FLOPs/Byte}
$$
:::details{summary="补充：怎样快速估算一个 kernel 的算术强度"}
以逐元素 c[i] = a[i] * b[i] + c[i] 为例：一次乘加按 2 FLOPs 计，读取 a、b、c 并写回 c，若都是 FP32，则最朴素的估算是：
$$
I\approx\frac{2}{4+4+4+4}=0.125\ \text{FLOPs/Byte}
$$
缓存命中、写分配和编译器变换会改变实际流量，所以这个数用于判断数量级，最终仍要由性能计数器验证。
:::
在途请求与延迟隐藏
Little’s Law：带宽、延迟与在途字节
一次显存访问的延迟很长，不代表显存带宽一定低。尚未完成的访存请求称为在途请求。对稳定运行时的平均量，Little’s Law 给出：
$$
\text{Bytes in flight}=\text{Bandwidth}\times\text{Latency}
$$
*每两秒通过一人、全程四十秒，对应平均二十人同时在扶梯上。*
延迟不变时，要维持更高带宽，就需要更多请求同时等待完成。这些请求既可以来自多个 warp，也可以来自同一线程内彼此独立的访存指令。
循环展开与独立访存
如果先读完、算完并写回一个元素，再开始处理下一个元素，两组访存就难以重叠。循环展开把多个元素的操作放进同一段代码，为编译器调度独立指令提供机会，这称为指令级并行（Instruction-Level Parallelism，ILP）。
```cpp
int i0 = blockIdx.x * (2 * blockDim.x) + threadIdx.x;
int i1 = i0 + blockDim.x;
float a0 = a[i0];
float b0 = b[i0];
float c0 = c[i0];
float a1 = a[i1];
float b1 = b[i1];
float c1 = c[i1];
c[i0] = a0 * b0 + c0;
c[i1] = a1 * b1 + c1;
```
*左侧时间线提前发出两组 load，右侧标出 24 bytes 的读取和 8 bytes 的写回。*
两组输入的等待时间可以重叠，所以在简化调度示例中，工作量翻倍而总周期数仍接近原先水平。
展开会延长中间值的存活时间，增加寄存器需求（因为访存时目的寄存器会占位锁定）。<下划线起>资源占用过高可能减少可驻留的 warp 数，或发生 spill，把部分值保存到 local memory！<下划线末>
Occupancy 与可发射的 warp
隐藏延迟有两种主要来源：
- 线程级并发：一个 warp 等待时，调度器切换到另一个 ready warp。
- 指令级并发：同一线程提前发出多条互不依赖的请求。
Occupancy 衡量一个 SM 上实际驻留的 warp 数相对硬件上限的比例。它受 block 大小、每线程寄存器数和每 block shared memory 用量共同限制。
Resident warp 指已经驻留在 SM 上的 warp，eligible warp 指当前满足指令发射条件的 warp。高 occupancy 只说明驻留得多；如果这些 warp 都在等待数据或同一队列，就仍然缺少可发射的工作。反过来，单个 warp 内有足够多独立指令时，较低 occupancy 也可能达到较高吞吐。
Coalescing 与数据布局
请求并发解决等待能否重叠的问题，访存合并（coalescing）则关系到每次搬运的数据有多少被实际使用。
连续访问、偏移访问与 sector 数量
NVIDIA GPU 中，warp 搬运全局内存数据的最小物理单位是 32-byte Sector。
对于 32 个 lane，每个 lane 连续读取一个 4-byte float，总有效数据恰好是 128 bytes。若起始地址对齐，它们覆盖 4 个 sector。下图的绿色区域就是这 4 个连续的 32-byte sector。
只偏移一个元素，就可能覆盖 5 个 sector。下图的访问箭头整体向右偏移，绿色区域随之多覆盖一段。逻辑上仍然是连续数组，物理传输量却增加了 25%。
AoS（Array of Structures，结构体数组）把同一对象的字段放在一起，访问相同字段时会产生跨步，有效数据只有 $32\times4=128$ bytes，却覆盖了 32 个 sector。
float4 向量化的条件
float4 把四个 FP32 元素组织成 16 bytes 的向量。编译器生成相应的向量访存指令时，它能减少 load/store 指令数量（因为指令也会花费时间，太多的时候形成拥塞）。使用时要分别检查对齐、边界和地址布局：
- 起始地址满足 16-byte 对齐；
- 长度尾部仍有正确的边界处理；
- 相邻 lane 的向量块最好连续。
:::details{summary="补充：一个带尾部处理的 float4 访存骨架"}
下面的加法示例要求 x、y、out 均满足 16-byte 对齐，三个数组互不重叠，n 非负且下标计算不溢出；启动线程数至少覆盖 $\lceil n/4\rceil$ 组元素。
```cpp
__global__ void add_vec4(
const float* __restrict__ x,
const float* __restrict__ y,
float* __restrict__ out,
int n) {
int vec = blockIdx.x * blockDim.x + threadIdx.x;
int i = vec * 4;
if (i + 3 < n) {
float4 xv = reinterpret_cast<const float4*>(x)[vec];
float4 yv = reinterpret_cast<const float4*>(y)[vec];
float4 ov{xv.x + yv.x, xv.y + yv.y,
xv.z + yv.z, xv.w + yv.w};
reinterpret_cast<float4*>(out)[vec] = ov;
} else {
for (int j = i; j < n; ++j) {
out[j] = x[j] + y[j];
}
}
}
```
完整的四元素组走向量路径，最后不足四个元素时逐个处理。设备分配的基地址通常满足对齐要求，但传入切片指针前仍需检查偏移。
:::
AoS 与 SoA 的字段布局
若 kernel 总是让相邻线程读取不同对象的同一字段，可以把前面的 AoS 改成 SoA（Structure of Arrays，数组结构体），为各字段分别存放连续数组：
```cpp
struct SoA {
float* u;
float* v;
float* w;
};
```
此时 u[tid] 在 warp 内天然连续，访存好。
Reduce 的分层归约
Reduce 把多个输入合成一个结果，例如数组求和 $s=\sum_{i=0}^{N-1}x_i$。
问题和解决方法
全局原子操作互锁
若每个线程都执行 atomicAdd(out, x[i])，就会有大量线程争用同一个地址。
可以先在 block 内算出部分和，再由一个线程提交结果：全局原子操作从“每个输入一次”减少到“每个 block 一次”。
Block 内可以用 shared memory 做树形归约。例如 8 个部分和，每轮两两相加，数量依次变为 $8\rightarrow4\rightarrow2\rightarrow1$。
Shared memory 的 bank conflict
在 32 个 bank、连续 32-bit 字映射到连续 bank 的 shared memory 模型下，对于同一 warp 中的不同线程：
- 访问不同 bank 可以并行
- **访问同一 bank 的*不同地址*会产生 bank conflict**
- **访问同一 bank 的*同一地址*无冲突，可以广播。**
*从左到右：不同 bank、两路冲突、同一地址广播。*
因此，树形归约除了每轮减少多少个部分和，还要看线程怎样访问这些部分和。
> bank conflict 与[coalescing](#连续访问偏移访问与-sector-数量)分别检查两种映射：前者是 Shared memory，看 bank；后者是 global memory 看 sector。
>
> 二维 float tile[32][32] 按列访问也会遇到 bank conflict 冲突，增加一列 padding 成为 [32][33] 可以改变 bank 映射，但会多占 shared memory。
先在线程内累加，再交换部分和
让每个线程先遍历多个输入，把和留在寄存器中，最后才参与一次 block 归约。
Grid-stride 循环：线程从各自的起始下标出发，每轮前进整个 grid 的线程数。同一轮中相邻线程仍读相邻元素，而每个线程可以处理多轮数据。这样能用较少的 block 覆盖整个数组，减少 block 归约和全局结果提交的次数（block 也不能少到让 GPU 缺少并行工作）。
Warp shuffle 的数据交换
一个 warp 内可以用 shuffle 直接交换寄存器值，省去 shared memory 中转。
block 归约分两级：各 warp 先算出一个和，把这些和写入 shared memory；同步后，第一个 warp 再归约它们。
*上方经 shared memory 中转，下方在 warp 内直接交换寄存器值。*
下面的循环让 lane 0 得到整个 warp 的和。前提是 32 个 lane 都执行每次 shuffle；输入不足时以零参与，而不提前退出。
```cpp
constexpr unsigned mask = 0xffffffffu;
for (int offset = 16; offset > 0; offset >>= 1) {
value += __shfl_down_sync(mask, value, offset);
}
```
*图中用 8 个 lane 演示：右侧经过 4、2、1 三级交换，lane 0 得到总和；左侧 Butterfly 则把总和传播到所有 lane。*
Mask 只声明参与线程，不能自动处理缺失的源 lane。这里限定 block 线程数为 32 的倍数，无效输入补零。
下面保留 block 级骨架，采用每轮处理两个元素的 grid-stride 变体。使用一维 block，线程数为 32 的倍数且不超过 1024。
:::details{summary="补充：block 级 Reduce 的核心骨架"}
```cpp
__inline__ __device__ float warp_reduce_sum(float value) {
constexpr unsigned mask = 0xffffffffu;
for (int offset = 16; offset > 0; offset >>= 1) {
value += __shfl_down_sync(mask, value, offset);
}
return value;
}
__global__ void reduce_sum(const float* x, float* out, int n) {
__shared__ float warp_sums[32];
int tid = threadIdx.x;
int i = blockIdx.x * blockDim.x * 2 + tid;
int stride = blockDim.x * 2 * gridDim.x;
float value = 0.0f;
for (; i < n; i += stride) {
value += x[i];
if (i + blockDim.x < n) {
value += x[i + blockDim.x];
}
}
value = warp_reduce_sum(value);
if ((tid & 31) == 0) {
warp_sums[tid >> 5] = value;
}
__syncthreads();
int warp_count = (blockDim.x + 31) / 32;
value = tid < warp_count ? warp_sums[tid] : 0.0f;
if (tid < 32) {
value = warp_reduce_sum(value);
if (tid == 0) {
atomicAdd(out, value);
}
}
}
```
代码先在 warp 内用 shuffle 得到部分和，再用 shared memory 把各 warp 的结果交给第一个 warp。grid 大小需要按数据规模选择。
:::
> 浮点加法不满足结合律，归约顺序和 block 提交顺序可能改变末位结果；验证时应使用合理的误差容限。
>
课件中的 Reduce 性能对比
课件按“全局原子加 → shared-memory 归约与访问调整 → grid-stride 局部累加 → warp shuffle → float4”逐步优化：
GEMM 的分块与数据复用
GEMM 计算 $C=AB$，其中 $A$ 为 $M\times K$，$B$ 为 $K\times N$，$C$ 为 $M\times N$。
Naive GEMM 的重复读取
若一个线程只计算一个输出：
$$
C_{ij}=\sum_{k=0}^{K-1}A_{ik}B_{kj}
$$
那么计算同一行输出的线程会重复读取 $A$ 的同一行，计算同一列输出的线程会重复读取 $B$ 的同一列。
*四个线程计算相邻的四个输出，共同依赖 $A$ 的两行和 $B$ 的两列。*
忽略缓存时，这种实现的算术强度约为[0.25 FLOPs/Byte](#gemm-的理想复用与实际访存)。
Shared Tiling 与 block 内复用
让一个 block 负责 $C$ 中大小为 $B_M\times B_N$ 的输出 tile。沿求和维度 $K$ 每次处理 $B_K$ 个元素：合作加载一个 $B_M\times B_K$ 的 $A$ 子块和一个 $B_K\times B_N$ 的 $B$ 子块，在 shared memory 中复用，累加到同一个输出 tile，直到遍历完 $K$。
*A、B 的灰块沿 K 维前进，C 的灰块位置不变，持续累加各轮结果。*
每个 $A$ 元素可以服务于 $B_N$ 列输出，每个 $B$ 元素可以服务于 $B_M$ 行输出。若每轮输入只加载一次，FP32 元素占 4 bytes，一次 FMA 按 2 FLOPs 计，忽略输出写回，则：
$$
I_{\text{CTA}}\approx
\frac{2B_MB_NB_K}{4(B_MB_K+B_KB_N)}
=\frac{B_MB_N}{2(B_M+B_N)}
$$
增大 $B_M,B_N$ 能提高复用，但也增加 shared memory 和输出累加器的需求，可能减少可驻留的 block 数。
核心流程如下：
```cpp
// As、Bs 是共享数组，ty、tx 是当前线程在输出 tile 中的坐标；acc 在外层循环前置零。
for (int k0 = 0; k0 < K; k0 += BK) {
load_a_tile_to_shared(As, A, k0);
load_b_tile_to_shared(Bs, B, k0);
__syncthreads();
for (int k = 0; k < BK; ++k) {
acc += As[ty][k] * Bs[k][tx];
}
__syncthreads();
}
```
第一次同步保证大家加载完，第二次保证大家读完，才能覆盖共享数组。越界输入应补零，所有线程仍需到达同步点
作加载还要沿行优先布局的连续列分配相邻线程，满足[coalescing](#连续访问偏移访问与-sector-数量)；采用 float4 时另外检查对齐和尾部。
Register Tiling 与线程内复用
Shared Tiling 减少了 global-memory 重复读取，但上面的内层循环每做一次 FMA，仍要从 shared memory 取一个 $A$ 和一个 $B$。如果这些加载指令限制了发射速度，就需要让每次 shared load 支持更多计算。
*课件将 shared-memory 指令归到 MIO 路径；统计中 Load/Store 指令最多。*
Register Tiling 让一个线程计算 $T_M\times T_N$ 个输出，并用寄存器保存这些累加器。每个 $k$ 只读 $T_M$ 个 $A$ 值和 $T_N$ 个 $B$ 值，再两两组合，完成 $T_MT_N$ 次 FMA。
以 $T_M=T_N=2$ 为例：四个输出原本各自需要两次读取，共八次；共享输入后，只需读两个 $A$ 和两个 $B$，就能更新四个输出。
*左侧小框是线程负责的 micro-tile；课件用 $16\times16$ 个线程、每线程 $2\times2$ 个输出，覆盖$32\times32$ 的 block tile。*
下面只展示一个；假定共享输入已就绪且访问不越界。
```cpp
// k 的更新，r、c 是线程输出小块在共享 tile 中的行列起点。四个累加器在遍历 K 前置零，并保留到所有分块计算结束
float a0 = As[r][k], a1 = As[r + 1][k];
float b0 = Bs[k][c], b1 = Bs[k][c + 1];
acc00 += a0 * b0;
acc01 += a0 * b1;
acc10 += a1 * b0;
acc11 += a1 * b1;
```
shared tile 供 block 内多个线程复用，register micro-tile 让同一线程的多个输出复用输入。 代价是更多寄存器；$T_M,T_N$ 过大可能降低 occupancy 或产生 spill。需把它们与 $B_M,B_N,B_K$ 和 block 线程数一起选择。
数据搬运与计算流水线
两级分块减少了搬运量，但若仍按“加载一块 → 计算一块 → 再加载”串行执行，计算下一块之前还要等待数据。流水线使用多个缓冲阶段（stage）：计算当前 tile 时，提前搬运下一块；使用前等待加载完成，复用缓冲区前确认旧数据已读完。
Kernel 优化检查顺序
| 观察 | 更可能的含义 | 回到代码先看什么 |
| ------------------------------------ | ----------------------------------------------------------- | ------------------------------------------------- |
| DRAM 吞吐低，long scoreboard 高 | 在途请求不足或依赖链太长 | grid 是否够大、是否能展开独立 load |
| DRAM 吞吐高，FMA 管线低 | 算术强度低，接近 memory-bound | 是否存在可利用的数据复用 |
| sectors/request 偏高 | warp 地址不连续或未对齐 | lane 到地址的映射、AoS/SoA、leading dimension |
| shared bank conflict 高 | 多 lane 命中同一 bank 的不同地址 | stride、转置访问、padding |
| resident warps 多，eligible warps 少 | warp 都在等待同一依赖或队列 | shared/MIO 压力、同步、数据依赖 |
| local memory 流量出现 | 可能发生寄存器 spill，也可能是局部数组等被放入 local memory | 编译报告、局部变量寻址、展开程度、micro-tile 大小 |
| FMA 利用率高但仍未到目标 | 计算管线已接近主瓶颈 | 指令调度、流水线、是否已接近合理上界 |
分析一个新的 CUDA kernel，可以按这个顺序走：
1. 写出数据流：每个输出读取什么，哪些输入会被多个线程复用。
2. 估算算术强度：先算算法上限，再算当前实现真正穿过 DRAM 的字节。
3. 定位 Roofline 区域：先判断 memory-bound 还是 compute-bound。
4. 检查请求并发：grid、resident warp、eligible warp 和单线程独立 load 是否足够。
5. 检查事务效率：lane 地址是否连续、对齐，sector 利用率是否合理。
6. 逐级放置数据：DRAM 负责容量，shared memory 负责 CTA 内复用，寄存器负责线程内复用。
7. 检查片上瓶颈：bank conflict、MIO 队列、同步和寄存器 spill。
8. 一次只改一个变量：保留每个版本的时间、带宽和关键计数器。
9. 评估流水线收益：确认加载等待仍影响关键路径，再比较异步搬运与额外缓冲的收益和开销。