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

1、602009，45（3）ComputerEngineeringandApplications計算機工程與應用塊三對角線性方程組的并行直接解法樊艷紅，呂全義，聶玉峰FANYan-hong，LVQuan-yi，NIEYu-feng西北工業(yè)大學應用數(shù)學系，西安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)境下求解塊三對角線性方程組的一種并行算法，該算法充分利用系數(shù)矩陣結構的特殊性，通過對系數(shù)矩摘要：陣進行適當分解及近似處理，使算法只在相鄰處理機間通信兩次。并從理論上給出了算法有效的一個充分條件。最后，在HP結果表明，實算與理論是一致的，并行性也很好。rx2600集群上進行了數(shù)值實驗，塊三對角線性方程組；矩陣分解；并行算法；并

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

6、更高的并行迭代算法。2算法考慮塊三對角線性方程組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ù)為p，記第i臺處理機為p（，p），設ii=1，m=pk（k2，kZ），若將系數(shù)矩陣A分解成兩矩陣的乘積：Dp=mmmmmmmmmmmmmmmmmA（i-1）k+1C（i-1）k+2Cm-1Am-1CmBm-1Am基金項目：陜西省自然科

8、學基金（theNaturalScienceFoundationofShaanxiProvinceofChinaunderGrantNo.2006A05）。作者簡介：樊艷紅（1982-），女，碩士研究生，主要研究方向：信息處理中的快速與并行算法；呂全義（1963-），女，碩士生導師，副教授，主要研究方向：并行算法；聶玉峰（1968-），男，碩士生導師，副教授，主要研究方向：有限元及其應用。樊艷紅，呂全義，聶玉峰：塊三對角線性方程組的并行直接解法設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的矩陣，設為，Z，ZnG為

10、×n的矩陣，設為H為Z，Z；n×n，i=1，p-1的矩陣，設為；且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可通過求解方程組：）（x+x）=u（T+T（9）來近似代替式（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性，只需相鄰處理機之間進行兩次數(shù)據(jù)傳遞，有效減小了通訊次數(shù)。下面介紹一下具體的計算步驟：（3）（1）求解方程組（6）Diui=fi，i=1

12、，p。同時求解方程組（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）的右端項和i+1，所求解向量亦作相應劃分。則求解方程組（1）等價于求解線性方程組：Fu=fTx=u（6）（7）-1（2）求解（Aik-CikYk-1，xk=fk-CikUk-1-BikU1i-BikZ1，i+1）（3）計算xj=Uj-Yj，1xk

15、（i）（i）（1）（1）（1）但在求解式（7）時因各處理機之間要相互傳遞數(shù)據(jù)而使并行性降低，為了只在相鄰處理機之間有數(shù)據(jù)傳遞，提高其并行性，現(xiàn)給T一個擾動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計算機工程與應用

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

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

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

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

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

21、3，則對任意的>0，存在正整數(shù)K，當k>K時，有：max（Qi+Si）<（i=1，p-1）i-1-15數(shù)值算例應用此算法在HPrx2600集群上進行數(shù)據(jù)試驗，以下算例為了與文獻5算法的結果是在HPrx2600集群上計算得到的，進行比較，文獻5算法的結果也是在同一個并行環(huán)境下重新計算得到的。例1對于形如式（1）的線性方程組：1111111111111111軍B<l1。證明由文獻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，右端項fi=1111111111111111111111111111111111ni×1樊艷紅，呂全義，聶玉峰：塊三對角線性方程組的并行直接解法，m），Bm=C1=O，當T=max（Qi+Si）=1.0e-i2009，45（3）632222m=500，ni=10時的計算結果見表2。表1例1m=20000時兩種方法的計算結果處理機臺數(shù)時間/s本文方法加速比效率時間/s文獻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時的計算結果見表5。表5處理機臺數(shù)時間/s本文方法加速比效率時間/s文獻5方法加速比效率98.3867例3，m=20000時兩種方法的計算結果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ù)時間/s本文方法加速比效率10.527320.27341.92880.964440.15623.37580.844080.07037.50070.9376Ax-f2

25、.1316e-141.4921e-131.4921e-132.6290e-13算例結果分析：（1）從三個例子看出，該算法適合在分布式存儲環(huán)境下并行求解塊三對角線性方程組，并行效率均高于80%，特別地，系數(shù)矩陣的階數(shù)越大時，效率會越高，精度也會更高。（2）當系數(shù)矩陣A不滿足定理的條件時，算法仍然有效；算例驗證了定理的條件只是一個充分條件。（3）從算例1和算例2可以看出，解的誤差與基本是一致的。（4）該算法優(yōu)于文獻5的算法，與前面的分析是一致的。Ax-f3.5527e-143.5527e-143.5527e-144.0856e-7例2對于形如式（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時的計算結果見表3，在m=100時的誤差見表4。表3例2，m=500000時兩種方法的計算結果處理機臺數(shù)時間/s本文方法加速比效率時間/s文獻5方法加速比效率10.039117.695324.09381.87970.93995.76951.74000.870042.08823.68