海光 DCU 进阶：从真实项目看高性能计算优化全流程

本文为先导杯创作者激励计划投稿文章，以一个在 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% 分位数的计算：

__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 字节）存储中间数据：

// 读取阶段：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 直接定位数据块

// 启动时自动探测数据偏移量
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）调度器让所有线程通过原子操作竞争任务：

std::atomic<int> 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 数据全部就绪时，即可立即启动：

for (int g = 0; g < gpus_per_node; g++) {
    // 阻塞等待：该 GPU 所需的 DOY 范围全部加载完成
    std::unique_lock<std::mutex> 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 传输的流水线重叠：

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<<<grid0, block, 0, s0>>>(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<<<grid1, block, 0, s1>>>(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 是必不可少的工具。它从硬件计数器层面告诉你程序到底瓶颈在哪里：

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 项目编译与提交模板：

#!/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 实战：基本指令与操作

上一篇：海光 DCU 基础认知与核心优势