并行PCGMRES算法的研究与实现修改(计算机学报)

1、 Krylov子空间上并行PCGMRES算法的研究及在边界元问题中的应用本文相关研究得到国家自然科学基金(70773035),河北省自然基金(F200600377)(一般不超过30个字)完成人：*队张三，李四，王五摘要：通过研究基于主从模式的并行计算模型牙Krylov子空间Gmres算法的基本理论，提出了一种Krylov子空间上带预校型的并行PCGMRES新算法，给出了求解线性方程组的算例和求解大型弹性边界元问题的算例，与并行GMRES(m)算法的运行结果进行比较之后表明，所设计的并行算法在保证计算精度的前提下，可以减少迭代次数，缩短计算时间，有很好的加速比和计算效率(。8行左右)关键词：Kr

2、ylov子空间,PCGMRES算法，并行算法边界元法中图分类号：O246，O241.6TheResearchofParallelPCGMRESAlgorithminKrylovSubspaceanditsApplicationinBEMAbstract:ThroughtheresearchoftheparallelcomputationalmodelbasedontheprincipalandsubordinatemodeandthebasictheoryofGmresAlgorithminKrylovsubspace,thisessayraisesanewparallelPCGMRESalg

3、orithmwhichpossessespredict-correctpattern,andshowsthecomputingexamplesforlinearequationsandlarge-scaleelasticproblem.AfterthecomparisonwiththeresultfromtheparallelGMRES(m)algorithm,itshowsthatthisdesignedparallelalgorithmcanreducetheiterationfrequency,shortenthecomputingtimeandobtainbetterspeedupra

4、tioandcomputingefficiencyatthepremiseofassuringthecomputationprecision.Keywords:KrylovSubspace;PCGMRESAlgorithm;ParallelAlgorithm;BEM.1引言对于大部分工程实际问题来说，线性方程组的求解是其中计算的核心部分，如有限元和边界元问题，最后都归结为求解大型线性方程组，但是目前大部分算法都局限于研究稀疏线性方程组的求解问题，而对于大型稠密线性方程组的求解却不是很完善，特别是在并行技术产生之后，研究计算大型稠密线性方程组的并行算法已经成为了一个科学计算中的重要问题。边界元法

5、作为一种强有力的数值计算方法，由于其高精度和降维等优点，被大量的应用于许多学科的分析计算中，尤其是在计算数学和计算力学领域被公认为是有限元法的重要补充。但是边界元法最后所得方程组的系数矩阵一般具有非对称性和稠密的特点，并且对于边界元求大规模问题时，若划分的单元较多时，其运算量和存储量将剧增，而已限制了边界元在大规模问题中的应用。随着高速网络技术的快速发展，并行计算已经成为大型计算的主要技术之一，机群系统已经成为并行计算的主要平台，但是一般并行算法只适用于中等粒度以上的并行，这样就有必要设计适合于网络并行的粗粒度并行算法。GMRES方法(GeneralizedMinimalResidualalg

6、orithm)是一种求解大型非对称线性方程组的Krylov子空间投影法，1986年由Y.Saad和M.H.chltz提出.该法基于Krylov向量的完全正交化，由于所需计算量和内存量较少等方面的优势，目前广泛应用于机械、计算力学、计算数学等工程领域。为了提高算法的计算效率，采用预条件子技术是一种非常有效的方法。刘晓明等将GMRES方法应用于油藏数值模拟计算问题，利用矩阵分块技术和拟消去法(PE)对系数矩阵进行了预处理。于春肖、杨爱民等1,2提出改进GMRES(m)算法收敛性的一种预条件方法，并论证方法的正确性。陈一鸣、杨爱民等34利用分治策略对GMRES方法进行了并行处理，用于快速解决矩阵和向

