29、地球物理代码在多核计算架构上的效率优化研究

地球物理代码在多核计算架构上的效率优化研究

1. 引言

在解决实际数值模拟问题时，需要开发能高效利用现代计算架构的并行程序。如今，不仅有多种数学方法和算法可供选择，计算架构也多种多样。因此，研究人员需综合考虑可用计算架构的能力，确定最佳选择。

近年来，超级计算机架构的全球趋势是采用带有计算加速器、多核处理器和协处理器的异构系统。这些系统性能强大且节能，例如俄罗斯科学院西伯利亚分院计算数学与数学地球物理研究所的西伯利亚超级计算机中心就安装了相关系统，如配备 NVIDIA 图形加速器的 NKS - 30T + GPU 集群和配备 Intel Xeon Phi 处理器的 NKS - 1P 集群，总峰值性能达 167 TFLOPS。开发适用于此类混合系统的高效软件虽需额外知识和时间，但能显著提升性能。

为简化数值建模问题中计算架构的选择和高效代码的开发，可使用基于本体方法开发的智能支持系统。该系统基于协同设计方法，从高效并行计算的角度出发，开发物理和数学模型、数值方法、并行算法及其软件实现。

2. 智能支持系统

智能支持系统用于解决超级计算机上计算密集型的数学物理问题，主要由以下部分组成：
- 知识库 ：包含计算方法和算法的本体、并行架构和技术的本体以及推理规则。用户可通过信息分析互联网资源研究知识库中的对象、查看它们之间的联系，并补充新对象。
- 专家系统 ：用户将待解决问题的规格提交给专家系统，推理引擎根据知识库中的本体对象和专家制定的推理规则构建问题解决方案。然后，利用软件组件库中的现有模块，结合计算算法和所选计算系统的架构生成并行代码。
- 仿真模块 ：帮助用户确定解决问题的最佳计算核心数量。

推理引擎用于处理本体模型，可检查本体的一致性，操作类、属性和实体的名称，并根据专家规则推断本体中未明确包含的信息。

3. 问题陈述与数值方法

模拟复杂三维非均匀介质中弹性波的传播，需要有效的并行和缩放算法。这一类问题包括研究岩浆火山特征介质中地震波的传播特性。

地震波在复杂弹性非均匀介质中的传播可通过求解弹性理论方程组及相应的初始和边界条件来描述。以位移向量 $\vec{U} = (U, V, W)^T$ 表示：
[
\rho \frac{\partial^2 \vec{U}}{\partial t^2} =
\begin{bmatrix}
A_1 \
A_2 \
A_3
\end{bmatrix}
\vec{U} + \vec{F}(t, x, y, z)
]
其中：
[
A_1 =
\left(
(\lambda + 2\mu) \frac{\partial^2}{\partial x^2} + \mu \left( \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} \right)
\right)
(\lambda + \mu) \frac{\partial^2}{\partial x \partial y}
(\lambda + \mu) \frac{\partial^2}{\partial x \partial z}
]
[
A_2 =
(\lambda + \mu) \frac{\partial^2}{\partial y \partial x}
\left(
(\lambda + 2\mu) \frac{\partial^2}{\partial y^2} + \mu \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial z^2} \right)
\right)
(\lambda + \mu) \frac{\partial^2}{\partial y \partial z}
]
[
A_3 =
(\lambda + \mu) \frac{\partial^2}{\partial z \partial x}
(\lambda + \mu) \frac{\partial^2}{\partial z \partial y}
\left(
(\lambda + 2\mu) \frac{\partial^2}{\partial y^2} + \mu \left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right)
\right)
]
这里，$t$ 是时间，$\rho(x, y, z)$ 是密度，$\lambda(x, y, z)$ 和 $\mu(x, y, z)$ 是拉梅系数。模拟的是一个各向同性的三维非均匀弹性复杂介质，其为一个平行六面体区域，其中一个面为自由表面（平面 $z = 0$），自由表面上的初始和边界条件为零。

