为什么“先调库”是性能工程原则

成熟 CUDA 库通常同时处理：

• 不同 GPU 架构与数据类型的 kernel dispatch；
• 小/中/大尺寸和特殊 shape；
• 对齐、尾块、数值稳定与特殊值；
• 多 stream、workspace、heuristic 与算法选择；
• 多年积累的性能调优和回归测试。

自写 kernel 的合理理由是库没有对应算法、需要跨算子融合、数据布局特殊，或 profile 已证明库调用边界是瓶颈。即便如此，也应把库结果作为正确性和性能基线。

生态地图

库名相近不表示抽象相同：例如 cuFFT 是 host 端计划并启动 kernel；cuFFTDx 把 FFT primitive 嵌入自定义 kernel，以便融合。

CUB：device-wide primitive 的标准基线

DeviceReduce 采用“两次调用”模式：第一次查询临时空间，第二次执行。

#include <cub/cub.cuh>

void* temp = nullptr;
std::size_t temp_bytes = 0;

CUB_CHECK(cub::DeviceReduce::Sum(
    temp, temp_bytes, d_in, d_out, n, stream));
CUDA_CHECK(cudaMallocAsync(&temp, temp_bytes, stream));
CUB_CHECK(cub::DeviceReduce::Sum(
    temp, temp_bytes, d_in, d_out, n, stream));
CUDA_CHECK(cudaFreeAsync(temp, stream));

生产代码不要每次调用都查询并分配。把 temp size 纳入 plan/cache，workspace 由调用方或 memory pool 复用。还要检查当前 CUB overload 的参数顺序与返回类型，因为 CCCL 接口会演进。

CUB 的层次：

Thread primitives
Warp primitives：WarpReduce / WarpScan / exchange
Block primitives：BlockLoad / BlockReduce / BlockScan / radix sort
Device algorithms：DeviceReduce / Scan / Select / Sort

写自定义 kernel 时优先组合 block/warp primitive，而不是重写易错的 shuffle、bank conflict 和临时存储逻辑。

Thrust：算法组合与迭代器

#include <thrust/execution_policy.h>
#include <thrust/transform_reduce.h>

struct square {
    __host__ __device__ float operator()(float x) const { return x * x; }
};

float sum = thrust::transform_reduce(
    thrust::cuda::par.on(stream),
    begin, end, square{}, 0.0f, thrust::plus<float>{});

显式 execution policy 和 stream，避免算法落到错误后端或默认 stream。注意返回 host scalar 的算法可能需要同步；若流水线要保持异步，优先选择输出到 device 的接口或 CUB primitive。

Fancy iterators 可避免中间数组：counting、transform、permutation、zip iterator 能把 map/gather/多数组视图融合进算法。但过深 iterator 表达式会让编译错误和调试困难，关键路径应保留小型测试。

cuBLAS 基本 GEMM

数学目标：

$$C = \alpha op(A)op(B) + \beta C$$

cublasHandle_t handle{};
CUBLAS_CHECK(cublasCreate(&handle));
CUBLAS_CHECK(cublasSetStream(handle, stream));

CUBLAS_CHECK(cublasSgemm(
    handle,
    CUBLAS_OP_N, CUBLAS_OP_N,
    m, n, k,
    &alpha,
    A, lda,
    B, ldb,
    &beta,
    C, ldc));

经典 cuBLAS API 按 column-major 解释矩阵。C/C++ row-major 数据可通过交换操作数/转置标志进行等价映射，但必须从下标公式推导，不能只“把 m、n 对调试试”。cuBLASLt 对 layout、epilogue 和算法选择提供更灵活描述。

cuBLASLt 的对象模型

Matmul descriptor
  ├─ compute type / scale type
  ├─ transpose
  └─ epilogue（bias、activation 等，按版本支持）

Matrix layout A/B/C/D
  ├─ dtype, rows, cols, leading dimension
  └─ batch/stride/order

Preference
  └─ max workspace bytes

heuristic → algorithm candidates → benchmark representative shapes → cache choice

Heuristic 返回候选，不保证在所有实际尺寸上最优。生产系统应按 (m,n,k,dtype,layout,transpose,epilogue,arch,workspace limit) 建 tuning cache，并为无候选或 workspace 不足保留 fallback。

