一种IDRstab算法高可扩展并行处理方法及系统

一种idrstab算法高可扩展并行处理方法及系统
技术领域
1.本发明属于高性能计算技术领域，尤其涉及一种idrstab算法高可扩展并行处理方法及系统。


背景技术：

2.本部分的陈述仅仅是提供了与本发明相关的背景技术信息，不必然构成在先技术。
3.随着超级计算机的发展和应用，复杂物理系统的高置信度数值模拟已经成为加速科学研究越来越重要的手段；在许多数值模拟计算中，大规模稀疏线性代数方程组求解时间所占的比重越来越大，在有的应用领域中（如惯性约束聚变、石油勘探、天气预报），其比重已达到80%以上，是整个高置信度数值模拟的关键点；对于大型稀疏线性方程组的求解，并行迭代法是唯一可行的方法；然而随着高性能计算机体系规模的扩大和节点计算能力的提高，全局通信已成为制约迭代算法有效求解大规模稀疏线性方程组的重要瓶颈；尤其是随着物理问题的日益复杂、数值模拟精度需求的不断提升，实际计算规模的数量级增长，使得迭代算法的并行实现及性能优化遇到前所未有的挑战。
4.idrstab(s,l)算法是最近提出的一种求解大型稀疏线性方程组的krylov子空间迭代法,以占用内存小、收敛速度快、实用性强的特点，引起了广泛关注；然而当前超级计算机的异构多核架构环境下，idrstab(s,l)算法的并行性能常常受限于内积计算的全局通信和同步问题；在 idrstab(s,l)算法中，每次循环需要 (3s+1)l+1全局通信；显然，当处理器数量的增加时，这些全局通信开销变得越来越重要，使得算法在上千核上难以发挥大规模并行计算机的性能，很难满足实际应用的需求。


技术实现要素：

