多色SSOR-PCG的MPI编程实现

林绍忠，许合伟

（1.长江科学院院长办公室，武汉 430010;2.黄河勘测规划设计有限公司，郑州 450003）

多色SSOR-PCG的MPI编程实现

林绍忠1，许合伟2

（1.长江科学院院长办公室，武汉 430010;2.黄河勘测规划设计有限公司，郑州 450003）

对称逐步超松驰预处理共轭梯度法（SSOR-PCG）是一种求解大型稀疏对称正定线性方程组的非常有效的迭代法。SSOR-PCG并行化的难点在于每步迭代都要求解2个三角方程组。采用一种改进的SSOR-PCG并行求解有限元方程组，并采用多色排序技术提高并行度。基于MPI模型开发了并行程序，通过测试，选择了有效的MPI通信函数。

SSOR-PCG;并行计算;多色排序;有限元方程组;MPI

1 研究背景

预处理共轭梯度法是求解大型稀疏线性方程组的极为有效的迭代法。预处理矩阵不同，便构成不同的预处理共轭梯度法。重要的预处理共轭梯度法有不完全因子分解预处理共轭梯度法（ICCG）和对称逐步超松驰预处理共轭梯度法（SSOR-PCG）［1］。

SSOR-PCG的预处理矩阵由方程组系数矩阵的上、下三角矩阵（和一个对角矩阵）的乘积构成，构造简单，存贮量少，算法程序简单。其每步迭代的主要工作是计算系数矩阵与向量的乘积以及求解2个稀疏三角方程组。SSOR-PCG并行化的难点在于三角方程组求解过程的递推特性。对稀疏线性方程组的未知量进行多色排序，使得同色未知量相互独立，是提高并行化程度的有效措施之一，已广泛应用于经典SOR，SSOR迭代法及以其作为预处理器的共轭梯度法等算法的并行化［2-5］。

改进的SSOR-PCG［6］避免了系数矩阵与向量的乘积运算，计算量比原迭代方法少，已得到推广应用［7-12］。本文采用改进的SSOR-PCG并行求解有限元方程组，采用图的顶点着色方法对有限元网格结点进行多色排序，基于MPI（Message Passing Interface）并行编程工具［13］开发了分布内存并行程序，通过测试选择了有效的MPI通信方式。

2 改进的SSOR-PCG

设线性方程组Ax=b为n阶对称正定稀疏线性方程组。记D为A的对角矩阵，L为A的严格下三角矩阵，ω（0＜ω＜2）为松弛因子，W=D/ω+L，V=（2-ω）D/ω，则SSOR-PCG的预处理矩阵（即SSOR的分裂矩阵）为M=WV-1WT。改进的SSORPCG迭代格式为［6］:

从式（1）可见，矩阵与向量的乘积运算为W-1（z-Vd），W-Tz，Vy及Vd。由于V为对角矩阵，后2个乘积的计算量少，且容易并行化。计算内积和标量与向量的乘积也容易并行化。W为下三角矩阵，所以求解2个三角方程组是每步迭代的主要计算量，也是算法并行化的难点所在。

3 有限元网格结点多色排序

图的顶点着色［14］是指图中任意2个相邻顶点都分配到不同的颜色，即同色顶点互不相邻。对有限元网格结点着色后，同色结点连续编号，不同色的结点顺序排列。这样，同色结点对应的刚度矩阵的主块为块对角矩阵。例如，图1（a）所示的有限元网格结点按4色排序后，结点编号见图1（b），下三角刚度矩阵的稀疏形式见图2。对于这样的三角方程组，按色号依次求解，但同色结点并行求解，并行实现容易，并行度高。

图1 有限元网格多色结点编号Fig.1Multicolor ordering of finite element mesh nodes

图2 多色排序下的下三角刚度矩阵Fig.2Lower triangular stiffness matrix after multicolor ordering

需要指出的是:①单元的所有结点都是相邻的，即单元的对角结点也是相邻的，在结点着色时应加以考虑;②一个结点有多个未知量时，图中的“×”表示1个子矩阵，其中对角线上的“×”是一个下三角矩阵;③线或块多色排序同理。

