海光 DCU 深度技术报告 · 中篇：HIP 编程实战与性能分析

基于 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<<<grid, 256>>>  → 统计均值和P90
  │     └── D2H: hipMemcpy(GPU结果 → CPU内存)
  └── writer线程串行写 NC 输出（HDF5 非线程安全）

核心计算全在 compute_stats_kernel 这一个 kernel 里完成。下面从最小的 HIP 程序讲起，逐步展开这个 kernel 的设计细节。每个结论都来自项目实测。

---

一、HIP Kernel 编程：从零到最优

1.1 一个最小的 HIP 程序

在讲复杂的 kernel 之前，先看一个完整的 HIP 程序长什么样：

#include <hip/hip_runtime.h>
#include <cstdio>

// 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<<<gridSize, blockSize>>>(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<<<grid, block>>>(...) 在 GPU 上执行并行计算
4. D2H：hipMemcpy(..., hipMemcpyDeviceToHost) 将结果从 GPU 拷回主机
5. 清理：hipFree 释放 GPU 显存

1.2 我们的实际 Kernel：compute_stats_kernel

实际的 kernel 比上面的例子复杂得多。它负责计算每个海洋格点的气候态均值和 90% 分位数：

__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 维度设计：

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：

// 每个线程处理一个像素，遍历 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 上解码还原：

// 存储侧（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。

基本用法：

rocprof --stats -o result.csv ./your_hip_app

在我们的项目中，对 kernel 做性能分析：

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 分为前半和后半：

hipStream_t s0, s1;
hipStreamCreate(&s0);
hipStreamCreate(&s1);

// Stream 0 处理前半 DOY
compute_stats_kernel<<<grid0, block, 0, s0>>>(..., d_mean, d_p90);
hipMemcpyAsync(h_mean, d_mean, half0_bytes, hipMemcpyDeviceToHost, s0);

// Stream 1 处理后半 DOY
compute_stats_kernel<<<grid1, block, 0, s1>>>(..., 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 重新映射，开销爆炸

我们的实现时序：

// 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 路径）：

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 路径）：

// 完全绕过 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 极限分析、何时停止优化。