7、量之间、矩阵和矩阵之间的运算。除上述方法外，本文将给出一种加快收敛速度、减小存储空间的并行带预校GMRES(PCGMRES)算法，此算法是基于B.K.Schmidt等给出的度量通信开销的方法5,6，是引入并行技术的Krylov子空间投影类方法7仙期间使用预测一校正策略，可以减少计算时间和存储量，提高计算效率。2PCGMRES算法2.1Krylov子空间上的Galerkin原理设所求方程组为Ax=b，其中A为一非奇异的大型矩阵，beRn为一给定向量。后面用到的范数都是2-范数。记K和L是Rn中的两个m维子空间,分别由5h和whmmii=1ii=1所张成。取XeRn为任意向量。令x=x+z，则原方

8、程组等价为Az=r，其中000r=b-Ax。00矩阵符号表示的Galerkin原理为：$.mm令V=C,v，v)和W=(w,w，w),其中m和m分别是K和Lm12mm12mii=1ii=1的基底。因此，可以把z表示成z=Vy，其中的y,eRm。可得WTAV人=WTr，mmmmmmmmm0假设WtAV为非奇异阵，则可得到近似解z=VWtAV-1WTr。mmmmmmm02.2GMRES算法若选取L=K，则称之为Arnoldi算法，若选取L=AK，则称之为GMRES算法。GMRES算法经过许多相关人员的研究，已经取得了切艮大的改进和完善，同时将它和各种预测一校正系统结合起来，已经成为了当前求解大型非

9、对称线性方程组的主要手段。若取K=span,Ar，Am-1r可得K的一组标准正交基：m000m-Vm+1r-AVy|=llr-AVHye1-Hmy10Azll就等Az|=|peHy|，于是在Rn中极小化r1m价于在K”中极小化|Pe1-Hj。即可以把它最终归结为求解最小二乘方程minpe-Hy|。GMRES算法的计算步骤为：已知VTV=I,则有|rm+1m+10min0m0m+1m任取X，计算r=f一Ax和v=r/|r|；000100迭代Forj=1,2,k,(直到满足条件)doh=(Av,v)(i=1,2,j)ijjiV=Av一工hvj+1jijii=1hj+1,j屮j+1llv=V/hj+

10、1j+1j+1,j(3)构造近似解x=x+Vy，其中y满足minJ(y)(J(y)=0e-HyII)0k0kkk111kk112.3预测一校正系统下的Gyres書法(PCGMRES)通过理论分析可知，如果，厂-1线性无关，当m=n时，GMRES算法可以给出方0i=0程的解析解。但是，对于实际数值计算来说，当m很大时，计算中需要保存所有的，ii=1这将引起存储空间过多和计算时间过长，并且当kTs时，矩阵V的各列之间的正交性k也会变得很差，此时，解会引起振荡而并不收敛。而对该算法引入预测校正系统后，可以通过改变迭代的初始值的重新开始技术来克服这个困难，这样就是预测校正系统下的Gmres算法(PCG

11、MRES)。PCGMRES算法的具体实现步骤为：yI取x=0,r=b-Ax,0=IIrII,v=r/0,V=lv1；000010迭代：Forj=1,2,mdoh=(Av,v)(i=1,2,j),ijjih=j+1,j11jVj+1Vj+1,v=V/hj+1j+1j+1,j(Hj-10Vj+1=(Vj,vj+1),Hj=Av一工hvjijii=1hijh.丿丿+1,j(j+1)Xj得到的H.为上Hessenberg矩阵，当j=1时第一列省略，并且有AV=VH。解最小二乘问题|r|=min0e-Hy获得y；1mmmywRm构造近似解x=x+Vmym0计算残余向量的模r预测-校正：若Irj3)mm+

12、1myeRmm-；mm;II=Ib一AxII;，则令x=mx即为所求；若|r|8，则不满足误差(5)(6)要求，此时令兀=X，重修预测初始值再返回(1)进行校正计算7可多次进行)。其中mAm为预测一校正的迭代步数。3并行PCGMRES算法Arnoldi算法和GMRES算法在构造Hessenberg矩阵元素和Krylov向量时都要用到长递推算式，引起存储量大和计算时间长，预校系统下的并行GMRES算法是为了克服这些缺点。PCGMRES算法主要计算包括：计算向量内积、矩阵和向量相乘、矩阵和矩阵相乘用QR分解求解最小二乘问题、预测校正重新启动等。对于大型线性问题来说，算法中的这几个部分的并行都是必要

