开放科学（资源服务）标志码（OSID）：

0 概述

全球区域同化预报系统（Global and Regional Assimilation Prediction System，GRAPES）是我国自主研发的新一代数值天气预报系统，能够改善天气预报模式的性能并有效提高天气预报精度。GRAPES动力框架的计算核心是赫姆霍兹方程的求解，该方程是原始大气方程组经过一系列离散化处理之后形成的大规模稀疏线性方程组。目前对线性方程组的求解方法主要有直接法和迭代法^[1-2]两种。由于计算机硬件资源的限制，使用直接法并不能完成有效求解，而迭代法具有硬件存储空间要求较小、原始矩阵在计算过程中保持不变等优点。因此，目前求解大规模稀疏线性方程组主要选择迭代法。

GRAPES数值天气预报系统采用的迭代法为广义共轭余差（Generalized Conjugate Residual，GCR）法。文献[3]利用基准测试程序对GRAPES并行应用进行分析，发现在迭代法计算过程中全局通信对整个系统的性能影响较大，指出在进行性能优化时需要考虑系统通信所带来的性能下降问题。文献[4]针对赫姆霍兹方程求解过程中通信占比高的现象，引入短向量计算替代长向量计算，以减少通信开销。文献[5]在使用重启动和截断两种方式进行性能优化的同时，还使用一种新的稀疏近似逆预条件子来加快算法收敛速度。

研究发现，增加重启动和截断处理的广义共轭余差法能够减少数据存储量，同时改善计算的局部性。文献[6]为改善GRAPES模式动力框架中存在的极点问题，利用阴阳网格设计新的动力框架，与原GRAPES模式框架相比，新的动力框架表现出较好的数值稳定性及计算性能。文献[7]利用MPI和OpenMP混合并行的方式提升GRAPES模式的并行度，以适应主流硬件架构既有分布式内存又有共享内存的情况。文献[8]结合国产高性能平台实现了基于MPI+共享变量编程模型的众核线程级并行的多级并行方案，提升了方程求解的运行效率。文献[9]基于PETSc科学计算库和Hypre预条件库实现了广义最小残差（Generalized Minimal Residual，GMRES）法，与GRAPES模式现有的GCR算法相比，GMRES算法在更高分辨率的条件下具有较好的性能。文献[10]则利用CPU设备结合线程级并行、矢量化和数据重用技术提升了算法计算性能。当前GRAPES模式动力框架中广义共轭余差法的优化工作主要集中在通用CPU及国产高性能机器领域，部分学者讨论了算法本身的改进方案，但大多是与CPU设备进行兼容，而结合GPU的性能优化相对较少^[11]。

本文针对GRAPES数值天气预报系统中的赫姆霍兹方程求解问题，分别通过MPI、MPI+OpenMP、CUDA三种并行方式实现求解大规模稀疏线性方程组的广义共轭余差法，并对测试结果进行对比分析，评估优化性能。

1 GRAPES模式中大规模稀疏线性方程组求解

GRAPES数值天气预报系统主要包括动力框架和物理过程参数化方案两部分^[12-13]，其中动力框架为核心部分，将动力框架方程组进行简化和离散后可得到GRAPES模式的基本预报方程^[14-16]。对该方程进行处理后，动力框架的主要计算内容就变成对一个包含大型稀疏矩阵的赫姆霍兹方程的求解^[17-19]。如图 1所示，赫姆霍兹方程的求解过程计算量庞大，占据整个过程求解计算量的30%以上。随着模式分辨率的提高，方程的计算量还会增加。因此，对赫姆霍兹方程进行高效求解成为动力框架性能提升的关键^[20]。

如图 2所示，GRAPES模式动力框架部分的赫姆霍兹方程中计算每一个格点的方程需要选取空间中与其相关的19个格点作为系数^[21]。

空间格点的计算公式如下：

