基于Jacobi-Davidson算法的大規(guī)模模態(tài)分析并行計算研究

范宣華, 陳 璞, 吳瑞安, 肖世富

(1.中國工程物理研究院 總體工程研究所，四川 綿陽 621900； 2.北京大學 力學與工程科學系，北京 100871)

模態(tài)分析是工程結(jié)構(gòu)動力學數(shù)值模擬過程中極其重要的環(huán)節(jié)之一，工程結(jié)構(gòu)的日趨復雜以及精度可靠性要求的不斷提高對數(shù)值模擬提出了新的挑戰(zhàn)，往往要求結(jié)構(gòu)有限元規(guī)模達到數(shù)百萬甚至千萬自由度以上。傳統(tǒng)的商業(yè)分析軟件如ANSYS、Nastran等在計算規(guī)模和效率等方面受到了很大限制，自主研究基于先進算法的并行模態(tài)分析程序顯得尤為重要。

模態(tài)分析在數(shù)學上可以歸結(jié)為大型稀疏矩陣特征值問題。對于該類問題的研究大多基于子空間投影技術(shù)，即借助一系列低維數(shù)的子空間，將一個大型稀疏矩陣A向這些子空間進行投影，將大型矩陣A的部分特征值問題轉(zhuǎn)化為中小規(guī)模矩陣的特征值問題，利用中小規(guī)模矩陣的一些成熟算法求解其特征對，并形成對A矩陣部分特征對的一個近似[1-3]。根據(jù)子空間構(gòu)造方式不同主要有兩條思路，一是Krylov子空間類算法；二是Davidson類算法。

Krylov子空間類算法在大型稀疏矩陣特征值問題的研究過程中占據(jù)了重要地位，其中以Lanczos方法[4]最為著名。后來Lanczos方法發(fā)展成很多的變種形式，其中Sorensen等[5]提出的隱式重啟動Arnoldi/Lanczos方法和Wu等[6]提出的thick-restart Lanczos方法最為著名，至今影響深遠。然而對于模態(tài)分析，為求解最小特征值部分，需進行譜變換處理，由此導致每步子空間迭代過程中都要求解一個大型線性病態(tài)方程組，迭代求解難以收斂，大多采用直接法進行矩陣分解[1-3]。然而分解過程導致了內(nèi)存的大量消耗，并在一定程度上破壞了原稀疏矩陣的結(jié)構(gòu)特性，而且分解還帶來了各并行進程間的通訊量劇增，不利于并行計算的擴展。

Davidson類算法屬于非Krylov類子空間投影算法，它們和Krylov子空間投影算法都采用了瑞利-里茲過程，但是構(gòu)造正交基的方式不同。實際上，Davidson類算法可以看成是Krylov子空間類算法的預處理版本[7]，即相當于在Krylov類算法每步子空間擴展時添加了一個預處理矩陣。該類算法中最值得關(guān)注的一種算法是J-D方法[8]，該方法通過在某一近似特征向量的正交補空間上求解“修正方程”擴展子空間向量，對修正方程的求解大多采用迭代法，即J-D算法通過內(nèi)-外層迭代求解特征對，外層迭代擴展子空間，內(nèi)存迭代求解修正方程。這樣以來算法主要涉及矩陣-向量乘、向量積等運算，不但特征矩陣的稀疏結(jié)構(gòu)得以保留，各并行進程間通訊量較Krylov子空間類算法也大大降低，非常有利于并行計算可擴展性能的提升。

本文的研究就是基于J-D算法，通過添加譜變換、預處理以及收縮和重啟動技術(shù)將其改造成適合大規(guī)模模態(tài)分析的算法；并基于PANDA框架[9]和特征值求解包SLEPc[10]，對一簡化飛行器有限元模型進行了大規(guī)模模態(tài)分析并行計算研究。

1 模態(tài)分析中的J-D算法設計

模態(tài)分析的數(shù)學本質(zhì)就是求解如下廣義矩陣特征值問題：

Kx=λMx

(1)

