1. CUDA科学计算加速导论

科学计算的核心矛盾，往往不是“算法是否正确”，而是“在可接受时间内能否算完”。以三维热传导显式差分、N 体引力模拟、稀疏矩阵迭代求解、蒙特卡洛积分等典型问题为例，计算量随网格规模、粒子数或样本数呈线性甚至超线性增长。一个 512³ 网格的显式热传导迭代，单步需要更新约 1.34 亿个格点；若迭代 1000 步，就是 1340 亿次格点更新。每次更新包含 6 次邻点读取、1 次中心点读取、若干次乘加与 1 次写回，实际访存与浮点运算交织，属于典型的计算密集型 + 访存密集型混合负载。CPU 单核即使以 3 GHz 主频、每周期 2 次浮点运算估算，理论峰值也只有约 6 GFLOP/s，且受内存带宽限制，实际有效吞吐往往更低。这正是 GPU 加速的切入点。

GPU 的收益并非无条件。它来自“大量同构线程 + 高带宽显存 + 高吞吐浮点单元”的组合：一块主流数据中心 GPU 可提供数千个 CUDA Core、数百 GB/s 到数 TB/s 的显存带宽，以及远高于 CPU 的 FP32/FP64 峰值。但收益边界由三个条件共同决定：并行度是否足够、数据是否可规则划分、数据传输与计算比是否合理。若一个问题只有几十个独立任务，或存在强数据依赖、频繁分支、随机访存，GPU 的利用率会迅速下降。更关键的是 PCIe/NVLink 传输开销：若每次计算只做很少的运算却要搬运大量数据，加速比可能被传输时间吞掉。因此，判断是否值得用 CUDA，应先估算“计算时间 / 传输时间”的比值，再决定是否把数据常驻显存、是否融合多个 kernel。

CUDA 在数值模拟与数据分析中的定位，可以概括为“把规则、密集、可并行的内层循环搬到 GPU，把控制流、I/O 和复杂逻辑留在 CPU”。典型模式是：CPU 负责网格生成、边界条件、参数扫描与结果落盘；GPU 负责 stencil 更新、稀疏矩阵向量乘、FFT、归约、直方图、距离矩阵等。以稀疏矩阵向量乘 y = A·x 为例，若 A 为 100 万阶、每行平均 20 个非零元，则单次 SpMV 约 2000 万次乘加，但需要读取约 2000 万个非零元及其列索引，访存远大于计算。此时优化重点不是提高浮点峰值，而是合并访存、使用共享内存缓存 x 的分段、减少索引开销。这说明 CUDA 优化必须回到“算术强度”这一基本量：每字节访存对应的浮点运算次数越低，越容易受带宽限制。

课程工具链围绕“可编译、可测量、可复现”组织。基础环境包括 NVIDIA 驱动、CUDA Toolkit（nvcc、cuda-gdb、Nsight Compute、Nsight Systems）、以及一个可用的 GPU（建议计算能力 7.0 及以上，以支持较完整的 Tensor Core 与异步拷贝特性）。编译侧使用 nvcc 将 .cu 文件分离为主机代码与设备代码，主机代码交给 gcc/clang，设备代码生成 PTX 与 SASS。性能度量基线建议固定为四类指标：加速比 S = T_cpu / T_gpu、有效显存带宽、浮点吞吐 GFLOP/s、传输时间占比。其中加速比必须说明 CPU 基线版本（单线程还是多线程、是否开启 -O3 -march=native），否则数字没有意义。

下面给出一个最小可运行的热传导 stencil 示例，用于建立第一条性能基线。示例参数为 512×512 网格、迭代 200 步、单精度浮点，边界固定为 0，初始中心区域为 100。CPU 版本用三重循环实现，GPU 版本用二维线程块映射。代码仅作教学示例，实际数值会随硬件变化。

// heat_cpu.c —— CPU 基线（单线程）
for (int t = 0; t < 200; ++t) {
  for (int i = 1; i < N-1; ++i)
    for (int j = 1; j < N-1; ++j)
      u_new[i][j] = 0.25f * (u[i-1][j] + u[i+1][j]
                           + u[i][j-1] + u[i][j+1]);
  swap(u, u_new);
}

// heat_gpu.cu —— GPU 版本核心 kernel
__global__ void heat_kernel(const float* u, float* un, int n) {
  int i = blockIdx.y * blockDim.y + threadIdx.y;
  int j = blockIdx.x * blockDim.x + threadIdx.x;
  if (i > 0 && i < n-1 && j > 0 && j < n-1) {
    int idx = i * n + j;
    un[idx] = 0.25f * (u[idx-n] + u[idx+n]
                     + u[idx-1] + u[idx+1]);
  }
}
// 启动配置：dim3 block(16,16); dim3 grid((n+15)/16,(n+15)/16);

运行后应记录：CPU 单线程耗时、GPU kernel 耗时、H2D/D2H 传输耗时、总耗时。若 GPU 总耗时没有明显低于 CPU，先检查传输占比是否过高；若 kernel 本身很快但总时间仍长，说明数据搬运是瓶颈，应考虑把多次迭代合并到一次 kernel 内、或使用 pinned memory 与多流重叠。若 kernel 耗时远高于预期，则用 Nsight Compute 查看显存吞吐、L2 命中率、占用率与分支发散情况。判断是否达成目标的依据不是“跑通了”，而是：加速比是否稳定、有效带宽是否接近硬件峰值的合理比例、传输占比是否被压到可接受范围。