数学模式与复现性

库可能为了吞吐选择 Tensor Core、TF32、reduced-precision reduction 或不同 reduction tree。明确配置：

• storage type、compute type、scale type；
• 是否允许 TF32 或低精度 reduction；
• workspace 与算法选择；
• 是否要求 deterministic/reproducible；
• NaN/Inf 与舍入行为要求。

“同一个 API”在不同 GPU、workspace 和版本下可能选择不同算法，导致性能和末位结果变化。严格复现需要阅读每个库的 reproducibility 说明并固定更多环境。

cuFFT：Plan 是性能对象

cufftHandle plan{};
CUFFT_CHECK(cufftPlan1d(&plan, n, CUFFT_C2C, batch));
CUFFT_CHECK(cufftSetStream(plan, stream));
CUFFT_CHECK(cufftExecC2C(plan, data, data, CUFFT_FORWARD));
CUFFT_CHECK(cufftDestroy(plan));

真实服务中应缓存 plan，而不是每次创建。高级 layout 用 cufftPlanMany/Xt 接口描述 embed、stride、distance 和 batch。

FFT pipeline 常见结构：

window/pack → R2C FFT → frequency-domain multiply → C2R FFT → scale/unpack

要检查：in-place 布局需要的 padding、实数变换的 Hermitian 存储、inverse 是否自动归一化、plan workspace、callback/graph capture 支持，以及 prime size 是否导致慢算法。

cuSPARSE：格式比 API 更重要

同一个稀疏矩阵可表示为 CSR、CSC、COO、BSR、SELL 等。性能取决于行长度分布、块结构、排序和操作方向。

CSR
row_offsets: [0, 2, 3, 6]
col_indices: [0, 4, 2, 0, 1, 4]
values:      [a, b, c, d, e, f]

row 0 → entries [0,2)
row 1 → entries [2,3)
row 2 → entries [3,6)

通用 API 通常使用 descriptor：

cusparseSpMatDescr_t matA{};
cusparseDnVecDescr_t vecX{}, vecY{};
// 创建 descriptor，查询 bufferSize，分配 workspace，再执行 cusparseSpMV。

与 CUB 相同，查询 workspace 后应复用。输入是否排序、index 32/64 位、index base、transpose、算法枚举和预处理阶段都会影响性能与正确性。

Sparse 选型流程

矩阵真的稀疏吗？
├─ 否 → dense GEMM 往往更快
└─ 是
   ├─ 规则小块非零 → BSR / block sparse
   ├─ 行长度相近   → SELL/ELL 类布局可评估
   ├─ 通用行稀疏   → CSR 基线
   └─ 结构化 2:4   → 检查 cuSPARSELt 与硬件支持

计算 FLOP/s 时使用实际非零运算量，同时报告预处理和格式转换是否摊销。只测 SpMV kernel 而忽略每轮重排，会夸大端到端收益。

cuSOLVER

常见任务包括 LU、QR、Cholesky、SVD、特征值和稀疏求解。求解器通常有：

1. 参数/handle 与 stream；
2. workspace size query；
3. factorization；
4. solve 或提取结果；
5. device info 检查。

info 是算法结果的一部分，不能只检查 API status。矩阵奇异、非正定或未收敛可能通过 info 报告，而 API 调用本身成功。

cuRAND：并行随机数的状态设计

两种路线：

• Host API 批量生成到 device buffer，易用且吞吐高；
• Device API 在自定义 kernel 内维护 generator state，便于融合。

可复现性由 generator、seed、sequence/subsequence、offset 与线程映射共同决定：

global seed
  ├─ logical stream 0 → subsequence 0
  ├─ logical stream 1 → subsequence 1
  └─ logical stream k → subsequence k

不要用线程的瞬时时间作为 seed，也不要让不同线程复用同一 state。改变 block size 后若逻辑样本映射改变，随机序列也可能变化；需要与 launch 几何无关的复现时，用样本 ID 定义 counter/subsequence。

NCCL 与 NVSHMEM

NCCL 提供针对 GPU 拓扑优化的 collective 和点对点通信，例如 AllReduce、AllGather、ReduceScatter。NVSHMEM 提供 GPU 可发起的 PGAS 通信、put/get、原子和同步。