式中：K為剛度矩陣，M為質(zhì)量矩陣。K和M可以通過對工程結(jié)構(gòu)進行有限元離散和積分得到，二者都是大型稀疏、對稱(半)正定矩陣。模態(tài)分析就是求解式(1)的多個低階特征對。

標準的J-D算法[8]是通過迭代求解一個“修正方程”來擴張子空間，通過設定一個目標值τ來控制特征值收斂范圍。對于大規(guī)模模態(tài)分析，需要圍繞“修正方程”求解、收縮和重啟動、譜變換等幾個方面對算法進行優(yōu)化或改進。

首先，對于廣義Hermitian特征值問題，修正方程表達式變?yōu)閇1]：

(2)

式中，*表示共軛轉(zhuǎn)置，θk和uk分別是里茲值和里茲向量，t是需要求解的正交補向量，rk是迭代產(chǎn)生的留數(shù)向量rk=Kuk-θkMuk。采用M-正交技術(shù)，可以將式(1)縮減成Schur形式[1-3]：

KQk=ZkDk

(3)

式中，矩陣Dk是k×k階對角陣，對角線元素為廣義特征問題(1)的k個特征值；Qk是關(guān)于M正交的n×k階矩陣，存放著k個相應的特征向量，Zk=MQk。通過式(3)可以將修正方程(2)轉(zhuǎn)換成如下形式：

(4)

對于該修正方程的求解，需要采取預處理技術(shù)，假定P是K-θkM的一個預條件子，即滿足P-1(K-θkM)≈I，則式(4)的預條件子可以取為：

(5)

施加預條件后，修正方程(4)變?yōu)椋?/p>

(6)

在采用J-D算法進行求解時，隨著外層迭代子空間維數(shù)的增長，存儲和計算量將會顯著增加，不利于大規(guī)模模態(tài)分析。為此，需要限制外部子空間的最大維數(shù)，當子空間達到最大維數(shù)mmax時，重新迭代擴張子空間，即重啟動。在重啟動時，為充分利用上次子空間擴張的信息，通常取最靠近目標值τ的mmin個里茲值對應的里茲向量作為重啟動初始迭代子空間。另一方面，隨著迭代進行，部分特征對會趨向收斂，需要將相應的里茲向量從迭代子空間中“剔除”，以便于算法繼續(xù)搜索下一個特征對。這時采用的是一種收縮策略，即將新的擴張子空間向量與已經(jīng)收斂的特征向量做正交化處理。

圖1 模態(tài)分析問題的J-D算法

綜合以上分析，給出模態(tài)分析的J-D算法如圖1所示。理論上，J-D算法當目標值設置τ高于最大特征值時，算法會收斂于最大特征值部分；當目標值τ低于最小特征值時，算法將收斂于最小特征值部分。而實際上，由于J-D算法的根基還是冪法和瑞利-里茲過程，里茲值總是趨向于最大特征值部分，雖然通過小的目標值設定和“修正方程”的反復修正會將最大特征值求解逐步轉(zhuǎn)到最小特征部分，但收斂速度會非常慢甚至不收斂。為此，我們考慮在J-D算法中添加一個常用的Shift-Invert譜變換[12]:

(7)

式中，σ為移位值，取值一般略低于最低階特征值。這樣以來，分別以矩陣M和K-σM代替原問題的K和M，以θ=1/(λ-σ)代替λ，將會得到θ的最大特征值部分，再反算得到最小的λ。若對應工程問題沒有剛體模態(tài)，可以取σ=0。這樣只需在J-D算法中將K和M調(diào)換即可，這時設置目標值τ略大于最低階特征值的倒數(shù)，將會得到很好的收斂效果。

2 并行計算實現(xiàn)