在协同设计框架下，对比位移和应力速度的公式，发现位移计算具有优势，速度更快且节省大量内存，这在解决大规模问题时尤为重要。

对于三维弹性理论的动态问题，有限差分法是最灵活且广泛使用的方法。它允许使用规则的平行六面体网格，且不损失可接受的精度，便于通过计算域分解和重叠子域的方法并行化算法。显式有限差分格式适合图形加速器架构，因为它们构建的计算网格可直接映射到 GPU 架构，且每个计算单元的计算相互独立。该方法也适用于向量化。

采用交错网格上的显式有限差分格式进行数值求解，网格系数的计算基于积分守恒定律，该格式在时间和空间上具有二阶近似。计算网格上每个节点的位移仅需相邻节点的值。

4. 多核系统

研究的计算系统涵盖了多种多核但架构差异较大的系统。一方面是经典处理器，如具有向量扩展和不同核心数的 Intel 处理器，以及具有高同步多线程（SMT）的 IBM 处理器；另一方面是拥有数百甚至数千个轻量级核心和复杂快速内存系统的 GPU。

Intel Broadwell 和 IBM POWER9 架构的优势在于它们具有较大的 L3 缓存，分别为每芯片 40 MB 和 120 MB，比 MCDRAM KNL 快得多。但 KNL 可使用多达 72 个核心（Intel Broadwell 为 16 个核心，IBM POWER9 为 12 个核心），每个核心每周期可处理 32 个双精度数（Broadwell 和 POWER9 为 16 个），因此使用 Intel KNL 进行计算向量化有望显著提高生产力。

以下是各多核系统的主要参数：
| 主架构 | 参数 | 值 |
| ---- | ---- | ---- |
| Intel Broadwell | 处理器 | 2 × Intel Xeon E5 - 2697A v4，2.6 GHz，16 核，SMT2 |
| | 内存 | 128 GB DDR4 RAM |
| | 峰值性能 | 1331.2 GFLOPS |
| Intel KNL | 处理器 | Intel Xeon Phi 7290 KNL，1.5 GHz，72 核，SMT4 |
| | 内存 | 16 GB MCDRAM，96 GB DDR4 RAM |
| | 峰值性能 | 3456 GFLOPS |
| IBM POWER9 | 处理器 | IBM POWER9 处理器，3.8 GHz，2 × 12 核，SMT8 |
| | 内存 | 32x32 GB DDR4 RAM |
| | 峰值性能 | 2918.4 GFLOPS |
| NVIDIA Fermi | 处理器 | 2 × Intel Xeon X5670 |
| | 内存 | 96 GB DDR4 RAM |
| | 加速器 | NVIDIA Tesla M2090，1300 MHz，512 核，6 GB GDDR5 |
| | GPU 峰值性能 | 1331 GFLOPS |
| NVIDIA Kepler | 处理器 | 2 × Intel Xeon E5 - 2650 v2 |
| | 内存 | 64 GB DDR4 RAM |
| | 加速器 | NVIDIA Tesla K40，745 MHz，2880 核，12 GB GDDR5 |
| | GPU 峰值性能 | 4291 GFLOPS |
| NVIDIA Pascal | 处理器 | IBM POWER8 |
| | 内存 | 256 GB DDR4 RAM |
| | 加速器 | NVIDIA Tesla P100，1480 MHz，3584 核，16 GB HBM2 |
| | GPU 峰值性能 | 10608 GFLOPS |

除基于 NVIDIA Fermi 的系统外，其他架构的测试计算均在 581 × 581 × 581 的网格上进行（总内存约 11 Gb）。对于 NVIDIA Fermi 系统，使用 581×581×193 的较小网格（总内存约 3.6 Gb）以适应加速器内存。效率和性能的计算基于执行时间。

4.1 Intel 和 IBM 处理器的优化

