Gauss—Seidel迭代法求解線(xiàn)性方程組的并行化研究

米瑞琪　于朋

【摘 要】Gauss-Seidel迭代法在工程計(jì)算中具有廣泛應(yīng)用。高維矩陣對(duì)線(xiàn)性方程組的求解效率提出了新的要求。本文提出了兩種模式下的Gauss-Seidel并行計(jì)算方法，并對(duì)比了在不同矩陣維數(shù)下的加速效率，得到了具有較高加速比的并行迭代方法。

【關(guān)鍵詞】Gauss-Seidel；迭代法；并行；MPI；矩陣運(yùn)算

1 問(wèn)題背景與提出

在實(shí)際工程應(yīng)用中經(jīng)常需要用到求解維數(shù)較高的線(xiàn)性代數(shù)方程組，對(duì)于線(xiàn)性代數(shù)方程組的求解方法可以分成兩類(lèi)：直接法和迭代法，其中直接法得到的是方程組的真解，其求解過(guò)程為通過(guò)有限步四則運(yùn)算，常用的實(shí)現(xiàn)方式是Gauss消去方法和矩陣的三角分解方法，這種計(jì)算方式在求解過(guò)程中會(huì)產(chǎn)生大量的非零元素，存儲(chǔ)量和計(jì)算量均比較大，因此常用于求解低階、稠密線(xiàn)性方程組。對(duì)于高維線(xiàn)性方程組，可以采用并行運(yùn)算的方式進(jìn)一步提高計(jì)算效率，Gauss-Seidel迭代法具有收斂速度快、計(jì)算穩(wěn)定性好等優(yōu)點(diǎn)，是求解高維線(xiàn)性方程組的常用方法，因此本文將對(duì)Gauss-Seidel迭代法的并行化進(jìn)行研究，期望獲取較高的運(yùn)行效率和加速比。

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

任務(wù)分配模式一：連續(xù)等分矩陣法

在《并行算法的設(shè)計(jì)與分析》一書(shū)中給出了一種矩陣分配的算法，設(shè)矩陣的階為N，采用m個(gè)線(xiàn)程進(jìn)行計(jì)算，則將矩陣等分成m塊，每一塊分到的任務(wù)量為N/m，每一個(gè)線(xiàn)程計(jì)算完自己的分量之后立刻向其他的線(xiàn)程進(jìn)行廣播。

由于Gauss-Seidel迭代法中，對(duì)于右端項(xiàng)的計(jì)算和中間項(xiàng)的計(jì)算是不需要用到下一次迭代值的，因此在這一階段各個(gè)線(xiàn)程工作量相同，同時(shí)完成，問(wèn)題的關(guān)鍵在于首項(xiàng)的計(jì)算，順序等分矩陣時(shí)，以四個(gè)線(xiàn)程為例，在計(jì)算左端項(xiàng)的時(shí)候，各個(gè)線(xiàn)程的工作次序如下圖所示。其中不同的顏色代表不同線(xiàn)程的任務(wù)，每個(gè)線(xiàn)程依次計(jì)算分給自己的一個(gè)子塊，每計(jì)算完一個(gè)分量后立刻沿上方箭頭所指方向進(jìn)行廣播。

任務(wù)分配模式二：離散等分矩陣法

為了改進(jìn)任務(wù)分配模式一中的缺陷，進(jìn)一步提高計(jì)算的并行程度，同時(shí)保持任務(wù)的均衡性，筆者對(duì)模式一中的任務(wù)分配方式作出了一定程度的改進(jìn)，設(shè)線(xiàn)程數(shù)為m，則將矩陣以每m行作為一個(gè)單位，順序分配給m個(gè)線(xiàn)程，前一個(gè)線(xiàn)程完成自己的計(jì)算任務(wù)之后，立刻將計(jì)算結(jié)果廣播給其余的m個(gè)線(xiàn)程，同時(shí)開(kāi)始下一個(gè)矩陣塊中自己應(yīng)完成的任務(wù)，直到最后將自己所應(yīng)完成的N/m個(gè)分量計(jì)算完畢。任務(wù)分配模式見(jiàn)左圖，其中不同的顏色代表不同的線(xiàn)程的計(jì)算任務(wù)，每個(gè)塊右端的箭頭表示變量的廣播方向。

