1、 Krylov子空間上并行PCGMRES算法的研究及在邊界元問題中的應(yīng)用本文相關(guān)研究得到國家自然科學(xué)基金(70773035),河北省自然基金(F200600377)(一般不超過30個(gè)字)完成人：*隊(duì)張三，李四，王五摘要：通過研究基于主從模式的并行計(jì)算模型牙Krylov子空間Gmres算法的基本理論，提出了一種Krylov子空間上帶預(yù)校型的并行PCGMRES新算法，給出了求解線性方程組的算例和求解大型彈性邊界元問題的算例，與并行GMRES(m)算法的運(yùn)行結(jié)果進(jìn)行比較之后表明，所設(shè)計(jì)的并行算法在保證計(jì)算精度的前提下，可以減少迭代次數(shù)，縮短計(jì)算時(shí)間，有很好的加速比和計(jì)算效率(。8行左右)關(guān)鍵詞：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引言對于大部分工程實(shí)際問題來說，線性方程組的求解是其中計(jì)算的核心部分，如有限元和邊界元問題，最后都?xì)w結(jié)為求解大型線性方程組，但是目前大部分算法都局限于研究稀疏線性方程組的求解問題，而對于大型稠密線性方程組的求解卻不是很完善，特別是在并行技術(shù)產(chǎn)生之后，研究計(jì)算大型稠密線性方程組的并行算法已經(jīng)成為了一個(gè)科學(xué)計(jì)算中的重要問題。邊界元法

5、作為一種強(qiáng)有力的數(shù)值計(jì)算方法，由于其高精度和降維等優(yōu)點(diǎn)，被大量的應(yīng)用于許多學(xué)科的分析計(jì)算中，尤其是在計(jì)算數(shù)學(xué)和計(jì)算力學(xué)領(lǐng)域被公認(rèn)為是有限元法的重要補(bǔ)充。但是邊界元法最后所得方程組的系數(shù)矩陣一般具有非對稱性和稠密的特點(diǎn)，并且對于邊界元求大規(guī)模問題時(shí)，若劃分的單元較多時(shí)，其運(yùn)算量和存儲量將劇增，而已限制了邊界元在大規(guī)模問題中的應(yīng)用。隨著高速網(wǎng)絡(luò)技術(shù)的快速發(fā)展，并行計(jì)算已經(jīng)成為大型計(jì)算的主要技術(shù)之一，機(jī)群系統(tǒng)已經(jīng)成為并行計(jì)算的主要平臺，但是一般并行算法只適用于中等粒度以上的并行，這樣就有必要設(shè)計(jì)適合于網(wǎng)絡(luò)并行的粗粒度并行算法。GMRES方法(GeneralizedMinimalResidualalg

6、orithm)是一種求解大型非對稱線性方程組的Krylov子空間投影法，1986年由Y.Saad和M.H.chltz提出.該法基于Krylov向量的完全正交化，由于所需計(jì)算量和內(nèi)存量較少等方面的優(yōu)勢，目前廣泛應(yīng)用于機(jī)械、計(jì)算力學(xué)、計(jì)算數(shù)學(xué)等工程領(lǐng)域。為了提高算法的計(jì)算效率，采用預(yù)條件子技術(shù)是一種非常有效的方法。劉曉明等將GMRES方法應(yīng)用于油藏?cái)?shù)值模擬計(jì)算問題，利用矩陣分塊技術(shù)和擬消去法(PE)對系數(shù)矩陣進(jìn)行了預(yù)處理。于春肖、楊愛民等1,2提出改進(jìn)GMRES(m)算法收斂性的一種預(yù)條件方法，并論證方法的正確性。陳一鳴、楊愛民等34利用分治策略對GMRES方法進(jìn)行了并行處理，用于快速解決矩陣和向

7、量之間、矩陣和矩陣之間的運(yùn)算。除上述方法外，本文將給出一種加快收斂速度、減小存儲空間的并行帶預(yù)校GMRES(PCGMRES)算法，此算法是基于B.K.Schmidt等給出的度量通信開銷的方法5,6，是引入并行技術(shù)的Krylov子空間投影類方法7仙期間使用預(yù)測一校正策略，可以減少計(jì)算時(shí)間和存儲量，提高計(jì)算效率。2PCGMRES算法2.1Krylov子空間上的Galerkin原理設(shè)所求方程組為Ax=b，其中A為一非奇異的大型矩陣，beRn為一給定向量。后面用到的范數(shù)都是2-范數(shù)。記K和L是Rn中的兩個(gè)m維子空間,分別由5h和whmmii=1ii=1所張成。取XeRn為任意向量。令x=x+z，則原方

