CUDA：线程、内存与 kernel

在 PyTorch 里写一句 c = a + b，GPU 上到底发生了什么？CPU 把一个叫 kernel 的函数交给 GPU，GPU 一次启动成千上万个线程，每个线程只算一两个元素。加法本身很便宜；一个 kernel 快不快，多半看每个 warp 搬来的字节有没有用上，以及同一份数据有没有被反复从显存读出来。 这一篇用一个 vector add 和一个 4×4 的矩阵乘法把这两件事算清楚，再接到 cuBLAS、CUTLASS 和 Triton。

待审核草稿。 这一篇还没有经过站主审阅，事实、取舍和措辞都可能修改；引用前请回到文中链接的一手资料核对。

文中的规则核对自 NVIDIA CUDA Programming Guide v13.4.2（2026-09-10 更新，旧名 CUDA C++ Programming Guide，下称 Guide）和 CUDA C++ Best Practices Guide 13.4（下称 BPG）。GPU 为什么做成这样，见 GPU 架构；怎样判断一个 kernel 被算力还是带宽卡住，见下一篇性能。

CPU 发起，GPU 执行

CUDA 把系统分成两边：CPU 和它的内存叫 host，GPU 和它的显存叫 device。程序总是从 CPU 开始运行。host 代码负责分配显存、在两边之间拷数据、启动 GPU 上的函数，再等结果。在 GPU 上执行、由 host 启动的函数叫 kernel，启动这个动作叫 launch（Guide §1.2.1）。

Launch 对 CPU 是异步的：kernel<<<...>>>(...) 返回时，kernel 可能还没开始跑（Guide §2.1.2）。所以不先同步就计时，量到的可能只是把任务放进队列的时间。PyTorch 训练循环里的同步与计时问题，见训练循环。

线程怎样分组：thread、warp、block、grid

一次 launch 创建一个 grid；grid 由许多同样大小的 thread block 组成，block 里是 thread。执行时，block 内每 32 个线程编成一个 warp，一起执行同一条指令（SIMT）。

这些规定见 Guide 的编程模型和 thread block cluster 两节。Compute capability 9.0 对应 Hopper（例如 H100）；更早的 GPU 没有 cluster 这一层。SM（streaming multiprocessor）是 GPU 上执行 block 的单元；按 H100 架构白皮书 V1.04，H100 SXM5 有 132 个。

一维时，每个线程用同一条公式算出自己负责的下标。记 $b$ 为 blockIdx.x，$D$ 为 blockDim.x，$t$ 为 threadIdx.x：

$$ i=b\cdot D+t. $$

取 $n=1000$ 个元素、每个 block 256 个线程，需要 $\lceil 1000/256\rceil=4$ 个 block，一共 1024 个线程。第 3 个 block（$b=2$）里 $t=5$ 的线程负责 $2\times256+5=517$。最后一个 block 覆盖 768–1023，其中 1000–1023 这 24 个线程越界，要靠 if (i < n) 跳过；它们恰好都在这个 block 最后一个 warp（992–1023）里，所以这个 warp 只有 8 条 lane 在干活。Block 大小取 32 的倍数，就是为了不在每个 block 末尾留下半空的 warp（Guide §1.2.2.2）。

一个完整的 vector add

下面是可以单独编译的最小程序（nvcc vec_add.cu -o vec_add），结构与 Guide §2.1.3 的显式内存管理示例一致。本页没有在 GPU 上编译运行它。

#include <cstdio>
#include <cuda_runtime.h>

__global__ void vecAdd(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];
    }
}

int main() {
    const int n = 1000;
    const size_t bytes = n * sizeof(float);

    float *h_a, *h_b, *h_c;
    cudaMallocHost(&h_a, bytes);
    cudaMallocHost(&h_b, bytes);
    cudaMallocHost(&h_c, bytes);
    for (int i = 0; i < n; ++i) { h_a[i] = i; h_b[i] = 2.0f * i; }

    float *d_a, *d_b, *d_c;
    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);

    int threads = 256;
    int blocks = (n + threads - 1) / threads;
    vecAdd<<<blocks, threads>>>(d_a, d_b, d_c, n);
    cudaError_t err = cudaGetLastError();
    if (err == cudaSuccess) err = cudaDeviceSynchronize();
    if (err != cudaSuccess) {
        printf("CUDA error: %s\n", cudaGetErrorString(err));
        return 1;
    }

    cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);
    printf("c[999] = %f\n", h_c[999]);

    cudaFree(d_a); cudaFree(d_b); cudaFree(d_c);
    cudaFreeHost(h_a); cudaFreeHost(h_b); cudaFreeHost(h_c);
    return 0;
}