這樣劃分任務(wù)，每個(gè)線(xiàn)程的工作量沒(méi)有發(fā)生變化，但是第一個(gè)達(dá)到終點(diǎn)的線(xiàn)程和最后一個(gè)達(dá)到終點(diǎn)的線(xiàn)程的時(shí)間相差僅為m個(gè)分量的計(jì)算時(shí)間。假設(shè)矩陣的階為10000，采用4個(gè)線(xiàn)程進(jìn)行計(jì)算（m=4），那么第一個(gè)線(xiàn)程和最后一個(gè)線(xiàn)程到達(dá)終點(diǎn)的時(shí)間差為m=4個(gè)分量的計(jì)算時(shí)間。反觀(guān)第一種任務(wù)分配模式，假設(shè)第一個(gè)線(xiàn)程計(jì)算完自己的N/m=10000/4=2500個(gè)分量，則需要再經(jīng)過(guò)計(jì)算2500、5000、7500個(gè)分量的時(shí)間之后二、三、四號(hào)線(xiàn)程才能到達(dá)終點(diǎn)，并且隨著矩陣階的增大，這個(gè)時(shí)間差異還會(huì)繼續(xù)增大。

根據(jù)上面的分析，在理論上模式二具有更優(yōu)的時(shí)間性能。筆者也進(jìn)行了實(shí)驗(yàn)，從實(shí)際操作上證明了模式二具有更優(yōu)的時(shí)間性能。因此選擇模式二作為并行化的基本算法，基本并行算法步驟可以描述如下：

（一）末項(xiàng)/中間項(xiàng)計(jì)算的并行化

（二）首項(xiàng)的并行計(jì)算

首項(xiàng)是下三角矩陣與本次迭代的結(jié)果相乘，因此不易實(shí)現(xiàn)并行化，但是借鑒并行計(jì)算中的流水線(xiàn)作業(yè)思想，可以采用步進(jìn)和廣播的方法進(jìn)行并行化。假設(shè)核數(shù)為4，則首項(xiàng)的并行化可以設(shè)計(jì)如下：

2.2 基于MPI的多核并行算法的設(shè)計(jì)

采用任務(wù)分配模式二，采用MPI的并行算法可以描述如下：

int N：矩陣的階數(shù)

int numProcs：用于計(jì)算的進(jìn)程數(shù)量

int m：每個(gè)進(jìn)程所分得的任務(wù)量，m=N/numProcs

int myID：標(biāo)識(shí)本進(jìn)程的進(jìn)程號(hào)

int task：記錄本次循環(huán)中本線(xiàn)程計(jì)算的未知量下標(biāo)

double A[m][numProcs]：每個(gè)進(jìn)程單獨(dú)

保存屬于自己的矩陣分塊

double x_old[N]：double類(lèi)型的N維數(shù)組，用于保存本次迭代分量的數(shù)值

double sum[N]：保存迭代分量計(jì)算過(guò)程中的累加結(jié)果；

double x_diff：相鄰兩次迭代運(yùn)算中的誤差，定義為：x_diff=（_1≤i≤N）|x_i^（（k+1） ）-x_i^（（k） ） |

bool x_status[N]：標(biāo)識(shí)每個(gè)迭代分量是否完成本次迭代

算法描述：

//進(jìn)入并行域

//0號(hào)進(jìn)程讀取矩陣A，右端項(xiàng)b，迭代初始值x0并向其他進(jìn)程進(jìn)行廣播

Do

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

For（i=0；i