最近10多年來，大批學者和工程師們基于各種大型稀疏矩陣特征值問題的算法，開發(fā)了很多的特征值并行求解軟件包，這些軟件包基本都是開源的，可以在因特網(wǎng)上進行下載。文獻[13]對這些數(shù)值軟件包進行了較為全面的概括。在這些軟件包中， SLEPc軟件包[10]得到了廣泛的關(guān)注，SLEPc基于MPI并行通訊標準，已經(jīng)集成了很多軟件庫和算法，并且具備與其他一些軟件庫的接口功能，用戶可以在其基礎上添加特征值解法器。2011年，SLEPc軟件包的3.1版本已經(jīng)對Davidson類算法進行了代碼實現(xiàn)。

雖然目前關(guān)于大規(guī)模稀疏矩陣特征值并行計算問題的算法理論和代碼開發(fā)非常多，但應用于具體工程問題的研究卻非常少。究其原因，一個數(shù)值算法應用于實際工程，還需要很多其他學科知識和輔助技術(shù)手段，以大規(guī)模模態(tài)分析并行計算為例，需要力學、計算機、數(shù)學等多學科知識和豐富的工程數(shù)值計算能力。目前找到模態(tài)分析并行實踐的見文獻[14]，該文獻采用AMG預處理技術(shù)，分別采用Lanzos方法、Davidson方法以及LOBPCG等方法實現(xiàn)了航母189萬自由度有限元模型的模態(tài)分析并行計算，被認為是當時最具挑戰(zhàn)性的工作。

對于模態(tài)分析并行計算的實現(xiàn)采用流程如下：

(1)通過商業(yè)有限元軟件如MSC.Patran進行前處理建模，并根據(jù)不同單元類型、不同材料和邊界約束等對有限元模型進行分組處理。然后將這些單元節(jié)點信息以MSC.Patran 的中性文件(neutral)形式進行輸出。目前采用MSC.Patran可以在普通工作站上建立千萬自由度以上的有限元模型。

(2)利用并行計算框架PANDA[9]生成質(zhì)量矩陣M和剛度矩陣K。PANDA框架是由中國工程物理研究院開發(fā)的一個有限元并行框架，包含了與前處理商業(yè)軟件的接口以及一些解法器接口。通過PANDA的接口功能，MSC.Patran中的有限元網(wǎng)格信息被轉(zhuǎn)換成PANDA可以識別的網(wǎng)格信息，再結(jié)合PANDA集成的Metis[15]圖剖分軟件對有限元網(wǎng)格進行區(qū)域分解，利用不同子區(qū)域并行生成質(zhì)量和剛度矩陣。

(3)通過PANDA框架調(diào)用特征值算法進行求解。將SLEPc及相關(guān)輔助軟件包集成到PANDA框架中，并根據(jù)以上算法改進對SLEPc軟件包的具體實現(xiàn)進行相應改進和添加，由PANDA框架調(diào)用改進后的J-D算法并行求解特征對。求解過程中各進程的并行通信由PANDA框架和SLEPc軟件包通過MPI進行管理。

3 模態(tài)分析算例

通過一個具體工程算例對Jacobi-Davidson算法進行一個詳細測評。采用的工程模型為一簡化飛行器有限元模型，如圖2所示。模型底部固定，采用八節(jié)點六面體單元，在MSC.Patran中分別建立了數(shù)十萬至上千萬自由度規(guī)模的測試算例，所有算例均求解前20階模態(tài)頻率。并行計算測試在某曙光5000A機群上進行，該機群由16個刀片節(jié)點構(gòu)成，每個節(jié)點由雙路四核Intel Xeon E5550 2.66 GHz處理器共享16 GB內(nèi)存，各節(jié)點間通過千兆網(wǎng)連接。

3.1 J-D算法參數(shù)設置與結(jié)果對比

從圖1中的Jacobi-Davidson算法可以看出，決定算法收斂速度的因素很多。為更合理選取各參數(shù)，首先對飛行器71萬(n=710 526)自由度規(guī)模進行了一系列測試。

圖2 飛行器有限元模型