8、程組等價(jià)為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?？傻肳TAV人=WTr，mmmmmmmmm0假設(shè)WtAV為非奇異陣，則可得到近似解z=VWtAV-1WTr。mmmmmmm02.2GMRES算法若選取L=K，則稱之為Arnoldi算法，若選取L=AK，則稱之為GMRES算法。GMRES算法經(jīng)過許多相關(guān)人員的研究，已經(jīng)取得了切艮大的改進(jìn)和完善，同時(shí)將它和各種預(yù)測一校正系統(tǒng)結(jié)合起來，已經(jīng)成為了當(dāng)前求解大型非

9、對稱線性方程組的主要手段。若取K=span,Ar，Am-1r可得K的一組標(biāo)準(zhǔn)正交基：m000m-Vm+1r-AVy|=llr-AVHye1-Hmy10Azll就等Az|=|peHy|，于是在Rn中極小化r1m價(jià)于在K”中極小化|Pe1-Hj。即可以把它最終歸結(jié)為求解最小二乘方程minpe-Hy|。GMRES算法的計(jì)算步驟為：已知VTV=I,則有|rm+1m+10min0m0m+1m任取X，計(jì)算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)構(gòu)造近似解x=x+Vy，其中y滿足minJ(y)(J(y)=0e-HyII)0k0kkk111kk112.3預(yù)測一校正系統(tǒng)下的Gyres書法(PCGMRES)通過理論分析可知，如果，廠-1線性無關(guān)，當(dāng)m=n時(shí)，GMRES算法可以給出方0i=0程的解析解。但是，對于實(shí)際數(shù)值計(jì)算來說，當(dāng)m很大時(shí)，計(jì)算中需要保存所有的，ii=1這將引起存儲空間過多和計(jì)算時(shí)間過長，并且當(dāng)kTs時(shí)，矩陣V的各列之間的正交性k也會變得很差，此時(shí)，解會引起振蕩而并不收斂。而對該算法引入預(yù)測校正系統(tǒng)后，可以通過改變迭代的初始值的重新開始技術(shù)來克服這個(gè)困難，這樣就是預(yù)測校正系統(tǒng)下的Gmres算法(PCG

11、MRES)。PCGMRES算法的具體實(shí)現(xiàn)步驟為：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矩陣，當(dāng)j=1時(shí)第一列省略，并且有AV=VH。解最小二乘問題|r|=min0e-Hy獲得y；1mmmywRm構(gòu)造近似解x=x+Vmym0計(jì)算殘余向量的模r預(yù)測-校正：若Irj3)mm+

12、1myeRmm-；mm;II=Ib一AxII;，則令x=mx即為所求；若|r|8，則不滿足誤差(5)(6)要求，此時(shí)令兀=X，重修預(yù)測初始值再返回(1)進(jìn)行校正計(jì)算7可多次進(jìn)行)。其中mAm為預(yù)測一校正的迭代步數(shù)。3并行PCGMRES算法Arnoldi算法和GMRES算法在構(gòu)造Hessenberg矩陣元素和Krylov向量時(shí)都要用到長遞推算式，引起存儲量大和計(jì)算時(shí)間長，預(yù)校系統(tǒng)下的并行GMRES算法是為了克服這些缺點(diǎn)。PCGMRES算法主要計(jì)算包括：計(jì)算向量內(nèi)積、矩陣和向量相乘、矩陣和矩陣相乘用QR分解求解最小二乘問題、預(yù)測校正重新啟動(dòng)等。對于大型線性問題來說，算法中的這幾個(gè)部分的并行都是必要

13、的和可行的。期間基本是根據(jù)分治策略，依據(jù)主從模式，設(shè)計(jì)粗粒度的并行算法。這對節(jié)點(diǎn)數(shù)量不多的機(jī)群系統(tǒng)來說是一個(gè)很好的辦法。對于大型計(jì)算問題來說，算法中的正交化、條件矩陣的建立、條件矩陣的相關(guān)計(jì)算、QR分解求最小二乘問題、預(yù)測校正重新啟動(dòng)是算法中計(jì)算的主要部分，因此在設(shè)計(jì)的并行算法是對這幾個(gè)部分的并行算法的有機(jī)結(jié)合。若令A(yù)二(AT,AT，,AT)t,b二(fT,fT，,fT)T，上式為方程組的分塊形式，各個(gè)塊可以分布到不同的處理機(jī)上，利用預(yù)校系統(tǒng)下并行GMRES算法來求解線性方程組，完成式(上1)的并行迭代求解。為了能快速得到式(1-1)的收斂解，在每一步迭代中都重新形成矩陣H。k并行PCGMRE

