米瑞琪 于 朋

（中国石油大学〈华东〉，山东 青岛266580）

Gauss-Seidel迭代法求解线性方程组的并行化研究

米瑞琪 于 朋

（中国石油大学〈华东〉，山东 青岛266580）

Gauss-Seidel迭代法在工程计算中具有广泛应用。高维矩阵对线性方程组的求解效率提出了新的要求。本文提出了两种模式下的Gauss-Seidel并行计算方法，并对比了在不同矩阵维数下的加速效率，得到了具有较高加速比的并行迭代方法。

Gauss-Seidel；迭代法；并行；MPI;矩阵运算

1 问题背景与提出

在实际工程应用中经常需要用到求解维数较高的线性代数方程组，对于线性代数方程组的求解方法可以分成两类：直接法和迭代法，其中直接法得到的是方程组的真解，其求解过程为通过有限步四则运算，常用的实现方式是Gauss消去方法和矩阵的三角分解方法，这种计算方式在求解过程中会产生大量的非零元素，存储量和计算量均比较大，因此常用于求解低阶、稠密线性方程组。对于高维线性方程组，可以采用并行运算的方式进一步提高计算效率，Gauss-Seidel迭代法具有收敛速度快、计算稳定性好等优点，是求解高维线性方程组的常用方法，因此本文将对Gauss-Seidel迭代法的并行化进行研究，期望获取较高的运行效率和加速比。

2 Gauss-Seidel迭代法的并行算法描述

任务分配模式一：连续等分矩阵法

在《并行算法的设计与分析》一书中给出了一种矩阵分配的算法，设矩阵的阶为N，采用m个线程进行计算，则将矩阵等分成m块，每一块分到的任务量为N/m,每一个线程计算完自己的分量之后立刻向其他的线程进行广播。

由于Gauss-Seidel迭代法中，对于右端项的计算和中间项的计算是不需要用到下一次迭代值的，因此在这一阶段各个线程工作量相同，同时完成，问题的关键在于首项的计算，顺序等分矩阵时，以四个线程为例，在计算左端项的时候，各个线程的工作次序如下图所示。其中不同的颜色代表不同线程的任务，每个线程依次计算分给自己的一个子块，每计算完一个分量后立刻沿上方箭头所指方向进行广播。

任务分配模式二：离散等分矩阵法

为了改进任务分配模式一中的缺陷，进一步提高计算的并行程度，同时保持任务的均衡性，笔者对模式一中的任务分配方式作出了一定程度的改进，设线程数为m，则将矩阵以每m行作为一个单位，顺序分配给m个线程，前一个线程完成自己的计算任务之后，立刻将计算结果广播给其余的m个线程，同时开始下一个矩阵块中自己应完成的任务，直到最后将自己所应完成的N/m个分量计算完毕。任务分配模式见左图，其中不同的颜色代表不同的线程的计算任务，每个块右端的箭头表示变量的广播方向。

这样划分任务，每个线程的工作量没有发生变化，但是第一个达到终点的线程和最后一个达到终点的线程的时间相差仅为m个分量的计算时间。假设矩阵的阶为10000，采用4个线程进行计算(m=4)，那么第一个线程和最后一个线程到达终点的时间差为m=4个分量的计算时间。反观第一种任务分配模式，假设第一个线程计算完自己的N/m= 10000/4=2500个分量，则需要再经过计算2500、5000、7500个分量的时间之后二、三、四号线程才能到达终点，并且随着矩阵阶的增大，这个时间差异还会继续增大。

根据上面的分析，在理论上模式二具有更优的时间性能。笔者也进行了实验，从实际操作上证明了模式二具有更优的时间性能。因此选择模式二作为并行化的基本算法，基本并行算法步骤可以描述如下：

（一）末项/中间项计算的并行化

末项与中间项在计算过程中不涉及未知量的运算，可以按行进行划分，将运算分配到每一个核上。由于系数矩阵A和常数项b保持不变，因此该项无论是在串行算法还是并行算法中均只计算一次即可。

假设线性方程组是n阶的，在m个核上进行并行计算，则将末项与中间项的矩阵划分为块，将m个划分好的矩阵顺序分配给m个核（核1,2,…,m分到的计算任务为

（二）首项的并行计算

首项是下三角矩阵与本次迭代的结果相乘，因此不易实现并行化，但是借鉴并行计算中的流水线作业思想，可以采用步进和广播的方法进行并行化。假设核数为4，则首项的并行化可以设计如下：

2.2 基于MPI的多核并行算法的设计

