在微观粒子世界中起主要作用的是强相互作用、电磁相互作用和弱相互作用, 描述强相互作用的量子色动力学(QCD)和描述电磁、弱相互作用的电弱统一理论(EW)统称为粒子物理中的标准模型. 就强相互作用而言, 其低能特性如夸克禁闭问题, 仍是需要解决的世纪难题. 实际上, 对于夸克禁闭以及与之密切相关的强子结构和强子性质等低能区强相互作用现象需要有效的非微扰方法. 格点QCD^[1]是目前公认有效的从第一原理出发的求解QCD的非微扰方法, 除标准模型基本参数之外没有任何额外模型参数, 因而其计算结果被认为是对强相互作用现象的可靠描述. 该格点QCD理论对不少重大物理研究有着不可替代的作用, 主要体现在用高能物理实验结果对现有理论的验证、用理论分析结果指导实验设计、以及进一步地寻找标准模型以外的新物理等方面.

格点QCD的研究手段是通过计算机进行大规模的矩阵运算和蒙特卡罗数值模拟, 其算法特点具有内在的高并行度和高可扩展性. 然而, 随着计算机硬件体系结构的日益复杂, 研发计算高效的格点QCD应用程序也变得愈加挑战性. 大规模格点QCD计算过程对高性能计算资源需求巨大, 其计算过程体现出计算密集型特征, 而其数据分析体现出数据密集型特征. 计算精度和数据分析效率直接反映了格点QCD计算方法的应用水平. 格点QCD已经逐渐成为与高能物理实验研究和理论解析研究并列的高能物理第三大分支, 研究范围已经拓展到强相互作用现象的各个方面. 美国、日本、欧洲等国家无一例外将格点QCD应用作为高性能计算的重要应用. 美国更是将其列为未来E级计算重要的应用之一, 凸显了格点QCD数值模拟内在的科学计算应用的价值^[2]. 近年来, 格点QCD的数值模拟方法发展出了大量理论和算法, 形成了美国USQCD^[3]和MILC^[4]、英国Grid^[5]和日本Bridge++软件^[6]等诸多优秀的格点QCD模拟软件.

从格点QCD数值计算的特点来分析, 其主要过程为矩阵向量间的运算, 并行化时涉及格点体系的区域分解并行划分, 边界格子进程间通信还有组态存储分析等任务^[7]. 如何实现高效的大规模格点QCD数值计算一直以来都是一个现实问题. 但我们研究发现, 格点QCD抽象分层软件架构的设计方法得到了重要应用, 软件框架的这种分层抽象的设计方法极大地提高了代码的开发效率和运行效率. 比如美国主导的USQCD软件框架分为4层: 最上层应用层包括Chroma、CPS和MILC等, 下面3、4层为支撑库层, 分别实现消息通讯、数据结构定义和基本算法实现, 而其中的QDP++支撑库是用来实现数据并行应用接口的, 提供逻辑上的数据并行操作, 隐藏底层数据依赖体系结构的操作, 提供上层统一的数据通信函数接口. 同样在这方面, 英国Peter Boyle等人2015年研发出了Grid软件^[5], 它抽象出了一个数据并行计算模式, 针对格点QCD的热点计算抽象了统一的数据并行接口, 支持多种平台的SIMD矢量化计算, 如AVX2、AVX512、NEON等形式, 同时也支持组合OpenMP和MPI等多种并行计算编程模式. 在实现线程化、SIMD操作可以实现比QDP++更有效.

Grid软件包近年来受到越来越多的应用, 利用它开展的物理研究包括K介子衰变中直接电荷共轭-宇称(CP)联合破坏参数的格点QCD计算^[8], 这能帮助物理学家理解宇宙中观测到的正反物质不对称性. 文献[8] 利用Grid软件改进了他们以前的物理测量代码, 性能得到显著提升. 另外核子结构的研究^[9]以及超出标准模型新物理的研究^[10]等也使用了Grid软件. 文献[11]中作者们初步尝试将Grid软件移植到GPU上, 他们比较了几种策略: OpenACC、OpenMP 4.5、Just-In-Time编译和CUDA并讨论了各自的优缺点.

