塊三對(duì)角線性方程組的并行直接解法

1、602009，45（3）ComputerEngineeringandApplications計(jì)算機(jī)工程與應(yīng)用塊三對(duì)角線性方程組的并行直接解法樊艷紅，呂全義，聶玉峰FANYan-hong，LVQuan-yi，NIEYu-feng西北工業(yè)大學(xué)應(yīng)用數(shù)學(xué)系，西安710072DepartmentofAppliedMathematics，NorthwesternPolytechnicalUniversity，Xian710072，ChinaE-mail：fanyanhongFANYan-hong，LVQuan-yi，NIEYu-feng.Improvedparallelalgorithmforsolvin

2、gblock-tridiagonallinearequations.ComputerEngineeringandApplications，2009，45（3）：60-63.Abstract：Aparallelalgorithmforblock-tridiagonallinearequationsondistributed-memorymulti-computersispresented.Makingfulluseofthespecialstructureofthecoefficientmatrix，thealgorithmisbasedondecomposingthecoefficientma

3、trixproperlyandapproximatelydisposingthematrix.Thecommunicationonlyneedstwicebetweentheadjacentprocessors.Intheory，thispapergivesasufficientconditionabouteffectivityofthisalgorithm.Finally，somenumericalresultsonHPrx2600clustershowthatpracticecomputingisconsistentwiththeory.Thealgorithmsparallelismis

4、preferable.Keywords：block-tridiagonallinearequations；decompositionofthematrix；parallelalgorithm；parallelefficiency；HPrx2600cluster提出了分布式環(huán)境下求解塊三對(duì)角線性方程組的一種并行算法，該算法充分利用系數(shù)矩陣結(jié)構(gòu)的特殊性，通過(guò)對(duì)系數(shù)矩摘要：陣進(jìn)行適當(dāng)分解及近似處理，使算法只在相鄰處理機(jī)間通信兩次。并從理論上給出了算法有效的一個(gè)充分條件。最后，在HP結(jié)果表明，實(shí)算與理論是一致的，并行性也很好。rx2600集群上進(jìn)行了數(shù)值實(shí)驗(yàn)，塊三對(duì)角線性方程組；矩陣分解；并行算法；并

5、行效率；關(guān)鍵詞：HPrx2600集群DOI：10.3778/j.issn.1002-8331.2009.03.017文章編號(hào)：10028331（2009）03-0060-04文獻(xiàn)標(biāo)識(shí)碼：A中圖分類號(hào)：TP301（2）mmmmmmmmmmmmmmmmmmmmmmmm1引言在科學(xué)與工程問(wèn)題中經(jīng)常遇到的許多微分方程，經(jīng)過(guò)適當(dāng)差分或有限元離散而形成系數(shù)矩陣是塊三對(duì)角的線性方程組，它們的求解是高性能并行計(jì)算的重要課題之一。目前針對(duì)求解塊三對(duì)角線性方程組的并行算法的研究已經(jīng)有了一些成果1-5，文獻(xiàn)5通過(guò)對(duì)系數(shù)矩陣進(jìn)行了分解與近似處理，構(gòu)造了具有良好的并行性的算法。借鑒文獻(xiàn)5的求解思想，構(gòu)造出了一個(gè)并行效率

6、更高的并行迭代算法。2算法考慮塊三對(duì)角線性方程組Ax=f，即：mmmmA1B1x1mf1mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm2mmmmmmmm-1mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmA=FT其中：m（）mD1mmCkAkBkmmm（D2）mmmC2kA2kB2kmF=mm（Dp-1）CA（p-1）k（p-1）kA（i-1）k+1C（i-1）k+2B（i-1）k+1A（i-1）k+2B（i-1）k+2B（p-1）k（Dp）mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmC

