病態(tài)線性方程組的并行迭代求解

黃麗嫦，林 結(jié)，黃 潤(rùn)

（1.佛山職業(yè)技術(shù)學(xué)院基礎(chǔ)教學(xué)部，廣東佛山528137；2.佛山職業(yè)技術(shù)學(xué)院電子信息系，廣東佛山528137）

病態(tài)線性方程組的并行迭代求解

黃麗嫦1，林 結(jié)1，黃 潤(rùn)2

（1.佛山職業(yè)技術(shù)學(xué)院基礎(chǔ)教學(xué)部，廣東佛山528137；2.佛山職業(yè)技術(shù)學(xué)院電子信息系，廣東佛山528137）

分析了病態(tài)線性方程組的相關(guān)概念及判別方法，給出了一種病態(tài)線性方程組并行迭代的求解算法。算法首先對(duì)病態(tài)線性方程組的系數(shù)矩陣進(jìn)行嚴(yán)格對(duì)角占優(yōu)預(yù)處理，在此基礎(chǔ)上，用并行的Jacobi迭代法進(jìn)行多步迭代求解。新算法易于在多核架構(gòu)的微機(jī)中實(shí)現(xiàn)，且數(shù)值實(shí)驗(yàn)也驗(yàn)證了算法具有良好的收斂性和并行性。

病態(tài)線性方程組；并行計(jì)算；雅可比迭代法

線性代數(shù)是數(shù)值計(jì)算方法的基礎(chǔ)，許多科學(xué)研究和工程技術(shù)等實(shí)際問(wèn)題的求解最終都可轉(zhuǎn)化為線性代數(shù)問(wèn)題來(lái)處理，如計(jì)算插值函數(shù)與擬合函數(shù)，最小二乘數(shù)據(jù)擬合，構(gòu)造求解微分方程的差分格式等［1-2］。當(dāng)方程組的系數(shù)行列式不為零時(shí)，根據(jù)克萊姆（Cramer）法則便可以求出方程組的惟一解，然而這一理論完善的方法在高階線性方程組中因計(jì)算量龐大而變得難以應(yīng)用［3］，因此如何建立有效而實(shí)用的線性方程組求解方法，具有重要的現(xiàn)實(shí)意義。目前，線性方程組的求解方法大致可以分為兩類：直接法和迭代法［4-5］。直接法是指假設(shè)在沒(méi)有舍入誤差的條件下，能在預(yù)定的運(yùn)算次數(shù)內(nèi)求得精確解的方法，由于在實(shí)際計(jì)算過(guò)程中總存在誤差，因此直接法存在數(shù)值計(jì)算的穩(wěn)定性問(wèn)題。而迭代法則是從解的某個(gè)近似值出發(fā)，通過(guò)構(gòu)造一個(gè)無(wú)窮序列去逼近精確解的方法，迭代法因?yàn)榫哂蟹椒ê?jiǎn)單且易于實(shí)現(xiàn)的優(yōu)點(diǎn)而更具應(yīng)用價(jià)值。本文在多核架構(gòu)的微機(jī)中，基于雅可比（Jacobi）迭代法研究了病態(tài)線性方程組的并行求解問(wèn)題，并設(shè)計(jì)實(shí)現(xiàn)相應(yīng)的迭代求解算法，數(shù)值實(shí)驗(yàn)也驗(yàn)證了算法具有良好的收斂性和并行性。

1　病態(tài)線性方程組的判別方法

1.1線性方程組的病態(tài)問(wèn)題

考慮n階非奇異線性方程組

其中，A=（aij）∈Rn×n且det（A）≠0，x=（x1，x2，…，xn）∈Rn，b=（b1，b2，…，bn）∈Rn，若式（1）中的系數(shù)矩陣A和右端項(xiàng)b的微小擾動(dòng)引起了解x的巨大變化，則稱該線性方程組為病態(tài)線性方程組，否則稱其為良性線性方程組。

線性方程組的病態(tài)性強(qiáng)弱程度一般可用下式（2）所示的條件數(shù)進(jìn)行定量描述［6］

1.2病態(tài)方程組的判別方法

利用式（2）對(duì)線性方程組的病態(tài)性進(jìn)行描述，其主要困難在于對(duì)‖A-1‖V的計(jì)算，當(dāng)然，可以通過(guò)求解方程組Ayi=I（i=1，2，…，n，I為n階單位矩陣）而得到A-1=［y1，y2，…，yn］，然后再求出‖A-1‖V，然而這樣所需的計(jì)算量卻不低于求解x的計(jì)算量，因此有必要引入簡(jiǎn)便而實(shí)用的線性方程組條件數(shù)估計(jì)方法。事實(shí)上，對(duì)于式（1）所示的非奇異線性方程組而言，由于det（A）≠0，則方程組有惟一解，即