从现有工作来看, 格点QCD计算软件是以分层软件架构为主, 以数据并行计算模式为核心, 如何提升格点QCD中数据并行计算效率和传输效率成为关键点, 这是现有 USQCD和Grid 等软件包重要的研发目标. 本文基于格点QCD计算Grid软件包, 开展格点QCD数据并行计算的理论分析, 研究数据并行计算模式在实现上的关键问题, 分析典型应用软件Grid和QDP++的软件架构实现, 最后结合具体的计算平台, 测试相关的应用软件性能, 针对实验数据开展分析, 验证相应的理论分析的有效性.

1 格点QCD模拟数据并行计算分析 1.1 格点QCD计算算法原理

格点QCD数值计算的基本思想为: 1)将四维连续时空离散化为一个四维超立方时空格子L_x×L_y×L_z×L_t, 这样体系的自由度变为有限可数; 2)在该时空格子上定义基本自由度夸克场和胶子场, 并在此基础上定义QCD 的格点作用量形式; 3)采用路径积分量子化的方式; 4)物理观测量的测量类似统计物理中的系综平均值. 于是问题就转化为一个包含近邻相互作用(因果律要求物理体系的相互作用必须是局域的)、多自由度体系的统计物理问题, 可以采用蒙特卡罗数值模拟方法对物理观测量进行计算. 格点QCD数值模拟计算的简单流程如下:

1)按照Boltzmann分布因子进行蒙特卡罗随机抽样, 产生QCD物质状态系综(规范场组态), 这个过程中会涉及大量的线性系统求解 ${\bf{Dx}} = {\bf{b}}$ 来用于组态更新, 其中涉及到维数达到 ${L_x} \times {L_y} \times {L_z} \times {L_t} \times 4 \times 3$ (一般超过千万数量级)的稀疏矩阵 ${D^{f,ab}_{\alpha \beta }}{(n|m)_{\begin{array}{*{20}{c}}\end{array}}}$ (以典型的威尔逊格点作用量为例):

$\begin{split} {D^{f,ab}_{\alpha \beta }}{(n|m)_{\begin{array}{*{20}{c}} \end{array}}} =& \left( {{m^{(f)}} + \frac{4}{a}} \right){\delta _{\alpha \beta }}{\delta _{ab}}{\delta _{n,m}} \\ &-\frac{1}{{2a}}\displaystyle\sum\limits_{\mu = \pm 1}^{ \pm 4} {{{\left( {{\bf{1}} - {\gamma _\mu }} \right)}_{\alpha \beta }}{U_\mu }{{(n)}_{ab}}{\delta _{n + \hat \mu ,m}}} \end{split}$

(1)

这里n和m的取值为整个格子 ${L_x} \times {L_y} \times {L_z} \times {L_t}$ , α和β为夸克自旋指标分别有4个取值, 下标a和b为夸克的颜色自由度分别有3个取值.

2)物理量测量, 这个过程涉及在第1)步生成的每一个组态上做大量的线性系统求解 ${\bf{Dx}} = {\bf{b}}$ . 与1) 的不同在于向量b的设置不一样以及解向量x这里是用于物理量的计算.

3)对物理量做组态平均, 通过数据拟合与分析获得物理结果.

总体而言, 格点QCD的主要计算热点就是求解大规模线性方程组 ${\bf{Dx}} = {\bf{b}}$ . 人们一般采用CG、BiCGstab、MinRes、GCR等迭代算法, 并且配合Even-odd、Multi-grid等适当的预处理算法, 在{b, Db, D2b, D3b,···} 张开的Krylov子空间里迭代求解所需的大规模线性方程组 ${\bf{Dx}} = {\bf{b}}$ . 在求解过程中, 这些算法会涉及大规模数据处理操作, 因而“矩阵D乘以向量”操作是计算热点的热点. 另外, 在超大规模计算中, “计算两个向量内积”由于涉及到全局归约, 也会产生额外开销.