7、2A2B2x2xm-1ffDi=（1）=Cm-1Am-1Bm-1mmmmmmmmmmmmmmmm，i=1，p-1Cik-2B（i-1）k+1A（i-1）k+2B（i-1）k+2Aik-2Cik-1Bik-2Aik-1CmAmxmf這里Ai是nm階方陣，Bi和Ci分別是ni×ni+1和ni×ni-1矩陣，xi，fi均為ni維列向量。假設(shè)處理機(jī)臺(tái)數(shù)為p，記第i臺(tái)處理機(jī)為p（，p），設(shè)ii=1，m=pk（k2，kZ），若將系數(shù)矩陣A分解成兩矩陣的乘積：Dp=mmmmmmmmmmmmmmmmmA（i-1）k+1C（i-1）k+2Cm-1Am-1CmBm-1Am基金項(xiàng)目：陜西省自然科

8、學(xué)基金（theNaturalScienceFoundationofShaanxiProvinceofChinaunderGrantNo.2006A05）。作者簡(jiǎn)介：樊艷紅（1982-），女，碩士研究生，主要研究方向：信息處理中的快速與并行算法；呂全義（1963-），女，碩士生導(dǎo)師，副教授，主要研究方向：并行算法；聶玉峰（1968-），男，碩士生導(dǎo)師，副教授，主要研究方向：有限元及其應(yīng)用。樊艷紅，呂全義，聶玉峰：塊三對(duì)角線性方程組的并行直接解法設(shè)F可逆，則T=FA，即：pppppppppppppppppppppppppppppppppp2009，45（3）pppppppppppppppppppp

9、ppppppppppp61-1p（）OpOpppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppppI1H1R1G2-Q2I2-S1H2R2-S2OQ2S1OS2O（O）OT=T=Qp-2OOQp-1（O）Sp-2OOOO-Qp-2Rp-2Gp-1-Qp-1Ip-1-Sp-2Hp-1Rp-1GpIp則：ppppppppppppppppppppppppppppppppI1H1R1G2I2H2R2×，i=1，p-1的單位矩nnI為×的單位矩陣。G為陣，nn×n，i=2，p-1的矩陣，設(shè)為，Z，ZnG為

10、×n的矩陣，設(shè)為H為Z，Z；n×n，i=1，p-1的矩陣，設(shè)為；且Y，Ynk-1k-1Ii為s=1（i-1）k+ss=1（i-1）k+sk-1k-1s=1s=1T+T=k-1TTTRp-2Gp-1Ip-1Hp-1Rp-1GpIps=1（i-1）k+sk-1（i-1）k1，ik-1，iTTTps=1（p-1）k+sp-11，pk，pik-1TTT可通過(guò)求解方程組：）（x+x）=u（T+T（9）來(lái)近似代替式（7），其中x滿足Tx=u，而式（9）具有較好的并行s=1（i-1）k+siki1，k-1，i它們滿足下列方程組：mmmmmmmmmmmmmmmmmmmmmmmmmmmmmm

11、mmmmmmmmmmmmA（i-1）k+1C（i-1）k+2B（i-1）k+1Bik-2Aik-1Cik-1B（i-1）k+1mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmY1，iY2，iYk-1，iZ1，iZ2，immmmmmmmmmmmmmmmmmmmmmmmmmmm=mmmmmmmmmmmmmOOBik-1mmmmmmmmmmmmm性，只需相鄰處理機(jī)之間進(jìn)行兩次數(shù)據(jù)傳遞，有效減小了通訊次數(shù)。下面介紹一下具體的計(jì)算步驟：（3）（1）求解方程組（6）Diui=fi，i=1

12、，p。同時(shí)求解方程組（3）（5），得解分別為：mmmmmmmmmmmmmA（i-1）k+1C（i-1）k+2Bik-2Aik-1=Cik-1B（p-1）k+1Zk-1，immmmmmmmmmmmmmmmmmmmmmmmmmmC（i-1）k+1OOYj，（-1）i=（4）Zj，（-1）i=k-1-j軍（A儀s=jjk-1-1（i-1）k+sBj=k-1，1）（i-1）k+s）（p（p-1）k+s（p-1）k+sipppppppppppppppppppppppppppppppp（10）j-1贊（A儀s=jj-1（i-1）k+sCj=1，k-1）（i-1）k+s）（11）Cj=1，k）（p-1）k+