错误检查的写法来自 Guide §2.1.7：kernel launch 本身不返回错误码，执行期错误要等下一次同步或 API 调用才报出来。正式代码应检查每个 CUDA API 的返回值，这里为了短，只检查了 kernel。

这个 kernel 每个元素读 8 字节、写 4 字节，只做 1 次加法。它的速度几乎完全由显存带宽决定，性能篇会把这句话换成数字。

内存：先分清放在哪里、谁能看见

前四列来自 Guide 的内存空间一节。最后一列的寄存器、shared memory 和 constant 上限取 Guide 附录 5.1 中 compute capability 9.0 一栏；显存容量与带宽取 H100 架构白皮书 V1.04 表 3 的 H100 SXM5。每个 block 用超过 48 KB 的 shared memory 时，必须改用动态分配并显式 opt-in。

各级存储的带宽和容量怎样对比，见 GPU 架构。此外，每个 SM 有 L1 cache，它和 shared memory 共用同一块 unified data cache，比例可以按 kernel 调整；全 GPU 共享一个 L2（H100 SXM5 为 50 MB）。Cache 由硬件管理，shared memory 则由程序员显式读写，可以当成一块“自己管的 cache”。

写 kernel 时最常见的一步，就是把会被反复用到的数据从 global 搬进 shared memory 或寄存器，在片上复用。下面两节分别回答：一次搬运有没有浪费；同一份数据被搬了几次。

Coalescing：一个 warp 一次读多少字节

Global memory 以 32 字节为单位传输。一个 warp 执行同一条 load 时，硬件把 32 个线程的地址合并（coalesce），需要几个 32 字节段，就发几个 transaction（compute capability 6.0 及以上，见 BPG §10.2.1 和 Guide §2.3.4.1）。下表假设每个线程读一个元素，数组起点按 32 字节对齐（cudaMalloc 至少保证 256 字节对齐）：

Stride 32 时，相邻线程相隔 128 字节，每个线程独占一个 32 字节段：搬了 1024 字节，只用上 128 字节。按列去读一个行优先（row-major）、每行至少 32 个 float 的矩阵，就是这种情况；Guide 的矩阵转置例子里，未优化版本的写入正是如此。

Coalescing 不要求“第 k 个线程读第 k 个元素”。只要一个 warp 的地址落在尽量少的 32 字节段里，段内顺序打乱也还是 4 个 transaction（BPG）。另外，表里数的是 transaction，不是耗时：相邻 warp 可能复用 cache 里多拉进来的数据。BPG 在 V100 上测“起点错开”的拷贝，带宽降到不错开时的约 9/10，而不是按 4/5 推算的那么低。

Shared memory tiling：同一份数据别从 global 读很多遍

矩阵乘 $C=AB$，三个矩阵都是 $n\times n$。最直接的 kernel 让每个线程算一个 $C_{ij}$：读 $A$ 的第 $i$ 行和 $B$ 的第 $j$ 列，各 $n$ 个元素。

取 $n=4$。16 个输出各读 8 个元素，共 128 次 global 读取；但 $A$ 和 $B$ 一共只有 32 个不同的元素，每个都被读了 4 遍。比如 $A_{ik}$ 会被第 $i$ 行的 4 个输出各读一次。

分块（tiling）把 $C$ 切成 $T\times T$ 的块，每个 block 负责一块。取 $T=2$：4 个 block，每个 block 有 $2\times2$ 个线程。沿求和方向分 $n/T=2$ 个阶段；每个阶段，block 先把 $A$ 的一个 $2\times2$ 子块和 $B$ 的一个 $2\times2$ 子块搬进 shared memory（各读 4 次），同步之后在片上做乘加。

一般地，分块把 global 读取从 $2n^3$ 降到 $2n^3/T$：每个元素被读的次数从 $n$ 遍降到 $n/T$ 遍。$n=4096$、$T=32$、FP32 时，不分块要发出约 $1.37\times10^{11}$ 次读取（约 550 GB），分块后约 $4.3\times10^{9}$ 次（约 17 GB），而三个矩阵本身只有约 0.2 GB。

#define TILE 32

