从“会写 Kernel”到“能交付系统”

一个完整 CUDA 项目不是单个快 kernel，而是一条可证明的链：

需求与误差预算
      ↓
CPU/高精度 reference + 数据集
      ↓
算法分解与数据布局
      ↓
库 / SIMT / Tile / 多 GPU 选型
      ↓
正确性测试 + Sanitizer
      ↓
端到端 benchmark + profiler
      ↓
瓶颈驱动的优化
      ↓
能力查询、fallback、部署矩阵与回归监控

下面四个用例覆盖规则 stencil、随机模拟、稀疏图和大模型推理，展示如何把前 17 章的知识连起来。

用例一：实时图像处理流水线

目标：批量 RGB 图像经过灰度化、Gaussian blur、Sobel 和 threshold，输出边缘图。要求吞吐稳定，支持任意宽高与 batch。

数据流

CPU decode
   │ pinned ring buffers
   ▼
H2D ─→ RGB→Gray ─→ Gaussian X/Y ─→ Sobel+threshold ─→ D2H/encoder
      stream 0 image k
H2D ─→ RGB→Gray ─→ Gaussian X/Y ─→ Sobel+threshold ─→ D2H/encoder
      stream 1 image k+1

先用 NPP/nvJPEG 等库建立基线。若格式或融合需求不匹配，再写自定义 kernel。

映射与布局

• block 用二维 shape，例如 32×8 或 16×16；
• 行主序图像让相邻 threadIdx.x 访问相邻像素；
• pitch/stride 是字节数，不可假定 width*sizeof(pixel)；
• blur/stencil 用 shared tile 加 halo；
• border mode 明确为 clamp、reflect、wrap 或 zero；
• RGB packed 与 planar 布局要分别 benchmark。

shared tile for radius=1

  halo halo halo halo halo halo
  halo  x    x    x    x  halo
  halo  x  output output x  halo
  halo  x  output output x  halo
  halo  x    x    x    x  halo
  halo halo halo halo halo halo

Kernel 骨架

template<int BX, int BY>
__global__ void sobel(const std::uint8_t* in, std::size_t in_pitch,
                      std::uint8_t* out, std::size_t out_pitch,
                      int width, int height) {
    __shared__ std::uint8_t tile[BY + 2][BX + 2];
    int x = blockIdx.x * BX + threadIdx.x;
    int y = blockIdx.y * BY + threadIdx.y;

    // 协作加载中心和 halo；所有线程必须到达 barrier。
    load_clamped_tile(tile, in, in_pitch, x, y, width, height);
    __syncthreads();

    if (x < width && y < height) {
        int gx = gradient_x(tile, threadIdx.x + 1, threadIdx.y + 1);
        int gy = gradient_y(tile, threadIdx.x + 1, threadIdx.y + 1);
        row_ptr(out, out_pitch, y)[x] = saturate(abs(gx) + abs(gy));
    }
}

重点不是复制辅助函数，而是 barrier 必须位于边界 if 外；否则边缘 block 的一部分线程提前退出会导致未定义行为。

优化次序

1. 先把 decode/transfer/kernel/encode 分段计时；
2. pinned buffers + 多 stream 重叠多张图；
3. 消除 RGB→Gray 与 blur 间的中间 round trip；
4. shared halo 降低重复 global load；
5. 批量稳定后用 CUDA Graph 降低 launch overhead；
6. 只有 profiler 显示算术热点，再考虑 vectorized load 或 fast math。

测试矩阵

用例二：Monte Carlo 期权定价

目标：模拟几何布朗运动终值，计算 discounted payoff 的均值和置信区间。

$$S_T=S_0\exp\left((r-\tfrac12\sigma^2)T+\sigma\sqrt{T}Z\right), \quad Z\sim\mathcal N(0,1)$$

$$\hat V=e^{-rT}\frac1N\sum_{i=1}^N \max(S_T^{(i)}-K,0)$$

数据流与算法

sample id
   ↓ counter/subsequence RNG
normal Z
   ↓ path/payoff（可与 RNG 融合）
payoff
   ├─ sum(payoff)
   └─ sum(payoff²)
          ↓
mean + variance + confidence interval

大数组不是必需：每个线程可生成多个样本并做局部累加，再用 CUB reduction 合并 {sum, sumsq}。这样减少显存流量。

数值与统计正确性