13、s）（A（p-1）k+1C（p-1）k+2Z1，pZ2，pZk，pBm-1Am-1=Cm-1mmmmmmmmmmmmmC（p-1）k+1OOmmmmmmmmmmmmmZj，（-1）p=其中：（5）軍軍軍軍軍軍軍軍軍軍軍軍軍軍軍軍軍軍軍軍軍軍軍軍軍軍軍j-1贊（A儀s=j-1（p-1）k+s軍A（i-1）k+1=A（i-1）k+1（j=2，k-1）-1軍軍AA）k+j-1B（i-1）k+j=A（i-1）k+j-C（i-1）k+j（i-1（i-1）k+j-1贊=AAikik（j=k-1，1）-1贊贊A（i-1）k+j=A（i-1）k+j-B（i-1）k+jA（i-1）k+j-1C（i-1）k+j-

14、1贊=AApkpk（j=1，k）-1贊贊A（p-1）k+j=A（p-1）k+j-B（p-1）k+jA（p-1）k+j-1C（p-1）k+j-1（i）（i）（i）（i+1）于是Qi=AikCikZk-1，i=2，p-1，Ri=Aik（Aik-CikYik-1，i=i，i-BikZ1，i+1）1，p-1，Si=AikBikY1，i=1，p-2，方程組（1）的右端項(xiàng)和i+1，所求解向量亦作相應(yīng)劃分。則求解方程組（1）等價(jià)于求解線性方程組：Fu=fTx=u（6）（7）-1（2）求解（Aik-CikYk-1，xk=fk-CikUk-1-BikU1i-BikZ1，i+1）（3）計(jì)算xj=Uj-Yj，1xk

15、（i）（i）（1）（1）（1）但在求解式（7）時(shí)因各處理機(jī)之間要相互傳遞數(shù)據(jù)而使并行性降低，為了只在相鄰處理機(jī)之間有數(shù)據(jù)傳遞，提高其并行性，現(xiàn)給T一個(gè)擾動(dòng)T：T=（Tij）m×mp其它子塊均為零矩陣。即：（8）其中Tij為ni階方陣Tik，i=1，p-2，T（i-1）k，（i+1）k=Si，ik=Qi，（12）（j=1，k-1）（i-1）（13）（14）（15）xj=Uj-Yj，ixk-Zj，ixkxj=Uj-Zj，pxk（p）（p）（p-1）（i）（j=1，k-1）（j=1，622009，45（3）ComputerEngineeringandApplications計(jì)算機(jī)工程與應(yīng)用

16、l1AikBikY1，<j+1l1+l3-13并行實(shí)現(xiàn)3.1存儲(chǔ)方法將Di，Bik，C（i-1）k+1，f（i-1）k+（，k）分配到第i臺(tái)處理機(jī)jj=1，P（，p）中。ii=1，111k于是，對(duì)任意給定的>0，存在K1=意的k>K1時(shí)，有：Si</2，i=1，p-2ln-ln2N，對(duì)任lnl-ln（l+l）3.2計(jì)算過(guò)程（1）第（ii=1，p-1）臺(tái)處理機(jī)Pi求解方程組Diui=fi及（i）T（i）TTT111，式（3）、（4），得到（，U1），（Uk-1）1Y1，Yk-1，i，iTT11；第1臺(tái)處理機(jī)P1求解方程組D1u1=f1及方程Z1，Zk-1，i，iTTT同理

