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

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

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

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

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

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

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

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

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

Kx=λMx

(1)

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

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

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

(2)

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

KQk=ZkDk

(3)

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

(4)

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

(5)

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

(6)

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

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

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

(7)

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

2 并行計(jì)算實(shí)現(xiàn)

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

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

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

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

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

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

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

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

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

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

圖2 飛行器有限元模型

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

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

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

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

表2 重啟動(dòng)向量個(gè)數(shù)與外層迭代及計(jì)算時(shí)間的關(guān)系

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

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

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

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

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

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

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

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

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

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

Sp=Ts/Tp

(8)

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

4 結(jié) 論

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

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

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

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

參 考 文 獻(xiàn)

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