{

task=i*numProcs+myID；//本次循環(huán)中需要計(jì)算的分量下標(biāo)

temp=x_old[task]；//保存該該分量的原始值，以便計(jì)算誤差

//計(jì)算右端項(xiàng)

sum[task]=b[i]/A[i][task]；//這里的A[i][task]對(duì)應(yīng)原迭代矩陣中的A[task][task]，這里因?yàn)閷?duì)矩陣重新分割，下標(biāo)也相應(yīng)發(fā)生變化

//計(jì)算中間項(xiàng)

For（j=task+1；j

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

For（int j=0；j

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

//全體進(jìn)程依次計(jì)算自己本次迭代中的左端項(xiàng)

For（int j=0；j

{

if（myID==j）//輪到當(dāng)前進(jìn)程計(jì)算自己的task對(duì)應(yīng)的分量并進(jìn)行廣播

{

for（int k=0；k

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

//計(jì)算出新的解向量

x_old[task]=sum[task]；

//將解進(jìn)行廣播

MPI_Bcast（&x_old[task]，1，MPI_DOUBLE，myID，MPI_COMM_WORLD）；

}

else//循環(huán)到的分量下標(biāo)不是當(dāng)前進(jìn)程的計(jì)算任務(wù)，則當(dāng)前進(jìn)程接收新解向量

{

//接收解向量

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]）； //計(jì)算本次迭代的誤差

}

//全體進(jìn)程向其他進(jìn)程廣播自己本次誤差，通過(guò)全體進(jìn)程誤差最大值判斷是否需要繼續(xù)迭代

MPI_Allreduce（&thread_diff，&x_diff，1，MPI_DOUBLE，MPI_MAX，MPI_COMM_WORLD）；

}While（x_diff

//若達(dá)到迭代精度則輸出運(yùn)算結(jié)果

x[n]=x_new[n]

3 數(shù)值實(shí)驗(yàn)和結(jié)論

可以測(cè)試串行算法的正確性，測(cè)試矩陣取為：

由于A(yíng)是嚴(yán)格對(duì)角占優(yōu)矩陣，因此G-S迭代法一定收斂，并且該方法的精確解為全1向量。

從表1中不難看出，隨著矩陣維數(shù)的增大，MPI算法沒(méi)有體現(xiàn)出并行的優(yōu)勢(shì)，這說(shuō)明對(duì)矩陣采用連續(xù)分塊的策略還是應(yīng)該改進(jìn)的。采用筆者的改進(jìn)算法之后的加速比與MPI對(duì)比如下：

表1 任務(wù)分配模式二下MPI與openMP加速比對(duì)比

與采用任務(wù)分配模式一和openMP相比，筆者改良的迭代法均有更高的加速比，加速比穩(wěn)定在1.7-1.8左右。理論上采用四個(gè)核加速比上限值為4，但是考慮到進(jìn)程之間通信需要一定的時(shí)間開(kāi)銷(xiāo)，以及并行工具本身也需要時(shí)間開(kāi)銷(xiāo)，因此加速比不會(huì)隨著核數(shù)增加而線(xiàn)性增長(zhǎng)。

【參考文獻(xiàn)】

[1]陳國(guó)良.并行算法的設(shè)計(jì)與分析[M].高等教育出版社，1994.

[2]黃麗嫦.Gauss-Seidel迭代法的多核并行運(yùn)算研究[J].科學(xué)技術(shù)與工程，2012（12）：2674-2692.

[3]周偉明.多核計(jì)算與程序設(shè)計(jì)[M].武漢：華中科技大學(xué)出版社，2009.

[4]劉其成，胡佳男，孫雪姣，等.并行計(jì)算與程序設(shè)計(jì)[M].北京：中國(guó)鐵道出版社，2014.

[5]張武生，薛巍，李建江，等.MPI并行程序設(shè)計(jì)實(shí)例教程[M].北京：清華大學(xué)出版社，2009.

[責(zé)任編輯：李書(shū)培]