> 基于 MCC2026 OSTIA SST 气候态计算项目（V0 → V10.24，13,214× 加速比）实战经验 > 硬件环境：Hygon C86 7185 + 4×DCU Z100L (gfx906) × 2节点 --- 上篇聊了 DCU 硬件长什么样、软件栈怎么搭。这一篇聊怎么写代码让它跑起来——从第一个 HIP 程序，到 kernel 设计、内存布局、rocprof 性能诊断，再到我们实际踩过的三个差点放弃的坑。 ### 题目背景 计算 **1991-2020 年全球海洋表面温度（SST）的气候态均值和 P90 分位数**。输入数据是 3060 个 NetCDF 文件（102 天 × 30 年），每个文件包含一张 1440×721 的 SST 格点图。输出是 184 个 NC 文件（92 天的均值 + 92 天的 P90）。 计算的核心逻辑：对每个海洋格点（1,038,240 个），取该日 ±5 天共 11 天 × 30 年 = 330 个样本，统计均值和 P90。像素之间完全独立——典型的"尴尬并行"（embarrassingly parallel，任务间几乎不需要通信或同步，各算各的，简单到"令人尴尬"的程度），很适合 GPU。 ### 代码整体结构 ``` main() ├── 读取坐标/元数据 ├── 预加载 3060 个 NC 文件 → 6.4 GB 内存池 pool[DOY][year][pixel] ├── 启动 N 个 gpu_worker 线程（每 GPU 一个） │ ├── H2D: hipMemcpy(pool子集 → GPU显存) │ ├── Kernel: compute_stats_kernel<<>> → 统计均值和P90 │ └── D2H: hipMemcpy(GPU结果 → CPU内存) └── writer线程串行写 NC 输出（HDF5 非线程安全） ``` 核心计算全在 `compute_stats_kernel` 这一个 kernel 里完成。下面从最小的 HIP 程序讲起，逐步展开这个 kernel 的设计细节。每个结论都来自项目实测。 --- ## 一、HIP Kernel 编程：从零到最优 ### 1.1 一个最小的 HIP 程序 在讲复杂的 kernel 之前，先看一个完整的 HIP 程序长什么样： ```cpp #include #include // GPU kernel：每个线程计算 output[i] = input[i] * 2.0 __global__ void double_kernel(const float *input, float *output, int N) { int i = blockIdx.x * blockDim.x + threadIdx.x; // 全局线程ID if (i < N) output[i] = input[i] * 2.0f; } int main() { int N = 1024; size_t bytes = N * sizeof(float); // 1. 分配 host 内存 float *h_in = new float[N]; float *h_out = new float[N]; for (int i = 0; i < N; i++) h_in[i] = (float)i; // 2. 分配 device 内存 float *d_in, *d_out; hipMalloc(&d_in, bytes); hipMalloc(&d_out, bytes); // 3. Host → Device hipMemcpy(d_in, h_in, bytes, hipMemcpyHostToDevice); // 4. 启动 kernel int blockSize = 256; int gridSize = (N + blockSize - 1) / blockSize; double_kernel<<>>(d_in, d_out, N); // 5. Device → Host hipMemcpy(h_out, d_out, bytes, hipMemcpyDeviceToHost); // 6. 清理 hipFree(d_in); hipFree(d_out); delete[] h_in; delete[] h_out; return 0; } ``` 这个简单的例子包含了所有 HIP 程序的五个步骤： 1. **分配**：`hipMalloc` 在 GPU 显存上分配空间 2. **H2D**：`hipMemcpy(..., hipMemcpyHostToDevice)` 将数据从主机内存拷到 GPU 3. **Kernel Launch**：`kernel<<>>(...)` 在 GPU 上执行并行计算 4. **D2H**：`hipMemcpy(..., hipMemcpyDeviceToHost)` 将结果从 GPU 拷回主机 5. **清理**：`hipFree` 释放 GPU 显存 ### 1.2 我们的实际 Kernel：compute_stats_kernel 实际的 kernel 比上面的例子复杂得多。它负责计算每个海洋格点的气候态均值和 90% 分位数： ```cpp __global__ void compute_stats_kernel( const short *__restrict__ pool, // 输入：pool[DOY][year][pixel] (short量化) const int *__restrict__ doy_to_pool, // 映射表：日期→pool索引 int doy_count, // 本kernel处理的DOY数 int doy_offset, // DOY偏移 float scale_factor, // 解码参数：scale float add_offset, // 解码参数：offset short fill_value, // 缺失值标记 double *__restrict__ out_mean, // 输出：气候态均值 [DOY][pixel] double *__restrict__ out_p90) // 输出：P90分位数 [DOY][pixel] { int pixel = blockIdx.x * blockDim.x + threadIdx.x; // 我的像素 int d_idx = blockIdx.y; // 我的DOY if (pixel >= GRID_SIZE || d_idx >= doy_count) return; // 为这个像素构建滑动窗口（±5天，共11个DOY） int tdoy = DOY_START + doy_offset + d_idx; // ... 11天 × 30年 = 330个样本的统计 ... } ``` **Grid 维度设计**： ```cpp dim3 block(256); // 每block 256线程 = 4个warp(wavefront) dim3 grid(GRID_SIZE / 256, // x方向：覆盖全部1,038,240像素 doy_count); // y方向：覆盖该GPU分配的DOY数 ``` 选择 256 线程/block 是一个经验性平衡： - **太小**（如 64）：调度开销大，CU 利用率低 - **太大**（如 1024）：寄存器压力过大，occupancy 进一步降低 - **256**：在 gfx906 上通常是最优的起点值 我们使用了 **2D grid**（不是常见的 1D grid），这允许一次 kernel launch 处理多个 DOY——避免 92 次独立 kernel launch 的累积开销（每次 ≈10 μs × 92 = ≈1 ms）。 **内核核心逻辑（伪代码）**——每个线程处理一个像素，统计 330 个样本的均值和 P90： ```cpp // 每个线程处理一个像素，遍历 11天×30年=330个样本 for (int w = 0; w < 11; w++) { // 滑动窗口 ±5天 const short *blk = pool[doy_index[w] * 30 * GRID_SIZE]; for (int yr = 0; yr < 30; yr++) { short raw = blk[yr * GRID_SIZE + pixel]; // 合并访问：相邻线程读相邻地址 if (raw == FILL_VALUE) continue; float v = raw * SCALE + OFFSET; nv++; sum += v; // Top-34 降序有序数组（插入排序维护），用于计算 P90 if (v > sorted[33]) { for (int i = 33; i >= 0; i--) { if (v > sorted[i]) { sorted[i+1] = sorted[i]; sorted[i] = v; } } } } } // P90: Hyndman & Fan Type=5, r = nv*0.9 + 0.5 double r = nv * 0.9 + 0.5; int k = (int)floor(r); double f = r - k; p90 = sorted[needed-1] + f * (sorted[needed-2] - sorted[needed-1]); ``` **为什么用有序数组而不是最小堆？** 事实上我们在初始版本使用了最小堆维护的top-k算法来计算P90分位，但在 MemUnitBusy 99.7% 的内存瓶颈场景下，ALU 有大量空闲周期。O(34) 次比较可以被内存等待完全掩盖（"免费的计算"），而数组的 warp 执行一致性更好，因为其所有线程执行相同的比较和移动操作，无分支发散。 ### 1.3 内存布局的学问：为什么 pool[DOY][year][pixel]？ 这是我个人认为整个项目中最重要、也最容易忽视的架构决策。 数据在内存中的排列方式直接决定了 GPU 的内存带宽利用率。在我们的 kernel 中，每个线程固定处理一个 pixel，遍历 330 个（DOY, year）组合。所以访问模式是： ``` 线程0：pool[doy0][yr0][p0], pool[doy0][yr1][p0], ..., pool[doy10][yr29][p0] 线程1：pool[doy0][yr0][p1], pool[doy0][yr1][p1], ..., pool[doy10][yr29][p1] ``` **关键**：当所有线程同时执行同一条指令（`short raw = blk[yr * GRID_SIZE + pixel]`）时，线程 0~63 访问的是： - 地址 A, A+2, A+4, ..., A+126 这 128 bytes 落在 2 个 cache line 内 → 只需要 2 次内存事务就完成整个 warp 的数据获取。这就是**合并访问（Coalesced Access）**。 **如果换成 pool[pixel][DOY][year]**： ``` 线程0：pool[p0][doy0][yr0], pool[p0][doy0][yr1], ... 线程1：pool[p1][doy0][yr0], pool[p1][doy0][yr1], ... ``` 线程 0 和线程 1 的地址相距 102×30×2 = 6,120 bytes → 不在同一 cache line → 每个线程触发独立的内存事务。64 线程需要 64 次事务，带宽利用率下降 **32 倍**。 **推而广之的规则**：GPU 数据结构的**最内层维度**（内存地址变化最快的那个维度）最好对应 warp 内线程 ID 的变化方向。在 HPC 术语中，这是 "Structure of Arrays (SoA)" 优于 "Array of Structures (AoS)" 的 GPU 版本。 ### 1.4 数据量化：float → short SST 原始数据是 float(4 bytes)，但为了提速，我们把数据量化为 short(2 bytes)，在 GPU 上解码还原： ```cpp // 存储侧（CPU预加载）：float → short const double inv_scale = 1.0 / 0.01; // scale_factor = 0.01 from NC metadata for (int i = 0; i < GRID_SIZE; i++) { double v = raw_double[i]; double sv = v * inv_scale; dest[i] = (sv == sv) ? (short)lround(sv) : FILL_VALUE; // NaN → 填充值 } // 计算侧（GPU kernel）：short → float 解码 float v = (float)((double)raw * 0.01 + 20.0); // scale_factor, add_offset ``` 效果：pool 从 12.7 GB 降到 6.4 GB，H2D 传输量减半。精度代价：气候态均值差 0.0002°C，完全可忽略。 --- ## 二、rocprof：GPU 性能分析的显微镜 ### 2.1 rocprof 是什么？ rocprof 是 ROCm 平台提供的**硬件性能计数器 Profiler**，可以直接读取 GPU 芯片上的硬件计数器——内存控制器忙碌时间、ALU 利用率、缓存命中率等。类比 CPU 世界的 `perf stat`。 基本用法： ```bash rocprof --stats -o result.csv ./your_hip_app ``` 在我们的项目中，对 kernel 做性能分析： ```bash rocprof --stats -o profile.csv ./get_climatology_hip \ --nc-path /public/home/achwjznh4b/Newdata/ \ --save-path /tmp/ --gpus-per-node 1 ``` `--stats` 输出每个 kernel 的硬件计数器汇总，`-o` 指定输出 CSV 文件。 ### 2.2 关键硬件计数器详解 我们从该项目 kernel 的 rocprof 报告中选取最重要的几个计数器： | 计数器 | 本项目值 | 含义 | 健康范围 | | ---------------- | :--------: | ------------------------------- | :--------------: | | **MemUnitBusy** | **99.71%** | HBM 内存控制器忙碌的时间比例 | <60%=计算瓶颈 | | VALUUtilization | 58.81% | 向量 ALU 执行指令的时间比例 | >80%=计算瓶颈 | | WriteUnitStalled | 23.49% | 写操作等待 HBM 完成的时间比例 | <10%=健康 | | L2CacheHit | 67.21% | L2 缓存命中的比例 | >80%=数据友好 | | FETCH_SIZE | 161 MB | 每 kernel 从 HBM 读取的总数据量 | 对比理论带宽 | | SQ_WAVES | 1,485,250 | 总共调度的 wavefront 数量 | 对比 grid 大小 | | arch_vgpr | 160 | 每线程使用的 VGPR 数量 | <64=高 occupancy | ### 2.3 三句话诊断法 拿到 rocprof 输出后，用这个三步流程判断 kernel 的瓶颈类型： **第一步：看 MemUnitBusy** 如果 >80%：**内存瓶颈**。GPU 的 HBM 控制器跟不上计算速度，ALU 在等待数据。这种情况下，优化 kernel 代码（减少除法、合并运算、循环展开）不会带来任何性能提升——因为指令执行时间被内存等待完全掩盖了。 如果 <60%：**计算瓶颈**。GPU 的算术单元是瓶颈。这时才值得优化指令流水、减少寄存器使用、提升 occupancy。 本项目：MemUnitBusy = 99.71% → **毫无疑问的内存瓶颈**。 **第二步：看 VALUUtilization** VALUUtilization 是向量 ALU 真正在执行计算的时间比例。如果 VALUUtilization 远小于 MemUnitBusy，说明 ALU 大部分时间在空转等数据——**证实了内存瓶颈诊断**。 本项目：VALUUtilization = 58.81%，而 MemUnitBusy = 99.71%。差值 41% 就是 ALU 在等待 HBM 数据的时间。**指令是免费的——你在内存等待期间执行再多指令都不会让程序变慢。** **第三步：看 L2CacheHit** 如果 >80%：数据访问模式很好地适应了 L2 缓存容量。 如果 <60%：数据量超过了 L2 容量（4 MB），大量访问直接穿透到 HBM。 本项目：L2CacheHit = 67.21%——处于中等水平。33% 的访问直接到 HBM 是因为 `pool[DOY][year][pixel]` 布局中，跨年的 stride 为 GRID_SIZE × 2 bytes = 2 MB，远超 L2（4 MB）的容量，L2 无法缓存跨年访问。 ### 2.4 一个常见误解：Occupancy 越高越好 GPU 厂商的宣传材料常常强调 "occupancy"（每 SIMD 驻留的 wavefront 数）的重要性。理论依据是：更多 wavefront 可以更好地隐藏内存延迟——当一个 wavefront 在等待内存时，SIMD 可以切换到另一个。 但这有一个前提：**内存瓶颈不是绝对瓶颈**。 当 MemUnitBusy 已经接近 100% 时，意味着 HBM 控制器已经从早忙到晚，没有多余的带宽给任何新 wavefront。此时提升 occupancy 只是在 HBM 控制器前加更长的排队——不会让数据传输变快。 事实上，**较低的 occupancy 在特定情况下反而更好**： - 更少的 wavefront → 更少的 L2 竞争 → 更高的单 wavefront L2 命中率 - 更少的 wavefront → 更少的 HBM bank 冲突 → 更高的有效带宽 我们 25% 的 occupancy 在 MemUnitBusy 99.7% 的背景下，不仅不是问题，甚至可能是最优的。 --- ## 三、H2D、Kernel、D2H 三级管线 ### 3.1 基本概念 GPU 计算的完整流程可以分为三个阶段： ``` Host 内存 ----H2D----> GPU 显存 ----Kernel----> GPU 显存 ----D2H----> Host 内存 (输入数据) (数据传输) (准备就绪) (并行计算) (计算结果) (回传) (拿到结果) ``` - **H2D (Host to Device)**：把输入数据从 CPU 内存通过 PCIe 总线传到 GPU 显存 - **Kernel**：GPU 执行并行计算 - **D2H (Device to Host)**：把计算结果从 GPU 显存传回 CPU 内存 这三个阶段的耗时分布（本项目，单 GPU）： - H2D：~100ms（1.35 GB short 数据 / PCIe 3.0 ×16 ≈ 14 GB/s） - Kernel：~900ms（46 DOY × 1M pixel × 330 sample） - D2H：~100ms（46 DOY × 1M pixel × 8 bytes double ≈ 370 MB） ### 3.2 HIP Stream 异步流水线 HIP Stream 允许我们**并行执行**这三个阶段中的不同部分： ``` 时间轴 → Stream 0: [H2D_chunk0][Kernel_chunk0][D2H_chunk0] Stream 1: [H2D_chunk1][Kernel_chunk1][D2H_chunk1] ↑ H2D_chunk1 和 Kernel_chunk0 重叠！ ``` 我们当前的实现使用 2 个 stream，将 DOY 分为前半和后半： ```cpp hipStream_t s0, s1; hipStreamCreate(&s0); hipStreamCreate(&s1); // Stream 0 处理前半 DOY compute_stats_kernel<<>>(..., d_mean, d_p90); hipMemcpyAsync(h_mean, d_mean, half0_bytes, hipMemcpyDeviceToHost, s0); // Stream 1 处理后半 DOY compute_stats_kernel<<>>(..., d_mean_half1, d_p90_half1); hipMemcpyAsync(h_mean_half1, d_mean_half1, half1_bytes, hipMemcpyDeviceToHost, s1); hipStreamSynchronize(s0); hipStreamSynchronize(s1); ``` **不过，当前 H2D 是同步的**（没有用 stream 包装），排在所有 kernel 之前。这是因为 H2D 和 kernel 用的是不同的硬件资源（PCIe 控制器 vs HBM 控制器），理论上可以通过三缓存完全重叠——这在我们下一阶段的优化计划中。 ### 3.3 hipHostRegister：让 GPU 直接看 Host 内存 默认情况下，`hipMemcpy` 的 H2D 传输需要两步： 1. 驱动将数据从 host 页拷贝到内部 **staging buffer**（一段提前分配的 GPU 可访问内存） 2. GPU DMA 引擎从 staging buffer 传输到显存 `hipHostRegister`（或 `hipHostMalloc`）将 host 内存页**锁定**（pin）在物理内存中，并注册到 IOMMU（I/O 内存管理单元），使 GPU DMA 引擎可以直接访问这些物理页面。这样就省掉了第一步的 staging buffer 拷贝。 **收益**：H2D 带宽约翻倍（在我们的场景中，~6 GB/s → ~12 GB/s）。 **代价**： - 注册操作本身有开销（我们实测 ~0.5-1s） - 被锁定的页不能换出（减少系统可用虚拟内存） - **必须在页面已被访问后注册**——如果注册时页面还没有物理内存（未 fault-in），后续每个 page fault 都触发 IOMMU 重新映射，开销爆炸 我们的实现时序： ```cpp // 1. mmap 匿名映射 pool 空间 short *pool = (short*)mmap(nullptr, pool_bytes, PROT_READ | PROT_WRITE, MAP_PRIVATE | MAP_ANONYMOUS, -1, 0); // 2. 预加载线程写入数据（触发 page fault，分配物理页） for (auto &t : preload_threads) t.join(); // 3. 页面全部 fault-in 后，注册到 IOMMU hipHostRegister(pool, pool_bytes, hipHostRegisterDefault); // 4. GPU 通过 DMA 直接读取 pool（无需 staging buffer 拷贝） // ... GPU 计算 ... // 5. 释放 hipHostUnregister(pool); munmap(pool, pool_bytes); ``` **关键**：步骤 2 必须在步骤 3 之前。如果注册时页面还没有物理内存（未 fault-in），后续每个 page fault 都触发 IOMMU 重新映射，开销爆炸。 --- ## 四、Debug 实战：三个差点让我们放弃的 Bug ### 4.1 P90 偏差 0.0126°C（Hyndman & Fan 公式） 这是整个项目中最让人抓狂的 bug——排查了四天。 **症状**：气候态均值 RMSE = 0.0000°C（完全匹配），P90 RMSE = 0.0126°C（系统性偏差）。这个量级非常诡异——远大于浮点舍入误差（~1e-12），又远小于典型的编程错误（1°C+）。 **排除过程**： 第一步，怀疑维度转置。检查发现 NC 输出的 `grid_dims` 确实是 `{lon_dim, lat_dim}` 写反了——C 行优先与 NetCDF Fortran 列优先不匹配。修复了这个 bug，但 **P90 RMSE 纹丝不动**。 第二步，怀疑堆选择算法有遗漏。把 GPU kernel 的 top-34 堆选择换成全排序（std::sort 330 个元素），期待 P90 变化。完全不变化——说明堆选择没有问题，选择的 top-34 值是正确的。 第三步，怀疑数据解码精度。把 short 解码从 float 提升到 double，仍然无变化。 第四步，**终于想到检查公式本身**。MATLAB `prctile` 函数有多个 variant： - MATLAB R2018a 默认使用 **Hyndman & Fan Type=5**：`r = N × 0.9 + 0.5` - 统计教材常见 **Type=7**：`r = (N - 1) × 0.9 + 1.0` 对于 N = 330： - Type=5：`r = 297.5` → k = 297, f = 0.5 → P90 = v₃₀₀ + 0.5 × (v₂₉₉ - v₃₀₀)（注意降序排列，v₃₀₀ 是第 300 大的值） - Type=7：`r = 297.1` → k = 297, f = 0.1 → P90 = v₃₀₀ + 0.1 × (v₂₉₉ - v₃₀₀) 两者选择相同的两个顺序统计量，但插值权重相差 0.4 → 恰好产生 0.0126°C 的全局均匀偏差。 **教训**：跨语言/跨库数值验证时，**第一步就应该核对数学公式**，而不是假设同名函数算出来的结果一样。 ### 4.2 HDF5 1.8 线程不安全（随机崩溃） NetCDF 4.3.0 底层依赖 HDF5 1.8.x，它有一个全局 `H5FL` (HDF5 Free List) 内存池。当多个线程同时调用 `nc_open()` / `nc_close()` 时，底层 HDF5 会并发操作这个全局池→ double-free 或野指针→随机 SIGSEGV。 **症状**：有时候运行正常，有时候在 `H5FL_fac_free` 崩溃，完全随机。 **两种绕过方案**： 方案 A（Fork 路径）： ```cpp for (int w = 0; w < 28; w++) { pid_t pid = fork(); if (pid == 0) { // 子进程：独立地址空间 → 独立的 H5FL 池 → 无竞争 read_nc_files(...); _exit(0); } } ``` fork 后父子进程有独立的地址空间（copy-on-write），天然隔离了 HDF5 全局状态。 方案 B（pread 路径）： ```cpp // 完全绕过 NetCDF/HDF5 库，直接用系统调用 int fd = openat(dir_fd, filename, O_RDONLY); pread(fd, buffer, GRID_SIZE * sizeof(double), data_offset); close(fd); ``` pread 是 POSIX 标准系统调用，100% 线程安全。不需要 fork，不需要锁。 方案 B 在性能上远超方案 A（sys 降 35%，page-fault 降 98%），但它只能用于原始数据块在文件中连续存储的场景——需要提前用 HDF5 API 的 `H5Dget_offset()` 确认数据偏移量。 ### 4.3 Fortran 列优先 vs C 行优先（16°C RMSE） MATLAB 使用 `mexcdf` 工具箱写 NC 文件，数据按 Fortran 列优先（column-major）顺序排列。C++ 的 `nc_put_var_double` 按 C 行优先（row-major）解释。 对于同一个 `[Lon=1440, Lat=721]` 矩阵： ``` C 行优先： element(lon, lat) = lon × 721 + lat → lat 变化最快 Fortran 列优先：element(lon, lat) = lon + lat × 1440 → lon 变化最快 ``` 如果搞错，经纬度互换 → **16°C 的 RMSE**。 **正确做法**：`grid_dims[2] = {lat_dim, lon_dim}`（NetCDF 的 dims 从左到右对应 C 的最慢变化到最快变化维度），并将数据按 `C_lon_lat[lon * LAT_NUM + lat]` 写入。 --- *下一篇：下篇《系统级优化工程》—— 13,214× 加速比的完整路线图、I/O 优化、Lustre 极限分析、何时停止优化。*