13、的和可行的。期间基本是根据分治策略，依据主从模式，设计粗粒度的并行算法。这对节点数量不多的机群系统来说是一个很好的办法。对于大型计算问题来说，算法中的正交化、条件矩阵的建立、条件矩阵的相关计算、QR分解求最小二乘问题、预测校正重新启动是算法中计算的主要部分，因此在设计的并行算法是对这几个部分的并行算法的有机结合。若令A二(AT,AT，,AT)t,b二(fT,fT，,fT)T，上式为方程组的分块形式，各个块可以分布到不同的处理机上，利用预校系统下并行GMRES算法来求解线性方程组，完成式(上1)的并行迭代求解。为了能快速得到式(1-1)的收敛解，在每一步迭代中都重新形成矩阵H。k并行PCGMRE

14、S算法的详细迭代步骤为(假设有P台处理器)：任取X，给定参数值g,a,卩，m。0在各个处理机P(i二1,2,P)上计算r(i)二b-AX，利用主从模式，经过通讯i0i0得到仃二昱r0(订和|匕,然后再发送r和|匕|给P，并计算得到气二r/r|。i=1迭代：DOk=1,n在P上计算Av，经过通讯得到Av=才Av。iikkiki=1在P上进行下列计算：ih=(Av,v),i=1,2,k；ikkiv=Av-迓hvk+1ki=1ikih=v|；v=v/h；k+1,kk+1k+1k+rk+1,k取a=max勺v1,v/i=1,2，,k0ik+1if=eTg(对PCGMRES算法)kmmIF(f)THEN

15、kX=X+Vyk0kkGOTO(4)ENDIFIF(k二mandaa)THEN(预测)0X二X+Vyk0kk令X二X0kGOTO(2)ENDIF取0二f-minf(i=1,2,k)0kiiIF(00)，THEN0取1满足minf=fiilX(i)=X+Vy(i)0lX二X(校正)0iGOTO(2)ENDIFENDDO在P上独立的完成有关的其他计算。i(保证V的正交性)，max11r|r|0iki保证过其中，除以外，大写字母表示矩阵，小写字母表示向量。其中nfg，|rj|(保证精度要求)，maX|v|a1im1程的稳定性）。该并行算法总体可以概括为：（1）在形成V和H矩阵的正交化过程中，将调用并

16、行内积算法、并行矩阵向量乘积（2）在求解最小二乘问题的过程中，将调用并行QR分解、并行矩阵向量乘积和并行矩阵乘积等算法。（3）在整个计算过程中始终遵循分治策略和主从模式，在并行计算的迭代过程中要进行先预测，后校正的过程。以上就是并行PCGMRES算法的基本计算格式和设计思路，它利用并行中的分块技术和内外存的交互技术，增加了解题的规模，加快了计算速度，降低了分析和通讯时间，因此该并行GMRES算法更适合于大规模工程问题的计算，为线性方程组的大规模应用提供了一种良好的方法。4并行PCGMRES算法在求解大型线性方程组中的应用在1000Mbps以太网机群中对算法进行模拟，程序米用王从模式结构，用MP

17、I+Fortran语言将其实现。我们在8节点的机群系统中分别用2个节点、4个节点和8个节点对并行PCGMRES算法进行模拟，并与已知串行和并行GMRES(m)算法的运行时间比较，其中的节点都是P42.6GHzx2。假设线性方程组的求解划分为2、4或8个子任务。在n=800，k=3,p=4时，给出预测-校正的迭代次数m的一些值，利用本文的并行算法加以计算，m对预测一校正次数l的影响如图1所示。JIJ28060050100150200250300M(迭代次数)40200080604020图1迭代次数对预测-校正次数的影响试设定预测-校正的迭代次数m=200，部分计算结果如表1所示。表1GMRES(

18、m)算法与PCGMRES算法的计算结果比较n算法串行时间(s)K=2P=2K=4P=4K=8P=8并行时间(s)加速比效率并行时间(s)加速比效率并行时间(s)加速比效率600GMRES(m)90.6658.491.550.7843.562.130.5344.222.050.26预测-校正GMRES90.6650.031.810.9140.682.230.5635.792.530.32800GMRES(m)130.5881.611.600.8053.082.460.6162.182.100.26预测-校正GMRES130.5871.661.820.9144.722.920.7350.462.5

19、90.321200GMRES(m)190.53119.081.600.8076.212.500.6291.602.080.26预测-校正GMRES190.53102.351.860.9360.333.160.7972.342.630.33表中的n表示矩阵的阶数，P为处理机台数，K为划分的任务数。从表1中可以看出提出的并行算法在1000M机群环境下，在P=2时获得了一定的加速比和效率；当n固定，p=K,K增大时，算法的加速比也应该随之增大，但实验中K=4时的加速比与K=8时的加速比相比较，反而有所降低，这是因为K增大时，传输的数据量强增大，而计算的通信开销增大引起的。这与理论的分析结论是吻合的。

20、5并行PCGmres算法在边界元中的应用边界元法计算弹性问题时，首先需要确定物体的边界曲面，然后将曲面离散成若干四边形线性单元或其它类型单元，并同时给出单元的节点组成及节点坐标。单元的节点编号为逆时针方向，且单元的外法矢指向边界外侧。根据边界上已知的力和位移条件，确定各节点和单元的边界条件。然后进入迭代求解，使用的方法是并行GMRES(m)算法。最终构造出新的近似解。如果近似解满足精度要求，迭代终止，否则以得到的近似解为初解继续迭代，直到满足设定精度为止。设有单个带空弹性体A(20mX20mX2.5m,3mX2.5m),其一端固定的情况下受拉力的计算模型及和边界离散结构模型分别如图2和图3。计

21、算模型的离散数据见表1。A物体的弹性模量E=210GPa,泊松比v=0.3,所受到的均布载荷P=104MPa。带空弹性体A的中间一列节点上的节点位移(Z方向)和面力的在基于上述并行PCGmres算法的边界元法以及与有限元法的计算结果之比较见图4、图5。计算中取PCGmres中m=200，计算的误差RK=2.843213900325107E-006时终止，并行PCGmres算法下的边界元求解总的计算时间为：4时15分27秒13毫秒，而并行GMRES(m)算法下的边界元求解总的计算时间为：5时30分50秒10毫秒。根据算例的结果与有限元法的结果进行比较可以看出，弹性体的位移和面力分布合理，证明了该

22、模型的正确性和求解的稳定性，即该数值算例的结果表明，使用并行PCGMRES算法求解边界元问题，模型具有较高的计算效率、精度和数值稳定性，与串行方法的结果比较，缩短了计算时间，提高了计算效率。图3离散结构模型(局部)图2计算模型SSERTS-Z5024681012141618结点数（从左到右X103）图5A物体的Z向位移压力（有限元法和边界元法比较）表2离散数据总结点数总单元数外边界结点数孔内结点数9.36X1049.36X1041.8X1043.2X4X1046结论与未来工作本文提出的求解线性方程组的并行PCGMRES算法，具有通信开销小，并行度较高的特点。通过理论分析和计算结果与已有算法的结果进行比较，然后更新了传统边界元法的结构，建立了基于PCGmres算法的并行边界元法的计算格式。然后以单个带孔弹性体为模型的受力问题，数值算例表明在边界元方法中使用并行PCGmres算法，具有较高的计算效率、精度和稳定性。本文的研究结果在计算数学、计算力学等领域具有十分广阔的应用前景。参考文献AiminYang,ChunfengLiu.TheResearchandApplicationoftheParallelAlgorithmforQRDecompositionofMatri