$ \begin{array}{l}{\left({\xi }_{{\mathit{\Pi}} 0}\right)}_{i, j, k}={B}_{1}{\left({\mathit{\Pi}} \right)}_{i, j, k}+{B}_{2}{\left({\mathit{\Pi}} \right)}_{i-1, j, k}+{B}_{3}{\left({\mathit{\Pi}} \right)}_{i+1, j, k}+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {B}_{4}{\left({\mathit{\Pi}} \right)}_{i, j-1, k}+{B}_{5}{\left({\mathit{\Pi}} \right)}_{i, j+1, k}+{B}_{6}{\left({\mathit{\Pi}} \right)}_{i+1, j+1, k}+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {B}_{7}{\left({\mathit{\Pi}} \right)}_{i+1, j-1, k}+{B}_{8}{\left({\mathit{\Pi}} \right)}_{i-1, j-1, k}+{B}_{9}{\left({\mathit{\Pi}} \right)}_{i-1, j+1, k}+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {B}_{10}{\left({\mathit{\Pi}} \right)}_{i, j, k-1}+{B}_{11}{\left({\mathit{\Pi}} \right)}_{j-1, j, k-1}+{B}_{12}{\left({\mathit{\Pi}} \right)}_{i+1, j, k-1}+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {B}_{{}_{13}}{\left({\mathit{\Pi}} \right)}_{i, j-1, k-1}+{B}_{14}{\left({\mathit{\Pi}} \right)}_{i, j+1, k-1}+{B}_{15}{\left({\mathit{\Pi}} \right)}_{i, j, k+1}+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {B}_{16}{\left({\mathit{\Pi}} \right)}_{i-1, j, k+1}+{B}_{17}{\left({\mathit{\Pi}} \right)}_{i+1, j, k+1}+{B}_{18}{\left({\mathit{\Pi}} \right)}_{i, j-1, k+1}+\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; {B}_{19}{\left({\mathit{\Pi}} \right)}_{i, j+1, k+1}\left(\mathrm{ }1\mathrm{ }\right)\end{array} $

将每一个空间格点的计算进行组合后可以得到一个大规模的稀疏线性方程组，该方程组可简化为：

$ \mathit{\boldsymbol{A}}\mathit{\boldsymbol{x}}=\mathit{\boldsymbol{b}} $

(2)

其中：A是一个大规模稀疏矩阵；x是规模为空间总格点数的解向量；b是与x相对应的方程右端向量。因此，GRAPES模式动力框架中赫姆霍兹方程的求解就是对一个大规模稀疏线性方程组的求解。在当前系统中，方程求解所采用的算法为广义共轭余差法，在满足精度要求的前提下，广义共轭余差法算法收敛较快，整体算法易于实现，综合表现较好^[22]。广义共轭余差法的具体内容如算法1所示。

算法1    广义共轭余差法

输入    矩阵A，初始解x₀，右端向量b

输出    近似解x

1.r₀ =b-Ax₀，p₀=r₀

2.for i=1：maximum number of iterations do

3.α_i-1=（r_i-1，Ap_i-1）/（Ap_i-1，Ap_i-1）

4.x_i=x_i-1+α_ip_i-1

5.r_i=r_i-1–α_iAp_i-1

6.if converged then exit；

7.for j=int[（i-1）/k]k，…，i-1 do

8.α_ij=-（Ar，Ap_j）/（Ap_j，Ap_j）

9.end for

10.$ {\mathrm{p}}_{\mathrm{i}}={\mathrm{r}}_{\mathrm{i}}+\sum _{\mathrm{j}=\mathrm{i}\mathrm{n}\mathrm{t}\left[\right(\mathrm{i}-1)/\mathrm{k}]\mathrm{k}}^{\mathrm{i}-1}{\mathrm{\alpha }}_{\mathrm{i}\mathrm{j}}{\mathrm{p}}_{\mathrm{j}} $

11.$ \mathrm{A}{\mathrm{p}}_{\mathrm{i}}=\mathrm{A}{\mathrm{r}}_{\mathrm{i}}+\sum _{\mathrm{j}=\mathrm{i}\mathrm{n}\mathrm{t}\left[\right(\mathrm{i}-1)/\mathrm{k}]\mathrm{k}}^{\mathrm{i}-1}{\mathrm{\alpha }}_{\mathrm{i}\mathrm{j}}\mathrm{A}{\mathrm{p}}_{\mathrm{j}} $

12.end for

2 广义共轭余差法的并行实现 2.1 问题模型

本文主要考虑问题模型Ax=b，等式中的A为稀疏矩阵，x为需要求解的向量，b为右端向量。系数矩阵具有明显的稀疏性，在矩阵存储格式的选择上，统一使用CSR（Compressed Sparse Row）格式^[23-24]。方程求解算法选择广义共轭余差法，算法主要计算内容包括稀疏矩阵向量乘、向量内积、向量数乘。

2.2 系数矩阵预处理优化

迭代法具有依赖系数矩阵条件数的特点，可能存在收敛速度慢甚至不收敛的情况，因此，通常使用预处理技术对算法进行优化。目前的预处理技术主要通过构建预条件子将原方程Ax=b转化为$ {\mathit{\boldsymbol{M}}}^{(-1)}\mathit{\boldsymbol{A}}\mathit{\boldsymbol{x}}={\mathit{\boldsymbol{M}}}^{(-1)}\mathit{\boldsymbol{b}} $等形式求解。常用的预条件子构造方法有稀疏近似逆预条件、不完全分解预条件等。本实验所采用的预条件技术为不完全分解LU预条件子（ILU）。预处理过程将原始矩阵A分为上三角矩阵U、下三角矩阵L及残差矩阵R，分解之后的矩阵可以满足不同条件P，例如与原系数矩阵的稀疏结构保持一致。构造ILU预条件子的过程如算法2所示。