其中，Ai（i=1，2，…，n）是把系數(shù)矩陣A中第i列元素用方程組右端項(xiàng)b代替后所得的n階矩陣。若式（1）中某些元素的擾動(dòng)使得det（A）、det（Ai）的增量分別為ΔT和ΔTi，當(dāng)（det（A）+ΔT）≠0時(shí)，則此時(shí)方程組對(duì)應(yīng)的解為

設(shè)系數(shù)矩陣A=［a1，a2，…，an］，ai∈Rn，當(dāng)向量組a1，a2，…，an的線性相關(guān)程度較小時(shí)，|det（A）|將取較大的數(shù)值，所以有|det（A）|>>|det（Ai）|，從而有xi≈x'i，即擾動(dòng)對(duì)方程組的解的影響較小，故此時(shí)的線性方程組是良態(tài)的；反之，當(dāng)向量組a1，a2，…，an的線性相關(guān)程度較大時(shí)，|det（A）|將取較小的數(shù)值，此時(shí)擾動(dòng)對(duì)方程組的解的影響較大，即線性方程組是病態(tài)的。

從上述分析可知，通過(guò)計(jì)算系數(shù)矩陣A的行列式便可判斷方程組的病態(tài)性，這里引入如下引理。

引理1［7］對(duì)于非奇異線性方程組Ax=b而言，設(shè)det（A）=Σ（-1）ta1p1a2p2…anpnK=max（a1p1a2p2…anpn），若K>>|det（A）|，則方程組是病態(tài)的。

引理1中，p1，p2，…，pn為系數(shù)矩陣A的元素的列標(biāo)號(hào)，t為排列p1，p2，…，pn的逆序數(shù)，Σ表示對(duì)n個(gè)數(shù)的所有排列p1，p2，…，pn取和［8］。

2　病態(tài)線性方程組的并行迭代求解算法設(shè)計(jì)

2.1線性方程組的并行迭代求解原理

利用迭代法求解方程組，就是構(gòu)造一個(gè)無(wú)限的向量序列，使它的極限是方程組的解向量。為了更好地發(fā)揮目前主流的微機(jī)多核計(jì)算功能，本文擬采用具有良好并發(fā)性能的Jacobi迭代法來(lái)研究病態(tài)線性方程組的求解問(wèn)題。

由于迭代法不能通過(guò)有限次算術(shù)運(yùn)算求得方程組的精確解，因此，收斂性與精度控制是迭代法的核心問(wèn)題。對(duì)于Jacobi迭代法，有如下引理。

（1）A為非奇異；

（2）Jacobi迭代法收斂。

為確保病態(tài)線性方程組Ax=b的Jacobi迭代收斂，根據(jù)引理2中的收斂性質(zhì)，可對(duì)病態(tài)線性方程組的系數(shù)矩陣A進(jìn)行如下分解

然后在方程組左右兩邊加上一個(gè)對(duì)角矩陣M，有

對(duì)式（6）進(jìn)行移項(xiàng)及整理，有

從而得到如下的迭代形式

于是，相應(yīng)的Jacobi迭代為

2.2算法的設(shè)計(jì)及實(shí)現(xiàn)

從式（9）可知，迭代表達(dá)式x1（k+1），x2（k+1），…，xn（k+1）之間沒(méi)有關(guān)聯(lián)，可將它們平均分布在CPU中的各個(gè)核中計(jì)算。因此，可設(shè)計(jì)如下病態(tài)線性方程組的Jacobi并行迭代算法。

算法名稱：病態(tài)線性方程組的Jacobi并行迭代算法。

輸入：病態(tài)線性方程組的矩陣系數(shù)An×n，常數(shù)向量bn×1，迭代精度e，初始解向量xn×1；

輸出：病態(tài)線性方程組的解向量xn×1。

算法描述：（設(shè)CPU的核數(shù)為P）

（1）初始化及預(yù)處理。

1）從系數(shù)矩陣A中分離出P，Q矩陣，其中P=aii｛｝，Q=aij｛ ｝，i≠j，1≤i，j≤n；

（2）并行迭代計(jì)算過(guò)程。

do｛

2）各計(jì)算核完成當(dāng)前迭代計(jì)算后，將自己所計(jì)算分量的新值廣播到所有計(jì)算核中；

