求解块三对角方程组的一种并行策略

求解块三对角方程组的一种并行策略

462011，47（13）ComputerEngineeringandApplications计算机工程与应用

求解块三对角方程组的一种并行策略

段治健1，杨永1，吕全义2，马欣荣3

DUANZhijian1，YANGYong1，LVQuanyi2，MAXinrong3

1.西北工业大学翼型叶栅空气动力学国防科技重点实验室，西安7100722.西北工业大学应用数学系，西安7100723.咸阳师范学院数学系，陕西咸阳712000

1.NationalKeyLabofAerodynamicDesignandResearch，NorthwesternPolytechnicalUniversity，Xi’an710072，China2.DepartmentofAppliedMathematics，NorthwesternPolytechnicalUniversity，Xi’an710072，China3.DepartmentofMathematics，XianyangNormalUniversity，Xianyang，Shaanxi712000，China

DUANZhijian，YANGYong，LVQuanyi，puterEn-gineeringandApplications，2011，47（13）：46-49.

Abstract：Thispaperfocusesonaparalleliterativemethodforsolvingblock-tridiagonallinearsystemsondistributed-memo-rymulti-computers.ThroughchoosingthebaseofsubspacebasedonGalerkintheory，thecommunicationonlyneedtwicebe-tweentheadjacentprocessorsperiterationstep.Furthermore，thesufficientconditionforconvergenceisgivenwhenthecoeffi-cientmatrixAisasymmetricpositivedefinitematrix.Finally，thenumericalexperimentsimplementedonHPrx2600clusterindicatethatthealgorithm’sparallelaccelerationratesandefficiencyarehigherthanthemulti-splittingmethod’s.Keywords：block-tridiagonallinearsystems；Galerkintheory；HPrx2600cluster；parallelism摘

要：提出了一种在MIMD分布式存储环境下求解块三对角线性方程组的并行算法。基于Galerkin原理适当取基构造算法，

使整个计算过程只在相邻处理机间通信两次，并给出了系数矩阵为对称正定矩阵时算法收敛的条件。在HPrx2600集群系统上进行的数值计算结果表明该算法与多分裂方法相比具有较高的加速比和并行效率。关键词：块三对角线性方程组；Galerkin原理；HPrx2600集群；并行性DOI：10.3778/j.issn.1002-8331.2011.13.014

文章编号：1002-8331（2011）13-0046-04

文献标识码：A

中图分类号：TP301

1引言

文献[6]克服了M-1r(s)

并行化处理的困难，提出了一种块二级科学计算与工程应用中的许多领域如计算流体力学、数

多分裂PE迭代算法，算法简单明了，收敛速度快；文献[7]充分值天气预报、石油地震数据处理等均离不开微分方程数值求利用系数矩阵的结构特点，提出了并行交替方向算法，并给出解，尤其是在计算流体力学中采用有限差分法、有限体积法或了系数矩阵为Hermite正定矩阵和M-矩阵收敛的充分条件；文有限元法离散求解偏微分方程时会遇到系数矩阵为块三对角献[8]基于二维重叠区域分解，对每个子区域上局部不完全LU线性方程组的求解问题。往往在实际计算过程中它们的求解分解所得到的上、下三角因子分别进行组合，给出一类全局并占整个问题的计算时间比重很大。因此研究块三对角线性方行不完全分解型预条件。

程组的高效算法显得尤为重要。

本文基于Galerkin原理，适当取基，给出了一种新的并行大规模稀疏线性方程组的求解一般采用迭代法，这是由算法，并与多分裂方法作了比较，结果表明算法简便可行且具于迭代法所需存储量和计算量比较少，且易于控制。1985年，有良好的并行性。

O’Leary和White提出了并行多分裂方法[1]并进行了深入研究，它被认为是一类重要的并行数值方法。目前各种迭代方2并行算法

法[2-8]不断涌现，使用最多的迭代法是Krylov子空间方法[2-3]，但设块三对角线性方程组Ax=b，可以表示为：

其收敛速度依赖于系数矩阵的特征值分布。特征值分布越集æçA1B1öççC÷æx1öæb1ö中收敛越快，相反则收敛很慢，甚至出现不收敛的情况。文ç2A2B2÷ç献[4]通过分裂系数矩阵，使算法具有很好的并行性；文献[5]将ç÷çxç÷çç2÷÷ççb2÷÷CA÷ç÷÷÷=ççç÷÷÷（1）

ç串行计算方法通过调用BLAS子程序达到处理机间负载平衡；

çn-1n-1Bn-1÷÷ççxn-1÷÷ççbè