• 累加采用 double 或经过验证的分层 float；
• variance = E[x²]-E[x]² 在方差很小时可能发生消减，必要时用更稳定算法；
• 固定 seed 只保证可复现实验，不证明随机质量；
• 每个逻辑样本拥有不重叠 sequence/counter；
• 与 Black-Scholes 闭式解对照，误差应随$1/\sqrt{N}$ 统计下降；
• 报告置信区间而不是只报一个浮点值。

并行粒度

thread i:
  local_sum = 0
  repeat K samples:
      z = rng(sample_id)
      payoff = simulate(z)
      local_sum += payoff
  block_reduce(local_sum)

grid result → second-stage reduce / CUB

每线程多个样本能摊销 RNG state 和 launch 成本，但增加寄存器及串行工作。扫描 K 并看 occupancy、instruction throughput 和总时间。

性能实验

用例三：稀疏图 BFS

目标：在 CSR 图上从源点计算无权最短层级。它展示 GPU 不擅长的负载不均、随机内存与动态 frontier。

CSR 与 Frontier

row_offsets[v]..row_offsets[v+1] → v 的邻边

frontier level L: [v2, v8, v9]
      │ expand neighbors
      ▼
candidate vertices
      │ visited atomic / duplicate removal / compact
      ▼
frontier level L+1

基线 Kernel

__global__ void expand(const int* frontier, int frontier_size,
                       const int* row, const int* col,
                       int* distance, int level,
                       int* next, int* next_size) {
    for (int i = blockIdx.x * blockDim.x + threadIdx.x;
         i < frontier_size;
         i += blockDim.x * gridDim.x) {
        int v = frontier[i];
        for (int e = row[v]; e < row[v + 1]; ++e) {
            int u = col[e];
            if (atomicCAS(&distance[u], -1, level + 1) == -1) {
                int pos = atomicAdd(next_size, 1);
                next[pos] = u;
            }
        }
    }
}

它正确直观，但高 degree vertex 由一个线程串行处理，next_size 又形成全局原子热点。

逐步优化

1. Warp-centric：一个 warp 合作处理高 degree vertex；
2. 分级调度：小 degree 每线程、中 degree 每 warp、大 degree 每 block；
3. Block queue：先写 shared queue，再批量 reserve global slots；
4. Prefix-sum compaction：用 scan 替代逐元素 global atomic；
5. Direction optimization：frontier 很大时从 unvisited 顶点检查 incoming 邻居；
6. Persistent/work stealing：degree 极不均时动态分配边范围。

degree distribution
small vertices ─→ thread lanes
medium vertices ─→ whole warp
large vertices ─→ whole block / split edge ranges

指标

BFS 的“有效工作”随图结构变化，不能只报 kernel 时间。至少报告：

• traversed edges per second；
• 每层 frontier size 与时间；
• 重复 candidate 比例；
• 原子次数与成功率；
• global load efficiency、cache hit、warp divergence；
• 与 CPU reference 的 distance 全量比较。

道路网、均匀随机图和幂律社交图表现完全不同，必须覆盖多种 degree 分布。

用例四：LLM Attention 推理片段

目标：解释 prefill 与 decode 下 attention 为什么需要不同 kernel，并把 GEMM、online softmax、KV cache 和多 GPU 串联起来。

$$O=\operatorname{softmax}\left(\frac{QK^T}{\sqrt d}+M\right)V$$

朴素数据流

Q,K ─→ GEMM ─→ score[N,N] ─→ mask/softmax ─→ P[N,N]
                                                │
P,V ───────────────────────────────────────────→ GEMM → O

中间 score 与 P 的$O(N^2)$ 显存流量是主要问题。FlashAttention 类算法按块计算 score，维护在线 softmax 统计量，不把完整矩阵写回 global memory。

Online Softmax 合并

对已经处理部分的最大值$m$、指数和$l$ 和输出累加$o$，加入新 tile 后：

$$m' = \max(m, m_{tile})$$