__global__ void tiledMatMul(const float* A, const float* B, float* C, int n) {
    __shared__ float As[TILE][TILE];
    __shared__ float Bs[TILE][TILE];
    int row = blockIdx.y * TILE + threadIdx.y;
    int col = blockIdx.x * TILE + threadIdx.x;
    float sum = 0.0f;
    for (int t = 0; t < n / TILE; ++t) {
        As[threadIdx.y][threadIdx.x] = A[row * n + t * TILE + threadIdx.x];
        Bs[threadIdx.y][threadIdx.x] = B[(t * TILE + threadIdx.y) * n + col];
        __syncthreads();
        for (int k = 0; k < TILE; ++k)
            sum += As[threadIdx.y][k] * Bs[k][threadIdx.x];
        __syncthreads();
    }
    C[row * n + col] = sum;
}

这是教学版：假设 $n$ 是 TILE 的倍数，用 dim3 block(TILE, TILE) 和 dim3 grid(n / TILE, n / TILE) 启动。第一个 __syncthreads() 保证整块搬完才开始用；第二个保证所有线程都用完了，下一阶段才覆盖。读 A 和 B 时，threadIdx.x 都落在地址的最低维，相邻线程读相邻地址，所以搬运本身也是 coalesced 的。

BPG 用同样的思路优化 $C=AB$（§10.2.3.2），在 Tesla V100 上把有效带宽从 119.9 GB/s 提到 195.5 GB/s，并指出提升来自去掉重复读取，而不是合并访问变好。读取次数降为原来的 1/32，速度却远没有快 32 倍：一部分重复读取本来就被 cache 挡住了；shared memory 也有自己的问题，它分成 32 个 bank，同一个 warp 访问同一 bank 的不同地址会被拆开串行执行。读取次数说明的是机制，不是耗时预测。

真正的 GEMM 库会继续分层：block 级的 tile 放 shared memory，warp 级的 tile 放寄存器，最内层交给 Tensor Core 的 MMA 指令，并把下一块的搬运藏到当前块的计算后面（CUTLASS 的 GEMM 层次说明）。FlashAttention 把同样的分块思想用到 attention 上，再配合 online softmax，见 FlashAttention 和 PagedAttention。

Occupancy：一个 SM 上同时能放几个 warp

一个 warp 等 global memory 的数据时，SM 可以切换去执行别的 warp。同时驻留的 warp 越多，越有机会把等待藏起来；这套延迟隐藏的硬件原理见 GPU 架构。Occupancy = SM 上活跃的 warp 数 ÷ 这个 SM 最多能驻留的 warp 数（Guide §2.3.7）。

一个 SM 还能不能再放一个 block，看几项资源里谁先用完。H100 属于 compute capability 9.0，Guide 附录 5.1 给出的每 SM 上限是：

设每个 block 256 个线程（8 个 warp），逐项去除：

可以用几行 Python 复算这张表：