1.2 格点QCD计算的数据并行计算特征

格点QCD数值模拟时, 四维超立方格子上的点被分配到多个处理核做并行计算, 三维中的示意图如图1所示.

由于QCD具有近邻相互作用特征, 每个处理器核心在涉及小立方边界处就会向其他相关核心发起通信以获得相邻立方边界的数据. 因此切分格子的方法会对计算效率产生明显的影响, 需考虑计算与通信的平衡.

稀疏矩阵 ${D^{f,ab}_{\alpha \beta }}{(n|m)_{\begin{array}{*{20}{c}}\end{array}}}$ 中包含的U矩阵为颜色空间3×3的复矩阵, U和颜色空间向量的乘积涉及到矩阵元相乘然后求和. 单个这样的矩阵乘向量操作通常不适合直接做向量化, 因为矩阵元的个数和寄存器的宽度通常不会恰好配合, 而且最后的求和操作无法向量化. 在Grid软件中, 多个这样的矩阵乘向量操作被组合在一起(数目根据硬件而定), 不同矩阵的矩阵元与向量元的乘积交织在一块并行, 之后的多个求和也可以同时进行, 从而实现高效的向量化计算.

1.3 数据并行计算的典型设计

Grid 最重要也是其突出的特点就体现在数据处理的矢量化设计上. 前面提到过, 在格点QCD计算中, 为了能够适应多节点并行计算, 通常是将格子在时空维度进行划分, 然后在单个节点上进行各种矩阵的计算, 这是很自然的想法. 在当时处理器没有向量化或向量化能力比较弱的情况下, 这种处理方法没有什么问题的. 早期的格点QCD程序如QDP++ 在考虑处理器矢量化时, 一般考虑的就是在单个矩阵计算时的矢量化(如图2, 其中等号左边列向量中实点的值等于等号右边矩阵中实点所在行与列向量的内积(图引自文献[5]), 其并行方法是在单个矩阵计算内, 将多个元素的内积及规约进行并行), 这需要 MPI进程之间进行很多的归约, 形成通讯上的瓶颈, 对性能造成影响.

但是数据的划分方式除了时空划分这种方式, 还可以通过每个格点上的自旋和洛仑兹指标进行划分, Grid 在设计时, 就部分采用了这种划分方式, 采用了与矩阵内矢量化不同的方式, 如图3, 其考虑两个或4个矩阵一起对每个元素进行计算, 在 SIMD 线上, 没有 MPI归约, 减少了通讯开销. 图3中, SIMD可实现多个矩阵向量乘法并行操作, 多个向量对等位上的数存储在一个SIMD向量去计算, 避免了SIMD位长限制和规约操作(图3引自文献[5]).

这种方式很适合格点 QCD计算, 考虑到 AVX512 有8个 double 宽度, 由于要计算的矩阵有4维旋量指标, 自然就适合4个矩阵一起进行计算, 而且矢量化也非常简便. Grid软件中采用的这种SIMD计算方式充分发挥了数据并行计算的潜在能力. 此外, 这种SIMD 矢量化方式对各种数值运算都进行了矢量化抽象, 只需要完成相应的数据矢量化代码, Grid就能够很容易地扩展到新的硬件架构中. 我们在开发格点QCD代码时, 也可以借鉴 Grid 优秀的设计思想.

2 实验结果与分析 2.1 AVX512与AVX2性能对比

选择Grid软件Wilson实例作为测试对象, 我们比较了在Intel Xeon Gold 5120 机器在采用AVX2和AVX512向量化版本的性能差异, 性能结果数据如表1和图4.

从图4中, 我们可以看出浮点计算性能在AVX512版本下优于AVX2, 但随着线程数(Threads)增加, 越来越接近总核数(28), 会发现其性能对比差异开始下降.

2.2 OpenMP共享内存测试

接下来分析线程绑定对数据并行计算的影响. 我们选择在2路2核的Intel Xeon Silver 4114机器上做线程绑定测试, 采用Intel编译器来编译Grid软件, 通过设置环境变量KMP_AFFINITY实现线程绑定, 主要分为两种模式, 即compcat和scatter. compact模式下, 线程将优先集中在一个CPU上进行计算; scatter模式下, 线程将均匀分布在多个CPU上进行计算, 访存变大.

