Content
初识 CUDA
> 本文为 LCPU Infra Seminars 第一部分课程笔记。回放：[1.1-Cuda Programming Model（主讲：郑熠）](https://www.bilibili.com/video/BV1g83w6DETm/)，[1.2-Triton/TileLang Tile Level Programming（主讲：王艺霏）](https://www.bilibili.com/video/BV1g83w6DETm/)
>
> 注：AI 参与部分笔记整理工作
GPU 与 CPU
CPU 追求低延迟。操作系统、业务逻辑和分支密集的代码不太规整，所以 CPU 会花很多晶体管做复杂控制、分支预测和多级缓存，尽量让当前任务早点结束。
GPU 追求吞吐量。它拿出更多面积放计算单元，适合让大量数据执行相似的操作。
*CPU 偏向复杂控制，GPU 偏向堆叠并行计算资源。*
| | CPU | GPU |
| --- | --- | --- |
| 更关心 | 单个任务多久结束 | 单位时间完成多少任务 |
| 擅长 | 不规则控制、串行依赖 | 规则的数据并行 |
| 等待内存时 | 依靠缓存和乱序执行 | 换一个就绪的 warp 继续跑 |
GPU 需要足够多的 thread 才有工作可调度；具体怎样用 warp 填住访存空档，放到[SM 调度与延迟隐藏](#sm-调度与-occupancy)再看。
CUDA 的层级
| 视角 | 层级 |
| --- | --- |
| 程序组织 | Grid → Block / CTA → Thread |
| 硬件执行 | SM → Warp → Thread |
grid、block、thread 用来描述工作；SM 是真正执行工作的硬件。warp 很少直接出现在启动配置里，却是理解性能最重要的一层。
Thread
thread 是一份 kernel 程序实例。所有 thread 跑同一段代码，再用自己的索引找到不同的数据。
向量加法里，第 $i$ 个 thread 只负责：
$$
C_i = A_i + B_i.
$$
Warp
NVIDIA GPU 把 32 个 thread 编成一个 warp。
一个 warp 通常发射同一条指令，只是每个 thread 操作的数据不同。
因此，下面这些问题都要按 32 个 thread 一组来看：
- 分支是否走在一起，也就是[warp divergence](#warp-divergence)；
- 地址是否连续，也就是[coalescing](#coalescing)；
- 寄存器数据如何交换，例如[warp shuffle](#warp-shuffle)；
- SM 还有多少 warp 可以切换。
Block / CTA
block 是 thread 的协作边界。
同一 block 内可以共享 shared memory，也可以用 __syncthreads() 同步。
不同 block 通常不应该依赖执行顺序，否则程序很难随 GPU 规模扩展。
Grid 和 SM
一次 kernel 启动会产生一个 grid，其中可以有远多于 SM 数量的 block。
硬件把 block 分批派给 SM；某个 block 一旦落到一个 SM 上，就会在那里跑到结束。
*block 多于 SM 没有问题，多出来的会排队。*
CUDA Kernel
CUDA 源文件通常以 .cu 结尾，由 nvcc 拆分 host 和 device 两部分：
| 源码部分 | 编译路径 | 运行位置 |
| --- | --- | --- |
| Host code | C++ 编译器 → CPU 目标代码 | CPU |
| Device code | 面向 CUDA 的中间指令表示 PTX → 面向具体 GPU 的 CUBIN | GPU |
| 修饰符 | 执行位置 | 用途 |
| --- | --- | --- |
| __host__ | CPU | 普通 CPU 函数，通常可以省略 |
| __global__ | GPU | 从 host 启动的 kernel |
| __device__ | GPU | 由 GPU 代码调用的辅助函数 |
host 负责分配、复制和启动；kernel 只在 device 上执行。
:::details{summary="基础语法："}
一个最小的向量加法 kernel 如下：
```cpp
__global__ void vec_add(
const float* a,
const float* b,
float* c,
int n
) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < n) {
c[i] = a[i] + b[i];
}
}
```
threadIdx.x 是 thread 在 block 内的位置，blockIdx.x 是 block 在 grid 内的位置，blockDim.x 是每个 block 的 thread 数。全局下标就是：
$$
i = \text{blockIdx.x} \times \text{blockDim.x} + \text{threadIdx.x}.
$$
memory 管理：
```cpp
int threads = 256;
int blocks = (n + threads - 1) / threads;
// 整数向上取整。它会多开几个 thread，所以 kernel 里的 i < n 不能删。
size_t bytes = static_cast<size_t>(n) * sizeof(float);
float *d_a = nullptr, *d_b = nullptr, *d_c = nullptr;
cudaMalloc(&d_a, bytes);
cudaMalloc(&d_b, bytes);
cudaMalloc(&d_c, bytes);
cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, h_b, bytes, cudaMemcpyHostToDevice);
vec_add<<<blocks, threads>>>(d_a, d_b, d_c, n);
cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c);
```
- cudaMalloc 分配 device memory，返回的是 device pointer；
- cudaMemcpyHostToDevice 和 cudaMemcpyDeviceToHost 指定复制方向；
- <<<blocks, threads>>> 决定 grid 中有多少 block、每个 block 有多少 thread；
- 所有权仍由 host 管理，成功分配的 device memory 最后要用 cudaFree 释放。
:::
:::details{summary="基础语法：二/三维 block"}
二维 block：
```cpp
int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
if (x < cols && y < rows) {
int i = y * cols + x;
c[i] = a[i] + b[i];
}
```
通常会配一个 dim3 block(16, 16)，再分别对行列向上取整：
```cpp
dim3 block(16, 16);
dim3 grid(
(cols + block.x - 1) / block.x,
(rows + block.y - 1) / block.y
);
```
.x 方向先连续，然后每 32 个 thread 组成一个 warp。设计二维 block 时，让相邻的 .x thread 对应相邻数据，通常比较省心。
对于三维 block，thread 的线性编号是：
$$
\text{linearId} = \text{threadIdx.x}
+ \text{threadIdx.y} \times \text{blockDim.x}
+ \text{threadIdx.z} \times \text{blockDim.x} \times \text{blockDim.y}.
$$
:::
Warp、分支与数据交换
CUDA 使用 SIMT（Single Instruction, Multiple Thread）模型，但是硬件执行时仍然希望 SIMD。
Warp divergence
同一 warp 内，如果一部分 thread 走 if，另一部分走 else，硬件通常要依次执行两条路径，并暂时屏蔽不属于当前路径的 thread。
```cpp
if (condition[threadIdx.x]) {
value = path_a(value);
} else {
value = path_b(value);
}
```
if 本身不是问题，warp 内同时出现几条路径才是问题：
分支 A 与 B 依次执行，未参与当前分支的 lane 暂时失活。
* if (i < n) 边界检查一般只影响最外侧几个 warp，不必紧张。*
Warp shuffle
warp 内可以直接交换寄存器值，不必先绕到 shared memory。下面五轮 shuffle 就能完成 32 个值的求和：
```cpp
define FULL_MASK 0xffffffff
for (int offset = 16; offset > 0; offset /= 2) {
value += __shfl_down_sync(FULL_MASK, value, offset);
}
```
*偏移依次取 16、8、4、2、1，最后 lane 0 拿到完整结果。*
其它 lane 只累加了归约树上的一部分数据，不能当成完整结果使用。
内存与数据复用
GPU 上常见的内存层次：
| 内存 | 谁能看到 | 特点 |
| --- | --- | --- |
| Register | 单个 thread | 最快，但数量有限 |
| Local memory | 单个 thread | 不保证速度，数据可能落到显存 |
| Shared memory | 同一 block | 片上、显式管理，和 L1 cache 是同一块物理存储器 |
| Constant memory | 整个 grid | 有缓存快一点，适合很多 thread 读取同一小块常量 |
| Global memory | 整个 device | 容量大、延迟高，访问方式很重要 |
*越靠近计算单元，容量越小，可见范围也越窄。*
寄存器也不是越多越好。每个 thread 占用太多寄存器，SM 能同时留下的 warp 就会减少；本想少写几条指令，最后可能把延迟隐藏能力一并砍掉。
Shared memory 与同步
shared memory 最常见的用法，是让一个 block 合作把数据从 global memory 搬到片上，用很多次，再写回去。上面的[内存层次图](#内存与数据复用)给出了完整数据路径。
:::details{summary="基础语法：Shared memory 与 block 级同步"}
```cpp
__global__ void reverse_block(const int* input, int* output) {
__shared__ int tile[128];
tile[threadIdx.x] = input[threadIdx.x];
__syncthreads();
output[threadIdx.x] = tile[127 - threadIdx.x];
}
```
__shared__ 声明的数组由同一个 block 的 thread 共享。
__syncthreads() 是 block 级屏障。屏障前，每个 thread 只知道自己的写入完成了；所有 thread 通过 __syncthreads() 后，才可以安全读取别人写的数据。
:::
同一个 block 里，要参加同一次 __syncthreads() 的线程，不能分别停在控制流里的不同 barrier 位置：
```cpp
// assuming blockDim.x is 128
__global__ void syncthreads_valid_behavior(int* input_data, int* output_data) {
__shared__ int shared_data[128];
shared_data[threadIdx.x] = input_data[threadIdx.x];
if (blockIdx.x > 0) { // CORRECT, uniform condition across all block threads
__syncthreads();
output_data[threadIdx.x] = shared_data[127 - threadIdx.x];
}
}
__global__ void syncthreads_invalid_behavior(int* input_data, int* output_data) {
__shared__ int shared_data[128];
shared_data[threadIdx.x] = input_data[threadIdx.x];
for (int i = 0; i < blockDim.x; ++i) {
if (i == threadIdx.x) { // WRONG, non-uniform condition
__syncthreads(); // Undefined Behavior
}
}
output_data[threadIdx.x] = shared_data[127 - threadIdx.x];
}
```
Coalescing
一个[warp](#warp)访问 global memory 时，硬件会把地址合并成若干次内存事务。相邻 thread 读取相邻地址，通常需要的事务更少；跨步或随机访问会让带宽利用率迅速掉下去。
*上面是连续访问，下面的跨步访问碰到了更多内存区间。*
相邻 thread 访问 input[i] 容易合并；改成 input[i * stride] 后，stride 越大，32 个地址通常会散到越多内存区间。
Shared memory tiling
矩阵乘法是[shared memory](#shared-memory-与同步)分块的标准例子。一个 block 负责 $C$ 的一个小块，把 $A$ 和 $B$ 的对应小块搬进 shared memory，然后沿 $K$ 维不断累加。
$$
C_{m,n} = \sum_{k=0}^{K-1} A_{m,k} B_{k,n}.
$$
*A、B 的 tile 先搬进 shared memory，再复用来计算一个 C tile。*
一次搬进片上的 $A$、$B$ 元素会被同一个 block 内的多个 thread 反复使用，省下大量 global memory 访问。
SM 调度与 Occupancy
一个 block 被派到 SM 后，会占用 thread、warp、register 和 shared memory 配额。只要配额还有余量，同一个 SM 就能同时驻留多个 block 和很多 warp。
某个 warp 发出显存请求后，调度器会换一个就绪 warp：
*等待没有消失，SM 只是把空档交给其它 warp。*
occupancy 描述实际驻留 warp 数占硬件上限的比例：
$$
\text{Occupancy}
= \frac{\text{Active Warps per SM}}
{\text{Maximum Warps per SM}}.
$$
高 occupancy 往往更容易隐藏延迟，但它不是分数，没必要机械追求 100%。一个访存乱跳、分支发散的 kernel，occupancy 再高也不会突然变快。
真正要查的是限制来自哪里：
| 现象 | 先看什么 |
| --- | --- |
| SM 经常空闲 | grid 是否太小，驻留 warp 是否太少 |
| 显存带宽上不去 | warp 的地址分布、布局和对齐 |
| 分支效率低 | 一个 warp 内的条件分布 |
| occupancy 低 | 每 thread 寄存器、每 block shared memory、block 大小 |
| kernel 很快，程序还是慢 | 数据复制、同步和启动次数 |
表里的检查项对应着前面的[warp divergence](#warp-divergence)、[coalescing](#coalescing)和[shared memory](#shared-memory-与同步)。
cudaOccupancyMaxPotentialBlockSize 可以给 block 大小一个起点，最终配置仍然要拿真实输入计时。
从 Thread 到 Tile
手写 CUDA 时，[thread 映射](#cuda-的层级)、[shared memory](#shared-memory-与同步)、同步和边界都由程序员安排。规则张量计算会反复出现同一种[分块](#shared-memory-tiling)，于是 Triton、TileLang 选择让程序员直接描述 tile：
1. 当前程序实例负责哪块输出；
2. 输入 tile 从哪里加载；
3. tile 上做什么计算；
4. 结果写到哪里；
5. 边缘的不完整 tile 如何屏蔽。
*CUDA 显式安排 thread；Triton 保留 tile 决策，把 tile 内映射交给编译器。*
Triton
CUDA kernel 常问“我是哪个 thread”；Triton kernel 先问“我是第几个 program”。每个 program 处理一个 tile，tile 内怎样映射到 thread 和 warp，交给编译器。
向量加法
```python
import triton.language as tl
@triton.jit
def add_kernel(x, y, z, n, BLOCK_SIZE: tl.constexpr):
pid = tl.program_id(0)
offsets = pid * BLOCK_SIZE + tl.arange(0, BLOCK_SIZE)
mask = offsets < n
x_values = tl.load(x + offsets, mask=mask, other=0.0)
y_values = tl.load(y + offsets, mask=mask, other=0.0)
tl.store(z + offsets, x_values + y_values, mask=mask)
```
pid 选中当前 tile，mask 负责末尾越界。这里没有 threadIdx.x；tile 内怎样落到具体 thread，由 Triton 编译器处理。
:::details{summary="基础语法：Python 包装函数、grid 与 kernel 启动"}
```python
import torch
import triton
def add(x, y):
z = torch.empty_like(x)
n = x.numel()
grid = (triton.cdiv(n, 1024),)
add_kernel[grid](x, y, z, n, BLOCK_SIZE=1024)
return z
```
- triton.cdiv 做整数向上取整，保证最后一个不完整 tile 也会启动；
- grid 决定 program 实例数量，不直接指定 thread 数；
- add_kernel[grid](...) 把运行时参数和 tl.constexpr 编译期参数一起传入 kernel。
:::
Triton 矩阵乘法
和前面的[CUDA shared-memory tiling](#shared-memory-tiling)一样，一个 Triton program 通常负责 $C$ 的一个 BLOCK_M × BLOCK_N tile：
*沿 $K$ 维分段加载输入，每一段都累加进同一个输出 tile。*
```python
acc = tl.zeros((BLOCK_M, BLOCK_N), dtype=tl.float32)
for k0 in range(0, K, BLOCK_K):
a_ptrs = A + offs_m[:, None] * stride_am \
+ (k0 + offs_k[None, :]) * stride_ak
b_ptrs = B + (k0 + offs_k[:, None]) * stride_bk \
+ offs_n[None, :] * stride_bn
a = tl.load(a_ptrs, mask=a_mask, other=0.0)
b = tl.load(b_ptrs, mask=b_mask, other=0.0)
acc += tl.dot(a, b)
tl.store(c_ptrs, acc, mask=c_mask)
```
程序员还在决定 tile 尺寸、数据布局和 $K$ 维分块，只是不再逐个安排 thread 怎样完成 tl.dot。
Triton 很适合 elementwise、reduction、matmul 和 attention 这类规则张量计算。不过代码短不代表底层问题没了：布局不理想、寄存器压力过大或者 tile 参数不合适时，仍然要看生成代码和 profiler。特别复杂的 pipeline、warp specialization 或硬件新特性，也未必容易在高层抽象里表达。
TileLang
TileLang 也按 tile 写 kernel，但它比 Triton 更愿意把[内存层次](#内存与数据复用)摆在源码里。看代码时，能直接找到 shared tile、register fragment、数据搬运和流水线。
常见 primitive：
- T.Kernel 定义 grid，一个实例通常对应一个 CTA；
- T.alloc_shared 分配 shared memory tile；
- T.alloc_fragment 分配寄存器中的累加块；
- T.copy 搬运数据；
- T.gemm 做 tile 级矩阵乘加；
- T.Pipelined 尝试重叠不同阶段的搬运和计算。
TileLang 矩阵乘法
TileLang 的矩阵乘法仍然沿 $K$ 维推进：输入 tile 进入 shared memory，累加结果留在 register fragment，最后写回输出。
*TileLang 省掉了 thread 坐标和手写内层循环，但保留了 shared memory、寄存器与搬运关系。*
:::details{summary="TileLang 矩阵乘法的完整结构"}
```python
import tilelang
import tilelang.language as T
@tilelang.jit
def matmul(M, N, K, BLOCK_M, BLOCK_N, BLOCK_K,
dtype="float16", accum_dtype="float32"):
@T.prim_func
def main(
A: T.Tensor((M, K), dtype),
B: T.Tensor((K, N), dtype),
C: T.Tensor((M, N), accum_dtype),
):
with T.Kernel(
T.ceildiv(N, BLOCK_N),
T.ceildiv(M, BLOCK_M),
threads=128,
) as (bx, by):
a_shared = T.alloc_shared((BLOCK_M, BLOCK_K), dtype)
b_shared = T.alloc_shared((BLOCK_K, BLOCK_N), dtype)
c_local = T.alloc_fragment((BLOCK_M, BLOCK_N), accum_dtype)
T.clear(c_local)
for k in T.Pipelined(
T.ceildiv(K, BLOCK_K),
num_stages=3,
):
T.copy(A[by * BLOCK_M, k * BLOCK_K], a_shared)
T.copy(B[k * BLOCK_K, bx * BLOCK_N], b_shared)
T.gemm(a_shared, b_shared, c_local)
T.copy(c_local, C[by * BLOCK_M, bx * BLOCK_N])
return main
```
- T.Kernel 决定 CTA grid 与每个实例使用的 thread 数；
- a_shared、b_shared 是 shared memory 中的输入块；
- c_local 是寄存器里的累加结果；
- T.Pipelined 沿 $K$ 维推进，并尝试把不同阶段的搬运与计算重叠。
:::
和 Triton 相比，TileLang 多给了程序员一些结构控制；和 CUDA 相比，它又少了大量 thread 级样板。它适合那些想明确安排数据复用和 pipeline、但不想从每个 lane 写起的 kernel。
CUDA、Triton 与 TileLang 对比
| | [CUDA C++](#cuda-kernel) | [Triton](#triton) | [TileLang](#tilelang) |
| --- | --- | --- | --- |
| 主要抽象 | thread、warp、block | program、tensor tile | tile、内存层次、pipeline |
| thread 映射 | 手写 | 主要由编译器完成 | 由编译器完成，结构控制更显式 |
| shared memory | 显式分配、显式同步 | 多数细节隐藏 | 显式 tile 和搬运 |
| 适合 | 特殊机制、底层控制、极限调优 | 快速实现规则张量 kernel | 显式复用和流水线的高性能 kernel |
硬件视角
前面的 Grid → Block → Warp → Thread 是编程视角，换到硬件，对应关系更接近：
$$
\text{GPU} \rightarrow \text{GPC} \rightarrow \text{SM} \rightarrow \text{Warp Scheduler} \rightarrow \text{Execution Pipeline}.
$$
先看整块 GPU
*左侧是 CPU 与系统内存，右侧是 GPU 与显存；GPU 内的 L2 为多个 SM 提供共享数据路径。*
CPU 通过 PCIe 或 NVLink 与 GPU 通信；GPU DRAM 经过 memory controller 连接到片上 L2。再往里，多个 SM 被组织成 GPC。具体有多少 GPC、每组包含多少 SM，会随架构和芯片型号变化，写 kernel 时不应依赖这些固定数量。
这里最重要的边界是：
- L2 和显存属于整块 GPU 的共享数据路径；
- register file、L1/shared memory 和执行流水线属于单个 SM；
- block 一旦被派到某个 SM，通常会在那里驻留到执行结束；
- 不同 block 仍然不能假设先后顺序。
再看一个 SM
*左侧放大一个 SM 子分区，右侧展示完整 SM，底部保留 TMA 与 L1/shared memory 等模块。*
上图以 Blackwell Ultra 为例
| 模块 | 它负责什么 | 源码里通常对应什么 |
| --- | --- | --- |
| Warp scheduler / dispatch | 从就绪 warp 中选出下一条指令并发射 | 分支、依赖关系、eligible warp |
| Register file | 保存驻留 thread 的局部状态 | 局部变量、累加器、register pressure |
| CUDA Cores | 执行普通标量算术指令 | FP32/整数运算、普通 FMA |
| Tensor Cores | 执行矩阵形态的乘加 | GEMM、MMA、部分高层 tile 运算 |
| LD/ST units | 形成并执行 load/store 请求 | global、shared、local memory 访问 |
| SFU | 执行部分特殊函数 | 倒数、平方根倒数、三角函数等 |
| L1 / Shared memory | 缓存数据或保存 block 显式管理的 tile | cache、__shared__、数据复用 |
把一次 Kernel 启动顺着硬件走完
1. Host 启动 kernel，产生一个由许多 block 组成的 grid。
2. 硬件寻找资源足够的 SM，把一个 block 的 register、shared memory 和 warp 配额一次性分配进去。
3. Block 中的 thread 按线性 thread ID 被组成多个 warp；warp scheduler 从就绪 warp 中选择指令。
4. 算术指令进入 CUDA Core、Tensor Core 或 SFU 等流水线；访存指令进入 LD/ST 路径。
5. Global load 若未在片上命中，会继续经过 L2 到 GPU DRAM。当前 warp 等数据时，调度器切换到另一个 ready warp，这就是前面讲的[延迟隐藏](#sm-调度与-occupancy)。
6. Block 执行结束后，它占用的 register 和 shared memory 被释放，SM 才能接收新的 block。
以后读 kernel，可以在三个视角之间来回切换：
- 工作视角：grid 怎样切成 block，block 怎样切成 warp；
- 调度视角：一个 SM 上有多少 warp，其中多少真正 ready；
- 数据视角：数据经过 global、L2、L1/shared 和 register 的哪一条路径。
这三个视角把前面的[divergence](#warp-divergence)、[coalescing](#coalescing)、[shared memory](#shared-memory-与同步)和[occupancy](#sm-调度与-occupancy)放回同一张硬件图里。某个 kernel 为什么慢，通常就是其中一条路径没有被喂满，或某一种片上资源先耗尽了。