4 SSOR-PCG的MPI并行程序设计

4.1 数据存储方式与数据划分

SSOR-PCG的求解过程仅涉及A的非零元素。为减少存储量，可只存储其下三角或上三角矩阵中的非零元素。稀疏矩阵的存储方式很多，为与算法相适应，本文采用常用的CSR（Compressed Sparse Row）存储方式按行压缩存储下三角矩阵，即用3个一维数组存储稀疏对称矩阵A的下三角矩阵:实型数组A按行存储非零元素的值，整型数组INDEX存储每个非零元素对应的列号，整型数组MA存储每行主元在数组A中的位置。式（1）中的W与A的下三角矩阵的区别仅在于主对角线元素不同，所以W存储于数组A中，而无须另开辟数组。

并行计算时，为了提高并行效率，应保证各个处理器（进程）的负载平衡，即各个处理器的计算量大致相当。对于本文算法，在尽量使数组A中同色结点对应的非零元素个数大致均匀地划分到各处理器的同时，还应保证一个结点的所有自由度对应的各行非零元素完整地划分给某个处理器。由于一个结点的自由度个数少，其所在行的非零元素个数很有限，所以容易做到负载平衡。据此负载平衡原则确定划分给各处理器的各种颜色的始末方程号N1，N2和累积方程数N3后，按行块卷帘方式对各种数组进行数据划分。

设有k种颜色，p个处理器。以数组A为例，它被分为p×k个子块:A=（A1，1，A2，1，…，Ap，1，A1，2，A2，2，…，Ap，2，…，A1，k，A2，k，…，Ap，k），其中（A1，j，A2，j，…，Ap，j）对应于第j种颜色。通过调用MPI通信函数MPI_SEND和MPI_RECV，每个处理器接收从存储数组A的主进程依次发来的k个子块并顺序存储于数组ALOC中，例如第i个处理器的ALOC =（Ai，1，Ai，2…，Ai，k）。主进程的ALOC可通过赋值获取。

INDEX，MA数组的划分、分发和存储方式类似，其k个子块分别存储于INDEXLOC和MALOC数组中。但是，MALOC的数值应作适当修改，以指向ALOC的对应位置。

4.2 三角方程组的并行求解

经多色排序后，三角方程组的并行求解过程见表1。其中，IRANK为进程号;Y为n阶数组，求解前仅存本进程所需的右端项子块。可见，并行求解过程与串行前代/回代过程基本一样，不同的是，它是按色号依次求解，每求出一种颜色的解后，要在进程间进行通信。

4.3 通信方式的选择

前代过程中，各进程要将其计算得到的Y（N1: N2）传递给其他进程。有多种MPI通信方式可以实现这个过程。测试表明，全收集函数MPI_ALLGATHERV效率高。

表1中的回代过程利用了按行块卷帘存储的下三角矩阵。由于对称性，上三角矩阵是按列块卷帘方式存储于每个处理器。为方便寻址，回代过程中每求出一个解，就按列扫描算法消减对应的右端项。各进程归约所需的所有进程的消减值的子块后，继续求解前一色。也有多种方式可以实现这种归约过程。测试表明，归约散发函数MPI_REDUCE_SCATTER效率高。规模不大、进程数不多时，也可选用全归约函数MPI_ALLREDUCE。然而，回代过程的通信开销大于前代过程，因为前者传递数据量多，而且有归约运算。若同时分别按行存储上、下三角矩阵，回代过程与前代过程类似，可用MPI_ALLGATHERV节省通信时间。

表1 三角方程组并行求解过程Table 1Parallel algorithm solving triangular equations

5 算例

本文所用计算机为分布共享内存计算机，由2台曙光个人高性能计算机PHPC100构成，两者之间通过高速网络实现数据通信。每台由5个计算节点组成，每个计算节点含2个CPU（2.4 Hz），每个CPU由6个核心组成。每个计算节点的内存16 G，供节点内的12个核心共用。操作系统为SUSE Linux Enterprise Server 10（x86_64），编译器为Intel Visual Fortran。