对于配备多核处理器的系统，这些处理器核心数量相对较少（几十个），时钟频率高且相互独立。现代处理器的特点是具有向量处理单元，可在每个时钟周期处理多个浮点数。例如，KNL 支持每周期处理 512 位，Broadwell 为 256 位，POWER9 为 128 位。同时，这些系统都支持同步多线程，Broadwell 每个核心有 2 个硬件流，KNL 有 4 个，POWER9 有 8 个。

为并行化弹性介质中地震波传播的三维建模算法，对空间内循环进行向量化，使用 OpenMP 对空间外循环进行并行化，并对齐所有主三维数组以实现有效内存访问。
- 向量化和缓存 ：对于 Intel 处理器，比较了两种向量化技术：使用编译器额外指令进行无数据依赖的自动向量化，以及使用内在函数的底层技术。结果表明性能无显著差异，说明 Intel 编译器对初始代码的自动向量化效果良好。一般来说，使用向量运算可使程序加速数倍，如 KNL 在 72 个 OpenMP 流的缓存内存模式下可加速 2.75 倍。
研究不同循环顺序对缓存效率和缓存请求数量的影响，发现循环顺序不当会使内存访问速度减慢数十倍，最佳循环顺序（zyx）与三维数组组件的存储顺序一致。POWER9 和 KNL 的缓存比 Broadwell 更大更快，yzx 顺序也能获得与 zyx 相近的性能提升，POWER9 的 yzx 和 zyx 顺序分别可提升 15.8 倍和 17.2 倍。KNL 中连续内存值可实现最大提升，yzx 顺序为 19.4 倍，zyx 顺序为 26.9 倍。
- 负载均衡 ：由于许多应用程序的计算负载分布不均，为减少计算时间，需平衡线程间的任务。可通过修改计算算法和配置 OpenMP 并行化参数来实现。系统的多核程度越高，负载均衡越重要。
在测试代码中，将两个外循环合并为一个，内循环选择 X 方向，这种优化可显著增加外循环的迭代次数，减少每次迭代的工作量，从而实现更好的负载平衡。与最佳的 zyx 循环顺序相比，Broadwell 在 32 个线程下可加速 1.48 倍，KNL 在 72 个线程下可加速 1.18 倍，POWER9 在 192 个线程下可加速 1.16 倍。
还可通过选择 OpenMP schedule 指令的参数进一步优化负载均衡。研究发现，对于 Broadwell 和 POWER9 架构，默认参数（静态，块大小等于循环迭代次数除以线程数）性能最佳；对于 KNL，引导调度器可带来 1.2 倍的小幅度提升；对于 POWER9 架构，调度器类型对测试代码的负载均衡影响不大。
- 可扩展性和多线程 ：图 7 展示了在各架构最佳代码设置下，多核处理器相对于单线程的强可扩展性。考虑 SMT 时，Broadwell 和 KNL 系统在线程数超过物理核心数时可扩展性下降，使用超线程效率降低。Broadwell 系统在 32 个线程下最大加速比为 8.8 倍，KNL 在 72 个线程下为 56.6 倍。POWER9 充分发挥了超线程能力，在 192 个线程（24 个物理核心）的测试代码中实现了 48 倍的加速。原因之一是向量代码和向量寄存器越小，SMT 的效果越明显。同时，具有更先进缓存系统的架构（KNL 和 POWER9）由于测试代码受内存限制，可实现更好的加速。
- Intel KNL 的扁平内存模式 ：KNL 是全功能处理器，无需像其前身 KNC 协处理器那样进行预处理和数据复制。为优化内存使用，仅使用 16 GB 的快速 MCDRAM 内存，可通过多个加速器扩展容量。在本测试中，使用扁平内存模式将所有主数组放置在 MCDRAM 内存中，比缓存模式加速 1.3 倍，混合模式无优势。尽管 Intel 停止了 Xeon Phi 系列的进一步开发，但其部分架构创新已应用于新的 Intel 处理器。

4.2 NVIDIA 加速器的优化

