块三对角线性方程组的并行直接解法

1、602009，45（3）ComputerEngineeringandApplications计算机工程与应用块三对角线性方程组的并行直接解法樊艳红，吕全义，聂玉峰FANYan-hong，LVQuan-yi，NIEYu-feng西北工业大学应用数学系，西安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提出了分布式环境下求解块三对角线性方程组的一种并行算法，该算法充分利用系数矩阵结构的特殊性，通过对系数矩摘要：阵进行适当分解及近似处理，使算法只在相邻处理机间通信两次。并从理论上给出了算法有效的一个充分条件。最后，在HP结果表明，实算与理论是一致的，并行性也很好。rx2600集群上进行了数值实验，块三对角线性方程组；矩阵分解；并行算法；并

5、行效率；关键词：HPrx2600集群DOI：10.3778/j.issn.1002-8331.2009.03.017文章编号：10028331（2009）03-0060-04文献标识码：A中图分类号：TP301（2）mmmmmmmmmmmmmmmmmmmmmmmm1引言在科学与工程问题中经常遇到的许多微分方程，经过适当差分或有限元离散而形成系数矩阵是块三对角的线性方程组，它们的求解是高性能并行计算的重要课题之一。目前针对求解块三对角线性方程组的并行算法的研究已经有了一些成果1-5，文献5通过对系数矩阵进行了分解与近似处理，构造了具有良好的并行性的算法。借鉴文献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维列向量。假设处理机台数为p，记第i台处理机为p（，p），设ii=1，m=pk（k2，kZ），若将系数矩阵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性，只需相邻处理机之间进行两次数据传递，有效减小了通讯次数。下面介绍一下具体的计算步骤：（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）时因各处理机之间要相互传递数据而使并行性降低，为了只在相邻处理机之间有数据传递，提高其并行性，现给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并行实现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分布式存储环境下进行并行计算。注若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的算法，又两者的载荷平衡度是相同的，所以该算法优于文献5的算法。（i）（i）其中Q1=O，Sp-1=O。所以max（Qi+Si），越小，精度越imax（Qi+Si）高。首先考虑A满足什么条件时，能i尽量的小。定理1若Ai非奇异，存在常数li>0（i=1，2，3），使得AiBi<l1，AiCi<l2（i=1，m），且1-l1-l2l

21、3，则对任意的>0，存在正整数K，当k>K时，有：max（Qi+Si）<（i=1，p-1）i-1-15数值算例应用此算法在HPrx2600集群上进行数据试验，以下算例为了与文献5算法的结果是在HPrx2600集群上计算得到的，进行比较，文献5算法的结果也是在同一个并行环境下重新计算得到的。例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时两种方法的计算结果处理机台数时间/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处理机台数时间/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时用本文方法的计算结果处理机台数时间/s本文方法加速比效率10.527320.27341.92880.964440.15623.37580.844080.07037.50070.9376Ax-f2

25、.1316e-141.4921e-131.4921e-132.6290e-13算例结果分析：（1）从三个例子看出，该算法适合在分布式存储环境下并行求解块三对角线性方程组，并行效率均高于80%，特别地，系数矩阵的阶数越大时，效率会越高，精度也会更高。（2）当系数矩阵A不满足定理的条件时，算法仍然有效；算例验证了定理的条件只是一个充分条件。（3）从算例1和算例2可以看出，解的误差与基本是一致的。（4）该算法优于文献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时两种方法的计算结果处理机台数时间/s本文方法加速比效率时间/s文献5方法加速比效率10.039117.695324.09381.87970.93995.76951.74000.870042.08823.68