计算对象为一个三维柱体结构，尺寸119 cm× 119 cm×238 cm，底部全约束。分析自重作用下的变形及应力。用8结点6面体有限元网格，共120× 20×120=172.8万个结点，分为8种颜色。每个结点3个自由度。SSOR-PCG迭代收敛控制参数ε=10-8。单个处理器串行求解方程组所用时间为443.4 s。

采用2种并行计算测试方案。方案1:采用1-10个计算节点，每个节点只设1个进程。方案2:采用10个计算节点，每个计算节点设不同数量的进程数。单个处理器串行求解时间与p个处理器并行计算时间之比，即加速比见表2和表3。

表2 测试方案1的加速比Table 2Speedups of test plan no.1

从表2可见，10个进程时，只存储下三角刚度矩阵和同时存储上、下三角矩阵2种情况的加速比分别为5.93和8.00。后者并行效率高，达到80%。因此，在内存容量允许的情况下，可同时存储上、下三角矩阵。

从表3可见，每个节点设3个进程时加速比最高，2种存储方式下的加速比分别为6.56和10.57。同样，同时存储上、下三角矩阵的效率高。每个节点多于3个进程时，因通信开销增大，加速比反而下降。在分布共享内存计算机上，若采用MPI+ OpenMP混合编程，将可获得更好的并行性能。

6 结论

多色SSOR-PCG的并行化程度高。在MPI编程模型下，仅按行存储刚度矩阵的下三角矩阵时，前代过程可采用MPI_ALLGATHERV通信函数;回代过程可采用MPI_REDUCE_SCATTER通信函数，规模不大、进程数不多时，也可选用MPI_ALLREDUCE通信函数;回代过程的通信开销大于前代过程。在内存容量允许的情况下，可同时分别按行存储上、下三角矩阵，以获得更高的并行性能。

基于OpenMP编程模型和MPI+OpenMP混合编程模型的多色SSOR-PCG的并行实现将另文发表。

表3 测试方案2的加速比Table 3Speedups of test plan no.2

［1］吕涛，石济民，林振宝.区域分解算法-偏微分方程数值解法新技术［M］.北京:科学出版社，1994.（LV Tao，SHIH T M，LIEM C B.Domain Decomposition Methods:New Numerical Techniques for Solving PDE［M］.Beijing:Science Press，1994.（in Chinese））

［2］ADAMS L，ORTEGA J M.A Multi-color SOR Method for Parallel Computation［C］∥International Association of Computing and Communication（IACC）.Proceedings of the 1982 International Conference on Parallel Processing. Bellaire，Michigan，August 24-27，1982:53-58.

［3］张汝清.并行计算结构力学［M］.重庆:重庆大学出版社，1993.（ZHANG Ru-qing.Parallel Computing in Structural Mechanics［M］.Chongqing:Chongqing University Press，1993.（in Chinese））

［4］吴建平，王正华，李晓梅.稀疏线性方程组的高效求解与并行计算［M］.长沙:湖南科学技术出版社，2004. （WU Jian-ping，WANG Zheng-hua，LI Xiao-mei.Efficient Solution and Parallel Computing of Sparse Linear E-quations［M］.Changsha:Hunan Science＆Technology Press，2004.（in Chinese））

［5］POOLE E L，ORTEGA J M.Multicolor ICCG Methods for Vector Computers［J］.SIAM Journal on Numerical A-nalysis，1987，24（6）:1394-1418.

［6］林绍忠.对称逐步超松驰预处理共扼梯度法的改进迭代格式［J］.数值计算与计算机应用，1997，（4）:266-270.（LIN Shao-zhong.Improved Iterative Format of Symmetric Successive over Relaxation-Preconditioned Conjugated Gradient Method［J］.Journal of Numerical Methods and Computer Applications，1997，（4）:266-270.（in Chinese））

