> 本文为**先导杯创作者激励计划**投稿文章，以一个在 DCU 上实现 ~14,000× 加速的真实计算项目为线索，系统梳理海光 DCU 的高阶编程技巧与性能优化方法论。 --- ## 一、先导杯与 DCU 实战 **先导杯**是国内规模最大的国产异构计算赛事之一，汇聚了来自高校和科研院所的众多高性能计算爱好者。今年的赛题中有多项任务基于**海光 DCU** 平台，参赛者需要在实际的 DCU 集群上完成代码开发与性能调优。 本文将以一个典型的大规模数据处理任务为例——从 10,950 个 NetCDF 文件的海洋温度数据中计算气候态统计量——展示如何在海光 DCU 平台上完成从算法设计到硬件极限的全流程优化。该项目最终在 **2 节点 × 4 海光 DCU Z100L** 上实现了 **4.75 秒** 的墙钟成绩，相对于初始串行版本实现了近 **14,000 倍加速**。 无论你是在准备先导杯，还是正在学习 DCU 编程，本文中的技术思路和代码示例都值得参考。 --- ## 二、项目全景：大规模科学计算在 DCU 上的落地 ### 2.1 计算任务概述 - **输入**：10,950 个 NetCDF 格式的逐日海表温度文件（30 年 × 365 天），每个文件 1,166 × 721 网格点 - **计算**：对夏季（6–8 月）92 天，以 11 天滑动窗口（±5 天）统计 30 年间每个像素点的气候态均值（Clim）和 90% 分位数（P90） - **输出**：92 个 NetCDF 文件 - **验证**：与参考数据逐日交叉验证 RMSE 这是一个典型的**数据密集型 + 尴尬并行**任务——大量 I/O 加上天然可分割的计算负载。 ### 2.2 硬件环境 | 组件 | 配置 | | ------ | ----------------------------------------- | | 加速卡 | 8 × 海光 DCU Z100L（gfx906），每节点 4 卡 | | CPU | 海光 x86-64 处理器 | | 存储 | ParaStor/Lustre 并行文件系统 | | 互连 | 节点间无通信需求（独立计算各 DOY 段） | ### 2.3 软件栈 | 层级 | 软件 | 说明 | | ---------- | ------------------------ | ----------------------- | | 操作系统 | CentOS 7.6 | 国产化适配 | | 作业调度 | SLURM | 多节点任务分配 | | GPU 运行时 | DTK / ROCm 25.04.4 | 海光加速卡工具链 | | 编译器 | hipcc（Clang/LLVM） | `--offload-arch=gfx906` | | 并行模型 | std::thread + HIP Stream | C++14 | | 数据格式 | NetCDF 4.3.0 / HDF5 1.8 | 科学数据存储标准 | | 性能分析 | rocprof | 硬件计数器采集 | --- ## 三、核心代码深度剖析 ### 3.1 HIP GPU 内核设计：一核两用，均值与分位数并行计算 在海光 DCU 上编程，GPU kernel 的设计决定了计算效率的上限。下面的内核函数同时完成均值和 90% 分位数的计算： ```cpp __global__ void compute_stats_kernel( const short *__restrict__ pool, // 预加载到显存的数据池 const int *__restrict__ doy_to_pool, // DOY → 池索引映射表 int doy_count, int doy_offset, float scale_factor, // 数据解码：scale_factor float add_offset, // 数据解码：add_offset short fill_value, // 填充值标记 double *__restrict__ out_mean, // 输出：均值 double *__restrict__ out_p90) // 输出：P90分位数 { // 2D Grid 组织：X = 空间像素，Y = 时间维度（DOY） int pixel = blockIdx.x * blockDim.x + threadIdx.x; int d_idx = blockIdx.y; if (pixel >= GRID_SIZE || d_idx >= doy_count) return; double sum = 0.0; int nv = 0; // Top-34 降序数组：P90 只需最大的 34 个值（330 × 0.1 + 0.5 ≈ 33.5） // 34 × 4B = 136 字节，完全驻留 GPU 寄存器，零 LDS 通信开销 float sorted[35]; #pragma unroll for (int i = 0; i < 34; i++) sorted[i] = -1e9f; // 滑动窗口：每个 DOY 的 11 天窗口映射到 pool 索引 int pool_ix[11]; for (int w = 0; w < 11; w++) { int d = tdoy - 5 + w; while (d < 1) d += 365; while (d > 365) d -= 365; pool_ix[w] = doy_to_pool[d]; } // 遍历 330 个样本，同时累积均值和维护 top-34 for (int w = 0; w < 11; w++) { int pi = pool_ix[w]; if (pi < 0) continue; const short *blk = pool + (size_t)pi * 30 * GRID_SIZE; for (int yr = 0; yr < 30; yr++) { short raw = blk[(size_t)yr * GRID_SIZE + pixel]; if (raw == fill_value) continue; // 解码：short → float（量化存储的反向操作） float v = (float)((double)raw * scale_factor + add_offset); nv++; sum += (double)v; // 插入排序维护 top-34 降序数组 // 因为 GPU 是内存瓶颈，ALU 闲置，多算的 34 次比较"免费" if (v > sorted[33]) { #pragma unroll for (int i = 33; i >= 0; i--) { if (v > sorted[i]) { sorted[i + 1] = sorted[i]; sorted[i] = v; } } } } } // 输出均值 int idx = d_idx * GRID_SIZE + pixel; out_mean[idx] = (nv > 0) ? sum / (double)nv : NAN; // Hyndman & Fan type=5 线性插值 P90 if (nv == 0) { out_p90[idx] = NAN; return; } double r = nv * 0.9 + 0.5; // 对齐标准分位数公式 int k = (int)floor(r); double f = r - (double)k; int needed = nv - k + 1; out_p90[idx] = sorted[needed - 1] + f * (sorted[needed - 2] - sorted[needed - 1]); } ``` **设计要点解读：** - **2D Grid 布局**：`blockIdx.x` 遍历百万级像素，`blockIdx.y` 按 DOY 分组。这种布局天然适配科学计算中"同一算法对不同时空点独立执行"的尴尬并行模式，也是 DCU 上最推荐的 kernel 组织方式。 - **Top-34 代替全排序**：当需要 P90 分位数（即第 90 百分位）且样本量固定为 330 时，只需维护最大的 34 个值。这 136 字节的数据完全驻留在寄存器中，避免了昂贵的高带宽内存（HBM）访问。 - **有序数组 vs 堆的选择**：理论上课本会推荐最小堆（O(log N) vs O(N)），但 rocprof 硬件计数器显示 DCU 的内存单元占用率已达 99.7%——内核是纯内存瓶颈。ALU 闲置的算力让"多余"的比较操作实际零成本，而有序数组带来了更好的 warp 一致性和编译器优化空间。 - **分位数公式对齐**：`r = nv * 0.9 + 0.5` 是 Hyndman & Fan type=5 插值法的标准形式。在实际项目中，仅这一行公式的错误（误用 type=7 而非 type=5）就导致了 0.0126°C 的精度偏差——数值计算的细节往往决定了科学结果的正确性。 ### 3.2 内存优化：量化存储 从 float（4 字节）降至 short（2 字节）存储中间数据： ```cpp // 读取阶段：double → short 量化 static bool read_sst_short(int dir_fd, const char *fname, short *dest, float *tmp, double *raw_buf) { int fd = openat(dir_fd, fname, O_RDONLY); pread(fd, raw_buf, GRID_SIZE * sizeof(double), g_data_offset); close(fd); const double inv = 1.0 / g_scale; #pragma GCC ivdep for (int i = 0; i < GRID_SIZE; i++) { double sv = raw_buf[i] * inv; dest[i] = (sv == sv) ? (short)lround(sv) : g_fill; } return true; } // 计算阶段：在 GPU kernel 中解码回 float float v = (float)((double)raw * scale_factor + add_offset); ``` **收益远不止 50% 的内存节省：** - 数据池从 12.7 GB 降至 6.4 GB，减少内存压力 - H2D（Host to Device）传输量减半 - IOMMU 页表压力减半 - 精度损失仅 0.0002°C（均值）和 0.0021°C（P90），远小于典型观测误差 0.1–0.5°C ### 3.3 I/O 优化：绕过中间层的裸读 NetCDF 格式底层基于 HDF5，标准 API 路径 `nc_open → nc_get_vara_float` 经过多层抽象：NetCDF 层 → HDF5 层 → B-tree 元数据查询 → 实际磁盘读取。当数据文件是"创建后不再修改"的连续存储时，这些元数据层纯粹是开销。 **关键技术：pread 直接定位数据块** ```cpp // 启动时自动探测数据偏移量 static off_t detect_data_offset(const std::string &nc_path) { // 扫描前几个文件，用 HDF5 C API 探测 "data" 数据集的磁盘偏移 for (const auto &fn : sample_files) { hid_t fid = H5Fopen(fn.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); hid_t dsid = H5Dopen2(fid, "data", H5P_DEFAULT); haddr_t addr = H5Dget_offset(dsid); // 获取物理偏移 H5Dclose(dsid); H5Fclose(fid); // 验证所有文件偏移一致 if (addr != first_offset) return 0; // 不一致则回退 } return first_offset; // 如 23488 } // 使用探测到的偏移量，直接 pread 裸读 pread(fd, raw_buf, GRID_SIZE * sizeof(double), g_data_offset); // 跳过 NetCDF→HDF5→B-tree 全部元数据解析 ``` **效果**：单次文件读取从数十次函数调用降至 1 次系统调用，预加载耗时降低约 30%。 ### 3.4 并发调度：工作窃取预加载器 传统做法是将任务静态分区分配给线程（线程 0 负责第 1-3 天，线程 1 负责第 4-6 天…），但在网络文件系统（如 Lustre）环境下，I/O 延迟抖动很大——快的线程会空闲等待慢的线程。 **工作窃取（Work Stealing）调度器**让所有线程通过原子操作竞争任务： ```cpp std::atomic next_pool_idx{0}; // 36 个预加载线程，通过 atomic fetch_add 竞争 DOY 索引 for (int w = 0; w < 36; w++) { preload_threads.emplace_back([&]() { float *tmp = new float[GRID_SIZE]; double *raw = new double[GRID_SIZE]; int pi; // 原子竞争：快的线程自动多抢 DOY，慢的少抢 while ((pi = next_pool_idx.fetch_add(1, std::memory_order_relaxed)) < node_pool_doys) { int doy = node_pool_first + pi; short *blk = pool + (size_t)pi * N_YEARS * GRID_SIZE; for (int yr = 0; yr < N_YEARS; yr++) { short *dest = blk + (size_t)yr * GRID_SIZE; read_sst_short(dir_fd, ..., dest, tmp, raw); } // 标记该 DOY 就绪，通知 GPU 流水线可以启动 pool_doy_ready[pi].store(1, std::memory_order_release); preload_cv.notify_all(); } delete[] raw; delete[] tmp; }); } ``` 这种调度方式在异构 I/O 环境下自动实现了负载均衡，省去了手动调优 worker 分配策略的工作量。 ### 3.5 I/O 与计算流水线重叠 GPU 不需要等所有数据加载完才开始计算——当某个 GPU 所需的 DOY 数据全部就绪时，即可立即启动： ```cpp for (int g = 0; g < gpus_per_node; g++) { // 阻塞等待：该 GPU 所需的 DOY 范围全部加载完成 std::unique_lock lock(preload_mutex); preload_cv.wait(lock, [&]() { for (int pi = gpu_pool_first_idx[g]; pi <= gpu_pool_last_idx[g]; pi++) { if (pool_doy_ready[pi].load( std::memory_order_acquire) == 0) return false; } return true; }); // 数据就绪，立即启动该 GPU 的计算线程 gpu_threads.emplace_back(gpu_worker, g, ...); } ``` 这样 GPU 计算和后续数据预加载完全重叠，墙钟 = max(I/O, GPU) + 少量串行开销，而非 I/O + GPU。 ### 3.6 双 HIP Stream 流水线 在单张 DCU 内部，利用 HIP Stream 实现 H2D 传输 → Kernel 执行 → D2H 传输的流水线重叠： ```cpp hipStream_t s0, s1; hipStreamCreate(&s0); hipStreamCreate(&s1); int half0 = doy_count / 2; int half1 = doy_count - half0; // Stream 0：前半数据 → Kernel → 结果回传 hipMemcpyAsync(d_pool, h_pool, half0_bytes, hipMemcpyHostToDevice, s0); compute_stats_kernel<<>>(d_pool, ..., d_mean, d_p90); hipMemcpyAsync(h_mean, d_mean, half0_bytes, hipMemcpyDeviceToHost, s0); // Stream 1：后半数据 → Kernel → 结果回传（与 Stream 0 并发） hipMemcpyAsync(d_pool + half0, h_pool + half0, half1_bytes, hipMemcpyHostToDevice, s1); compute_stats_kernel<<>>(d_pool, ..., d_mean1, d_p90_1); hipMemcpyAsync(h_mean + half0, d_mean + half0, half1_bytes, hipMemcpyDeviceToHost, s1); hipStreamSynchronize(s0); hipStreamSynchronize(s1); ``` HIP Stream 的多流并发模式是 CUDA Stream 的国产化等价物——对于熟悉 CUDA 的开发者，迁移到 DCU 时编程模型几乎无需改变。 --- ## 四、硬件计数器驱动的性能分析 在海光 DCU 上进行性能调优，**rocprof** 是必不可少的工具。它从硬件计数器层面告诉你程序到底瓶颈在哪里： ```bash rocprof --stats ./your_dcu_program ``` 典型输出（gfx906 架构，单 GPU 模式）： | 指标 | 数值 | 解读 | | ---------------- | :--------: | ----------------------------------------------- | | **MemUnitBusy** | **99.71%** | 内存单元满负荷——HBM 带宽吃满 | | VALUUtilization | 58.81% | 向量 ALU 近半时间在等数据到来 | | L2CacheHit | 67.21% | L2 命中率尚可，但 33% 的 miss 仍在触发 HBM 访问 | | WriteUnitStalled | 23.49% | 写回单元约 1/4 时间在等待 | | FETCH_SIZE | 161 MB | 单次 kernel 从 HBM 读取的数据量 | | Occupancy | ~25% | 每 SIMD 仅 1 个 wave（VGPR=160 寄存器限制） | **核心结论：当 MemUnitBusy 达到 99.7% 时，继续优化 kernel 代码的 ALU 效率是徒劳的。** HBM 带宽已被榨干，此时优化策略应该转向： - 减少每次 kernel 的数据读取量（如量化存储、更紧凑的数据布局） - 改善 L2 缓存命中率（如调整访问模式、数据重排） - 减少 Host ↔ Device 数据传输（如使用固定内存 hipHostRegister） --- ## 五、优化历程中的关键方法论 从该项目 24 个版本的迭代中，提炼几条适用于 DCU 开发的通用原则： ### 5.1 Profiler 先行，数据驱动 不要凭直觉猜测瓶颈。项目中初始版本的 MATLAB Profiler 显示 NetCDF I/O 占 68%、日期格式化占 28%、真正的数学计算不到 1%。这个数据直接决定了后续所有优化都围绕 I/O 展开，而非去抠数值算法的细节。 **实战建议**：CPU 端用 perf stat，GPU 端用 rocprof，每一轮优化后对比计数器变化。 ### 5.2 不要迷信框架 任务天然是"尴尬并行"（各 DOY 独立计算，节点间零通信），引入 MPI 反而增加了 ~3.5s 的启动开销和大量代码复杂度。最终方案是两个独立进程通过 SLURM `srun -n 2` 启动，零 MPI 依赖。 **实战建议**：在选择分布式框架前，先问自己"节点间真正需要交换什么数据"。如果答案是"什么都不需要"，就不要引入通信框架。 ### 5.3 线程 vs 进程的选择取决于底层 API 的线程安全性 早期版本使用 `std::thread` 多线程并发读取 NetCDF 文件时频繁崩溃——根因是底层 HDF5 1.8 库的全局内存池非线程安全。解决方案是改用 `fork()` 创建子进程，每个子进程拥有独立的 HDF5 状态。 但在引入 pread 裸读（绕过 HDF5 库）后，又切回了 `std::thread`——因为 `pread()` 是线程安全的系统调用，fork 的进程隔离成为不必要的开销。 **实战建议**：理解你依赖的库的线程安全语义，据此决定并行模型。当底层库有全局状态竞争时，进程隔离是有效的兜底方案。 ### 5.4 知道何时停止 最终版的性能画像： - CPU：IPC 2.44，L1 miss 0.69%，Branch miss 0.22% —— CPU 端已无优化空间 - GPU：MemUnitBusy 99.7% —— HBM 带宽已吃满 - 墙钟占比：I/O 75%，GPU 25% 此时继续优化 kernel 代码的性价比为零。进一步提升需要升级存储（如本地 NVMe 缓存）或换用更新的 DCU 架构。 **实战建议**：用硬件计数器验证是否触达硬件极限。触达后，把精力留给架构级改进而非微优化。 --- ## 六、完整的 SLURM 编译运行脚本 以下是一个可直接参考的海光 DCU 项目编译与提交模板： ```bash #!/bin/bash #SBATCH -p kshdmcc2026 # DCU 专用计算分区 #SBATCH -N 2 # 申请 2 个节点 #SBATCH -n 2 # 2 个任务（每节点 1 个独立进程） #SBATCH --gres=dcu:4 # 每节点 4 张 DCU 加速卡 #SBATCH --exclusive # 节点独占 #SBATCH --mem=0 # 使用节点全部可用内存 #SBATCH -J dcu_project set -e # ── 1. 自动探测可用的 DTK 版本 ────────────────────────────── DTK_VER="" for try_ver in "dtk-25.04.4" "dtk-24.04.1" "dtk-22.10.1"; do candidate="/public/software/compiler/rocm/${try_ver}/hip/bin/hipcc" if [ -f "$candidate" ]; then DTK_VER="$try_ver" break fi done [ -z "$DTK_VER" ] && { echo "错误: 找不到可用的 DTK 工具链"; exit 1; } ROCM_PATH="/public/software/compiler/rocm/${DTK_VER}" HIPCC="${ROCM_PATH}/hip/bin/hipcc" export LD_LIBRARY_PATH="${ROCM_PATH}/lib:${ROCM_PATH}/hip/lib:$LD_LIBRARY_PATH" # ── 2. 编译（面向 gfx906 架构 ─────────────────────────── cp source.cpp source.hip $HIPCC -std=c++14 -O3 -march=native \ --offload-arch=gfx906 \ -I"${NC_INC}" -I"${H5_INC}" \ -o dcu_program source.hip \ -L"${NC_LIB}" -l:libnetcdf.so.7 -l:libhdf5.so.8 rm -f source.hip # ── 3. 双节点运行 ────────────────────────────────────────── time srun -N 2 -n 2 --ntasks-per-node=1 \ ./dcu_program --gpus-per-node 4 ``` --- ## 七、总结 海光 DCU 的 HIP 编程模型与 CUDA 高度相似，从 NVIDIA GPU 迁移的学习成本很低。但要在 DCU 上真正"用好算力"，核心挑战往往不在 kernel 代码本身，而在于： 1. **I/O 数据供给**：I/O 瓶颈往往是整个流水线的短板，需要从文件系统特性、数据格式、传输策略三个维度入手 2. **硬件计数器驱动决策**：rocprof 的 MemUnitBusy、perf 的 IPC 等硬件计数器是判断"是否还需要继续优化"的硬指标 3. **架构简洁性**：先跑通最简单的并行方案，再用 profiler 定位瓶颈，逐步迭代——不要一开始就引入复杂框架 希望本文能帮助你在海光 DCU 上的开发与优化少走弯路。也欢迎大家关注**先导杯**赛事，在真实的国产异构平台上锤炼自己的 HPC 技能。 --- > 上一篇：[海光 DCU 实战：基本指令与操作](./02-海光DCU实战-基本指令与操作.md) > > 上一篇：[海光 DCU 基础认知与核心优势](./01-海光DCU基础认知与核心优势.md)