超级会员免费看
大规模并行系统的优化相对论代码：数据结构对并行效率的影响
1. 引言
在计算过程中，为有效使用内存而进行数据结构化仍然是一个非常有趣的问题。目前，英特尔和 AMD 这两家 CPU 芯片巨头都在生产支持高级矢量扩展（AVX）指令的处理器，AMD 还为 PC 芯片添加了 AVX - 512 指令，我们也期待着支持 AVX - 512 的服务器芯片尽快问世。
尽管现代编译器功能强大，但对于大多数现代科学代码而言，在没有进行数据结构化的情况下进行自动矢量化，其性能提升并不显著。这意味着单输入多数据（SIMD）编程的问题仍然存在，我们需要运用现代技术来解决数据结构化问题。
在过去十年里，我们基于流体动力学方法开发了天体物理代码，使用 Fortran 和 C/C++ 语言进行代码开发。在性能研究中，我们尝试通过在代码中添加相应的 C++ 内联函数来使用 AVX - 512 SIMD 指令，并基于不同矢量寄存器的数据对齐程序进行了各种代码优化。然而，使用底层 AVX - 512 内联函数会使代码与计算节点架构绑定，即使芯片都支持 AVX - 512，由于软件接口实现不同，基于英特尔 AVX - 512 内联函数的代码也无法在 AMD 系统上使用。
为了兼顾代码的互操作性和良好性能，我们可以采用一些方法，比如从数组的结构（AoS）到结构的数组（SoA）的代码优化，或者简化代码循环等。以下是 AoS 和 SoA 的代码示例：
1.1 AoS 循环示例
struct point3D {
float x;
float y;
float z;
};
struct point3D points[N];
float get_point_x(int i) { return points[i].x; }
1.2 SoA 循环示例
struct pointlist3D {
float x[N];
float y[N];
float z[N];
};
struct pointlist3D points;
float get_point_x(int i) { return points.x[i]; }
AoS 表示法更直观易懂，但对于 SIMD 代码，使用 SoA 数据表示法更为合适。在我们的工作中，我们尝试将英特尔的 SDLT（SIMD 数据布局模板）库应用于最新的基于 C++ 的天体物理代码 HydroBox3d。SDLT 是一个 C++11 模板库，它提供的容器可以使用能够生成高效 SIMD 矢量代码的布局来表示“普通旧数据”对象数组。
2. 数学模型和数值方法
我们在性能测试中使用了 HydroBox3D 代码的相对论流体动力学版本。下面简要介绍其数学模型和数值方法。
2.1 物理变量定义
首先定义物理变量：密度 $\rho$、速度矢量 $\mathbf{v}$、压力 $p$，并设光速 $c \equiv 1$。洛伦兹因子 $\Gamma$ 定义为：
$\Gamma = \frac{1}{\sqrt{1 - (\mathbf{v}/c)^2}} = \frac{1}{\sqrt{1 - \mathbf{v}^2}}$
2.2 状态方程
我们在代码中使用理想气体模型，状态方程以特殊焓 $h$ 的方程形式表示：
$h = 1 + \frac{\gamma}{\gamma - 1} \frac{p}{\rho}$
其中 $\gamma$ 是绝热指数。此时，声速 $c_s$ 定义为：
$c_s^2 = \frac{\gamma p}{\rho h}$
2.3 相对论流体动力学方程
为了写出特殊相对论流体动力学方程，我们引入守恒变量：相对论密度 $D = \Gamma\rho$，相对论冲量 $M_j = \Gamma^2\rho h v_j$（速度矢量 $\mathbf{v}$ 的分量为 $v_j$，$j = x, y, z$），以及全相对论能量 $E = \Gamma^2\rho h - p$。方程组可以写成如下形式：
$\frac{\partial}{\partial t} \begin{pmatrix} \Gamma\rho \ \Gamma^2\rho h v_j \ \Gamma^2\rho h - p \end{pmatrix} + \sum_{k = 1}^{3} \frac{\partial}{\partial x_k} \begin{pmatrix} \rho\Gamma v_k \ \rho h\Gamma^2 v_j v_k + p\delta_{jk} \ (\Gamma^2\rho h - p)v_k + p v_k \end{pmatrix} = 0$
其中 $\delta_{jk}$ 是克罗内克符号。逆变换需要求解非线性方程：
$f(p) = \Gamma^2\rho h - p - E = 0$
我们使用迭代牛顿法求解，公式为：
$p_{m + 1} = p_m - \frac{f(p_m)}{f’(p_m)}$
其中函数 $f$ 的导数形式为：
$f’(p) = \left(\frac{\gamma}{\gamma - 1}\Gamma^2 - \frac{M^2\Gamma^3}{(E + p)^3}\right) \left(D + 2\frac{\gamma}{\gamma - 1}p\Gamma\right) - 1$
当达到所需精度时，我们停止迭代过程。在我们的例子中，方程组可以重写为矢量形式：
$\frac{\partial U}{\partial t} + \sum_{k = 1}^{3} \frac{\partial F_k}{\partial x_k} = 0$
对于任意单元格，戈杜诺夫格式的形式为：
$\frac{U_{i + \frac{1}{2}, k + \frac{1}{2}, l + \frac{1}{2}}^{n + 1} - U_{i + \frac{1}{2}, k + \frac{1}{2}, l + \frac{1}{2}}^n}{\tau} + \frac{F_{x, i + 1, k + \frac{1}{2}, l + \frac{1}{2}}^
- F_{x, i, k + \frac{1}{2}, l + \frac{1}{2}}^
}{h_x} + \frac{F_{y, i + \frac{1}{2}, k + 1, l + \frac{1}{2}}^
- F_{y, i + \frac{1}{2}, k, l + \frac{1}{2}}^
}{h_y} + \frac{F_{z, i + \frac{1}{2}, k + \frac{1}{2}, l + 1}^
- F_{z, i + \frac{1}{2}, k + \frac{1}{2}, l}^
}{h_z} = 0$
其中 $h_{x,y,z}$ 是空间网格的步长，$F^*$ 是相应变量通过单元格边界的通量。这些值可以通过求解黎曼问题得到。数值求解器基于局部模板上的分段抛物线方法（PPML）。
2.4 计算流程
下面是 HydroBox3D 代码实现的数值方法的 mermaid 流程图：
graph TD;
A[定义物理变量] --> B[计算洛伦兹因子];
B --> C[根据状态方程计算特殊焓和声速];
C --> D[构建相对论流体动力学方程];
D --> E[求解逆变换的非线性方程];
E --> F[使用戈杜诺夫格式计算];
F --> G[求解黎曼问题得到通量];
3. CPU 效率评估
我们使用了西伯利亚超级计算机中心 NKS - 1P 系统的计算节点，每个节点配备了两颗英特尔至强可扩展第二代 6248R（24 核，3 GHz）CPU 和 192 GB 的 DRAM4 内存，同时也在配备两颗英特尔至强可扩展第二代 8268（24 核，2.9 GHz）CPU 和 192 GB 的 DRAM4 内存的节点上进行了测试。为了对比服务器 CPU 和工作站的性能，我们还在配备英特尔酷睿 i9 10980XE CPU 和 64 GB 的 DRAM4 的 PC 上进行了测试，所有这些 CPU 都支持 AVX - 512 指令。
我们在 NKS - 1P 系统和工作站 PC 上都安装了相同的英特尔 oneAPI 基础 + HPC 工具。在本次研究中，我们使用了 HydroBox3D 代码的 C++ 11 版本，并使用英特尔的 SDLT SIMD 模板对代码进行了修改，以自动生成高效代码。以下是欧拉阶段数据初始化和循环修改的代码对比：
3.1 原始 MPI 版本数据初始化及部分主循环
build_parabola(rho, parrho);
build_parabola(ux, parux);
build_parabola(uy, paruy);
build_parabola(uz, paruz);
build_parabola(p, parp);
for (i = 1; i < NX - 1; i++)
for (k = 1; k < NY - 1; k++)
for (l = 1; l < NZ - 1; l++)
{
// Riemann solver
SRHD_Lamberts_Riemann(
rho[index(i, k, l)],
ux[index(i, k, l)], uy[index(i, k, l)],
uz[index(i, k, l)],
p[index(i, k, l)], ...
...
3.2 优化后的 SDLT + OpenMP 版本数据初始化及部分主循环
build_parabola_sdlt(rho, parrho_sdlt);
build_parabola_sdlt(ux, parux_sdlt);
build_parabola_sdlt(uy, paruy_sdlt);
build_parabola_sdlt(uz, paruz_sdlt);
build_parabola_sdlt(p, parp_sdlt);
auto parrho_access = parrho_sdlt.const_access();
auto parux_access = parux_sdlt.const_access();
auto paruy_access = paruy_sdlt.const_access();
auto paruz_access = paruz_sdlt.const_access();
auto parp_access = parp_sdlt.const_access();
#pragma omp for
for (i = 1; i < NX - 1; i++)
{
for (k = 1; k < NY - 1; k++)
{
#pragma ivdep
for (l = 1; l < NZ - 1; l++)
{
// Riemann solver
#pragma forceinline recursive
SRHD_Lamberts_Riemann(
parrho_access[i][k][l].u(),
parux_access[i][k][l].u(),
paruy_access[i][k][l].u(),
paruz_access[i][k][l].u(),
parp_access[i][k][l].u(), ...
...
3.3 性能评估工具及操作步骤
我们使用了 oneAPI 中的英特尔 Advisor 工具进行性能评估，以下是编译和 Advisor 命令行的操作步骤：
1.
编译 HydroBox3D 代码
：
mpiicc -o srhd64 main.cpp -std=c++17 -g -debug inline -debug-info -qopenmp -xcascadelake -axcascadelake -prec-div - -qopt-report=5 -O3 -qopt-zmm-usage=high -fpermissive -no-ipo
- 运行英特尔 Advisor 进行调查分析 ：
mpirun -n 1 advisor --collect=survey --data-limit=0 --trace-mpi --project-dir=./advresults --./srhd64
- 运行英特尔 Advisor 进行行程计数分析 ：
mpirun -n 1 advisor --collect trip_counts --data-limit=0 --project-dir ./advresults --trace-mpi --flop --no-trip-counts --./srhd64
3.4 不同 CPU 的峰值性能数据
| CPU 型号 | 标量加法峰值 (GFLOPS) | 双精度矢量加法峰值 (GFLOPS) | 双精度矢量 FMA 峰值 (GFLOPS) | 单精度矢量 FMA 峰值 (GFLOPS) |
|---|---|---|---|---|
| 英特尔至强 6248R (两颗) | 342 | 1974 | 3949 | 7899 |
| 英特尔至强 8268 (两颗) | 333 | - | 3948 | 7897 |
| 英特尔酷睿 i9 10980XE | 134 | - | 875 | 1710 |
3.5 性能测试结果
3.5.1 英特尔至强 6248R CPU
- 原始 MPI 版本 ：网格大小对性能影响不大，但为了与工作站 PC 进行性能比较，我们使用了 128³ 的网格。大多数函数受 DRAM 和标量加法峰值性能的限制，只有一些抛物线构建过程超出了这些限制，这表明原始 MPI 版本的 HydroBox3D 代码在没有优化的情况下无法进行自动矢量化。所有函数的算术强度都小于 1 FLOPs/Byte，该版本代码的总性能仅为 1 GFLOPS，欧拉阶段主循环的性能为 0.21 GFLOPS。
- 优化后的 SDLT + OpenMP 版本 ：添加 OpenMP 编译指令和英特尔 SDLT 模板对数据结构进行优化后，性能得到了显著提升。大多数函数受 CPU 不同级别的缓存限制，一些函数达到了芯片矢量指令的峰值性能。欧拉阶段计算量最大的过程达到了 25 GFLOPS，代码总性能达到了 47 GFLOPS。
3.5.2 英特尔至强 8268 CPU
优化后的代码性能结果与 6248R CPU 相近，但 8268 CPU 的零售价格是 6248R 的两倍。
3.5.3 英特尔酷睿 i9 10980XE CPU
未优化和优化代码测试的峰值性能略有不同，可能是因为在测试中无法取消 BIOS 和 Windows 操作系统的所有硬件性能优化，SDLT + OpenMP 版本的代码可能在 CPU 和内存频率提升的情况下运行。欧拉阶段的性能从 1 GFLOPS 提升到了 22 GFLOPS，总性能从 1 GFLOPS 提升到了 16 GFLOPS。
3.6 性能对比总结
在工作站 PC 上的所有测试都启用了超线程（36 线程），而在 NKS - 1P 计算节点上的测试是在 48 核（两颗 CPU）且未启用超线程的情况下进行的。未优化版本的代码在英特尔酷睿 i9 10980XE 和两颗英特尔至强 6248R 上的运行情况相似，工作站 CPU 上欧拉阶段的自动矢量化效果更好。优化后，两颗服务器 CPU 的性能提升了 50 倍，工作站 CPU 的性能提升了 20 倍。我们可以推测，在 NKS - 1P 节点上启用超线程后性能将进一步提升，并且在未优化的科学代码情况下，工作站 CPU 有时可以达到与服务器 CPU 相同的性能。
4. 结论
目前，高性能超级计算机上的科学代码要实现最大性能，只能通过针对处理器矢量运算单元的特殊优化来达成。对于英特尔和 AMD 公司的最新 CPU 芯片，有必要使用 AVX - 512 指令集。现代 C++ 编译器能够针对处理器每个周期处理八个双精度值的 FMA 指令优化代码中的循环，但软件开发人员仍需协助编译器。主要问题在于创建结构数组形式的数据结构，并对数据进行 AVX - 512 寄存器对齐。
虽然基于类似汇编的内联函数编写的代码速度最快，但这种代码对于非计算机科学专业人员来说难以理解和进行后续修改。英特尔公司为此开发了特殊的 C++ 11 模板（英特尔 SDLT 库），以实现易于使用的结构数组软件设计方法。
我们在自主研发的天体物理 HydroBox3D 代码上对该库进行了测试。在配备英特尔至强可扩展第二代 CPU 的高性能计算系统的计算节点上，我们取得了最佳的性能提升效果。优化后的 OpenMP + SDLT 版本代码的性能比简单的 MPI 版本程序高出五十倍。
值得注意的是，基于支持 AVX - 512 的芯片构建的现代工作站 PC 系统，在运行完全优化的科学程序时也能表现出出色的性能，这得益于其较高的频率和相当的核心数量。目前，我们正在等待支持 AVX - 512 的最新 AMD 服务器芯片，以便进行对比测试。
4.1 优化总结
以下是本次优化的关键要点总结：
1.
数据结构优化
：从数组的结构（AoS）转换为结构的数组（SoA），并使用英特尔 SDLT 库来辅助数据布局，有助于编译器生成高效的 SIMD 矢量代码。
2.
代码优化
：添加 OpenMP 编译指令实现并行化，提高代码的执行效率。
3.
性能评估
：使用英特尔 Advisor 工具进行性能分析，通过 Roofline 分析图表直观地了解代码各函数的性能与 CPU 和内存峰值性能的对比情况。
4.2 性能提升效果对比
| 代码版本 | 英特尔至强 6248R CPU 总性能 (GFLOPS) | 英特尔酷睿 i9 10980XE CPU 总性能 (GFLOPS) |
|---|---|---|
| 原始 MPI 版本 | 1 | 1 |
| SDLT + OpenMP 版本 | 47 | 16 |
4.3 未来展望
- 硬件适配 ：等待 AMD 支持 AVX - 512 的服务器芯片推出后，进行进一步的对比测试，以评估代码在不同硬件平台上的性能表现。
- 代码优化 ：继续探索其他优化方法，进一步提高代码的性能，例如针对不同的计算任务进行更细致的优化。
- 应用拓展 ：将优化后的代码应用于更广泛的天体物理问题，解决更复杂的科学计算任务。
4.4 优化流程总结 mermaid 图
graph TD;
A[选择合适的数据结构] --> B[使用英特尔 SDLT 库];
B --> C[添加 OpenMP 编译指令];
C --> D[使用英特尔 Advisor 工具进行性能评估];
D --> E[根据评估结果进行代码优化];
E --> F[获得性能提升的代码];
综上所述，通过合理的数据结构优化和使用英特尔 SDLT 库，我们成功地提升了 HydroBox3D 代码在大规模并行系统上的性能。未来，随着硬件技术的不断发展和优化方法的持续改进，我们有望在天体物理计算领域取得更加优异的成果。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