$$l' = e^{m-m'}l + e^{m_{tile}-m'}l_{tile}$$

$$o' = e^{m-m'}o + e^{m_{tile}-m'}o_{tile}$$

最终输出为$o/l$。这允许流式处理 K/V tiles，同时保持数值稳定。

Prefill 与 Decode

Decode 中 KV cache 的布局和分页元数据往往比单个矩阵指令更重要。

logical sequence pages
token 0..15  → physical page 7
token 16..31 → physical page 2
token 32..   → physical page 9

block table → gather K/V tiles → attention

工程组合

• QKV projection 和 output projection 优先用 cuBLASLt/CUTLASS；
• attention 核心可用优化框架、CUTLASS 或 Tile DSL，不从朴素 kernel 开始生产化；
• RoPE、bias、mask、dequant 可在带宽允许时融合；
• CUDA Graph 降低 decode 的重复 launch 开销；
• KV pages 用 memory pool/VMM 管理时要区分逻辑容量和物理驻留；
• tensor parallel 使用 NCCL collective，流水线并行还要处理 stage 间通信；
• 量化必须同时定义权重、activation、scale、accumulator 和输出 dtype。

评估指标

不要只报某个 attention kernel 的 TFLOP/s：

• time to first token；
• inter-token latency；
• tokens/s 与 batch/sequence 的关系；
• KV cache bytes/token 和最大并发序列；
• HBM bandwidth、Tensor Core utilization；
• 数值误差、困惑度或任务质量；
• 多 GPU scaling 与通信占比。

一个可复用的项目目录

cuda-project/
├─ CMakeLists.txt
├─ include/
│  ├─ error_check.hpp
│  ├─ device_buffer.hpp
│  └─ api.h                 # 稳定 C ABI
├─ src/
│  ├─ kernels.cu
│  ├─ dispatch.cpp
│  └─ reference.cpp
├─ tests/
│  ├─ correctness.cpp
│  ├─ edge_cases.cpp
│  └─ sanitizer_cases.cu
├─ benchmarks/
│  ├─ bench_kernel.cpp
│  └─ datasets.json
├─ tools/
│  ├─ profile_nsys.ps1
│  └─ profile_ncu.ps1
└─ docs/
   ├─ numerical_contract.md
   └─ performance_report.md

测试与 benchmark 分开：测试追求覆盖错误，benchmark 追求稳定测量；两者都不应把首次 context/JIT 开销意外混入或排除。

Feature Query 与 Dispatch

编译成功不等于当前 GPU 能执行全部路径。启动时建立 capability 表：

struct Capabilities {
    int cc_major;
    int cc_minor;
    bool managed;
    bool concurrent_managed;
    bool memory_pools;
    bool cooperative_launch;
    bool clusters;
    bool vmm;
};

input problem + numerical policy + device capabilities
                         │
                         ▼
dispatch
├─ optimized_sm100
├─ optimized_sm90 (cluster/TMA)
├─ generic_sm80 (async copy)
└─ portable fallback

不要仅用 GPU 商品名 dispatch；使用 Compute Capability、API attribute 和库 heuristic。MIG 配置还可能改变可用 cluster/SM 资源。

性能回归门禁

每个基准保存：

• GPU UUID/型号、CC、功耗/时钟策略；
• Driver、Toolkit、库、编译器与 git revision；
• 输入 shape、dtype、分布和随机 seed；
• 编译参数、算法 ID、workspace；
• warmup、repeat、统计量；
• kernel time、端到端 time、bytes/FLOPs 和误差。

门禁用区间而非单次值。例如以稳定机器的 median 与 P95 比较，只有超过噪声阈值才失败。共享开发机上的自动频率和其他进程会制造假回归。

生产故障分层排查

程序无法启动
└─ Driver/Toolkit/架构/动态库/容器挂载 → 01、16

结果错误或偶现
└─ 越界/竞态/同步/生命周期/数值 → 05、09、11、12

Kernel 慢
└─ 算法/bytes/occupancy/stall/布局 → 04、06、08

端到端慢但 kernel 快
└─ launch/分配/传输/同步/通信 → 05、10、14、17

新架构能力未生效
└─ capability、编译目标、dispatch、生成指令 → 13、15、16

CUDA 术语速查

API 速查

毕业检查表

学完后，应能对一个 CUDA 项目给出以下可验证答案：

1. 工作如何映射到 grid/block/warp 或 tile program？
2. 每个输出元素需要搬运多少字节、做多少运算？
3. 哪些数据在 register/shared/L2/HBM/host/fabric 中流动？
4. 所有跨线程、stream、device、进程依赖如何建立 happens-before？
5. 数值误差预算与 reference 是什么？
6. 为什么选择库、SIMT、Tile 或多 GPU？
7. profiler 证据指向哪个瓶颈？
8. 优化对端到端时间而非单 kernel 有何影响？
9. 不支持某项能力的 GPU 走哪条 fallback？
10. 构建、测试、部署与性能结果能否在另一台机器复现？

官方参考

• CUDA Programming Guide: Feature Survey
• CUDA C++ Best Practices Guide
• CUDA Samples
• CUDA GPU Compute Capability
• Nsight Systems
• Nsight Compute
• Compute Sanitizer