数据并行梯度同步：
GPU0 ─┐
GPU1 ─┼─ ReduceScatter → 每卡一段结果 → AllGather → 每卡完整结果
GPU2 ─┤
GPU3 ─┘

如果目标是标准 collective，不要先自己用 P2P/VMM 写 ring。只有特殊通信融合、设备侧细粒度生产消费或库不覆盖的拓扑协议，才值得下沉。

Library Handle 与 Stream 规则

多数库 handle 带有可变状态：当前 stream、math mode、pointer mode、workspace 等。不要在多个 host thread 中无锁修改并共享一个 handle。

常见模式：

每 device / 每 host worker
  └─ handle bundle
      ├─ cuBLAS/cuBLASLt
      ├─ cuSPARSE/cuSOLVER
      └─ associated stream + reusable workspace

若库支持传入用户 workspace，优先由上层 allocator 管理，减少隐藏分配和跨库内存峰值。Graph capture 前确认该 API、算法和内部分配路径是否 capture-safe。

Pointer Mode 与标量生命周期

cuBLAS 等库可让 alpha/beta 位于 host 或 device。Host pointer mode 下 API 可能在调用时读取标量；device pointer mode 可保持异步并让标量来自前序 GPU 工作。

host pointer mode:   CPU scalar ─→ library call
device pointer mode: GPU scalar ─stream dependency→ library kernel

配置和标量地址必须一致。栈上 host 标量、临时 device 标量和异步执行的生命周期都要清楚。

CUDA Graph 与库

重复调用库工作流时，可以尝试 stream capture：

cudaStreamBeginCapture(stream, cudaStreamCaptureModeGlobal);
// library calls + custom kernels + async copies
cudaStreamEndCapture(stream, &graph);
cudaGraphInstantiate(&exec, graph, 0);

但并非所有 API/算法/版本都支持 capture；隐藏同步、首次 lazy initialization、动态 workspace 分配可能失败。先在 capture 前创建 handle、plan、descriptor、workspace 并预热，再捕获稳态路径。

统一的 Library Call 模板

1 validate dtype/shape/layout/device
2 choose/create handle and set caller stream
3 build or lookup descriptor/plan
4 query or lookup workspace size
5 get reusable workspace from pool
6 choose algorithm under precision/determinism/workspace constraints
7 launch asynchronously
8 check immediate API status
9 at test/debug sync boundary check asynchronous errors and result
10 cache plan; record telemetry

自写 Kernel 与库的组合点

高价值优化通常发生在边界：

• 用自写 kernel pack/unpack 到库最喜欢的布局；
• 用 cuBLASLt epilogue 融合 bias/activation，省掉额外 launch；
• 在 CUB reduction 前融合输入 transform；
• 用 cuFFTDx 在自定义 kernel 内融合 FFT 前后处理；
• 用 CUTLASS/CuTe 构建库未覆盖的 GEMM-like fusion；
• 用 NCCL group/stream 安排通信与计算重叠。

融合是否值得，以减少的 global bytes 和 launch overhead 对增加的寄存器、workspace、编译复杂度做权衡。

基准设计

报告 plan/heuristic/preprocess 是否计入，明确冷启动与稳态。只展示最有利的单个尺寸不构成可靠选型。

练习

1. 比较自写 reduction、CUB DeviceReduce 与 Thrust，在小/中/大 n 下记录同步和 workspace 成本。
2. 用 cuBLASLt 实现 bias epilogue，与 cuBLAS GEMM + bias kernel 比较 global 流量和时间。
3. 对多种 FFT 长度做 sweep，观察质数、大质因数和 2 的幂尺寸差异。
4. 为同一稀疏矩阵生成 CSR/BSR/SELL，比较转换成本、SpMV 和多次复用后的盈亏平衡点。
5. 用 NCCL AllReduce 测 1/2/4 卡，并在时间线上验证通信计算 overlap。

官方参考

• CUDA C++ Core Libraries
• CUB
• Thrust
• cuBLAS
• cuFFT
• cuSPARSE
• cuSOLVER
• cuRAND
• NCCL User Guide
• CUTLASS