CnA÷çnøèx÷çnøèbn-1÷÷÷

nø基金项目：陕西省教育厅科研项目（theResearchFoundationofDepartmentofEducationofShaanxiunderGrantNo.09JK809）；咸阳师范学院重点课程项目（No.200812014）。

作者简介：段治健（1980—），男，博士研究生，研究方向：理论与计算流体力学、信息处理中的快速与并行算法。E-mail：zhijian_duan@收稿日期：2009-09-08；修回日期：2010-01-18

求解块三对角方程组的一种并行策略

段治健，杨永，吕全义，等：求解块三对角方程组的一种并行策略

2011，47（13）47

其中Ai是ni阶方阵，Bi和Ci分别是ni´ni+1和ni´ni-1矩为如下形式：

阵，xi和bi均为ni维列向量。假设处理机台数为p，记第i台

æ处理机为PçD1öæyi(i=12p)，为了问题讨论的方便，设ç1öæçrˉT

1rˉ1öT÷çn=pt(t³2tÎZ)。基于Galerkin原理[9]，适当选取子空间的基

çD2÷÷ç÷ççç÷çy2÷÷çrˉ÷÷=ç2rˉ2÷

底，进而推出并行算法。

èD÷ççp÷øçèy÷÷pøçèr

ˉT÷prˉp÷ø2.1Galerkin原理

其中

任意给定x0ÎRn，令x=x0+z，则方程组Ax=b等价于Di=

æAz=rrT（2）

çç(i-1)t+1A(i-1)t+1r(i-1)t+1rT(i-1)t+1B(i-1)t+1r(i-1)t+2

ö0

ç其中rnçrTç(i-1)t+2C(i-1)t+2r(i-1)t+1rT(i-1)t+2A(i-1)t+2r(i-1)t+2rT

(i-1)t+2B(i-1)t+2r÷÷(i-1)t+3÷0=b-Ax0。设R中的两个m维子空间Km和Lm的基ç÷

çm

ç分别为{vm

èrT÷÷T

÷÷it-1Cit-1rit-2rit-1Ait-1rit-1ø

i}

i=1

和{wi}

i=1

，寻找式（1）近似解zmÎKm，使得

(ræ0-Azm)^Lm，

即çy(i-1)t+1öær(i-1)t+1ö(ryçy÷ç÷i=ç0-Azm)^wi(i=12m)

（3）ç(i-1)t+2÷çr(i-1)t+2÷ç÷构造Vç÷÷rˉi=ççç÷(i=12p)m=(v1v2vm)和Wm=(w1w2wm)，则由

èy÷çr÷÷÷it-1øèit-1ø

zmÎKm可得：

则有

zm=η1v1+η2v2++ηmvìm=Vmym

x(k+1)=T

íi-1)t+jx(k)ï((i-1)t+j+y(i-1)t+jr(i-1)t+j其中yWT

ï

m=(η1η2ηm)ÎRm,则式（3）可写为m(r0-îx(k+1)(k)it=xit

AV=0,即(WT)

TT

mym)mAVmym=Wmr0，如果WmAVm可逆，则有

(i=12p;j=12t-1)yr

(k+1)

m=(

WT

AV1

TT)

m

m)

-Wmr0

，因此zm=Vm(

Wm

AVm)

-1

WTr。显然，

=b-Ax