采用任务分配模式二，采用MPI的并行算法可以描述如下：

int N:矩阵的阶数

int numProcs:用于计算的进程数量

int m:每个进程所分得的任务量，m=N/numProcs

int myID：标识本进程的进程号

int task:记录本次循环中本线程计算的未知量下标

double A[m][numProcs]:每个进程单独

保存属于自己的矩阵分块

double x_old[N]:double类型的N维数组，用于保存本次迭代分量的数值

double sum[N]:保存迭代分量计算过程中的累加结果；

double x_diff：相邻两次迭代运算中的误差，定义为：x_diff=(_1≤i≤N)|x_i^((k+1))-x_i^((k))|

bool x_status[N]:标识每个迭代分量是否完成本次迭代

算法描述：

//进入并行域

//0号进程读取矩阵A，右端项b，迭代初始值x0并向其他进程进行广播

Do

x_old[n]=x[n];x_diff=diff;

For(i=0;i＜m;i++)//每个进程完成自己的m个分量的计算任务

{

task=i*numProcs+myID;//本次循环中需要计算的分量下标

temp=x_old[task];//保存该该分量的原始值，以便计算误差

//计算右端项

sum[task]=b[i]/A[i][task];//这里的 A[i][task]对应原迭代矩阵中的 A

[task][task],这里因为对矩阵重新分割，下标也相应发生变化

//计算中间项

For(j=task+1;j＜N;j++)

sum[task]=sum[task]-A[i][j]/A[i][task]*x_old[j];

For(int j=0;j＜i*numProcs;j++)//用之前i次循环中已经更新的分量

计算左端项（一部分）

sum[task]=sum[task]-A[i][j]/A[i][task]*x_old[j];

//全体进程依次计算自己本次迭代中的左端项

For(int j=0;j＜numProcs;j++)

{

if(myID==j)//轮到当前进程计算自己的task对应的分量并进行广播

{

for(int k=0;k＜j;k++)

sum[task]=sum[task]-A[i][i*numProcs+k]/A[i][task]*x_old[i*numProcs+k];

//计算出新的解向量

x_old[task]=sum[task];

//将解进行广播

MPI_Bcast(&x_old[task],1,MPI_DOUBLE,myID,MPI_COMM_WORLD);

}

else//循环到的分量下标不是当前进程的计算任务，则当前进程接

收新解向量

{

//接收解向量

MPI_Bcast(&x_old[i*numProcs+j],1,MPI_DOUBLE,j,MPI_COMM_

WORLD);

}

}

if(abs(x_new[i]-x_old[i])＞x_diff)

x_diff=abs(x_new[i]-x_old[i]);//计算本次迭代的误差

}

//全体进程向其他进程广播自己本次误差，通过全体进程误差最

大值判断是否需要继续迭代

MPI_Allreduce (&thread_diff,&x_diff,1,MPI_DOUBLE,MPI_MAX,

MPI_COMM_WORLD);

}While（x_diff＜diff）//达到迭代精度则结束迭代，否则继续循环

//若达到迭代精度则输出运算结果

x[n]=x_new[n]

3 数值实验和结论

可以测试串行算法的正确性，测试矩阵取为：

由于A是严格对角占优矩阵，因此G-S迭代法一定收敛，并且该方法的精确解为全1向量。

从表1中不难看出，随着矩阵维数的增大，MPI算法没有体现出并行的优势，这说明对矩阵采用连续分块的策略还是应该改进的。采用笔者的改进算法之后的加速比与MPI对比如下：

表1 任务分配模式二下MPI与openMP加速比对比

与采用任务分配模式一和openMP相比，笔者改良的迭代法均有更高的加速比，加速比稳定在1.7-1.8左右。理论上采用四个核加速比上限值为4，但是考虑到进程之间通信需要一定的时间开销，以及并行工具本身也需要时间开销，因此加速比不会随着核数增加而线性增长。

［1］陈国良.并行算法的设计与分析[M].高等教育出版社，1994.

［2］黄丽嫦.Gauss-Seidel迭代法的多核并行运算研究[J].科学技术与工程，2012（12）：2674-2692.

［3］周伟明.多核计算与程序设计[M].武汉：华中科技大学出版社，2009.

［4］刘其成，胡佳男，孙雪姣，等.并行计算与程序设计[M].北京:中国铁道出版社，2014.

［5］张武生，薛巍，李建江，等.MPI并行程序设计实例教程[M].北京：清华大学出版社，2009.

［责任编辑：李书培］