對于目標值τ的選取，通過式(7)的譜變換，已經(jīng)將最小特征值問題轉(zhuǎn)換成了算法更容易收斂的最大特征值問題。故τ的選取只需要略大于最小特征值的倒數(shù)即可，對于本算例，最小特征值的倒數(shù)在10-7量級，故取τ=1e-6。需要注意的是，τ的取值不能小于最小特征值的倒數(shù)，否則會出現(xiàn)基頻附近漏根的情形。對于外層迭代子空間的最大維數(shù)mmax，由于施加了重啟動策略，一般選取求解特征值個數(shù)的兩倍左右即可，本算例取mmax=41。更大的子空間維數(shù)會帶來內(nèi)存需求的快速增加。

修正方程的求解屬于內(nèi)層迭代，采用了Jacobi預處理和CG迭代解法。內(nèi)層迭代如需要足夠精確，則需要足夠多的迭代步數(shù)，會嚴重影響算法的效率。故一般只采用一個近似迭代過程，主要控制迭代的步數(shù)，算法本身的精度通過外層迭代相對誤差控制，一般取1e-4就足夠了。內(nèi)層迭代和外層迭代之間存在一個平衡關(guān)系，一般說來，內(nèi)層迭代的增加會帶來外層迭代的減少，但內(nèi)層迭代占用的更多時間也會使得整個計算時間增多。表1給出了71萬模型在16進程下內(nèi)層迭代步數(shù)所致外層迭代步數(shù)及算法總時間的關(guān)系(重啟動向量個數(shù)固定為8個)。從表中可以看出，內(nèi)層迭代取30步左右達到最佳。需要說明的是，表1和后面表格中所指的計算時間均是調(diào)用J-D算法的運行時間，不包含并行區(qū)域分解和形成總體剛度、質(zhì)量矩陣的時間(事實上，并行區(qū)域分解和形成總體剛度、質(zhì)量矩陣的時間較J-D算法運行時間要少得多，基本可以忽略)。

當子空間維數(shù)達到最大維數(shù)時進行重啟動，重啟動向量的個數(shù)也會對收斂速度造成影響。重啟動向量一般選取最接近目標值的幾個里茲值對應的里茲向量作為初始向量空間，這些初始向量包含了大量上次子空間擴張的信息，其個數(shù)的選取也需要權(quán)衡，選的過少，則不能充分利用上次子空間擴張的信息，選的過多，則除繼承了上次子空間擴張的有用信息外，會摻雜進來很多不需要的信息。表2給出了71萬模型在16進程下重啟動向量個數(shù)選取所致外層迭代步數(shù)及算法總計算時間的關(guān)系(內(nèi)層迭代固定為30步)。在重啟動向量取12個左右達到最佳。

表1 內(nèi)層迭代與外層迭代及計算時間的關(guān)系

表2 重啟動向量個數(shù)與外層迭代及計算時間的關(guān)系

通過以上測試和比較，可以基本確定J-D算法求解前20階模態(tài)頻率時的各控制參數(shù)，目標值τ取1e-6，子空間最大維數(shù)取41，內(nèi)層迭代步數(shù)取30，重啟動向量維數(shù)取12，外層迭代控制誤差1e-4。

除J-D算法本身參數(shù)優(yōu)化設置外，還就算法求解的正確性與Krylov子空間類算法中的隱式重啟動Lanczos方法(IRLM)求解結(jié)果進行了對比，IRLM采用了Cholesky分解直接求解譜變換后的線性方程組，數(shù)值精度接近計算機浮點運算精度，通常被認為是“精確解”[1-2]。表3對比了兩種方法前10階的模態(tài)頻率結(jié)果。從中可以看出，二者計算結(jié)果幾乎完全一致，最大相對誤差在0.3%以內(nèi)。

表3 J-D算法與IRLM算法計算結(jié)果比較

3.2 大規(guī)模并行測試

在上節(jié)參數(shù)選取和算法正確性得以驗證的前提下，通過該飛行器簡化模型對J-D算法進行了并行可擴展性測試，主要包括計算規(guī)模和并行加速比。

