基于Jacobi-Davidson算法的大規模模態分析并行計算研究

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

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

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

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

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

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

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

1 模態分析中的J-D算法設計

模態分析的數學本質就是求解如下廣義矩陣特征值問題：

Kx=λMx

(1)

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

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

首先，對于廣義Hermitian特征值問題，修正方程表達式變為[1]：

(2)

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

KQk=ZkDk

(3)

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

(4)

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

(5)

施加預條件后，修正方程(4)變為：

(6)

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

圖1 模態分析問題的J-D算法

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

(7)

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

2 并行計算實現

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

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

對于模態分析并行計算的實現采用流程如下：

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

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

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

3 模態分析算例

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

3.1 J-D算法參數設置與結果對比

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

圖2 飛行器有限元模型

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

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

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

表1 內層迭代與外層迭代及計算時間的關系

表2 重啟動向量個數與外層迭代及計算時間的關系

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

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

表3 J-D算法與IRLM算法計算結果比較

3.2 大規模并行測試

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

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

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

圖4 三種不同規模模態分析的并行加速比曲線

表4 128進程下不同自由度規模的并行計算結果統計

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

Sp=Ts/Tp

(8)

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

4 結 論

本文針對大規模模態分析問題，對J-D算法開展了一系列適應性改進、程序實現和并行測試研究。通過這些研究可以發現：

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

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

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

參 考 文 獻

[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.