14、S算法的詳細(xì)迭代步驟為(假設(shè)有P臺處理器)：任取X，給定參數(shù)值g,a,卩，m。0在各個(gè)處理機(jī)P(i二1,2,P)上計(jì)算r(i)二b-AX，利用主從模式，經(jīng)過通訊i0i0得到仃二昱r0(訂和|匕,然后再發(fā)送r和|匕|給P，并計(jì)算得到氣二r/r|。i=1迭代：DOk=1,n在P上計(jì)算Av，經(jīng)過通訊得到Av=才Av。iikkiki=1在P上進(jìn)行下列計(jì)算：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(預(yù)測)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上獨(dú)立的完成有關(guān)的其他計(jì)算。i(保證V的正交性)，max11r|r|0iki保證過其中，除以外，大寫字母表示矩陣，小寫字母表示向量。其中nfg，|rj|(保證精度要求)，maX|v|a1im1程的穩(wěn)定性）。該并行算法總體可以概括為：（1）在形成V和H矩陣的正交化過程中，將調(diào)用并

16、行內(nèi)積算法、并行矩陣向量乘積（2）在求解最小二乘問題的過程中，將調(diào)用并行QR分解、并行矩陣向量乘積和并行矩陣乘積等算法。（3）在整個(gè)計(jì)算過程中始終遵循分治策略和主從模式，在并行計(jì)算的迭代過程中要進(jìn)行先預(yù)測，后校正的過程。以上就是并行PCGMRES算法的基本計(jì)算格式和設(shè)計(jì)思路，它利用并行中的分塊技術(shù)和內(nèi)外存的交互技術(shù)，增加了解題的規(guī)模，加快了計(jì)算速度，降低了分析和通訊時(shí)間，因此該并行GMRES算法更適合于大規(guī)模工程問題的計(jì)算，為線性方程組的大規(guī)模應(yīng)用提供了一種良好的方法。4并行PCGMRES算法在求解大型線性方程組中的應(yīng)用在1000Mbps以太網(wǎng)機(jī)群中對算法進(jìn)行模擬，程序米用王從模式結(jié)構(gòu)，用MP

17、I+Fortran語言將其實(shí)現(xiàn)。我們在8節(jié)點(diǎn)的機(jī)群系統(tǒng)中分別用2個(gè)節(jié)點(diǎn)、4個(gè)節(jié)點(diǎn)和8個(gè)節(jié)點(diǎn)對并行PCGMRES算法進(jìn)行模擬，并與已知串行和并行GMRES(m)算法的運(yùn)行時(shí)間比較，其中的節(jié)點(diǎn)都是P42.6GHzx2。假設(shè)線性方程組的求解劃分為2、4或8個(gè)子任務(wù)。在n=800，k=3,p=4時(shí)，給出預(yù)測-校正的迭代次數(shù)m的一些值，利用本文的并行算法加以計(jì)算，m對預(yù)測一校正次數(shù)l的影響如圖1所示。JIJ28060050100150200250300M(迭代次數(shù))40200080604020圖1迭代次數(shù)對預(yù)測-校正次數(shù)的影響試設(shè)定預(yù)測-校正的迭代次數(shù)m=200，部分計(jì)算結(jié)果如表1所示。表1GMRES(

18、m)算法與PCGMRES算法的計(jì)算結(jié)果比較n算法串行時(shí)間(s)K=2P=2K=4P=4K=8P=8并行時(shí)間(s)加速比效率并行時(shí)間(s)加速比效率并行時(shí)間(s)加速比效率600GMRES(m)90.6658.491.550.7843.562.130.5344.222.050.26預(yù)測-校正GMRES90.6650.031.810.9140.682.230.5635.792.530.32800GMRES(m)130.5881.611.600.8053.082.460.6162.182.100.26預(yù)測-校正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預(yù)測-校正GMRES190.53102.351.860.9360.333.160.7972.342.630.33表中的n表示矩陣的階數(shù)，P為處理機(jī)臺數(shù)，K為劃分的任務(wù)數(shù)。從表1中可以看出提出的并行算法在1000M機(jī)群環(huán)境下，在P=2時(shí)獲得了一定的加速比和效率；當(dāng)n固定，p=K,K增大時(shí)，算法的加速比也應(yīng)該隨之增大，但實(shí)驗(yàn)中K=4時(shí)的加速比與K=8時(shí)的加速比相比較，反而有所降低，這是因?yàn)镵增大時(shí)，傳輸?shù)臄?shù)據(jù)量強(qiáng)增大，而計(jì)算的通信開銷增大引起的。這與理論的分析結(jié)論是吻合的。

