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

米瑞琪 于 朋

（中國石油大學〈華東〉，山東 青島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.

［責任編輯：李書培］