算法2    ILU分解

输入    矩阵A，分解条件P

输出    预处理子M

1.for i=2，3，…，n do

2.for k=1，2，…，i-1 and（i，k）∉P do

3.A_ik=A_ik/A_kk

4.for j=k+1，…，n and（i，j）∉P do

5.A_ij=A_ij-A_ik*A_kj

6.end for

7.end for

8.end for

2.3 MPI并行及MPI+OpenMP混合并行

通信开销是MPI程序中不可忽视的部分，尤其是随着进程规模的增大，在计算时间减少的同时也会增加程序的通信开销^[25-26]。因此，在设计MPI并行算法时，需要平衡计算任务划分与通信函数使用之间的关系。在广义共轭余差法中，每一个迭代步骤之间都存在数据依赖，并行任务主要选择在每一个迭代步骤内部进行。MPI并行方式将计算任务根据分配的进程数量进行划分，每个进程处理不同的子任务，最后利用进程间通信进行数据同步。

本文实验所测试的MPI+OpenMP混合并行方式，在MPI并行方式的基础上，利用多线程实现了更细粒度的任务划分。OpenMP是共享内存并行编程模式，与MPI并行模式相比可减少数据通信的时间。在矩阵向量乘部分的并行处理中，MPI+OpenMP混合并行任务划分如图 3所示。矩阵数据按行子域划分为不同的row_part，每个MPI进程只计算所属row_part与向量的乘积，各进程并行执行计算任务。同时在MPI进程内部利用OpenMP模式将row_part划分为更细粒度的subrow_part，结合线程级并行完成子任务计算。

2.4 GPU并行优化 2.4.1 数据传输优化

利用CUDA并行求解方程时主要开销在于GPU与CPU之间的数据传输，为减少数据传输时间开销，避免PCI-E总线上的数据传输成为程序性能提升的瓶颈，将矩阵A以及各个中间向量等都保存在GPU上。在计算过程中利用显存进行通信，避免和主机端进行频繁的数据传输。无法避免的数据传输过程只出现在需要进行收敛判断及主机端计算的情况下，同时在数据传输时采用页锁定内存（pinned memory）进一步提高数据传输速度。如图 4所示，与可分页内存（pageable memory）相比，页锁定内存能够减少一次数据传输操作，提高数据传输效率。

2.4.2 访存优化

GPU对全局存储器的访问速度也会对程序性能造成一定影响，利用CUDA模型可以实现线程块间及线程块内两级并行，从而完成对全局存储器的合并访问。初始的数据存储格式由于将需要连续访问的数据离散存储，不满足线程合并访问的条件，因此需要对数据存储格式进行转换。将三维格式的矩阵转化为一维格式存储之后，每一个网格点在计算过程中同一线程束内的线程所访问的内存空间是连续的，满足合并访存的要求。如图 5所示，当同一线程束内的线程访存连续时，能够将多次访存操作减少为一次访存操作，从而减少计算过程中的访存开销。

2.4.3 存储器优化

在GPU计算过程中，存在一部分数据只需要线程块内部的线程访问。对于这一部分数据，利用共享存储器进行保存。共享存储器的访问延迟比全局存储器低，在计算过程中线程块内部的线程只需要访问共享存储器中，而不需要去访问全局存储器，从而进一步提高访存的速度。

在向量内积的计算过程中，并行方式如图 6所示。内积计算过程可以表示成$ \sum\limits _{i=1}^{n}{a}_{i}{b}_{i} $形式，向量a和向量b分别对应图中Vector a和Vector b。当核函数启动后，将向量子任务subVector分配到不同的线程块中。由于向量内积涉及到对同一内存地址的修改操作，因此先利用线程块访问各自私有的共享存储器进行局部规约操作，这一部分计算内容可以并行执行。当所有subVector都完成局部归约后，再进行最后的规约操作，从而提高整体的计算速度。

3 实验与分析 3.1 实验环境及测试数据

实验中共使用5个CPU节点，节点处理器为Intel^®Xeon^® CPU E5-2692 v2 @ 2.20 GHz。每个节点两块CPU（共24核心），节点内存为64 GB；实验使用GPU为Tesla T4，显存容量为16 GB；CUDA版本为v10.1。实验所使用的测试数据包含系数矩阵A、初始解x₀、右端向量b，其中系数矩阵A规模为360×180×38，共2 462 400个网格点。

3.2 结果分析