def resident_blocks(threads, regs_per_thread, smem_kb_per_block):
    limits = [32, 2048 // threads, 65536 // (regs_per_thread * threads)]
    if smem_kb_per_block:
        limits.append(228 // smem_kb_per_block)
    return min(limits)

for regs, smem in [(32, 0), (64, 0), (32, 48), (80, 0)]:
    blocks = resident_blocks(256, regs, smem)
    print(regs, smem, blocks, blocks * 256 // 32 / 64)

这个估算省掉了两处细节。第一，寄存器按 warp 以 256 个为粒度分配（BPG §11.1.1 在 compute capability 7.0 的例子里说明）；上面选的寄存器数让每个 warp 的用量正好是 256 的倍数。第二，每个 block 能用的 shared memory（227 KB）比每个 SM 的上限（228 KB）少 1 KB，CUDA 12.4 版 Guide 说明这 1 KB 留给系统使用；按每个 block 多占 1 KB 算（49 KB），仍是 4 个 block，结论不变。实际 kernel 先用 nvcc --resource-usage 看寄存器和 shared 用量，再用 Nsight Compute 的 occupancy calculator 或 cudaOccupancyMaxActiveBlocksPerMultiprocessor 核对。

Occupancy 不是越高越好。BPG 写明，超过某个点后提高 occupancy 不会再变快；少驻留几个 warp、让每个线程多用寄存器，有时能靠指令级并行把等待藏住（BPG §11.1）。Occupancy 回答的是“有没有足够的 warp 可以切换”，不是“这个 kernel 快不快”。

Streams：让拷贝和计算重叠

Stream 是一条按顺序执行的工作队列：同一个 stream 里的拷贝和 kernel 依次完成，不同 stream 之间可以交错，甚至同时进行（Guide §2.5.2）。要让 host 与 device 之间的拷贝真正和 kernel 重叠，需要同时满足：host 缓冲区是 pinned 内存；用 cudaMemcpyAsync；拷贝和 kernel 放在不同的非默认 stream 里。条件来自 BPG §10.1.2。用普通的 pageable 内存时，cudaMemcpyAsync 仍然正确，但会退化成同步行为，不再与别的工作重叠（Guide §2.5.2.3）。

BPG 的分段写法是把数据切成 nStreams 段，每段的“拷贝 → 计算”放进自己的 stream。下面假设 h_x 是 pinned 内存，stream[] 已经用 cudaStreamCreate 创建，n 能被 nStreams * 256 整除：

size_t chunk = n / nStreams;
for (int s = 0; s < nStreams; ++s) {
    size_t off = s * chunk;
    cudaMemcpyAsync(d_x + off, h_x + off, chunk * sizeof(float),
                    cudaMemcpyHostToDevice, stream[s]);
    process<<<chunk / 256, 256, 0, stream[s]>>>(d_x + off, chunk);
}
cudaDeviceSynchronize();

手算一次：整段拷贝要 $t_T=8$ ms，整段计算要 $t_E=10$ ms，顺序做是 18 ms。切成 4 段后，第 1 段拷完就能开始算，后面的拷贝藏在计算后面。BPG 给的粗略估计是 $t_E+t_T/n_{\rm streams}=10+2=12$ ms（适用于 $t_E>t_T$；反过来约为 $t_T+t_E/n_{\rm streams}$）。这些毫秒数是教学假设。能不能同时做两个方向的拷贝，还取决于设备有几个 copy engine（asyncEngineCount）。

在 PyTorch 里，对应的做法是 DataLoader 设 pin_memory=True，再用 .to(device, non_blocking=True)（CUDA semantics）。多 stream 下张量生命周期的坑，见训练循环。

“边搬边算”也发生在 kernel 内部：compute capability 8.0 起，global 到 shared 的拷贝可以异步发出（LDGSTS）；9.0 起有 Tensor Memory Accelerator（TMA），专门批量搬运多维数组的子块（Guide §4.12）。GEMM 和 FlashAttention 这类 kernel 靠它们把下一块数据的搬运藏到当前块的计算后面。

手写 kernel 之外：cuBLAS、cuDNN、CUTLASS 与 Triton

日常训练和推理很少需要手写 CUDA。cuBLAS 是建立在 CUDA runtime 上的 BLAS 实现（另有 cuBLASLt 等接口），PyTorch 在 CUDA 上做 BLAS 运算时默认调用 cuBLAS（见 preferred_blas_library 的说明）。cuDNN 提供 attention、卷积、矩阵乘、归一化等深度学习原语，并支持把多个运算按图融合。CUTLASS 是开源的 CUDA C++ 模板库，把上面那种分层 tiling 拆成可以组合的部件；CUTLASS 4 起还有 Python 的 CuTe DSL，截至 2026 年 10 月，README 列出的最新版本是 2026 年 9 月发布的 4.8.0。Guide 也建议通过这些库，而不是内联 PTX 去用 Tensor Core（附录 5.1）。

Triton（Tillet、Kung、Cox，MAPL 2019）换了一层抽象：程序员按 block 写程序、直接操作整块数据，coalescing、shared memory 的分配与同步、异步拷贝的调度和 Tensor Core 指令的选择交给编译器。PyTorch torch.compile 的默认后端 TorchInductor 在 NVIDIA、AMD 和 Intel GPU 上用 Triton 生成 kernel（PyTorch 2.14 文档）。CUDA 自己在 v13 的文档里也加入了类似的 tile 编程模型：Python 有 cuTile，CUDA Tile C++ 从 Toolkit 13.3 起提供（Guide §2.4）。不管用哪一层，本篇的三个问题都还在：每个 warp 搬的字节用上没有，数据复用了几次，SM 上有没有足够的 warp 可以切换。这些库怎样被框架串起来，见框架。

动手时先试什么

Sectors/request 的理想值引自 Nsight Compute Profiling Guide。本页没有在 GPU 上跑这些实验；第二列是根据上面引用的文档推出的预期，要在自己的硬件上验证。

接着读性能：roofline、算术强度与 MFU，把这里的“搬了多少字节”换算成时间上限。