計算規(guī)模方面，在128個進程下分別測試了飛行器有限元模型數(shù)十萬到上千萬自由度規(guī)模并行計算情況。圖3給出了1 025萬自由度規(guī)模前20階模態(tài)頻率計算過程中的迭代收斂歷史，從中可以看出迭代誤差隨外層迭代的變化情況，當?shù)`差低于1e-4，便得到一個特征對，算法繼續(xù)迭代并尋找下一個特征對。在繼續(xù)尋找下一個特征對的過程中，由于要不斷更新求解“修正方程”尋找新的特征方向，計算誤差在迭代過程中有所震蕩，但經(jīng)過多次修正后會在有限的迭代步內(nèi)逐漸減小直至收斂(誤差低于1e-4)。表4給出了不同規(guī)模下內(nèi)存需求、外層迭代步數(shù)、計算時間及一階固有頻率的變化情況。從中可以看出，J-D算法的內(nèi)存占用很小，并且隨計算規(guī)模增大呈現(xiàn)線性增長趨勢，這與Krylov子空間類算法相比是一個很大的優(yōu)勢，Krylov子空間類算法由于采用了直接法求解譜變換后的線性方程組，導致內(nèi)存占用非常大(經(jīng)測試IRLM方法在112萬自由度規(guī)模所占用的內(nèi)存便高達32 GB)，而且隨規(guī)模增大，所需內(nèi)存的增長也遠超線性。隨計算規(guī)模增大，外層迭代步數(shù)有所增加；1 025萬自由度規(guī)模和71萬自由度規(guī)模對應的一階固有頻率變化在8%以上，隨著規(guī)模增大，一階頻率逐漸趨于穩(wěn)定，從這點也可以看出大規(guī)模并行計算的必要性。

圖3 1 025萬自由度規(guī)模的迭代收斂歷史

圖4 三種不同規(guī)模模態(tài)分析的并行加速比曲線

表4 128進程下不同自由度規(guī)模的并行計算結(jié)果統(tǒng)計

測試過程中發(fā)現(xiàn)各進程所占用內(nèi)存基本一致，并行計算負載平衡良好。為進一步驗證J-D算法的并行可擴展性，還就幾種自由度規(guī)模進行了加速比測試。加速比是衡量并行計算程序性能和效率的基本評價方法，其定義如下：

Sp=Ts/Tp

(8)

式中，Sp為加速比，Ts為串行程序在并行計算機單處理器上的執(zhí)行時間，Tp為相應并行程序在p個處理器執(zhí)行所需的時間，Sp與p的比值即為并行效率。理想情況下，p個處理器上并行計算時間為串行計算時間的1/p，此時有Sp=p。但受到并行通信等因素的影響，在實際計算中Sp一般小于p，Sp越接近p，則說明并行程序性能越好，并行效率越高。圖4給出了71萬、688萬和1 025萬三種自由度規(guī)模在不同進程下測得的并行計算加速比曲線，經(jīng)計算三種情形在128核時的并行效率分別為67.8%、81.9%和88.1%。從圖4可以看出，三種規(guī)模在曙光5000A機群128個CPU核內(nèi)的并行加速性能非常優(yōu)異，三個加速比曲線接近線性加速，并且計算規(guī)模越大，曲線越靠近理想加速曲線，并行性能越好。究其原因，一方面J-D算法通過內(nèi)-外層迭代求解特征值問題，避開了Krylov子空間類算法的矩陣分解運算，所涉及的運算主要是矩陣-向量積或向量之間的一些內(nèi)積運算，不改變矩陣本身的結(jié)構(gòu)特性，通信量相對較少；另一方面，隨著計算規(guī)模的增大，各子區(qū)域邊界自由度數(shù)相對子區(qū)域內(nèi)部自由度數(shù)變少，從而各進程間的通訊計算量較分配到每個并行進程內(nèi)部的計算量也相應變小，并行性能得以更好的提升。

4 結(jié) 論

本文針對大規(guī)模模態(tài)分析問題，對J-D算法開展了一系列適應性改進、程序?qū)崿F(xiàn)和并行測試研究。通過這些研究可以發(fā)現(xiàn)：

(1)本文借助PANDA框架和特征值軟件包所建立的并行計算體系完全可以滿足大規(guī)模模態(tài)分析需求，對現(xiàn)有開源程序和數(shù)值計算軟件包進行集成可以為大規(guī)模并行程序開發(fā)節(jié)省大量人力和物力。

(2)改進后的J-D算法能夠很好地適應千萬自由度以上規(guī)模的模態(tài)分析，具有內(nèi)存占用少，并行可擴展性能優(yōu)異等優(yōu)點。通過選取不同的控制參數(shù)和施加譜變換技術(shù)，可以有效提高J-D算法外層迭代的收斂速度。

(3)本文在128核的并行機群上的測試結(jié)果表明，各種規(guī)模在最大CPU核數(shù)內(nèi)均未出現(xiàn)計算拐點，即J-D算法可以繼續(xù)向更多的CPU核數(shù)擴展，可以推斷利用該算法并基于更高性能并行計算機上在10余分鐘內(nèi)解決千萬自由度以上規(guī)模的模態(tài)分析問題是完全可能的。

參 考 文 獻

[1]Bai Z, Demmel J, Dongarra J, et al. Templates for the solution of algebraic eigenvalue problems：A practical guide[M].SIAM, Philadelphia, 2000.

[2]Saad Y, Numerical methods for large eigenvalue problems[M].Halsted Press, Div. of John Wiley &Sons, Inc., New York, 1992.

[3]Watkins D S. The matrix eigenvalue problem：GR and Krylov Subspace Methods[M]. SIAM, Philadelphia, 2007.

[4]Lanczos C, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators[J].Journal of Research of the National Bureau of Standards, 1950, 4(45):255-282.

[5]Calvetti D，Reichel L，Sorensen D C. An implicitly restarted lanczos method for large symmetric eigenvalue problems[J].ETNA, 1994,2:1-21.

[6]Wu K， Simon H. Thick-restart lanczos method for symmetric eigenvalue problems[R]. Report 41412, Lawrence Berkeley National Laboratory,1998

[7]Crouzeix M, Philippe B, Sadkane M. The davidson method[J].SIAM Journal on Scientific Computing,1994,15(1):62-76.

[8]Sleijpen G L G，Van der Vorst H A, A Jacobi-Davidson iteration method for linear eigenvalue problems[J]SIAM Journal on Matrix Analysis and Applications, 1996, 17(2):401-425.

[9]Shi G M，Wu R，Wang K Y. Discussion about the design for mesh data structure within the parallel framework[R].9th World Congress on Computational Mechanics and 4th Asian Pacific Congress on Computational Mechanics, Sydney, Australia, 2010.

[10]Hernandez V，Roman J E，Vidal V，et al. SLEPc：a scalable and flexible toolkit for the solution of eigenvalue problems[J].ACM Trans. Math. Softw, 2005, 31(3):351-362.

[11]Barrett R，Berry M，Chan T F，et al. Templates for the Solution of Linear Systems：building blocks for iterative methods[M].SIAM, Philadelphia, PA, USA, 1994.

[12]Meerbergen K， Spence A， Roose D. Shift-invert and cayley transforms for detection of rightmost eigenvalues of nonsymmetric matrices[J]BIT Numerical Mathematics, 1994, 34(3)：409-423.

[13]Hernandez V，Roman J E，Tomas A, et al. A survey of software for sparse eigenvalue problems. Tech. Rep[M]SLEPc Technical Report STR-6, Universidad Politecnica de Valencia, 2012.

[14]Arbenz P， Hetmaniuk U L，Lehoucq R B，et al. A comparison of eigensolvers for large-scale 3D modal analysis using AMG-preconditioned iterative methods, Internat[J]. Numerical Methods in Engineering，2005, 64(2)：204-236.

[15]Karypis G，Kumar V. Metis-a software package for partitioning unstructured graphs, partitioning meshes, and computing fill-reducing orderings of sparse matrices. version 4.0[R]Tech. Report, University of Minnesota, Minneapolis, 1998.