实验测试了数据在不同计算方式下的计算残差。结果表明，使用3种并行方式实现的广义共轭余差法与串行方式能够得到一致的计算结果，验证了不同并行方式求解稀疏线性方程组的正确性和有效性。

3.2.1 预处理优化结果

对原始GCR算法进行预处理优化后，方程的收敛速度获得大幅提升。如图 7所示，在误差精度要求为1E-10的条件下，预处理算法收敛所需的迭代次数从134次下降至21次。ILU预处理的引入能够改善方程的求解性能，同时也会增加新的计算内容。但是在并行算法中，计算量已经不是限制性能的最大因素，例如在MPI程序中，随着并行度的增加，进程之间的通信成为限制程序性能的最大因素。在利用MPI实现的并行GCR算法中，每一次迭代过程都会进行数据通信，所以，利用预处理技术加快算法收敛速度之后，整体的计算性能可以获得提升。

3.2.2 MPI并行结果

虽然使用MPI并行方式能够通过提高进程数目减少子任务的计算量，从而有效提升计算部分速度，但是也会增加数据通信的时间开销。MPI并行算法在实现过程中避免了冗余的通信开销，只对部分必要的向量更新操作进行通信。MPI并行算法在不同进程规模下的运行时间如图 8所示。可以看出，随着进程规模的增加，运行时间呈先下降后上升的趋势。这是因为MPI算法的进程数量会影响进程之间的通信开销，当进程数为32时，并行计算任务所带来的性能提升不足以抵消增加的通信开销，程序整体性能就会下降。由于MPI并行方式存在可扩展性问题，因此只通过增加进程数量的方式来提高性能并不能高效地利用计算资源。

3.2.3 MPI+OpenMP混合并行结果

与MPI并行方式相比，MPI+OpenMP混合并行方式充分结合了共享内存的优势，在不改变原始MPI并行方式粗粒度通信的情况下，提高了进程内部的线程级并行度。图 9为MPI并行与MPI+OpenMP混合并行在不同核数规模下的计算性能对比。可以看出，当计算所使用的核数过少时，利用OpenMP进行更细粒度的并行划分之后性能并没有明显改善，但随着核数的增加，OpenMP的共享内存优势得以发挥，计算性能相比MPI并行方式有一定提升。当核数为64时，MPI并行方式已经出现性能下降，但是MPI+OpenMP混合并行方式的性能依然较稳定。因此，在MPI并行方式由于可扩展性原因造成计算效率下降的情况下，使用MPI+OpenMP混合并行方式可以获得更好的性能。

3.2.4 GPU并行结果

基于GPU的CUDA并行方式在优化措施上采用了共享存储器优化、全局存储器访存优化以及主从设备数据传输优化。CUDA并行与MPI+OpenMP混合并行在算法求解过程中主要步骤的性能对比如图 10所示。与MPI+OpenMP混合并行相比，CUDA并行在矩阵向量乘部分（ar）有32%的性能提升，向量数乘部分（p & ap，x & r）有88%的性能提升，向量内积部分（beta，alpha，residual）有80%的性能提升。

3.2.5 优化性能对比

对不同优化方式的优化性能进行对比，结果如图 11所示。可以看出，通过预处理优化减少算法迭代次数对整体计算性能提升明显，加速比达到4.5。在预处理的基础上，使用MPI并行方式的性能获得提升，加速比达到11.7。结合OpenMP并行方式进行细粒度任务划分之后，混合并行计算性能相对于MPI并行方式有一定提升，加速比达到15.8。同时数据也显示CUDA并行方式相对于MPI+OpenMP混合并行方式能够获得约50%的性能提升，加速比达到23.4。基于GPU的并行方式在计算时间上的开销远小于基于CPU的并行方式，这得益于GPU数量众多的计算核心。在数据通信方面，GPU并行方式线程级的通信相比于CPU并行方式也具有优势，不需要进行各个进程之间的数据通信，所以基于GPU的CUDA并行方式能获得较好的计算性能。

4 结束语

大规模线性方程组的求解作为科学计算中的核心问题，长期以来都是研究者关注的重点，如何充分利用目前的高性能计算机资源对求解问题的计算性能进行优化是重要的研究方向。本文从GRAPES数值天气预报系统动力框架中的赫姆霍兹方程求解问题出发，使用不同并行方式对广义共轭余差法进行并行及优化，并对计算性能进行对比分析。实验结果表明，在MPI并行、MPI+OpenMP混合并行及CUDA并行3种优化方式中，基于GPU的CUDA并行方式能够获得更好的计算性能。后续将建立多机混合并行的优化模型，测试MPI+CUDA的多节点优化性能，同时还将分析多重网格预处理、稀疏近似逆预处理等其他预处理方式对迭代法收敛速度的影响。