(k+1m0

如何选择K当k为偶数时，VTmAVm是一个对角矩阵，此时式（4）形式

m和Lm使得WTm

AVm可逆是问题的关键。2.2并行算法推导

如下：

为了讨论问题方便，不妨取n为偶数，记

æçrTtAtrt

ö÷æytöærTtrtö

֍T

çrTçy2t÷÷çræ(k)(k)

(ç2tA2tr2t

çrT÷2i=çråii=12ptnriö

)

÷çç÷÷=çtr2t÷èj-1+1j=1

ån=1j÷jøçè

rT

ççptAptrpt÷ø

èy÷pt÷øç÷èrTptrpt÷ø令n0=1，并记

则

æìx(k+1)Q=(q)

çr100ö

ç0rí=x(k)ïitit+yitrit

1q2qpt=çç20÷÷

kç÷=ï(k+1)(k)

(i=12p;j=12t-1)ç÷îx(i-1)t+j=x

(i-1)t+jè

00r÷pt÷

ør

(k+1)

=b-Ax

(k+1)

æçr(k)100ö以上便是求解方程的并行迭代格式。

çç÷()÷k2.3

并行实现

çrn100÷2.3.1

存储方法

ç()÷

kç0rn(0)

1+1

0÷将初始向量x+jCç(i-t)+j及A(i-1)t+jB(i-1)t(i-1)t+jD(i-1)t+j，精

ç÷

度要求ε，分别存入到第i台处理机Pi(i=12p)中。ç0r(÷k)nn20÷1+ç2.3.2

循环计算过程

ç÷(k)÷

ç00rn

÷（1）第i台处理机P(0)

(0)

i计算r=b-Ax

，k=0。

nçåj=1j-1+1÷ç÷

（2）P(k)

<ε，

则x(k)

i判断，若ri为近似解，则停机，否则。

ç(k)

÷

（3）Pç00ri并行计算

èån÷n÷j=1j

ø

æ显然QçD1öæyˉT

k是列正交的。因此采用以下基底：

çç÷qqçD2÷ççy1öæçr1rˉ1ö2÷÷T÷t-1qçç÷÷÷çç÷=çrˉ2rˉ2÷Vìï(1

t+1q(p-1)t-1q(p-1)t+1qpt-1)èDç÷÷çp÷çøèy÷m=Wm=í

p÷øçïèr

ˉTprˉp÷øî(qtq2t

qpt)ìíx(k+1)t+j=x(k)

则方程组的解为x

(k+1)

=x

(k)

+V(k)

ï(i-1)(i-1)t+j+y(i-1)t+jr(i-1)t+jmy

，其中y

(k)

满足方程

ï

(V

T

)

k)

x(k+1)(kî)it=xit

mAVm)

y