20、5并行PCGmres算法在邊界元中的應(yīng)用邊界元法計(jì)算彈性問題時(shí)，首先需要確定物體的邊界曲面，然后將曲面離散成若干四邊形線性單元或其它類型單元，并同時(shí)給出單元的節(jié)點(diǎn)組成及節(jié)點(diǎn)坐標(biāo)。單元的節(jié)點(diǎn)編號為逆時(shí)針方向，且單元的外法矢指向邊界外側(cè)。根據(jù)邊界上已知的力和位移條件，確定各節(jié)點(diǎn)和單元的邊界條件。然后進(jìn)入迭代求解，使用的方法是并行GMRES(m)算法。最終構(gòu)造出新的近似解。如果近似解滿足精度要求，迭代終止，否則以得到的近似解為初解繼續(xù)迭代，直到滿足設(shè)定精度為止。設(shè)有單個(gè)帶空彈性體A(20mX20mX2.5m,3mX2.5m),其一端固定的情況下受拉力的計(jì)算模型及和邊界離散結(jié)構(gòu)模型分別如圖2和圖3。計(jì)

21、算模型的離散數(shù)據(jù)見表1。A物體的彈性模量E=210GPa,泊松比v=0.3,所受到的均布載荷P=104MPa。帶空彈性體A的中間一列節(jié)點(diǎn)上的節(jié)點(diǎn)位移(Z方向)和面力的在基于上述并行PCGmres算法的邊界元法以及與有限元法的計(jì)算結(jié)果之比較見圖4、圖5。計(jì)算中取PCGmres中m=200，計(jì)算的誤差RK=2.843213900325107E-006時(shí)終止，并行PCGmres算法下的邊界元求解總的計(jì)算時(shí)間為：4時(shí)15分27秒13毫秒，而并行GMRES(m)算法下的邊界元求解總的計(jì)算時(shí)間為：5時(shí)30分50秒10毫秒。根據(jù)算例的結(jié)果與有限元法的結(jié)果進(jìn)行比較可以看出，彈性體的位移和面力分布合理，證明了該

22、模型的正確性和求解的穩(wěn)定性，即該數(shù)值算例的結(jié)果表明，使用并行PCGMRES算法求解邊界元問題，模型具有較高的計(jì)算效率、精度和數(shù)值穩(wěn)定性，與串行方法的結(jié)果比較，縮短了計(jì)算時(shí)間，提高了計(jì)算效率。圖3離散結(jié)構(gòu)模型(局部)圖2計(jì)算模型SSERTS-Z5024681012141618結(jié)點(diǎn)數(shù)（從左到右X103）圖5A物體的Z向位移壓力（有限元法和邊界元法比較）表2離散數(shù)據(jù)總結(jié)點(diǎn)數(shù)總單元數(shù)外邊界結(jié)點(diǎn)數(shù)孔內(nèi)結(jié)點(diǎn)數(shù)9.36X1049.36X1041.8X1043.2X4X1046結(jié)論與未來工作本文提出的求解線性方程組的并行PCGMRES算法，具有通信開銷小，并行度較高的特點(diǎn)。通過理論分析和計(jì)算結(jié)果與已有算法的結(jié)果進(jìn)行比較，然后更新了傳統(tǒng)邊界元法的結(jié)構(gòu)，建立了基于PCGmres算法的并行邊界元法的計(jì)算格式。然后以單個(gè)帶孔彈性體為模型的受力問題，數(shù)值算例表明在邊界元方法中使用并行PCGmres算法，具有較高的計(jì)算效率、精度和穩(wěn)定性。本文的研究結(jié)果在計(jì)算數(shù)學(xué)、計(jì)算力學(xué)等領(lǐng)域具有十分廣闊的應(yīng)用前景。參考文獻(xiàn)AiminYang,ChunfengLiu.TheResearchandApplicationoftheParallelAlgorithmforQRDecompositionofMatri