海光 DCU 编程入门——以气候态海温百分位计算为教学案例
·
摘要
本文面向海光 DCU 的编程初学者，全程以一套真实的气候态海温第 90 百分位计算程序为主线，从零开始介绍 DCU 上的 HIP 编程。每引入一个概念（线程、显存、stream、核函数），都立即用 SST 程序中的对应代码来演示。读者完成本文后，应能理解 SST 程序的完整 GPU 计算流程，并具备将类似气候数据计算任务移植到 DCU 的能力。
关键词：海光 DCU；HIP 入门；SST 气候态；百分位计算；实例教学
第 1 章 计算问题与 DCU 概览
1.1 SST 百分位计算问题
本文全程围绕一个具体的科学计算问题展开：
给定 1991 年至 2020 年（30 年）全球逐日海表温度（SST）数据，网格分辨率 0.25 度（经度 1440 格点 × 纬度 721 格点），对每个海洋格点、每个日历日（DOY 152~243，共 92 天），计算其第 90 百分位 SST。
这是海洋热浪监测的基础算法。对每个格点，需要从 30 年 × 11 天滑动窗口 = 330 个有效温度值中找出第 90 百分位。
数据规模：
-
每文件：1440 × 721 × 8 字节 ≈ 8.3 MB
-
总文件数：30 年 × 约 100 天 × 8 个窗口天数 ≈ 24,000 个文件
-
总数据量：约 200 GB
计算特征：
-
每个格点完全独立（完美并行）
-
陆地格点为 NaN（约占 30%，可跳过）
-
需要大量浮点运算（累加、比较、排序）
1.2 为什么用 DCU
DCU 擅长的正是这种"对大量独立数据重复执行相同运算"的任务：
CPU 的方式（串行，慢）： for 每个格点(104 万个): for 每一年(30 年): for 每一天(11 天): 累加、累乘、判断 ... DCU 的方式（并行，快）： [格点 0 的线程] ← 同时运行 [格点 1 的线程] ← 同时运行 [格点 2 的线程] ← 同时运行 ... ← 4096 个流处理器同时工作
1.3 硬件规格
| 参数 | 值 | 在本项目中的意义 |
|---|---|---|
| 计算单元 (CU) | 64 | 至少需要 64 × 256 = 16384 线程填满 |
| 流处理器 | 4096 | 每个 CU 64 个 |
| 显存 | 16 GB HBM2/卡 | 可存放约 2 年的全球 SST 数据 |
| Wavefront | 64 线程 | 线程块大小应设 64、128 或 256 |
| 每 CU 寄存器 | 64 KB | 每线程最多约 64 个寄存器 |
| 节点配置 | 4 块 DCU | 4 块分工处理全球格点 |
1.4 SST 程序在 DCU 上的计算模式
CPU (主控) ┌──────────┐ │ MPI rank │ ← 8 个 MPI 进程 │ 负责 12 天 │ ← 每个进程处理约 12 个 DOY └─────┬────┘ │ ┌───────────────┼───────────────┐ │ │ │ ┌───┴───┐ ┌───┴───┐ ┌───┴───┐ │ DCU 0 │ │ DCU 1 │ ... │ DCU 3 │ │ 259560 │ │ 259560 │ │ 259560 │ │ 个格点 │ │ 个格点 │ │ 个格点 │ └───────┘ └───────┘ └───────┘
每个 MPI rank 处理约 12 个 DOY（约 22 个相关日历日），4 块 DCU 平分全球 104 万个格点。
第 2 章 HIP 程序结构和 SST 主循环
2.1 SST 程序的 HIP 骨架
打开 SST 程序的 main() 函数，可以看到完整的 HIP 四步法：
int main(int argc, char** argv)
{
// ── 初始化 MPI ──
MPI_Init(&argc, &argv);
int rk, nr; // rank 编号和总数
MPI_Comm_rank(MPI_COMM_WORLD, &rk);
MPI_Comm_size(MPI_COMM_WORLD, &nr);
// ── 第 1 步：分配 ──
// CPU 侧：host_ring 用于存放读进来的数据
double* host_ring[2];
hipHostMalloc(&host_ring[0], batch_bytes); // pinned!
hipHostMalloc(&host_ring[1], batch_bytes);
// DCU 侧：d_ring 存放传输到显存的数据
double* d_ring[4]; // 每块 DCU 一个
for(int g = 0; g < 4; g++){
hipSetDevice(g);
hipMalloc(&d_ring[g], d_ring_bytes);
}
// ── 第 2 步：传输 ──
for(int g = 0; g < 4; g++){
hipSetDevice(g);
hipMemcpy2DAsync(d_ring[g], ..., host_ring[cur], ..., st[g]);
}
// ── 第 3 步：计算 ──
for(int g = 0; g < 4; g++){
hipSetDevice(g);
ka_moment<<<blk, thr, 0, st[g]>>>(d_ring[g], ...);
}
// ── 第 4 步：传回 ──
for(int g = 0; g < 4; g++){
hipSetDevice(g);
hipMemcpy(h_result + offset, d_result[g], ..., hipMemcpyDeviceToHost);
}
MPI_Finalize();
return 0;
}
2.2 SST 数据的存储结构
SST 数据在内存中的排列方式是理解整个程序的关键：
host_ring[cur] 的内容（按行排列）： ┌──────────────────────────────────────────┐ │ 第 0 年第 0 天： [格点0][格点1]...[格点N] │ ← gs 个格点 │ 第 0 年第 1 天： [格点0][格点1]...[格点N] │ ← gs 个格点 │ ... │ │ 第 0 年第 21 天： [格点0][格点1]...[格点N] │ │ 第 1 年第 0 天： [格点0][格点1]...[格点N] │ │ ... │ │ 第 4 年第 21 天： [格点0][格点1]...[格点N] │ └──────────────────────────────────────────┘ 共 BATCH×ns 行（BATCH=5 年，ns≈22 天） 每行 gs=1,038,240 个格点，每个格点 8 字节
2.3 数据划分到 4 块 DCU
size_t gs = 1440 * 721; // 1,038,240 个格点
size_t ck = (gs + 3) / 4; // 259,560 格点/DCU
for(int g = 0; g < 4; g++){
size_t o = g * ck; // 本 DCU 的起始格点
size_t sz = (g == 3) ? (gs - o) : ck; // 本 DCU 的格点数
// hipMemcpy2DAsync 只拷贝本 DCU 负责的列
HIP_CHECK(hipMemcpy2DAsync(
d_ring[g], sz * 8, // 目标（DCU 端）
host_ring[cur] + o, gs * 8, // 源（CPU 端，偏移 o 列）
sz * 8, // 每行拷贝宽度
(size_t)by * ns, // 行数
hipMemcpyHostToDevice,
st[g]));
}
2.4 四种内存传输方向
// 1. CPU → DCU：输入数据 hipMemcpy(d_ring, host_ring, size, hipMemcpyHostToDevice); // 2. DCU → CPU：计算结果 hipMemcpy(h_result, d_result, size, hipMemcpyDeviceToHost); // 3. DCU → DCU：设备内部拷贝（本项目未使用） hipMemcpy(d_dst, d_src, size, hipMemcpyDeviceToDevice); // 4. 不传输：零拷贝（只在本项目 AVX2 预筛中使用） // CPU 直接读取 host_ring 做预分类，不涉及 DCU
第 3 章 SST 程序中的核函数
3.1 核函数的输入输出关系
SST 程序的核心计算在 ka_moment 核函数中完成。它的逻辑是：
输入: d_ring[g] — 一批 SST 数据（5 年 × 22 天 × 259560 格点） 输出: d_mom[g] — 4 阶矩（m1,m2,m3,m4）每个 DOY × 每个格点 d_cnt[g] — 有效值计数每个 DOY × 每个格点
3.2 ka_moment 核函数
// 计算每个格点的 4 阶统计矩
// 每个线程负责一个格点，遍历所有年份和天数
__global__ void ka_moment(
const double* __restrict__ r, // 输入：逐日 SST
float* mom, // 输出：矩 [doy][grid][4]
int* cnt, // 输出：计数 [doy][grid]
size_t sz, // 本 DCU 格点数
int ns, // 总天数（ns≈22）
int ndo, // 要算的 DOY 数（ndo≈12）
int nd, // 窗口大小（nd=11 天）
int ny, // 本批年数（ny=5 年）
size_t ls) // 行跨度（=sz）
{
// 当前线程的全局 ID → 对应一个格点
size_t ti = blockIdx.x * blockDim.x + threadIdx.x;
if(ti >= sz) return; // 边界检查
// 遍历每个 DOY
for(int di = 0; di < ndo; di++){
// 从显存加载已累加的矩到寄存器
float* m = &mom[((size_t)di * sz + ti) * 4];
float m1 = m[0], m2 = m[1], m3 = m[2], m4 = m[3];
int vc = cnt[di * sz + ti];
// 遍历所有年份和窗口天数
for(int y = 0; y < ny; y++)
for(int d = 0; d < nd; d++){
// 计算输入数据的位置
// 第 y 年第 d 天偏移 = (y*ns + di + d) 行 × ls 列 + ti
double v = r[((size_t)y * ns + di + d) * ls + ti];
if(!isnan(v)){ // 跳过陆地 NaN
float fv = (float)v;
vc++;
// 4 个 FMA：累加 4 阶矩
m1 += fv; // 一阶矩（和）
m2 += fv * fv; // 二阶矩（平方和）
m3 += fv * fv * fv; // 三阶矩（立方和）
m4 += fv * fv * fv * fv; // 四阶矩（四次方和）
}
}
// 写回显存
m[0] = m1; m[1] = m2; m[2] = m3; m[3] = m4;
cnt[di * sz + ti] = vc;
}
}
每个线程的工作量：
-
遍历 12 个 DOY × 5 年 × 11 天 = 660 个 SST 值
-
对每个值：1 次 NaN 检查 + 4 次 FMA
-
每格点总共：660 × 5 = 3,300 次浮点运算
3.3 线程块配置
size_t sz = 259560; // 每 DCU 的格点数 int thr = 64; // 每块 64 线程 = 1 wavefront int blk = (sz + thr - 1) / thr; // 需要 4056 个块 ka_moment<<<blk, thr, 0, st[g]>>>(d_ring[g], ...);
为什么设 64？ 因为 DCU 的最小调度单位是 wavefront（64 线程）。设 64 刚好一个 wavefront 完全填满，没有浪费。设 128 或 256 也没问题，但 64 在当前 kernel 中寄存器使用最优。
3.4 最终结果：从矩到 P90
计算完 4 阶矩后，用 Cornish-Fisher 展开求出 P90：
// 从 4 矩计算第 90 百分位 SST float n = (float)vc; // 有效值个数 float mu = m1 / n; // 均值 float s2 = m2/n - mu*mu; // 方差 float s = sqrtf(s2); // 标准差 float sk = (偏度公式); // 偏度 float ku = (峰度公式) - 3; // 超额峰度 // Cornish-Fisher 展开 float z = 1.2816f; // 标准正态 P90 分位数 float cf = z + (z*z-1)*sk/6 + (z*z*z-3*z)*ku/24 - (2*z*z*z-5*z)*sk*sk/36; float p90 = mu + s * cf; // 最终 P90 海温
第 4 章 SST 程序的多 DCU 协作
4.1 4 块 DCU 的分工
全球 SST 网格共 1,038,240 个格点。4 块 DCU 各负责约 259,560 个格点：
// DCU 0：格点 0000000~0259560 // DCU 1：格点 0259561~0519120 // DCU 2：格点 0519121~0778680 // DCU 3：格点 0778681~1038240
每个 batch 中，host_ring 存放了完整网格 × 若干年份。传输时每块 DCU 只取自己那部分：
for(int g = 0; g < 4; g++){
HIP_CHECK(hipSetDevice(g));
size_t o = g * ck; // 本 DCU 的起始列
size_t sz = gsz[g]; // 本 DCU 的列数
// 2D 异步拷贝：只取 o 到 o+sz 列
HIP_CHECK(hipMemcpy2DAsync(
d_ring[g], sz * 8,
host_ring[cur] + o, gs * 8,
sz * 8, height,
hipMemcpyHostToDevice, st[g]));
// 本 DCU 只算自己的格点
ka_moment<<<blk, thr, 0, st[g]>>>(d_ring[g], ...);
}
4.2 同步 4 块 DCU
每个 batch 结束时需要等待所有 DCU 完成：
// 顺序同步
for(int g = 0; g < 4; g++){
HIP_CHECK(hipSetDevice(g));
HIP_CHECK(hipStreamSynchronize(st[g]));
}
// 或用 std::async 并行同步（更快）
std::future<void> sf[4];
for(int g = 0; g < 4; g++)
sf[g] = std::async([&, g](){
HIP_CHECK(hipSetDevice(g));
HIP_CHECK(hipStreamSynchronize(st[g]));
});
for(int g = 0; g < 4; g++) sf[g].get();
4.3 合并 4 块 DCU 的计算结果
所有 batch 计算完成后，从 4 块 DCU 取回各自的 P90 结果，拼成完整全球网格：
for(int di = 0; di < ndo; di++){
int doy = ds + di; // 当前 DOY
// 从 4 块 DCU 各取一部分
for(int g = 0; g < 4; g++){
HIP_CHECK(hipSetDevice(g));
size_t o = g * ck;
size_t sz = gsz[g];
size_t so = di * sz; // 第 di 个 DOY 的偏移
// 取回该 DCU 负责的格点部分
HIP_CHECK(hipMemcpy(hc + o, dc[g] + so,
sz * 8, hipMemcpyDeviceToHost));
HIP_CHECK(hipMemcpy(hp + o, dp[g] + so,
sz * 8, hipMemcpyDeviceToHost));
}
// hc 和 hp 现在包含了完整的全球网格
// 写入 NetCDF 文件
write_netcdf(doy, hc, hp, lat, lon);
}
第 5 章 SST 程序的内存管理
5.1 三种内存在 SST 程序中的角色
// ── hipHostMalloc：H2D 传输缓冲区（高速）── // 存放从磁盘读出的 SST 数据 // 两个 ring 交替使用，实现流水线 double* host_ring[2]; hipHostMalloc(&host_ring[0], bB); // ~915 MB hipHostMalloc(&host_ring[1], bB); // ~915 MB // pinned memory，DMA 直传，带宽 ~12 GB/s // ── hipMalloc：DCU 显存（大容量）── // 存放传输后的 SST 数据和计算结果 double* d_ring[4]; // 输入数据环 ~228 MB/DCU float* d_mom[4]; // 4 阶矩 ~50 MB/DCU int* d_cnt[4]; // 计数 ~12 MB/DCU double* dc[4]; // 均值输出 ~2 MB/DCU double* dp[4]; // P90 输出 ~2 MB/DCU // ── malloc：CPU 内存（通用）── // 存放文件路径、临时变量等 char(*sa)[256]; // 文件路径表 double* hc; // 均值结果（CPU 暂存） double* hp; // P90 结果（CPU 暂存）
为什么 host_ring 必须用 hipHostMalloc？
因为 hipMemcpy2DAsync 要求源/目标内存是 page-locked 的。如果用 malloc，异步拷贝会退化为同步拷贝，流水线效果归零。
5.2 SST 程序退出前的清理
// 释放 pinned memory（释放 /dev/shm 空间）
for(int i = 0; i < 4; i++) hipHostFree(host_ring[i]);
// 释放 CPU 内存
free(hc); free(hp); free(sa);
// 释放 DCU stream（轻量）
for(int g = 0; g < 4; g++){
hipSetDevice(g);
hipStreamDestroy(st[g]);
// 注意：不调用 hipFree——OS 在进程退出后自动回收显存
// 跳过 hipFree 可节省约 1 秒
}
// MPI 清理
MPI_Finalize();
// 清扫 /dev/shm 残留
system("rm -rf /dev/shm/* 2>/dev/null");
第 6 章 SST 程序的异步流水线
6.1 为什么 SST 程序需要异步流水线
SST 程序的性能瓶颈分析：
一批数据的时间构成： 文件读取： ~2.3 秒 （从 Lustre 文件系统读取 20 个 NetCDF 文件） 内存传输： ~0.1 秒 （CPU → DCU，PCIe 3.0 x16） GPU 计算： ~0.2 秒 （ka_moment 核函数） ────────────────────────────── 合计： ~2.6 秒
如果串行执行，6 批需要 15.6 秒。通过流水线可以让读取和计算重叠，总时间降至约 8 秒。
6.2 双缓冲流水线的 SST 实现
// host_ring[0] 和 host_ring[1] 交替
// cur：当前 GPU 正在处理的 ring
// next：异步读正在填充的 ring
int cur = 0;
double tr0 = omp_get_wtime();
// 预先读第一批（同步）
#pragma omp parallel for
for each file in batch 0:
read_file(file, host_ring[0]);
tr += omp_get_wtime() - tr0;
for(int b = 0; b < nb; b++){
int next = 1 - cur;
// ★ 步骤 A：异步读取下一批（不阻塞 CPU）
std::future<void> rd;
if(b + 1 < nb){
rd = std::async([&](){
#pragma omp parallel for
for each file in batch b+1:
read_file(file, host_ring[next]);
});
}
// ★ 步骤 B：GPU 处理当前批（与步骤 A 同时进行）
// B1：H2D 传输
for(int g = 0; g < 4; g++)
hipMemcpy2DAsync(d_ring[g], host_ring[cur], ..., st[g]);
// B2：DCU 计算
for(int g = 0; g < 4; g++)
ka_moment<<<..., st[g]>>>(d_ring[g], ...);
// B3：同步
for(int g = 0; g < 4; g++)
hipStreamSynchronize(st[g]);
// ★ 步骤 C：等下一批读完
if(rd.valid()) rd.get();
cur = next;
}
6.3 时序图
无流水线（6 批 × 2.6s = 15.6s）： 读0 → 传0 → 算0 → 读1 → 传1 → 算1 → 读2 → ... 双缓冲流水线（约 8s）： 读0 → ─传0→算0── → ─传1→算1── → ─传2→算2── → ... ↘ 读1 ↗ ↘ 读2 ↗ ↘ 读3 ↗ （读1与传0+算0重叠，读2与传1+算1重叠，以此类推）
6.4 多环流水线（SST 项目进阶）
当 SST 数据量增大时（BATCH 更大或 ns 更长），读时间增长，2 个 ring 不够用：
2 环（读 2.3s，算 0.2s）： 读0 → 算0 → 等读1(等2.1s!) → 算1 → 等读2(等2.1s!) → ... GPU 利用率：~10% 4 环（提前 3 批启动异步读）： 读0 → ─算0→算1→算2→等读3(等1.5s)→算3→算4→算5→... GPU 利用率：~40%
环数选择原则：R ≥ T读 / T算 + 1。当 T读 = 2.3s、T算 = 0.2s 时，R ≥ 12.5，理论需要 13 环。实际中 4 环已有显著改善。
第 7 章 SST 程序的调试与性能分析
7.1 监控 SST 程序运行时 DCU 状态
# SST 程序运行时，另开一个终端 rocm-smi # 应该能观察到： # - 4 块 DCU 显存使用同时增长（每块 ~2GB） # - GPU 利用率在 ka_moment 执行期间达到 80-95% # - 温度逐渐上升，最终稳定在 60-70°C
7.2 检查 SST 数据是否正确
// 在核函数中输出调试信息（仅用于小数据调试）
__global__ void debug_kernel(const double* d, int n)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if(idx < 10){ // 只打印前 10 个格点
printf("d[%d] = %f\n", idx, d[idx]);
}
}
注意：printf 在核函数中极慢，只用于调试，不要留在正式代码中。
7.3 SST 程序常见错误
| 现象 | 在 SST 程序中的原因 | 解决方法 |
|---|---|---|
hipErrorInvalidValue |
数据量算错导致 dpitch 未 256 字节对齐 | gsz[g] * 8 必须是 256 的倍数 |
| 结果全是 NaN | 忘了从 DCU 拷回结果 | 加 hipMemcpy(DeviceToHost) |
| 部分格点结果不对 | 边界越界 | 核函数加 if(ti >= sz) return; |
| 速度比 CPU 还慢 | 忘了设 -O3 或没用 hipHostMalloc |
编译加 -O3，内存改 hipHostMalloc |
| 程序结束时卡住 | MPI_Finalize 等待所有 rank |
检查是否所有 rank 都调用了 |
| 第二次运行报 /dev/shm 不足 | 上次未清理 | 程序退出前加 rm -rf /dev/shm/* |
第 8 章 总结
8.1 从 SST 程序学到的 DCU 编程要点
| DCU 概念 | SST 程序中的对应 |
|---|---|
hipMalloc |
分配 d_ring、d_mom、d_cnt 等 DCU 显存 |
hipHostMalloc |
分配 host_ring，用于高速 H2D 传输 |
hipSetDevice |
在 4 块 DCU 之间切换 |
hipMemcpy2DAsync |
将 SST 数据按行传输到 DCU |
<<<blk, thr>>> |
每 DCU 用 4056×64 线程覆盖 259560 格点 |
hipStreamSynchronize |
每 batch 结束后等待所有 DCU |
| stream | 4 块 DCU 各用独立 stream，可并发执行 |
| wavefront (64线程) | thr=64 保证每个 wavefront 完全填满 |
8.2 SST 程序的核心性能数据
| 版本 | GPU 计算时间 | 总运行时间 | 加速比（相对基线） |
|---|---|---|---|
| 串行 CPU | (无法运行) | > 10 分钟 | — |
| 基线 GPU | ~8.0 秒 | 14.0 秒 | 1.0× |
| 统计矩优化 | 0.35 秒 | 10.8 秒 | 1.3× |
| +4 环流水线 | 0.6 秒 | 6.8 秒 | 2.1× |
| +滑动窗口矩 | 0.08 秒 | 4.2 秒 | 3.3× |
GPU 计算时间从 8.0 秒降到 0.08 秒（99% 降幅），总时间从 14.0 秒降到 4.2 秒（3.3× 加速）。
8.3 进阶学习路径
-
理解 SST 程序的完整流程 ← 你现在在这里
-
学习
__launch_bounds__寄存器优化 — 进一步提升 kernel 性能 -
学习 AVX2 CPU 预筛 — 让 CPU 在等待 I/O 时做有用的事
-
学习滑动窗口矩 — 利用数学性质减少 GPU 运算量
-
学习多环流水线调优 — 根据实际 T读/T算 调整环数
参考文献
[1] ROCm Documentation: HIP Programming Guide. https://rocm.docs.amd.com [2] Hobday, A.J., et al. "A hierarchical approach to defining marine heatwaves." Progress in Oceanography, 141:227-238, 2016. [3] Cornish, E.A. & Fisher, R.A. "Moments and cumulants in the specification of distributions." Revue de l'Institut International de Statistique, 5(4):307-320, 1937. [4] NetCDF Documentation. https://www.unidata.ucar.edu/software/netcdf/
免费领 150 小时云算力，进群参与显卡、AI PC 幸运抽奖
更多推荐
3
0- 0
已为社区贡献2条内容
所有评论(0)