5.为克服上述现有技术的不足，本发明提供了一种idrstab算法高可扩展并行处理方法及系统，引入tsqr分解优化正交化过程，通过循环移位将并行tsqr分解与内积计算移位到同个同步点上，在同一个同步点上同时执行并行tsqr分解与内积计算，大大减少idrstab算法的全局通信次数，克服现有idrstab算法的全局通信耗时太长的局限性，提高算法的可扩展能力，加快算法求解速度，使得用户能够快速使用算法求解实际应用问题。
6.为实现上述目的，本发明的一个或多个实施例提供了如下技术方案：本发明第一方面提供了一种idrstab算法高可扩展并行处理方法；一种idrstab算法高可扩展并行处理方法，用于高置信度数值模拟中大型稀疏线性方程组ax=b的求解，包括以下步骤：设置并行计算参数值，创建并行执行环境；将大型稀疏线性方程组ax=b的系数矩阵a与右端向量b划分到每个处理器上；计算矩阵a与初始向量x0的矩阵乘积ax0，计算每个处理器的局部残差，将所有处理器划分成s个虚拟区域，计算每个虚拟区域的局部块状正交对偶残差向量；
生成初始块状向量；循环迭代更新近似解、残差、块状正交向量与块状向量，直到残差向量范数值小于拟定的阈值，停止循环，输出解向量；在块状正交向量的生成和更新中，引入tsqr分解优化正交化过程，对原算法idrstab（s,l）中的步骤进行循环移位，将并行tsqr分解与内积计算调整在同一个同步点上，在同一个同步点上同时执行并行tsqr分解与内积计算。
7.进一步的，所述并行计算参数值,包括处理器总数p、双正交投影的次数s 、gmres多项式算子的个数l、算法停止阈值tol。
8.进一步的，将矩阵a与右端向量b按m行划分，包括将矩阵的非零局部矩阵块c
p
、t
p
、d
p
与b
p
按上标分布到标识相同的第p个处理器上。
9.进一步的，所述将所有处理器划分成s个虚拟区域，每块区域中的处理器总数小于等于，前s-1块区域中处理器的个数是，第s块区域中的处理器个数是，其中表示向上整数。
10.进一步的，所述计算每个虚拟区域的局部块状正交对偶残差向量，具体为：划分到第i块区域()的所有处理器，将各自的局部残差向量范数值发送到i区域中的处理器上，处理器将i区域内的所有处理器发送的局部残差向量范数值进行求和并将结果返回到该区域的各进程中，最后该区域中各进程计算局部残差向量乘以求和值的倒数，将所得结果作为局部块状正交对偶残差向量的第i列的非零部分。
11.进一步的，所述引入tsqr分解优化正交化过程，包括舍去原算法中的mgs正交过程，更换为并行tsqr分解。
12.进一步的，所述对原算法idrstab（s,l）中的步骤进行循环移位，包括将算法步骤中的某些块状向量之间或向量之间的内积计算，通过向量替换，移动到并行tsqr分解步骤，与并行tsqr分解一同进行归约操作。
13.进一步的，所述在同一个同步点上同时执行并行tsqr分解与内积计算，具体为：各个进程创建一个out数组，用于存放对块状向量执行简约qr分解后所得的一个上三角形与局部块状向量之间或局部向量之间的内积值；定义新的函数，将来自不同处理器的两个数组做简约qr分解与求和计算；基于新的函数，定义归约操作，将最后结果存储到out数组，并将out数组返回到根进程中；在根进程计算参数，并将所得参数传输到其他处理器内存中。
14.进一步的，所述循环迭代更新近似解、残差、块状正交向量与块状向量，具体步骤为：步骤（1）、基于块状向量u0与u1，更新近似解与残差向量，执行矩阵向量乘计算，生成1个gmres多项式算子；步骤（2）、执行idr步骤，生成1个块状正交向量vk；
步骤（3）、循环迭代步骤（1）与步骤（2）共l次，生成l个gmres多项式算子rk与l个块状正交向量vk；步骤（4）、基于l个gmres多项式算子，更新近似解与平滑修正残差，判断是否收敛，满足输出解向量，停止迭代，否则继续；步骤（5）、更新块状向量u0与u1，返回步骤（1）。
15.本发明第二方面提供了一种idrstab算法高可扩展并行处理系统。
16.一种idrstab算法高可扩展并行处理系统，用于高置信度数值模拟中大型稀疏线性方程组ax=b的求解，包括环境初始模块、数据划分模块、块状正交对偶残差向量生成模块、块状向量初始模块、循环求解模块：环境初始模块，被配置为：设置并行计算参数值，创建并行执行环境；数据划分模块，被配置为：将大型稀疏线性方程组ax=b的系数矩阵a与右端向量b划分到每个处理器上；块状正交对偶残差向量生成模块，被配置为：计算矩阵a与初始向量x0的矩阵乘积ax0，计算每个处理器的局部残差，将所有处理器划分成s个虚拟区域，计算每个虚拟区域的局部块状正交对偶残差向量；块状向量初始模块，被配置为：生成块状向量u0与块状向量u1；循环求解模块，被配置为：循环更新近似解、残差、块状正交向量与块状向量，直到残差向量范数值小于拟定的阈值，停止循环，输出解向量；在块状正交向量的生成和更新中，引入tsqr分解优化正交化过程，对原算法idrstab（s,l）中的步骤进行循环移位，将并行tsqr分解与内积计算调整在同一个同步点上，在同一个同步点上同时执行并行tsqr分解与内积计算。
17.以上一个或多个技术方案存在以下有益效果：本发明以减少全局通信次数为目的，合理设计整体方案的并行处理方法，充分考虑计算和通信开销的平衡等，提出了一种idrstab算法高可扩展并行处理方法，可以提高算法的可扩展性，使得算法性能不会因处理器的增加而下降，从而使稀疏线性方程组在高性能计算机上可快速求解，进而减少了实际问题的求解时间。
18.本发明引入tsqr分解优化正交化过程，循环移位消除迭代步之间的内积计算的数据依赖性以便于将tsqr分解与内积计算调整到同一个同步点上，在同一个同步点上同时执行tsqr分解与内积计算，大大减少了现有idrstab算法的全局通信次数，克服现有idrstab算法的全局通信耗时太长的局限性，加快算法求解速度，使得用户能够快速使用算法求解实际应用问题。
19.本发明附加方面的优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
20.构成本发明的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。
21.图1为第一个实施例的方法流程图。
22.图2为第二个实施例的系统结构图。
23.图3为第一个实施例中同一个同步点上的并行tsqr分解与内积计算同时执行示意图。
24.图4为第一个实施例中并行试验总执行时间与全局通信时间图。
25.图5为第一个实施例中并行试验全局通信占比与求解时间改善量图。
26.图6为第一个实施例中并行试验总执行时间图。
具体实施方式
27.应该指出，以下详细说明都是示例性的，旨在对本公开提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本公开所属技术领域的普通技术人员通常理解的相同含义。
28.需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本公开的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。
29.在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
30.实施例一本实施例公开了一种idrstab算法高可扩展并行处理方法；如图1所示，一种idrstab算法高可扩展并行处理方法，包括：s1.启动mpi并行执行环境、获取参数。
31.s11. 调用mpi_init函数初始化 mpi并行程序的执行环境，预先定义mpi_comm_world通信域(获取所有进程)；s12. 调用mpi_comm_rank函数来获得当前进程标识；s13. 调用mpi_comm_size函数来获取通信域包含的进程总数；s14. 根进程(0进程)读取处理器的总个数p，双正交投影的次数s，gmres多项式算子的个数l，算法停止阈值tol等相关参数，同时设置每个处理器中局部矩阵维数为m行m列，整个求解问题规模为；s15. 根进程调用mpi_bcast函数将p,m,s,l,n,tol等相关参数拷贝到其他处理器中。
32.s2.数据划分。
33.s21. 将大型稀疏线性方程组ax=b的系数矩阵a与右端项b按m行进行划分成p块，将局部非零矩阵块c
p
、t
p
、d
p
与局部向量b
p
按上标划分到相同标识的处理器上，即第p处理器；s3.生成n行s列的块状正交对偶残差向量。
34.各个进程执行局部矩阵向量乘，计算相应处理器中的局部初始残差与局部初始残差量范数值，将所有处理器虚拟划分为s块区域，在第i块()区域内的进程将相应处理器上的局部初始残差范数值发送给指定进程，指定进程求和并将结果返回到该区域的
(1)各个进程调用dcopy函数将处理器中局部残差赋值给局部块状向量第一列，令i=1；(2)参照步骤s31，执行局部矩阵向量乘计，将结果赋值给的第i+1列,赋值给的第列；(3)重复(2)步骤s次，具体的，从i等于2开始到i等于s结束；s42. 在同一个同步点上，对块状向量u0进行并行tsqr分解的同时执行与内积运算，这些浮点计算与全局通信按照以下关键路径执行：各个进程创建一个数组out，用于存放tsqr分解后所得的上三角形和内积值；定义pqrs函数，对来自不同处理器的两个数组out做简约qr分解与求和计算；基于pqrs函数，定义归约操作，将最后结果存储到out数组，并将out数组返回到根进程中；在根进程执行并将结果传播到其他处理器内存中。
38.参阅图3，其网络拓扑为二叉树结构；(1) 各个进程调用allocate函数动态创建临时数组out，调用degmqr函数对各自处理器的局部块状向量进行简约qr分解，生成一个上三角矩阵，存储到out的从第1列到s列out(:,1:s)；(2) 各个进程执行与内积运算，并将结果分别存储到的数组out(:,s+1:2s)位置与out(:，2s+1)位置；(3) 定义pqrs函数，将来自不同处理器的两个数组out数组的前部分上下排列成矩阵后，执行简约qr分解，获得的一个上三角形矩阵存储到out(:,1:s)；将两个数组的后部分求和存储到out(:,s+1:2s+1)的位置；(4) 再调用mpi_op_create函数将pqrs函数与op操作联合起来定义归约mpi_op_pqrs操作；(5) 根进程调用mpi_reduce函数，将所有进程的out进行mpi_op_pqrs归约操作，将out数组返回到根进程中；(6) 根进程计算out(:,1:s)的逆矩阵存储到数组t-1
(全局变量)中，将数组out(:,s+1:2s)位置的元素赋值给，计算；再将out(:,2s+1)位置的元素赋值给，计算；(7) 调用deallocate释放out数组存储空间；(8) 根进程调用mpi_bcast将相应处理器中的t-1
、、广播到其他处理器内存中；(9) 最后调用mpi_op_free函数将定义的归约操作撤销。
39.s43.各个进程调用omp parallel num_threads函数派生num_threads个线程，将局部块状向量与按m/num_threads行分配给每个线程块，定义t-1
为全局变量，每个线程执行、运算，当所有线程都完成了乘法操作后，主线程里将结果存储到处理器内存中。
40.s5.开始while循环迭代。
41.输入参数k等于1(具体的，“一次循环”是指从k等于1开始到k等于l结束)。
42.s51. 更新近似解与残差，生成1个gmres多项式算子rk。
43.(1)各个进程调用omp parallel num_threads函数将局部向量、局部块状向量
…
，按m/num_threads行分配给每个线程块，并将矩阵定义为全局变量，各个线程执行局部向量更新操作，当所有线程都完成了操作，从线程将终止，最后在主线程里将结果存储相应的处理器中；(2)各个进程参照步骤s31，执行矩阵向量乘运算，最后将结果赋值给。
44.s52. 执行idr步骤，完成s次基于正交投影，生成一个块状向量vk。(输入参数i等于1，具体一次“idr步骤”是指从i等于1开始到i等于s结束)。
45.(1)各个进程调用allocate函数动态创建临时数组，将赋值给临时数组v
p
；(2)各个进程调用ddot函数计算局部内积，调用mpi_reduce函数将所有进程的局部内积值求和归约，并向根进程输出计算结果，随后根进程计算，调用mpi_bcast广播将b到各个进程；(3)各个进程调用omp parallel num_threads函数将块状向量
…
按m/num_threads行分配给每个线程块，并将矩阵b定义为全局变量，各个线程执行向量更新操作，当所有线程都完成操作，从线程将终止，最后在主线程里将结果存储到相应的处理器中；(4)最后各个进程参照步骤s31，执行局部矩阵向量乘计算；(5)各进程调用allocate函数动态创建临时数组，各个进程调用dcopy函数将赋值给临时数组，再将赋值给，同时将赋值给的第i列；调用deallocate释放数组存储空间；(6)若i小于s，则，返回s53的第(1)步骤继续迭代；否则执行s54。
46.s53. 生成块状正交向量vk且同时执行内积计算。
47.(1) 参照步骤s42方式，在同一个同步点上对执行并行tsqr分解的同时执行内积计算，最后将根进程的t-1
,,广播到各个进程的处理器中；(2) 各个进程调用omp parallel num_threads函数将局部块状向量
按m/num_threads行分配给每个线程块，并将矩阵t-1
定义为全局变量，各个线程执行向量更新操作，当所有线程都完成操作后，再执行向量替换操作，从线程将终止，最后在主线程里将结果返回存储到处理器内存中。
48.s54. 循环移位：判断k的大小是否小于l，若满足则k=k+1,返回第 s52步并继续执行到第s54步；否则返回到第s52步并执行到第s53步后，跳转到第s56步。
49.s55. 生成块状正交向量v
l
且同时执行内积计算。
50.（1）参照步骤s42方式，各个进程调用allocate函数动态创建临时数组out，数组out1，数组out2，调用degmqr函数对各自处理器内的局部块状向量进行简约qr分解，生成一个上三角矩阵，存储到out(:,1:s)，与此同时做矩阵的局部部分的内积计算，输出结果分别存储到数组out(:,s+1: ls)位置，out1与out2位置；(2) 定义pqrs1函数，即输入来自不同处理器的数组out，前部分归约进行qr分解，即上下排列成矩阵，进行qr分解，获得的一个上三角形矩阵存储到out(:,1:s)，后续内积计算值分别求和存储到out(:,s+1:ls)位置，out1与out2位置；(3) 调用mpi_op_create函数将pqrs1函数与op操作联合起来定义归约mpi_op_pqrs1操作；(4) 根进程调用mpi_reduce函数，将所有进程的out，out1与out2数组进行mpi_op_pqrs1归约操作，将最后结果分别存储到out，out1与out2数组，并将这些数组返回到根进程中；(5) 根进程计算out(:,1:s)的逆矩阵存储到数组t-1
中，将数组out(:,s+1:2s)，out(:,2s+1:ls)位置的元素分别赋值给ε1与ε2，将out2(:，1)，out2(:，2：l+1)位置的元素分别赋给ε3与ε4，同时数组out2位置中的元素计算出；根进程调用mpi_bcast将根进程的ε1、ε2、ε3、ε
4、
g广播到各个进程的处理器中；(6) 调用deallocate释放out数组存储空间，调用mpi_op_free函数将定义的归约操作撤销。
51.s56. 更新近似解与平滑修正残差，判断算法是否收敛。
52.(1)各进程调用omp parallel num_threads函数将局部向量按m/num_threads行分配给每个线程块，并将矩阵定义为全局变量，各线程执行向量更新，当所有线程都完成操作，从线程将终止，最后在主线程里将结果返回到相应处理器中；(2) 根进程调用dnrm2函数来计算局部残差2-范数的平方，根进程调用
mpi_reduce函数，求残差向量2-范数，判定残差的2-范数是否小于拟定的阈值时，满足跳出while主循环，且调用deallocate释放v数组存储空间；否则继续。
53.s57. 生成块状向量u0与块状向量u1。
54.(1)各进程调用omp parallel num_threads函数将局部块状向量按m/num_threads行分配给每个线程块，并将矩阵g定义为全局变量，各线程执行块状向量更新，当所有线程都完成操作，从线程将终止，最后在主线程里将结果返回相应处理器中；(2)各进程计算参数，返回第s51步，继续迭代直到收敛。
55.本发明所提出的一种求解大型稀疏线性方程组的idrstab算法高可扩展并行处理方法，对idrstab算法的并行方案的每个技术细节提出了相应的解决方案，采用该并行方法，提高算法的可扩展性，使得算法性能不会因处理器的增加而下降，从而使稀疏线性方程组在高性能计算机上可快速求解，进而减少了实际问题的求解时间。
56.利用上述方法，采用天河二号高性能计算机系统和惯性约束聚变icf大规模数值模拟软件包中的测试算例二维扩散方程进行验证，算例规模从6.912e+06到1.728e+07；测试平台信息如下：采用7000刀片节点，每个节点2*12核(xeon e5-26992v2)，主频2.200ghz，64gb内存；采用mpi库，gcc 4.1.2进行编译；为了保证单个处理器的规模为，分别采用96~4800个处理器对本发明实现的idrstab高可扩展并行算法进行验证，其中参数s,l分别选为4，2与2，1。
57.在天河二号高性能计算机系统上的并行试验结果如图4、5、6所示，横轴是处理器的个数，对98~4800个处理器分别进行了测试，并记录了总执行时间、通信时间以及通信占比，其中，、分别表示idrstab算法与idrstab高可扩展并行算法的总执行时间，、分别表示idrstab算法与idrstab高可扩展并行算法的通信时间；表示idrstab算法的通信时间占总执行时间之比，表示idrstab高可扩展并行算法的通信时间占总执行时间之比，表示idrstab高可扩展并行算法比idrstab算法在总执行时间上的改进率，表示idrstab高可扩展并行算法比idrstab算法在通信时间上的改进率。
58.图4、5、6充分说明，随着问题规模逐渐扩大，处理器个数增加，通信时间占比是越来越重；问题规模为48004800时，idrstab算法的全局通信占比趋于72.2%，而idrstab高可扩展并行算法的全局通信占比趋于51.3%。与此同时，idrstab高可扩展并行算法比idrstab算法在总执行时间上省了45.5%，全局通信时间节省了61.32%；这说明本发明的idrstab高可扩展并行算法可大幅度降低全局通信时间，进而可快速求解实际应用问题；也就是说，在天河二号最高测试到4800个核心数，相比串行实现，该问题通信时间大幅度减少，最高加速比达到了60.3%，从而整体地大幅度降低了求解时间。
59.实施例二
本实施例公开了一种idrstab算法高可扩展并行处理系统；如图2所示，一种idrstab算法高可扩展并行处理系统，用于高置信度数值模拟中大型稀疏线性方程组ax=b的求解，包括环境初始模块、数据划分模块、块状正交对偶残差向量生成模块、块状向量初始模块、循环求解模块：环境初始模块，被配置为：设置并行计算参数值，创建并行执行环境；数据划分模块，被配置为：将大型稀疏线性方程组ax=b的系数矩阵a与右端向量b划分到每个处理器上；块状正交对偶残差向量生成模块，被配置为：计算矩阵a与初始向量x0的矩阵乘积ax0，计算每个处理器的局部残差，将所有处理器划分成s个虚拟区域，计算每个虚拟区域的局部块状正交对偶残差向量；块状向量初始模块，被配置为：生成块状向量u0与块状向量u1；循环求解模块，被配置为：循环更新近似解、残差、块状正交向量与块状向量，直到残差向量范数值小于拟定的阈值，停止循环，输出解向量；在块状正交向量的生成和更新中，引入tsqr分解优化正交化过程，对原算法idrstab（s,l）中的步骤进行循环移位，将并行tsqr分解与内积计算调整在同一个同步点上，在同一个同步点上同时执行并行tsqr分解与内积计算。
60.上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。