［7］林绍忠，苏海东.大体积混凝土结构仿真应力分析快速算法及应用［J］.长江科学院院报，2003，20（6）:19 -22.（LIN Shao-zhong，SU Hai-dong.Fast Algorithms for Stress Analysis Simulating Construction Process of Massive Concrete Structures and Applications［J］.Journal of Yangtze River Scientific Research Institute，2003，20（6）:19-22.（in Chinese））

［8］程昭，陈胜宏，田胜利，等.大规模P型有限元方程组的修正SSOR-PCG解法［J］.上海交通大学学报，2002，36（11）:1608-1611.（CHENG Zhao，CHEN Sheng-hong，TIAN Sheng-li，et al.Modified SSOR-PCG Solver for Large Scale P-Version FEM Equations［J］. Journal of Shanghai Jiaotong University，2002，36（11）: 1608-1611.（in Chinese））

［9］李静.陈健云，周晶.改进的SSOR-PCG迭代法在接触问题研究中应用［J］.大连理工大学学报，2006，46（4）:533-537.（LI Jing，CHEN Jian-yun，ZHOU Jing.Application of an Improved SSOR-PCG Iteration Method to Contact Problems［J］.Journal of Dalian University of Technology，2006，46（4）:533-537.（in Chinese））

［10］陈国荣，李皇胜，李红健.大体积混凝土温度场高速求解方法的改进［J］.河海大学学报，2009，37（4）:396-399.（CHEN Guo-rong，LI Huang-sheng，LI Hong-jian. Improvement of Fast Solution Method for Mass Concrete Temperature Field［J］.Journal of Hohai University（Natural Sciences，2009，37（4）:396-399.（in Chinese））

［11］张杨.高阶三维数值计算方法理论与应用研究［D］.南京:河海大学，2009.（ZHANG Yang.Theory and Application for High Order Three Dimensional Numerical Method［D］.Nanjing:Hohai University，2009.（in Chinese））

［12］李根.基于模拟的水岩耦合变形破坏过程及机理研究［D］.大连:大连理工大学，2011.（LI Gen.Simulation-Based Research on Rock Deformation and Failure Process under Water-Rock Coupling［D］.Dalian:Dalian University of Technology，2011.（in Chinese））

［13］都志辉.高性能计算并行编程技术:MPI并行程序设计［M］.北京:清华大学出版社，2001.（DU Zhi-hui. Parallel Programming Techniques in High Performance Computing:MPI Parallel Programming［M］.Beijing:Tsinghua University Press，2001.（in Chinese））

［14］孙惠泉.图论及其应用［M］.北京:科学出版社，2004. （SUN Hui-quan.Graph Theory and Its Application［M］Beijing:Science Press，2004.（in Chinese））

（编辑:刘运飞）

MPI-Based Implementation of Multicolor SSOR-PCG

LIN Shao-zhong1，XU He-wei2
（1.Yangtze River Scientific Research Institute，Wuhan430010，China; 2.Yellow River Engineering Consulting Co.Ltd.，Zhengzhou450003，China）

The method of symmetric successive over relaxation-preconditioned conjugate gradient（SSOR-PCG）is a very effective iterative method for solving large scale sparse symmetric positive-definite linear set of equations.The difficulty in the parallelization of the SSOR-PCG lies in solving two triangular equation systems in each iteration.In this research，an improved SSOR-PCG is applied to parallel solve finite element equations，and the multicolor ordering technique is used to increase the degree of parallelism.A MPI-based parallel program is coded and choice is made for efficient MPI communication routines by tests.

SSOR-PCG;parallel computing;multicolor ordering;finite element equations;MPI（Message Passing Interface）

TP311.52

A

1001-5485（2013）05-0082-04

10.3969/j.issn.1001-5485.2013.05.018

2013，30（05）:82-85

2012-07-26;

2012-09-27

中央级公益性科研院所基本科研业务费项目（CKSF2011016）

林绍忠（1960-），男，福建福安人，教授级高级工程师，工学博士，主要从事水工结构数值分析研究，（电话）027-82820007（电子信箱）Linsz@mail.crsri.cn。

book=85,ebook=180