17、可證，對(duì)任意給定的>0，存在K2=對(duì)任意的k>K2時(shí)，有：Qi</2，i=2，p-1ln-ln2N，lnl-ln（l+l）221）得到（3U1），（Uk-1）（1）TT（1）T111，；第p臺(tái)處理Y1，Yk-1，1，1TT（p）T（p）TTTT故，對(duì)任意給定的>0，存在K=maxK1，K2N，對(duì)任意的k>K時(shí)，有：max（Qi+Si）</2+/2=，i=1，p-1i1機(jī)Pp求解方程組Dpup=fp及方程（5）得到（，U1），（Uk-1）111。Z1，Zk，p，pTT（i）由于該定理的條件與文獻(xiàn)5中定理1的相同，所只要m足夠大，即m以由文獻(xiàn)5中的定理2可知，（

18、2）第i臺(tái)處理機(jī)Pi將Z1，U1傳給處理機(jī)Pi-1。由方程組i，（12）得到xk，再將xk的值傳給Pi+1；）在P1中計(jì)算式（13），同時(shí)在P（，p-1）中計(jì)算式（3ii=2，（14）且在Pp中計(jì)算式（15），得到方程組（1）的解。由上可知，該算法在計(jì)算過(guò)程中只需并行通信兩次，因而具有良好的并行性和可擴(kuò)展性，很適合在MIMD分布式存儲(chǔ)環(huán)境下進(jìn)行并行計(jì)算。注若m不能被p整除，則一些處理機(jī)按順序存m/p+1行塊，另一些存m/p個(gè)行塊，同時(shí)將xi、fi中相應(yīng)的向量存入對(duì)應(yīng)的處理機(jī)中，達(dá)到處理機(jī)的負(fù)載大體均衡，以減少等待時(shí)間。（i）（i）ln-ln（1+）Tp通過(guò)求解式（6）、（9）得到原+1時(shí)，ln

19、l1-ln（l1+l3）-1方程組的近似解x+x滿足x。因而得到與文獻(xiàn)5x相同的算法有效的充分條件。4.2運(yùn)算量比較在此，將該算法的計(jì)算量與文獻(xiàn)5的算法進(jìn)行比較。（1）由上面的算法可知，在處理機(jī)P（，p-1）計(jì)算ii=1，×1的矩陣，而文獻(xiàn)5中111nn×1的矩陣，相應(yīng)的u、f也比文獻(xiàn)11nk-1k-1s=1（i-1）k+ss=1（i-1）k+sk-1s=1（i-1）k+siiDiui=fi時(shí)，Di為Di為4誤差分析及運(yùn)算量比較4.1誤差分析xTTT由文獻(xiàn)5，而-1x1-TTT由式（10）可知Tmax（Qi+Si）（i=1，p-1）i-11nk-1s=1（i-1）k+s5少

20、一個(gè)nik階子塊；（2）前p-1臺(tái)處理機(jī)只要求解nik階方程組（12）就可得到xk的值，而文獻(xiàn)5需要求解一個(gè)nik+nik+1階方程組才能得到xk的值。（3）兩者都需要兩次并行通信，且該算法的通信量低于文獻(xiàn)5的算法。由此可見(jiàn)，該算法的計(jì)算量低于文獻(xiàn)5的算法，又兩者的載荷平衡度是相同的，所以該算法優(yōu)于文獻(xiàn)5的算法。（i）（i）其中Q1=O，Sp-1=O。所以max（Qi+Si），越小，精度越imax（Qi+Si）高。首先考慮A滿足什么條件時(shí)，能i盡量的小。定理1若Ai非奇異，存在常數(shù)li>0（i=1，2，3），使得AiBi<l1，AiCi<l2（i=1，m），且1-l1-l2l

21、3，則對(duì)任意的>0，存在正整數(shù)K，當(dāng)k>K時(shí)，有：max（Qi+Si）<（i=1，p-1）i-1-15數(shù)值算例應(yīng)用此算法在HPrx2600集群上進(jìn)行數(shù)據(jù)試驗(yàn)，以下算例為了與文獻(xiàn)5算法的結(jié)果是在HPrx2600集群上計(jì)算得到的，進(jìn)行比較，文獻(xiàn)5算法的結(jié)果也是在同一個(gè)并行環(huán)境下重新計(jì)算得到的。例1對(duì)于形如式（1）的線性方程組：1111111111111111軍B<l1。證明由文獻(xiàn)5中定理1的證明可知，Aiil1+l3-1因此：Y1，（-1）i+1=軍儀As=1-1k-1k-2-1ik+sBik+s<ll+l111k-1-14-1-1-14Si=-AikBikY1，=i

