# 海光DCU深度使用手记（中）：优化攻坚——kernel设计、显存博弈与踩坑实录 ## 前言 上篇讲了海光DCU的基本情况和第一次移植，最终端到端跑了21.5秒——DCU纯计算只要4.4秒，但I/O拖了后腿。这一篇我们深入DCU的内部优化：怎么设计kernel才能把DCU的计算单元喂饱、怎么管理显存避免OOM、以及整个过程中踩过的各种真实的坑。 先剧透一下：经过几轮迭代，我们最终把DCU方案的端到端时间压到了理论下限附近，虽然最后因为I/O策略的取舍没有选DCU作为最终提交方案（这个决策本身也值得讨论），但整个过程积累的DCU开发经验比最终跑出来的数字值钱得多。 ## 一、Kernel设计的核心：时序分块（Temporal Blocking） ### 1.1 问题的本质 我们的任务是对92个DOY，每个DOY收集330个样本（30年×11天窗口）做统计。最naive的GPU方案是每个DOY独立算——一次kernel launch处理一个DOY。这样要启动92次kernel，每次处理330个网格的数据，总共要读92×330=30360个网格。但相邻DOY的窗口高度重叠——DOY 152和DOY 153的窗口有10天是重复的。如果每个DOY独立读数据，同一个网格文件会被DCU从显存里读11次，显存带宽完全浪费在重复读取上。 这就是时序分块的动机：**把多个连续DOY的计算合并到一个kernel里，共享重叠的窗口数据，减少显存读取次数**。 ### 1.2 具体设计 我们取B=10个DOY为一个批次。10个连续DOY覆盖的窗口范围是：最早DOY的窗口起始到最晚DOY的窗口结束，跨度是B + WINDOW_DAYS - 1 = 10 + 11 - 1 = 20天。每个批次只需要20个唯一网格（per year），而不是之前的10×11=110个独立读取。 ``` 批次内10个DOY共享的网格窗口： DOY 152: [147, 148, ..., 157] ← 11天 DOY 153: [148, 149, ..., 158] ... DOY 161: [157, 158, ..., 167] └── 总共20个唯一网格（per year） ``` 显存读取量从每个DOY 330次独立采样，变为每个批次共享20×30=600个网格（一次性加载到共享内存），读取量减少约5.5倍。 ### 1.3 Kernel实现要点 ```cpp __global__ void temporal_batch_kernel( const float *__restrict__ pool, const unsigned char *__restrict__ land_mask, double *__restrict__ clim_out, double *__restrict__ p90_out, size_t grid_cells, int B, // 批次DOY数（通常=10） int batch_total_days, // 批次覆盖的唯一天数（=20） int n_unique, // 唯一网格数（=30年×20天=600） float fill_val, int has_fill ) { const int T = 16; size_t c = (size_t)blockIdx.x * T + threadIdx.x; if (c >= grid_cells) return; // 共享内存：协同加载本workgroup需要的所有网格数据 __shared__ float smem[600 * 16]; // 600 grids × 16 cells // Phase 1: 16个线程协同加载600个网格到共享内存 // 每个线程负责加载 ~600/16 ≈ 38个grid的(cell位置)数据 for (int i = threadIdx.x; i < n_unique; i += T) { smem[i * T + threadIdx.x] = pool[(size_t)i * grid_cells + c]; } __syncthreads(); // 确保所有网格数据加载完成 // Phase 2: 对批次内的每个DOY分别计算（共享内存中读取，极快） for (int j = 0; j < B; j++) { size_t out_idx = (size_t)j * grid_cells + c; if (land_mask[c]) { clim_out[out_idx] = -9999.0; p90_out[out_idx] = -9999.0; continue; } // 从共享内存收集330个样本 float samples[330]; int ns = 0; for (int yi = 0; yi < 30; yi++) { for (int di = 0; di < 11; di++) { int grid_idx = yi * batch_total_days + j + di; float v = smem[grid_idx * T + threadIdx.x]; if (isfinite(v) && (!has_fill || v != fill_val)) { samples[ns++] = v; } } } // 计算均值和P90（在寄存器/本地内存中操作samples数组） if (ns == 0) { clim_out[out_idx] = NAN; p90_out[out_idx] = NAN; } else { double sum = 0.0; for (int k = 0; k < ns; k++) sum += samples[k]; clim_out[out_idx] = sum / ns; // Quickselect for P90 // ... (下节详述) } } } ``` 这个kernel有几个设计要点： **第一，用共享内存做协同加载**。600个网格、每个网格取16个cell的数据 = 9600个float ≈ 37.5KB共享内存。海光DCU（gfx906架构）的LDS（Local Data Share）是64KB per CU，37.5KB刚好能塞进去还有余量。加载完成后用 `__syncthreads()`确保所有线程都拿到数据，然后Phase 2的读取全部在共享内存里完成——延迟从HBM2的几百个cycle降到LDS的十几个cycle。 **第二，land_mask提前判断**。陆地的格点没有有效SST数据，直接填NaN跳过。海洋格点才做真正的样本收集。这个分支在GPU上不会造成太大的warp divergence，因为大部分格点要么全是海洋要么全是陆地，一个warp内部的执行路径比较一致。 **第三，samples数组放寄存器/私有内存**。每个线程的 `samples[330]`约1320字节，DCU的寄存器文件一个CU有256KB，16个线程×1320字节≈21KB——寄存器装不下的时候会溢出到L1缓存，但L1也比HBM2快一个数量级。 ### 1.4 GPU上的Quickselect P90分位数的计算在CPU上用的是 `std::nth_element`——部分排序，O(n)复杂度。在DCU上，标准库没法用，得手写一个device端的quickselect： ```cpp __device__ float dev_quickselect_kth(float *arr, int n, int k) { int lo = 0, hi = n - 1; while (lo < hi) { float pivot = arr[(lo + hi) / 2]; int a = lo, b = hi; while (a <= b) { while (arr[a] < pivot) a++; while (arr[b] > pivot) b--; if (a <= b) { float t = arr[a]; arr[a] = arr[b]; arr[b] = t; a++; b--; } } if (k <= b) hi = b; else if (k >= a) lo = a; else break; } return arr[k]; } ``` 这个quickselect有几个需要注意的地方： 1. **不要用递归**。GPU的调用栈非常浅，递归的quickselect大概率栈溢出。老老实实写while循环。 2. **pivot选中间值**。虽然不是最优策略，但在GPU上做复杂的pivot选择（比如median-of-medians）得不偿失——额外的分支和内存访问比省下的比较次数更贵。 3. **arr是线程私有的**。每个线程独立操作自己的 `samples[330]`，不涉及线程间同步，所以不需要atomic操作。 计算P90（90%分位数）时，配对线性插值（精确匹配MATLAB的 `prctile`行为）： ```cpp // MATLAB prctile线性插值: r = (N-1) * P/100 double r = (ns - 1) * 0.9; int lo = (int)r; double w = r - lo; // 小数部分 // 第lo个元素 float v_lo = dev_quickselect_kth(samples, ns, lo); float v_hi; if (lo + 1 < ns) { v_hi = dev_quickselect_kth(samples, ns, lo + 1); } else { v_hi = v_lo; } p90_out[out_idx] = v_lo + w * (v_hi - v_lo); ``` 这里有个效率上的trade-off：每个格点要调两次quickselect（一次找第lo个、一次找第lo+1个）。如果ns很大（接近330），两次quickselect的重复搜索开销不小。但实测下来，因为samples数组在L1/寄存器里，330个元素的quickselect非常快（约几百ns），每个格点调两次完全可接受。 ## 二、显存管理：怎么在16GB里塞下一切 ### 2.1 显存压力分析 我们的数据规模： - 全量输入：3060个网格 × 1440×721格点 × 4字节(float) = **12.7 GB** - 输出Climmean：92 DOY × 1440×721 × 8字节(double) = **0.76 GB** - 输出P90：92 DOY × 1440×721 × 8字节(double) = **0.76 GB** - 陆地掩码：1440×721 × 1字节 = **1 MB**（可忽略） - **DCU显存总需求 ≈ 14.2 GB** 当时我们用的DCU大概率是16GB HBM2的版本，14.2GB已经非常接近上限了。如果再考虑DCU驱动和运行时的一些开销，基本是贴着极限在跑。 这就是为什么我们选择做**全量预读 + 单次H2D**策略的原因之一：如果分批处理，每次都要重新从磁盘读文件再做H2D，I/O压力更大。还不如一次性全读进来，DCU算完就释放。 ### 2.2 显存优化手段 如果你的数据规模更大（比如更多年份、更高分辨率），以下几种策略可以考虑： **（1）分批处理**。不一次性上传全部3060个网格，而是分2-3批。每批上传一批、计算对应的DOY、下载结果、释放显存。代价是H2D/D2H次数增多。 **（2）用float替代double存储中间结果**。虽然MATLAB精度验证要求高，但我们实测P90用float算，误差仍然在0.01°C量级，远低于赛题的1°C要求。输出时再转double，显存可以省一半。 **（3）land_mask压缩**。用1bit表示每个格点（bitset），而不是1字节的 `unsigned char`。1440×721个格点只需要约130KB，还可以进一步减少随机访问。不过DCU上没有原生的bit操作指令，需要用位运算模拟，代码会复杂一些。 **（4）即时计算陆地掩码**。不给land_mask分配独立显存，而是在kernel里通过判断样本是否全为NaN来动态判定陆地格点。但这会让每个线程多读好几次float，而且逻辑分支更复杂。 在我们的场景中，因为16GB刚好够用，没有进一步压显存。但如果你用的是8GB版本的海光DCU，上述手段就得认真考虑了。 ### 2.3 hipHostRegister的坑 上篇提到用 `hipHostRegister`pin住mmap内存可以加速H2D传输。这里有一个实际踩过的坑： `hipHostRegister`之后必须调用 `hipHostUnregister`才能安全释放内存，否则DCU驱动还在引用这块物理页。而且**先unregister再munmap**，顺序不能反过来： ```cpp // 正确顺序 HIP_CHECK(hipHostUnregister(full_data)); // 先解除DCU的pin munmap(full_data, full_buf_cells * sizeof(float)); // 再释放内存 // 错误顺序——可能导致segfault munmap(full_data, full_buf_cells * sizeof(float)); hipHostUnregister(full_data); // 操作已释放的内存！ ``` 这个bug的隐蔽之处在于它不一定会立即崩溃——有时候munmap后的虚拟地址恰巧还没被重新分配，`hipHostUnregister`可能"成功"返回，但在后续某个时间点突然segfault。我们调了快一个小时才发现是这个顺序问题。 ## 三、I/O策略的演进：从串行等待到全量预读 ### 3.1 Round G1：最初方案的问题 上篇提到第一版DCU方案跑了21.5秒，其中： - DCU计算：4.4s（很快） - Upfront I/O：~12s（瓶颈！） - Refill wait：14.65s（更大的瓶颈！） 等一下，为什么还有"Refill wait"？因为G1最初的设计是分批次处理的——不是一次性全读。每批10个DOY算完后，要等下一批的数据从磁盘读进来、上传到DCU。这个"等"的时间就有14.65秒，占了总时间的68%。 问题的根因是我们用的Lustre并行文件系统：同时打开300个文件做读取，MDS（元数据服务器）根本处理不过来，单个文件的打开延迟从几个毫秒飙升到几十甚至上百毫秒。**批量小文件的元数据操作是Lustre的死穴**，也是很多HPC应用的共同痛点。 ### 3.2 Round G3：全量预读策略 既然分批读慢，那就一次性全读。Round G3的方案是： ``` 启动时一次性读完所有3060个文件到CPU内存（mmap ~16GB） ↓ 单次H2D上传到DCU显存 ↓ DCU计算全部92个DOY（10个批次，每批10个DOY，最后一批2个） ↓ 单次D2H下载全部结果 ↓ 分批写出92个NetCDF结果文件 ``` 关键改动： - 初始I/O用6个进程（不是30个）并行读——减少Lustre MDS竞争 - Refill从"fork 30进程读300文件"变成"memcpy从CPU缓冲区"——耗时从1.6秒降到0.05秒 - 3060个文件一次性读完后，DCU计算期间零磁盘I/O 这个策略把I/O等待从14.65秒降到了约0秒（除了最开始的12秒upfront read），DCU计算时间4.4秒不变，整体在理论上可以达到16-17秒左右。虽然最终没有在超算上完整验证这个版本（后面会讲为什么），但从设计逻辑上是自洽的。 ### 3.3 为什么最终选了CPU方案而不是DCU方案 这是很多人会问的问题——DCU不是更快吗？为什么不提交DCU版本？ 说实话，这个决策当时队内也讨论了很久。最终选CPU方案，基于几个判断： 1. **稳定性优先**。CPU版本经过十几轮迭代，在超算上跑了无数次，稳定性已经充分验证。DCU版本在比赛截止前只做了有限的测试——而比赛规则明确说精度是第一位的，速度是加分项。万一DCU版本有隐藏的精度问题（比如HIP math库的浮点行为跟CPU不完全一致），风险太大。 2. **可复现性**。CPU版本只需要GCC+OpenMP+NetCDF三个依赖，任何超算节点都能跑。DCU版本需要特定的DTK版本和DCU硬件——如果评测节点的DCU型号不同、DTK版本不同，可能直接跑不起来。 3. **时间成本**。比赛截止时，CPU版本已经稳定跑9秒，足够拿满分了。继续调DCU版本的边际收益有限。但如果比赛允许更长的优化时间，DCU方案的潜力肯定更大——纯计算4.4秒 vs CPU的12秒，差距是实打实的。 这个决策在赛后复盘的时候我们仍然认为是正确的——**在有限时间和资源下，稳定可靠的方案优于潜力更大但不稳定的方案**。但这不代表DCU的探索没有价值——恰恰相反，它让我们积累了整套的DCU开发经验，这些经验在后续的项目中一定会用上。 ## 四、踩坑记录：DCU开发中的各种真实问题 ### 坑1：netCDF 4.4.1不是线程安全的 这个问题在CPU版本上就坑过我们一次，DCU版本又遇到了。netCDF底层依赖HDF5 1.8.x，而HDF5 1.8.x内部有全局状态，多线程并发 `nc_open`会导致随机segfault。 **解决方案**：永远用进程而不是线程做并行I/O。`fork()`创建子进程，每个子进程有独立的netCDF状态，互不干扰。子进程完成后用 `_exit(0)`直接退出（不用 `exit(0)`，避免调父进程的atexit处理函数）。 这条教训的价值超越了DCU本身——它适用于任何需要并行NetCDF/HDF5 I/O的场景。 ### 坑2：SLURM cgroup的进程数限制 我们在CPU版本上尝试过"预fork持久化工作进程池"——启动时fork 30个worker，用pipe信号触发读写。想法很好：避免热路径中反复fork的开销。结果程序一启动就被SIGKILL杀了。 原因：SLURM的cgroup会统计进程数。30个长驻子进程+父进程=31个进程，超过了cgroup的pids限制。但短命子进程（fork后几毫秒就_exit退出）不受影响，因为cgroup采样有一个时间窗口。 **教训**：在SLURM管理的超算环境里，长驻子进程的开销不仅仅是内存，还有cgroup限制。短命子进程虽然每次fork有开销，但更安全。 ### 坑3：DCU显存碎片 我们遇到了一个有意思的问题：`hipMalloc`返回成功，但后续的kernel launch报"out of memory"。排查发现是显存碎片导致的——虽然总空闲显存够，但没有足够大的连续块。 **解决方案**：把大块显存在程序启动时一次性分配完（先占坑），避免中间穿插小块分配和释放。这个做法有点粗暴但很有效。 ### 坑4：hipcc与GCC的ABI兼容 `hipcc`基于clang，而我们链接的NetCDF库是用GCC编译的。如果两个编译器的C++ ABI不一致（特别是 `_GLIBCXX_USE_CXX11_ABI`宏），链接时会出现 `undefined reference`错误。 在我们的环境中，DTK 25.04.4的 `hipcc`恰好跟GCC的ABI兼容，没有遇到这个问题。但如果你的环境不同，可能需要加 `-D_GLIBCXX_USE_CXX11_ABI=0`这样的编译选项。 ### 坑5：DCU温度墙与降频 DCU跟所有加速卡一样有TDP限制。我们在测试时发现，如果把kernel设计的计算密度太高（比如全部用寄存器、减少global memory访问），DCU温度会快速上升，触发降频保护，计算时间反而变长。 解决思路是适当增加一些memory-bound的操作（比如隔几个cycle读一次共享内存），让DCU的计算单元不要一直全速运转。这个优化方向比较玄学，没有精确公式，只能靠反复测试找平衡点。 --- 在下篇我会做一个系统性的总结：DCU和CPU在什么场景下各有什么优势、给DCU初学者的一份实战checklist、以及我对国产加速卡生态的整体感受和未来展望。