6
6
1
15
海光 DCU 编程入门——以气候态海温百分位计算为教学案例
2026-07-02 12:11:47
28
## 摘要
**本文面向海光 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<<>>(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<<>>(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<<>>(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 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 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** |
| `<<>>` | **每 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 进阶学习路径
1. **理解 SST 程序的完整流程** ← 你现在在这里
2. **学习 `__launch_bounds__` 寄存器优化** — 进一步提升 kernel 性能
3. **学习 AVX2 CPU 预筛** — 让 CPU 在等待 I/O 时做有用的事
4. **学习滑动窗口矩** — 利用数学性质减少 GPU 运算量
5. **学习多环流水线调优** — 根据实际 T*读/T*算 调整环数
---
## 参考文献
**[1] ROCm Documentation: HIP Programming Guide. **[https://rocm.docs.amd.com](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/](https://www.unidata.ucar.edu/software/netcdf/)
全部评论（0）
只看楼主
暂无数据--
楼主
0 点赞
0 帖子
楼主热帖
暂无数据常见问题新帖
暂无数据