｝while（diff＞e）

（3）輸出解向量xn×1。

3　算法的性能測(cè)試

為了檢驗(yàn)上述設(shè)計(jì)算法的有效性，本文在Intel Xeon E54504核3.0 GHz CPU（每個(gè)計(jì)算核的一級(jí)緩存各由32 KB數(shù)據(jù)緩存和32 KB指令緩存組成，二級(jí)緩存容量為12 MB，）、KingSton DDR3 1 333 MHZ 4 GB內(nèi)存及Red Hat Enterprise Linux 6.1操作系統(tǒng)的環(huán)境中進(jìn)行了相關(guān)的數(shù)值試驗(yàn)，程序采用OpenMP和C++語(yǔ)言進(jìn)行編寫。

分別用Jacobi迭代法和本文改進(jìn)的算法對(duì)不同階數(shù)的Hilbert經(jīng)典病態(tài)問(wèn)題進(jìn)行試驗(yàn)，其中，Hilbert病態(tài)問(wèn)題的系數(shù)矩陣如下式（10）所示［10］

而b=An［1 1…1］T，在計(jì)算精度e=1×10-12的條件下，得到如表1所示的計(jì)算結(jié)果。

表1　Jacobi迭代法與本文改進(jìn)算法的計(jì)算性能比較

實(shí)驗(yàn)數(shù)據(jù)表明，在同樣計(jì)算精度的條件下，由于加入有效的預(yù)處理，因此本文的改進(jìn)算法在計(jì)算效能上得到了有效改善。

4　結(jié)語(yǔ)

本文介紹了病態(tài)線性方程組的有關(guān)概念以及判別方法，在基于多核架構(gòu)的微機(jī)中提出了求解病態(tài)線性方程組的并行迭代算法，由數(shù)值試驗(yàn)得知，算法的改善效果是明顯的。

［1］丁麗娟，程杞元.數(shù)值計(jì)算方法［M］.2版.北京：北京理工大學(xué)出版社，2005：12-74.

［2］李大明.數(shù)值線性代數(shù)［M］.北京：清華大學(xué)出版社，2010：1-36.

［3］呂全義.線性代數(shù)［M］.北京：清華大學(xué)出版社，2015：1-21.

［4］陳長(zhǎng)軍.數(shù)值計(jì)算理論與實(shí)現(xiàn)［M］.武漢：華中科技大學(xué)出版社，2014：32-66.

［5］沈艷.高等數(shù)值計(jì)算［M］.北京：清華大學(xué)出版社，2014：20-60.

［6］徐樹(shù)方.矩陣計(jì)算的理論與方法［M］.北京：北京大學(xué)出版社，1995：54-75.

［7］鄭孝勇，張東俊.病態(tài)線性方程組的判定方法［J］.數(shù)學(xué)的實(shí)踐與認(rèn)識(shí)，2006，36（9）：166-169.

［8］居余馬.線性代數(shù)簡(jiǎn)明教程［M］.北京：清華大學(xué)出版社，2004：1-95.

［9］李慶揚(yáng).數(shù)值計(jì)算原理［M］.北京：清華大學(xué)出版社，2000：140-241.

［10］林勝良.病態(tài)線性方程組解法研究［D］.杭州：浙江大學(xué)，2005：51-54.

【責(zé)任編輯：王桂珍foshanwgzh@163.com】

Parallel iterative method for solving ill-conditioned linear equations

HUANGLi-chang1，LINJie1，HUANGRun2

（1.Department ofBasic Teaching，F(xiàn)oshan Polytechnic，F(xiàn)oshan 528137，China；2.Department ofElectronic Information，F(xiàn)oshan Polytechnic，F(xiàn)oshan 528137，China）

This paper analyzes the concepts and discriminatimg methods of ill-conditioned linear equations，presents a parallel iterative algorithm.The pretreatment algorithm based on ill conditioned coefficient matrix of the linear equations is strictlydiagonallydominant，based on multi step iterative method for solvingthe equations after pretreatment with parallel Jacobi iterative method.The new algorithm is easy to implement on computer multi-core architecture，and the numerical experiments show that the algorithm has good convergence and parallelism.

ill-conditioned linear equations；parallel computation；Jacobi iteration

O241.6

A

1008-0171（2016）04-0013-05

2016-01-20

佛山職業(yè)技術(shù)學(xué)院校級(jí)科研基金資助項(xiàng)目（KY2013Y07）

黃麗嫦（1962-），女，廣東中山人，佛山職業(yè)技術(shù)學(xué)院講師。