(k=VTmr

(（4）

显然，关键的问题是求解式（4），取基方式使得算法具有

(i=12p;j=12t-1)很好的并行性。

r

(k+1)

=b-Ax

(k+1)

当k为奇数时，VTmAVm是一个块对角矩阵，此时式（4）成

并进行一次通讯，将x(k+1)

(i-1)t+1传给Pi-1(i=23p)。

求解块三对角方程组的一种并行策略

482011，47（13）ComputerEngineeringandApplications计算机工程与应用

（4）P-1

i并行计算

FGGT

(GGT)-1

(FTAF)-1

((GT)T

GT)(GT)T

GTFTA=

ærTtAtrtöæyTççrT

÷töærtrtöF(FTAF)-1

FTA

ç2tA2tr2t÷ççy2t÷÷ççrTr÷2t2t÷并且有

ç÷ççç÷÷=ç÷T

èrTAç÷÷ptptrpt÷øèyptøçèrTptrpt÷ø

e

(k+1)

T

Ae

(k+1)

=e

(k)

T

Ae

(k)

-e

(k)

AV(TAVT

mVmm)

+

V(k)

mAe-

ìx(k+1)=x(k)ïíitit+yitrit

e

(k)

T

AVT+

k

ï(k+1)(k)

(i=12p;j=12t-1)m(VmAVm)

VT

()mAe

+îx(i-1)t+j=x

(i-1)t+jk

T

e()AVT

AVr

(k+1)

=b-Ax

(k+1)

)+

VTAV(VTAV)+

VTAe(k

)m

(V

m

mm

m

m

m

m

=k

T

kk

T

并进行一次通讯，将x(k+1)

it

传给Pi+1(i=12p-1)。

e()Ae()-e()AV(VTAV)+

VTAe(k

)m

m

m

m

=

（5）k+1Þk，转向（2）。

e

(k)

T

Ae

(k)

-e

(k)

T

AF(FTAF)-1

FTAe

(k)

从上面的计算过程可以看出，每次迭代只需相邻两台处理机间进行两次通信，因而使得算法具有很好的并行性。适又由FTF=I则一定存在矩阵SÎRn´(n-r)

r，使得(FS)

合在MIMD分布式存储环境下进行并行计算。

为正交矩阵，记(FS)=Q，则FTAe2£QTAe2=Ae2，联

立引理3.2得：

3收敛性分析

由上述推导过程知2

e

(k+1)

2

£e(k)

2

-e

(k)

T

T

-1

AA

AF

2λmin((F

AF))£e(k

)

2A

-

x

(k+1)

222

2

r(k

)=x

(k)

+V(k)



my

=()

+

(k)

+

e

(k)

T

AF

k

çæ2

èI-VVTTmmAVmVmAö÷øx£e()A

-=

2

2

VT+

T

2m(VmAVm)

VmAx*

k

2k

k

k

-2

其中x为方程组Ax=b的解，记S()()*

A

r(k

)k=Vm(

VT

AV+

Te()

(r,r)e()m

m)

Vm

A，则上,A-1

r(k

)

2

式为x(k+1)

=(I-S(k)

k)x+S*kx。若记e

(k)

=x

(k)

-x*，则有误差

由引理3.3得：

2方程e

(k+1)=(I-S(k)

k)e

，为证收敛性给出以下几个引理。

A

引理3.1设矩阵Vn´m2e(k+1

)£e(k

)

2e(k)

mÎR

r

，则一定存在一个满秩分解

A

A

-

A-12A=

2

Vn´rm=FG，其中GÎRr´m

r

æ2

k+1

rFÎR

且FTF=Ir。

证明由文献[10]知，V，ÎRr´mç1-1öe(k)e(0

)

2m=FG其中GrFÎRn´rr。èA-1A£æ1-1ö

2A÷ç2øè

A-1

÷A

2A2ø设F=(fn

(k+1)2

1f2fr)fiÎR(i=12r)，则由Vm的结

因此，

e(k+1

)A£æç1-1ö

e(0

)

构得：

è

A-1

÷A

，得证。

2A2øæçf12200ö

4数值算例

FTF=ç÷

ç0

f22÷20ç÷

例1已知块三对角线性方程组系数矩阵形如式（1），其中

ç÷

æ4-1è

00f2r÷2øA=ççç-14-1ö÷ç÷

i记F

ˉ=diagæç÷Bi=Ci=-IC1èf-1-1

-112f22,frö2øÎRr´r,F=FFˉ,G=Fˉ-1G得：çè

-=B1÷÷，

m=0-41-41÷

øt´tVm=FG=FFˉ-1FˉG=FG，且FTF=F

ˉTFTFFˉ=IrbTx(0)T

引理3.2[10]设AÎRn´n是对称正定矩阵，则

i=(111)t´1m=9600mi=20i=(000)t´1，在

maxHPrx2600集群上进行并行测试，计算结果见表1、表2。

x2

¹0R(x)=λmax(A)

引理3.3

[11]

设矩阵FÎR

n´r

且FT

F=IR

n´n

表1

本文算法计算结果表

rAÎ是对称

P1248正定矩阵，则有

T963.9996

492.8244251.7298132.3971λ-1

)³λ

S1.97463.86587.3501min((FTAF)

min

(A-1)

定理3.1设AÎR

n´n

是对称正定矩阵，则

E0.98730.96640.9188L

41144124

4126

4126

e

(k+1)



(k+1)2

A

£æç1-1öA-1e(0)



Δ9.8879E－119.9845E－119.8893E－119.8908E－11

è

2A÷2øA

表2例1采用多分裂方法计算结果表

证明由e(k+1)

=e

(k)

-VT+

T

(k)

m(

Vm

AVm)

VmAe

，及引理3.1得：

P1248V(T)

+

TT+

T

T134.2484

69.449739.688225.3379mVmAVmVmA=Vm(GFTAFG)VmA=S1.93303.38265.2983VG)+

(FTAF)+

(GT)+

VT

m(mA=

E0.9665

0.84560.6623-1

L

10531067

1067

1067

VTmGT(GG)-1(FTAF)

-1

((GT)T

GT

)

(GT)T

VTmA=

Δ1.0002E－101.4842E－101.4842E－101.4842E－10

求解块三对角方程组的一种并行策略

段治健，杨

永，吕全义，等：求解块三对角方程组的一种并行策略

2011，47（13）49

例2已知椭圆型偏微分方程（4）本文算法较文献[7]算法具有较高的加速比和并行效率。

C2

2

ux¶u+Cy¶+(C1sin2πx+C2)¶u¶x2¶y

2+5结束语

(D1sin2πy+D2)¶

u+Eu=0

提出的并行算法针对给定模型较多分裂算法具有高效性

0£x£10£y£2和边界条件

和较好的并行性，并给出了当系数矩阵为对称正定矩阵时的u|收敛性证明。算法实现简单，存储量较小，每次迭代只需相邻x=0=u|x=1=10+cosπy

u|处理机之间进行两次通信，具有良好的并行性，对在MIMD分y=1=u|y=1=10+cosπx

其中CE是常数。取两种模型

布式存储环境下求解大型块三对角线性方程组是一种较好的xCyC1C2D1D2（1）C并行算法。

x=Cy=E=1C1=C2=D1=D2=0致谢本文算法的测试得到西北工业大学高性能计算研究与（2）Cx=Cy=C1=D1=E=1C2=D2=0

发展中心的大力支持，在此表示衷心的感谢！

采用差分格式得到的线性代数方程组，

h=1/48。在HPrx2600集群上进行并行测试，计算结果见表3、表4。

参考文献：

表3

本文算法模型1计算结果表

[1]O’LearyDP，WhiteRE.Multi-splittingsofmatricesandparal-P1248lelsolutionoflinearsystems[J].SIAMJAlgDiscMeth，1985，6T365.6764

122.693567.450740.3048（4）：630-640.

S2.50444.55557.6237[2]SaadY.Iterativemethodsforsparselinearsystems[M].Boston：

E1.25221.13890.9530PWSPubCo，1996.

L

3687539720

39662

37742

Δ[3]李晓梅，吴建平.Krylov子空间迭代法及其并行计算[J].计算机科

9.9992E－119.9972E－119.9991E－119.9927E－11

学，2005，32（1）：19-20.

表4

本文算法模型2计算结果表

[4]CuiXining，LüQuanyi.Aparallelalgorithmforblock-tridiago-P1248nallinearsystems[J].AppliedMathematicsandComputation，T441.9931

143.120368.593550.84562006，173：1107-1114.

S2.64135.51117.4348[5]骆志刚，李晓梅.块三对角线性方程组的一种分布式并行算法[J].

E1.32071.37780.9293计算机学报，2000，23（10）：1028-1034.

L

4678947112

47253

43684

Δ[6]刘兴平，徐涛.块二级多分裂迭代方法[J].计算机工程与应用，

9.9917E－119.9921E－119.9993E－119.9986E－11

1999，35（1）：44-46.

注对于本文例2，多分裂方法失效。

[7]段治健，吕全义，马欣荣.带状线性方程组的并行交替方向算法[J].

其中，P为处理机台数，T为运行时间（单位：s），S为加计算机工程与应用，2009，45（20）：54-56.

速比，

E为效率，L为迭代次数，Δ为Ax-b[8]吴建平，宋君强.块三对角线性方程组的一类二维区域分解并行不

¥。完全分解预条件[J].计算物理，2009，26（2）：191-199.算例分析：

[9]张凯院，徐仲.数值代数[M].西安：西北工业大学出版社，2000.（1）本文算法较多分裂算法具有较高的加速比和并行效[10]程云鹏，张凯院，徐仲.矩阵论[M].西安：西北工业大学出版社，

率，且加速比成线性增长。

1999.

（2）算法取基方式比较简单，具有良好的并行性。[11]谷同祥，刘兴平，迟学斌.对称正定问题多搜索方向共轭梯度法的

（3）算法设计简单明了，易于实现。

收敛性理论[J].计算数学，2004，26（1）：117-128.

（上接29页）

systems[J].InformationSciences，1998，112：39-49.

[3]GuanYY，WangHK.Set-valuedinformationsystems[J].Informa-[10]KryszkiewiczM.Rulesincompleteinformationsystems[J].Infor-tionSciences，2006，176：2507-2525.

mationSciences，1999，113：271-292.

[4]ZadehLA.Fuzzysets[J].InformationandControl，1965，8（3）：

[11]DemriS，OrlowskaE.Incompleteinformation：Structure，infer-338-353.

ence，complexity[M].Heidelberg：Springer-Verlag，2002.

[5]TurkenB.Intervalvaluedfuzzysetsbasedonnormalforms[J].

[12]WuWZ.Attributereductionbasedonevidencetheoryinin-FuzzySetsandSystem，1986，20：191-210.

completedecisionsystems[J].InformationSciences，2008，178：[6]GorzafczaryB.Interval-valuedfuzzycontrollerbasedonverbal

1355-1371.

modalofobject[J].FuzzySetsandSystem，1988，28：45-53.[13]孟广武.区间值Fuzzy集的基本理论[J].应用数学，1993（2）：

[7]GongZengtai，SunBingzhen，ChenDegang.Roughsettheory

212-217.

fortheinterval-valuedfuzzyinformationsystems[J].Information[14]曾文艺，李洪兴，施煜.区间值模糊集合的分解定理[J].北京师范

Sciences，2008，178：1968-1985.

大学学报：自然科学版，2003，39（2）：171-177.

[8]SunBingzhen，GongZengtai，ChenDegang.Fuzzyroughsetthe-[15]曾文艺，李洪兴，施煜.区间值模糊集合的表现定理[J].北京师范

oryfortheinterval-valuedfuzzyinformationsystems[J].Informa-大学学报：自然科学版，2003，39（4）：444-447.

tionSciences，2008，178：2794-2815.

[16]张文修，吴伟志，梁吉业，等.粗糙集理论与方法[M].北京：科学出

[9]KryszkiewiczM.Roughsetapproachtoincompleteinformation

版社，2003.

本文来源：https://www.bwwdw.com/article/rxf4.html