GPU 拥有大量（数百或数千个）高度简化的计算核心，运行频率低，且配备复杂的内存系统，部分由软件控制。为高效使用 GPU，需使用比核心数量多得多的线程。

以下是几种 GPU 的主要特性对比：
| 显卡型号 | 核心数 | 内存大小（GB） | L2（KB） | 共享内存（KB） | 内存带宽（GB/s） |
| ---- | ---- | ---- | ---- | ---- | ---- |
| Tesla M2090 | 512 | 6 | 768 | 16/48 | 177 |
| Tesla K40 | 2880 | 12 | 1536 | 16/48 | 288 |
| Tesla P100 | 3584 | 16 | 4096 | 64 | 733 |

为减少主机和设备之间的传输，测试代码中的所有计算都在 GPU 上进行。首先在 CPU 上生成差分方案各点的介质系数数组，并复制到加速器内存中。计算过程中，在特定时间步将几个平面的波场图像传输回主机（平均每 500 - 10000 步传输一次图像即可）。这种方法虽将问题规模限制在设备内存大小范围内，但可避免因多次复制数据而造成的损失。后续可通过使用多个带图形加速器的节点来扩大问题规模。
- 线程块大小 ：在主 CUDA 核心中，每个线程计算一个差分网格中当前时间步的向量分量。增加线程的计算负载（如一个线程计算多个点或沿一个轴进行整个循环）会增加内核执行时间。例如，在 Tesla M2090 卡的 512 个线程上处理 360³ 的问题时，每个线程仅计算计算网格上一个点所需分量的 3D 网格版本比每个线程沿 z 轴循环的 2D 网格版本快近 2 倍。
研究了软件性能与线程块维度和大小的关系，将每个 GPU 的计算域划分为三维块网格，x 分量的块大小必须是 warp 长度的倍数，具体块大小需根据每个算法通过经验选择。每个 GPU 的最大线程数受其计算能力限制（M2090 为 512，K40 和 P100 为 1024）。以下是部分测量结果：

Tesla M2090	BSize	PGain	BSize	PGain	BSize	PGain	BSize	PGain
	512		256		128		64
	8×8×8	0.52	16×8×4	1.0	16×4×4	1.05	16×4×2	0.96
	16×2×2	0.99	32×4×4	1.47	32×4×2	1.49	32×2×2	1.48
	32×2×1	1.26	64×4×2	1.42	64×2×2	1.65	64×2×1	1.38
	64×1×1	1.35	128×2×2	1.45	128×2×1	1.44	128×1×1	1.41
	256×2×1	1.2	256×1×1	1.41	512×1×1	0.98

Tesla K40	BSize	PGain	BSize	PGain	BSize	PGain	BSize	PGain
	1024		512		256		128
	16×8×8	1.0	16×8×4	1.11	16×4×4	1.04	16×4×2	1.02
	32×8×4	1.53	32×4×4	1.67	32×4×2	1.63	32×2×2	1.55
	64×4×4	1.64	64×4×2	1.7	64×2×2	1.82	64×2×1	1.72
	128×4×2	1.63	128×2×2	1.7	128×2×1	1.78	128×1×1	1.73
	256×2×2	1.44	256×2×1	1.58	256×1×1	1.69	512×2×1	1.13
	512×1×1	1.43

Tesla P100	BSize	PGain	BSize	PGain	BSize	PGain	BSize	PGain
	1024		512		256		128
	16×8×8	1.0	16×8×4	1.02	16×4×4	1.02	16×4×2	0.98
	32×8×4	1.25	32×4×4	1.29	32×4×2	1.23	32×2×2	1.23
	4×64×4	0.49	32×8×2	1.23	4×4×64	0.48	32×16×1	1.12
	64×4×4	1.24	64×4×2	1.26	64×2×2	1.25	64×2×1	1.08
	128×4×2	1.25	128×2×2	1.27	128×2×1	1.12	128×1×1	1.08
	256×2×2	1.22	256×2×1	1.14	256×1×1	1.13	512×2×1	1.1
	512×1×1	1.14