具体测试方案: 分别在scatter模式下和compact模式下进行线程绑定的测试, 分别测试线程数为1, 2, 3, 4, 8, 10, 12, 16, 20, 30, 32, 36和40下的计算性能, 得到结果如表2, 同时图示化对比结果如图5所示.

由于程序存在大量的cache访问, 在线程数分别为2, 3, 4, 8, 10, 12, 16和20的情况下, 线程绑定导致共享相同cache, 造成冲突, 竞争使用cache, 造成cache频繁读写, 降低访存速度, 相比分散线程, 性能要差. 如果应用程序需要大量的访存操作, 需要大量的cache, 将OpenMP线程绑定核时, 尽量将线程分布在每个CPU, 可以很好的提升性能, 在本次测试中, 性能较compact下最高提升1.7倍. 如果程序的每个线程需要共享大量数据, 则将所有线程绑定在一个CPU时, 性能更好.

2.3 ARM 架构向量化性能对比

如今越来越多的服务器开始采用ARM架构的处理器, 比如华为的鲲鹏处理器, 亚马逊的Graviton处理器, 飞腾的FT系列处理器等等, “天河三号”原型系统便使用了 64 位 ARM 架构的64核心处理器 FT-2000plus和Matrix-2000加速卡, 其与 Intel 的 Xeon 架构类似. ARM架构的CPU向量化操作 (asimd) 有NEON技术. 我们在“天河三号”原型系统上, 测试了向量化前后Grid的性能对比情况.

GCC 编译器可以通过-march选项设置相关参数来启用向量化编译. 首先我们测试在没有开启向量化编译情况下Grid 的性能. 在每个节点上, 我们使用 1 个 MPI进程和64个 OpenMP 线程进行测试, 每个节点负责的子格子大小是固定的16⁴. 我们记录程序的计算性能和计算时间, 如表3所示, 其中, 第2列和第3列是计算 Dslash 和 Wilson fermion 时单个节点的总计算性能(Dslash 是格点QCD中的最重要的算符, Wilson fermion 是格点QCD中的一种费米子), 第4列和第5列则是计算所消耗的总时间和通讯时间.

在开启向量化(asimd)后, Grid的性能测试结果如表4, 各列含义与表3相同.

作为对比, 结合表3和表4, 在图6中画出了表示性能的Wilson列数值随节点数变化的曲线图, 其中w/o asmid代表没有开启asimd; 而w asimd表示开启了asimd.

从图6可以看出, 一方面, 随着节点数增加, 单个节点计算能力慢慢下降, 但总计算能力呈线性增加, 这表明Grid 有很好的弱扩展性; 另一方面, 目前我们只是采取了基本的自动化向量化编译, 可以看出矢量化后程序总体计算性能有一定提升, 而且加入矢量化之后, 同样具有非常好的弱扩展性, 可以预测在更大节点数时, 也能能够保持较好的计算效率.

矢量化前后的计算时间、通讯时间和总时间随节点数的变化情况如图7所示.

可以看到的是, 两种情况下计算时间基本一样, 减少明显的是通讯时间. 这是可以理解的, 在开启向量化之前, MPI间进行同步和归约时会有瓶颈, 每次迭代计算时有大量的同步; 而在向量化后, MPI会一次同步和归约多个元素, 因此同步和归约的次数会大大减少, 这一部分不再是瓶颈, 因此时间也相对减少了很多.

3 结论与展望

本文围绕格点QCD计算中采用数据并行模型提升上层应用计算效率的关键问题, 开展了计算算法的理论分析, 针对典型应用软件Grid开展了面向数据并行计算模型的软件架构分析, 并且结合具体计算平台, 测试了Grid应用软件的性能, 包括向量化特征和大规模可扩展性分析. 最终阐述了数据并行计算模式, 对格点QCD计算的有效性和重要性.