22、+1）（i=1，2，m）。AiBi，iC=0.5（i=1，2，-111s=1Ai=軍（A儀k-14-1-14-1-1ik+sBik+s）1111111111111111ni×ni，Bi=Ci=-I，右端項(xiàng)fi=1111111111111111111111111111111111ni×1樊艷紅，呂全義，聶玉峰：塊三對(duì)角線性方程組的并行直接解法，m），Bm=C1=O，當(dāng)T=max（Qi+Si）=1.0e-i2009，45（3）632222m=500，ni=10時(shí)的計(jì)算結(jié)果見(jiàn)表2。表1例1m=20000時(shí)兩種方法的計(jì)算結(jié)果處理機(jī)臺(tái)數(shù)時(shí)間/s本文方法加速比效率時(shí)間/s文獻(xiàn)5方法加速

23、比效率25.5527125.4844215.07421.69060.845316.17191.58010.790047.71093.30500.82628.32813.06830.767183.44537.39690.92463.96486.44490.8056Ci=-122-12其中Bm=C1=O，在m=20000，ni=15時(shí)的計(jì)算結(jié)果見(jiàn)表5。表5處理機(jī)臺(tái)數(shù)時(shí)間/s本文方法加速比效率時(shí)間/s文獻(xiàn)5方法加速比效率98.3867例3，m=20000時(shí)兩種方法的計(jì)算結(jié)果198.8224260.28121.63940.819760.51951.62570.8129425.14453.93020.9

24、82530.17193.26090.8152813.07817.55630.944515.01956.55060.8188Ax-f3.5527e-143.5527e-143.5527e-143.5527e-14Ax-f3.5527e-143.7303e-143.7303e-143.7303e-14Ax-f2.1316e-141.8474e-131.8474e-132.1316e-14表2例1m=500時(shí)用本文方法的計(jì)算結(jié)果處理機(jī)臺(tái)數(shù)時(shí)間/s本文方法加速比效率10.527320.27341.92880.964440.15623.37580.844080.07037.50070.9376Ax-f2

25、.1316e-141.4921e-131.4921e-132.6290e-13算例結(jié)果分析：（1）從三個(gè)例子看出，該算法適合在分布式存儲(chǔ)環(huán)境下并行求解塊三對(duì)角線性方程組，并行效率均高于80%，特別地，系數(shù)矩陣的階數(shù)越大時(shí)，效率會(huì)越高，精度也會(huì)更高。（2）當(dāng)系數(shù)矩陣A不滿足定理的條件時(shí)，算法仍然有效；算例驗(yàn)證了定理的條件只是一個(gè)充分條件。（3）從算例1和算例2可以看出，解的誤差與基本是一致的。（4）該算法優(yōu)于文獻(xiàn)5的算法，與前面的分析是一致的。Ax-f3.5527e-143.5527e-143.5527e-144.0856e-7例2對(duì)于形如式（1）的線性方程組：41-201Ai=，Bi=，Ci=-341-20T-3，fi=3，0，1，11，3，2m×1（i=1，2，m），其中，Bm=C1=O。在m=500000時(shí)的計(jì)算結(jié)果見(jiàn)表3，在m=100時(shí)的誤差見(jiàn)表4。表3例2，m=500000時(shí)兩種方法的計(jì)算結(jié)果處理機(jī)臺(tái)數(shù)時(shí)間/s本文方法加速比效率時(shí)間/s文獻(xiàn)5方法加速比效率10.039117.695324.09381.87970.93995.76951.74000.870042.08823.68