x 分量为 32（warp 大小）或 64（warp 大小 × 2）且在其他分量上均匀分布的块大小最为有利。最大可能的块大小并非最佳选择，因为 GPU 寄存器在块线程间分配，可能会不足。这种优化在代码修改方面并不复杂，却能显著提高性能。

地球物理代码在多核计算架构上的效率优化研究

5. 优化总结与规则制定

通过对不同多核系统的研究和优化实践，总结出以下针对不同架构的优化规则，这些规则有助于选择合适的并行化和优化方法，提高地球物理代码的性能。

5.1 Intel 和 IBM 处理器

• 向量化 ：使用向量运算可显著加速程序，Intel 处理器的自动向量化效果良好，可直接使用编译器的相关指令。对于 KNL，在 72 个 OpenMP 流的缓存内存模式下，向量化加速可达 2.75 倍。
• 循环顺序 ：选择与三维数组存储顺序一致的循环顺序（zyx）可提高缓存效率，减少内存访问延迟。POWER9 和 KNL 因缓存优势，yzx 顺序也能获得较好性能提升。
• 负载均衡 ：对于计算负载不均的应用，可通过合并外循环、选择合适的 OpenMP schedule 指令参数来平衡线程任务。如测试代码中，(zy)x 循环顺序结合负载均衡优化，Broadwell 在 32 个线程下比 zyx 顺序加速 1.48 倍。
• 可扩展性 ：Broadwell 和 KNL 在线程数超过物理核心数时，超线程效率降低；POWER9 可充分发挥超线程能力。具有先进缓存系统的架构（KNL 和 POWER9）在内存受限的代码中表现更佳。
• 内存模式 ：Intel KNL 使用扁平内存模式，将主数组放置在 MCDRAM 内存中，可加速 1.3 倍。

5.2 NVIDIA 加速器

• 计算位置 ：将所有计算放在 GPU 上进行，减少主机和设备之间的数据传输。先在 CPU 生成介质系数数组并复制到加速器内存，计算过程中按需将波场图像传输回主机。
• 线程块大小 ：选择 x 分量为 32 或 64 且在其他分量上均匀分布的线程块大小，避免使用最大可能的块大小，以提高性能。

6. 智能支持系统的应用

智能支持系统可简化计算架构的选择和高效代码的开发过程，其工作流程如下：

graph LR
    A[用户提交问题规格] --> B[专家系统]
    B --> C{推理引擎}
    C --> D[使用知识库和推理规则构建解决方案]
    D --> E[生成并行代码]
    E --> F[仿真模块确定最佳核心数]

1. 用户将待解决问题的规格提交给专家系统。
2. 推理引擎根据知识库中的本体对象（计算方法、并行架构等）和专家制定的推理规则，构建问题解决方案。
3. 利用软件组件库中的现有模块，结合计算算法和所选计算系统的架构，生成并行代码。
4. 仿真模块帮助用户确定解决问题的最佳计算核心数量。

7. 结论

通过对多种多核计算架构（Intel Broadwell、Knights Landing、IBM POWER9、NVIDIA Fermi、Kepler、Pascal）的研究和优化，得出以下结论：
- 不同架构具有不同的特点和优势，针对特定架构进行优化可显著提高地球物理代码的性能。
- 智能支持系统基于本体方法和协同设计理念，能有效帮助用户选择合适的计算架构和优化方法，简化开发过程。
- 制定的优化规则可用于指导不同架构下并行算法和程序的开发，提高计算效率。

未来，随着计算技术的不断发展，多核架构将更加多样化和复杂。进一步研究和优化地球物理代码在新架构上的性能，以及完善智能支持系统的功能，将是后续工作的重点。同时，可探索如何将这些优化方法应用到其他领域的计算密集型问题中，提高整体计算效率。

总之，通过对地球物理代码在多核计算架构上的效率优化研究，为解决复杂的数值模拟问题提供了有效的方法和工具，有助于推动相关领域的发展。