田 瑾，龔 利，史小衛(wèi)，徐 樂(lè)

(1.西安電子科技大學(xué)天線與微波技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西西安710071;2.華東師范大學(xué)科技創(chuàng)新與發(fā)展戰(zhàn)略研究中心，上海200062)

0 引言

矢量有限元(vector finite element method，F(xiàn)EM)是近40年來(lái)受到眾多學(xué)者關(guān)注的一種求解邊值問(wèn)題的電磁計(jì)算方法，已被廣泛應(yīng)用于求解各類電磁問(wèn)題［1－3］。金建銘［1］對(duì)該方法進(jìn)行了較為全面、系統(tǒng)的介紹和說(shuō)明，包括閉域和開(kāi)域問(wèn)題的幾種求解方法。由于現(xiàn)實(shí)中研究的問(wèn)題多是無(wú)限區(qū)域的，因此需要加入邊界條件將開(kāi)域問(wèn)題轉(zhuǎn)化為閉域問(wèn)題來(lái)求解。通常采用的邊界條件有3種:吸收邊界條件(absorbing boundary condition，ABC)、理想匹配層(perfectly matched layer，PML)［2］和邊界積分方程又叫邊界元方法(boundary integral method，BI)［3］。前2種方法結(jié)合FEM方法離散目標(biāo)得到的方程組的系數(shù)矩陣都是稀疏、對(duì)稱的復(fù)方陣，而第3種方法結(jié)合FEM方法計(jì)算得到的方程組的系數(shù)矩陣是部分稀疏部分滿秩的矩陣，針對(duì)后面的方法采用內(nèi)觀法分寫(xiě)方程組可以得到一個(gè)稠密系數(shù)復(fù)矩陣方程組和一個(gè)稀疏復(fù)矩陣方程組。

可以看出，不論采用何種邊界條件，離散后的線性系統(tǒng)都可以歸結(jié)為求解方程

(1)式中:A是稀疏/稠密、對(duì)稱的復(fù)矩陣;x是待求列向量;b是列向量。通常采用直接法或迭代法求解(1)式。直接法得到的解在理論上是準(zhǔn)確的，其缺點(diǎn)是高計(jì)算量和內(nèi)存占用空間;迭代法可減少計(jì)算量提高計(jì)算速度，但若直接將迭代算法應(yīng)用于大型有限元方程，收斂速度不能得到有效保證。由于迭代方法解的收斂速度主要由迭代矩陣譜特性決定，比較有效的求解方程組的方法是采用預(yù)處理方法改善迭代矩陣譜特性再結(jié)合迭代法進(jìn)行計(jì)算［4－6］。文獻(xiàn)［7］的數(shù)值結(jié)果表明，分解算法具有最快的收斂速度，也就是說(shuō)，可以考慮采用不完全分解對(duì)方程組進(jìn)行預(yù)處理，這樣既能保證迭代法收斂速度，同時(shí)又不用付出完全分解所消耗的大量?jī)?nèi)存和時(shí)間。

直接采用分解法不能很好地利用緩存，而基于混合拓展喬里斯基(expanded Cholesky，ECM)［8］的多波前算法(multifrontal algorithm，MF)［9］(ECM/MF)從轉(zhuǎn)化內(nèi)存間接尋址為直接尋址的思想出發(fā)設(shè)計(jì)了一個(gè)向右算法，計(jì)算對(duì)稱復(fù)系數(shù)矩陣方程組，該算法的分解操作是在2個(gè)稠密下三角矩陣內(nèi)完成的，可以利用基本線性代數(shù)系統(tǒng)庫(kù)函數(shù)(basic linear algebra subprograms，BLAS)［10］，達(dá)到很高的運(yùn)算峰值。并且，在計(jì)算中，符號(hào)分解后引入縮放矩陣改善矩陣條件數(shù)以保證數(shù)值分解數(shù)值穩(wěn)定性和之后的迭代計(jì)算的收斂速度。

預(yù)處理中采用的不完全分解法是在ECM/MF方法的基礎(chǔ)上應(yīng)用舍棄策略。文獻(xiàn)［11］提出的舍棄策略在第j列運(yùn)算中保留了分解因子的嚴(yán)格下三角部分的第j列中的coll個(gè)元素，coll為該列對(duì)應(yīng)系數(shù)矩陣A嚴(yán)格下三角部分該列的非零元個(gè)數(shù)，這種算法可以預(yù)測(cè)存貯空間并考慮了元素?cái)?shù)值的大小，但這種情況下可能舍棄元素過(guò)多，產(chǎn)生的預(yù)條件子顯得過(guò)于粗糙，可能會(huì)導(dǎo)致迭代多次才能收斂或者不收斂。文獻(xiàn)［12］提出的舍棄策略使分解因子擁有額外填充元，該算法保留了分解因子的嚴(yán)格下三角部分的第j列中的coll+p個(gè)元素，p是一個(gè)大于0的控制存貯空間的參數(shù)，然而這種情況下可能舍棄的元素又太少，從而使得預(yù)處理時(shí)間過(guò)長(zhǎng)。本文結(jié)合上述2種方法的優(yōu)點(diǎn)提出一種舍棄策略既能夠保留部分填充元又能夠控制預(yù)處理的時(shí)間。

針對(duì)(1)式中復(fù)系數(shù)矩陣對(duì)稱的特點(diǎn)，本文考慮基于ECM/MF法構(gòu)造一個(gè)高效率的不完全分解預(yù)處理方法。該算法是基于MF在超節(jié)點(diǎn)消去樹(shù)的指導(dǎo)下在稠密矩陣內(nèi)計(jì)算的原理，同時(shí)引入縮放矩陣以改善系數(shù)矩陣的條件數(shù)，并應(yīng)用提出的舍棄策略生成預(yù)條件子，然后結(jié)合迭代法計(jì)算采用矢量有限元技術(shù)分析電磁問(wèn)題產(chǎn)生的方程組。

1 預(yù)處理矩陣的構(gòu)造

考慮電磁問(wèn)題中(1)式的系數(shù)矩陣是復(fù)對(duì)稱矩陣，A=LLT－E是一個(gè)不完全分解，L是下三角矩陣，則預(yù)處理線性系統(tǒng)可以表示為

預(yù)處理的目的是使得L－1AL－T收斂比原始矩陣A快，即迭代矩陣特征譜IN－L－1AL－T集中在0附近。下面首先簡(jiǎn)要說(shuō)明ECM/MF分解算法，并給出舍棄策略，然后基于ECM/MF算法和舍棄策略提出不完全分解預(yù)處理方法并描述分解中用到的優(yōu)化算法，最后給出用于迭代的方程組。

1.1 ECM/MF分解算法

ECM/MF分解算法在第k步計(jì)算的過(guò)程可由圖1形象的給出，陰影部分為計(jì)算第k步時(shí)已被存貯在硬盤(pán)文件中的部分，黑色部分k為第k步正在分解的部分，黑色部分q為第k步正在更新的部分。采用ECM/MF分解法分解矩陣優(yōu)勢(shì)有2點(diǎn):1)下三角矩陣L按列存貯，便于存貯和處理;2)每次只對(duì)一列元素進(jìn)行分解操作，易于加入舍棄策略，從而構(gòu)造不完全分解矩陣。

1.2 舍棄策略

為了保證預(yù)條件子不過(guò)于粗糙，同時(shí)避免預(yù)條件時(shí)間過(guò)長(zhǎng)，本文以節(jié)點(diǎn)q為例提出的舍棄策略如下。

For j=1，…，m

dκ =τ‖low{A*j}‖F(xiàn)

判斷Sq的第j列

判斷1 如果該列有小于coll個(gè)非零元的

絕對(duì)值都≤dκ，則保留該列coll個(gè)最大非零元;

判斷2 否則取t=min(coll2，coll+p)，

如果有大于t個(gè)非零元的絕對(duì)值都滿足不小于某一

數(shù)值dκ，則保留該列t個(gè)最大非零元;

判斷3 否則

保留該列所有≥dκ的非零元;

判斷結(jié)束。

End for

其中j表示節(jié)點(diǎn)q相鄰子節(jié)點(diǎn)的局部列編碼;m是節(jié)點(diǎn)q相鄰子節(jié)點(diǎn)數(shù);* 表示A矩陣j列所有非零元的行號(hào);Sq表示節(jié)點(diǎn)q的波前陣;coll2表示包括填充元在內(nèi)的波前陣該列的非零元;p表示分解因子L中每列最大的填充元個(gè)數(shù)，一般取 p=nz(L)/N，nz(L)是分解因子L的非零元個(gè)數(shù)，N是矩陣K的階數(shù)，K表示正在分解的矩陣;τ是舍棄容忍因子。預(yù)條件子的第j列保留的非零元素的個(gè)數(shù)至少為coll個(gè)元素，最多為t個(gè)元素。最后每列具體保留多少個(gè)，通過(guò)參數(shù)τ來(lái)確定，這樣就不會(huì)浪費(fèi)一些不必要的計(jì)算時(shí)間和存儲(chǔ)空間。注意，如果選取的參數(shù)τ接近無(wú)限大的話，那么預(yù)條件子的第j列保留的元素個(gè)數(shù)就約等于coll個(gè)，這恰好是文獻(xiàn)［11］提出來(lái)的算法。相反，如果選取的參數(shù)τ接近無(wú)限小或者為0，且coll2≥coll+p時(shí)，那么保留的元素個(gè)數(shù)就等于coll+p個(gè)，這恰好是文獻(xiàn)［12］提出來(lái)的算法。

圖1 ECM/MF分解法計(jì)算模式Fig.1 Decomposition model of ECM/MF

1.3 ECM/MF不完全分解算法

在ECM/MF算法基礎(chǔ)上，引入縮放矩陣并應(yīng)用提出的舍棄策略，構(gòu)造出下列ECM/MF不完全分解算法。算法的流程如下。

1)根據(jù)系數(shù)矩陣的稀疏程度重排序方程組，采用反向卡特希爾－麥基(reverse Cuthill－Mckee，RCM)重排序稀疏矩陣系統(tǒng)，采用最小自由度(minimum degree，MMD)方法重排序中等稀疏和稠密矩陣系統(tǒng)，并令重排序后的矩陣為Aox=bo。

2)根據(jù)系數(shù)矩陣的正定程度計(jì)算對(duì)角矩陣P來(lái)縮放原系數(shù)矩陣，當(dāng)原系數(shù)矩陣是正定系統(tǒng)時(shí)，P的對(duì)角元素為原系數(shù)矩陣的對(duì)角元素的幅值，否則，采用文獻(xiàn)［13］的方法計(jì)算，使 ‖P－1AoP－1‖j∞∈［1 －μ，1+μ］，其中0μ1，‖‖j∞表示第j列的無(wú)窮范數(shù)，Ao表示A優(yōu)化重排序后的矩陣。令縮放后的方程為 Ky=f，其中 K=P－1AoP－1，y=Px，f=P－1bo，bo表示b優(yōu)化重排序后的向量。

3)采用ECM/MF方法計(jì)算方程組，以第j列為例，對(duì)波前陣F(j)用Submatrix－ECM分解，分解過(guò)程如下

(3)式中:Uj是分解進(jìn)行中的ir×ir矩陣;Low{}表示取括號(hào)中矩陣包含對(duì)角元的下三角矩陣。圖2給出了稠密矩陣空間Q1和Q2中的數(shù)據(jù)分解順序。Q1和Q2內(nèi)的箭頭表示在左邊節(jié)點(diǎn)波前陣計(jì)算結(jié)束后右邊節(jié)點(diǎn)未更新的波前陣在該矩陣內(nèi)替代左邊節(jié)點(diǎn)波前陣數(shù)據(jù)，兩個(gè)矩陣之間的箭頭表示子節(jié)點(diǎn)更新父節(jié)點(diǎn)。箭頭附近的數(shù)字表示計(jì)算操作順序。節(jié)點(diǎn)1更新陣一旦形成立刻集成到其父節(jié)點(diǎn)4上，分解因子向量存入硬盤(pán)文件中，然后用節(jié)點(diǎn)2數(shù)據(jù)替換節(jié)點(diǎn)1數(shù)據(jù)，以此類推完成其他節(jié)點(diǎn)分解。也就是說(shuō)，僅在2個(gè)稠密下三角矩陣內(nèi)就可以完成所有節(jié)點(diǎn)分解操作。

圖2 消去樹(shù)和兩個(gè)稠密矩陣中的分解順序Fig.2 Decomposition order of elimination tree and Execution order of two dense matrices

4)對(duì)波前陣執(zhí)行舍棄策略，則不完全分解后的L矩陣就是預(yù)條件子。

1.4 迭代求解

在預(yù)處理矩陣構(gòu)造完畢后，采用迭代法求解矩陣方程Koz=fo，其中Ko=L－1KoL－1，z=Ly，fo=Lf。

基于上述縮放矩陣和舍棄策略的應(yīng)用，本文實(shí)施的這種基于ECM/MF的不完全分解的預(yù)處理技術(shù)比傳統(tǒng)的不完全分解預(yù)處理技術(shù)具有更好的穩(wěn)定性和通用性。

2 數(shù)值算例

為了驗(yàn)證開(kāi)域有限元方程組采用MF不完全分解預(yù)處理迭代法計(jì)算的正確性以及所編程序的可靠性，本文分別基于FEM/ABC，F(xiàn)EM/PML和FEM/BI方法采用提出的不完全分解預(yù)處理共軛梯度法(conjugate gradient，CG)計(jì)算了多種目標(biāo)的雙站或后向雷達(dá)散射截面圖(radar cross section，RCS)。算法用Fortran編寫(xiě)，預(yù)處理CG法的收斂精度為1E－4，運(yùn)行環(huán)境為 CPU 2.4 GHz，1.0 GByte RAM。

2.1 基于FEM/ABC方法計(jì)算RCS

分別計(jì)算尺寸為1.5λ×1.0λ×1.0λ理想導(dǎo)電矩形腔體后向散射和邊長(zhǎng)為0.755λ金屬立方體的雙站散射截面(θinc=180°，φinc=90°)問(wèn)題。對(duì)于理想導(dǎo)電矩形腔體，虛擬邊界為與腔體同中心，且半徑r=1.35λ的球;對(duì)于金屬立方體，虛擬邊界為與腔體同中心，且半徑r=0.855λ的球。圖3給出理想導(dǎo)電腔體的后向散射截面(HH極化)和金屬立方體雙站雷達(dá)散射截面圖，圖3中實(shí)線和虛線是不完全分解預(yù)處理LG法的結(jié)果，空心方塊是直接法求解的結(jié)果。由圖3可以看出，本文方法與直接法計(jì)算結(jié)果一致。表1對(duì)比了本文不完全分解預(yù)處理CG法與直接法計(jì)算結(jié)果，表1中填充元個(gè)數(shù)比是不完全分解與完全分解的填充元個(gè)數(shù)比，可以看出，采用本文方法所需的內(nèi)存遠(yuǎn)小于直接法，且計(jì)算時(shí)間明顯優(yōu)化。

圖3 理想矩形導(dǎo)電腔體后向散射截面(HH極化)和金屬立方體雙站雷達(dá)散射截面圖Fig.3 Backscatter RCS of a perfectly conducting rectangular inlet and bistatic RCS of a metallic cube

表1 不完全分解預(yù)處理CG法與直接法計(jì)算結(jié)果對(duì)比Tab.1 Comparison of incomplete factorization preconditioned CG and direct results

2.2 基于FEM/PML方法計(jì)算方向圖

半波偶極天線PML模型如圖4所示。圖4中，PML層厚度為0.1λ，空氣層厚度為0.15λ，方程組有237 824個(gè)未知量，2 283 402個(gè)非零元。圖5給出了半波偶極天線的方向圖，可見(jiàn)，采用不完全分解預(yù)處理CG法和直接求解的結(jié)果一致。本文方法與直接法的時(shí)間比值約為1/10，由此說(shuō)明，在節(jié)省計(jì)算時(shí)間上，本文方法具有明顯的優(yōu)勢(shì)。

2.3 基于FEM/BI方法計(jì)算RCS

半徑r=0.5λ，介電常數(shù)εr=2.0的均勻介質(zhì)球的雙站雷達(dá)散射截面。采用內(nèi)觀法將FEM/BI方程組分寫(xiě)成兩部分，稀疏矩陣部分采用RCM重排序，稠密矩陣部分采用 MMD重排序，矩陣乘用BLAS 2里面的DGMM命令計(jì)算。有限元部分有28 394個(gè)未知量，邊界元部分有3 126個(gè)未知量。數(shù)值結(jié)果與Mie解析解比較見(jiàn)圖6。

圖6 均勻介質(zhì)球的雙站雷達(dá)散射截面圖(θθ極化)Fig.6 Bistatic RCS of homogeneous dielectric sphere

3 結(jié)束語(yǔ)

本文提出了一種基于ECM/MF法的不完全分解預(yù)處理方法，該算法能夠保證預(yù)條件子足夠精細(xì)同時(shí)節(jié)約了構(gòu)造時(shí)間。該算法改善了矩陣條件數(shù)和內(nèi)存占用性能。數(shù)據(jù)結(jié)果表明，本文不完全分解預(yù)處理方法結(jié)合迭代法能夠成功的應(yīng)用于求解矢量有限元方程組，并在很大程度上改善存貯性能且提高計(jì)算速度。

［1］JIN J M.The finite element method in electromagnetics［M］.USA:Wiley－IEEE Press，1993.

［2］DAVIDSON D B，BOTHA M M.Evaluation of a spherical pml for vector FEM applications［J］.IEEE Transactions on Antennas and Propagation，2007，55(2):494－498.

［3］LI P.Coupling of finite element and boundary integral methods for electromagnetic scattering in a two－layered medium［J］.Journal of Computational Physics，2010，229(2):481－497.

［4］LI S S，RUI P L，CHEN R S.An effective sparse approximate inverse preconditioner for vector finite element analysis of 3d em problems［C］//Antennas and Propagation Society International Symposium［s.l.］:IEEE Press，2006.

［5］PING X W，CUI T.The Factorized sparse approximate inverse preconditioned conjugate gradient algorithm for finite element analysis of scattering problems［J］.Progress In Electromagnetics Research，2009(98):15－31.

［6］LI L，HUANG T，JING Y，et al.Application of the incomplete Cholesky factorization preconditioned Krylov subspace method to the vector finite element method for 3－D electromagnetic scattering problems［J］.Computer Physics Communications，2010，181(2):271－276.

［7］SHENG xin－qing，PENG zhen.Further cognition of hybrid FE/BI/MLFMA ［J］.Acta Electronica Sinica.2006，34(1):93－98.

［8］TIAN Jin，LV Zhi－qing，SHI Xiao－wei，et al.An efficient approach for multifrontal algorithm to solve non－positive－definite finite element equations in electromagnetic problems［J］.Progress In Electromagnetics Research，2009(95):121－133.

［9］QU Y，F(xiàn)ISH J.Multifrontal incomplete factorization for indefinite and complex symmetric systems［J］.International journal for numerical methods in engineering，2002，53(6):1433－1459.

［10］LIU J W H.The multifrontal method for sparse matrix solution:Theory and practice［J］.Society for Industrial and Applied Mathematics，1992，34(1):82－109.

［11］JONES M T，PLASSMANN P E.An improved incomplete Cholesky factorization［J］.ACM Transactions on Mathematical Software(TOMS)，1995，21(1):5－17.

［12］LIN C J，MORé J J.Incomplete Cholesky factorizations with limited memory［J］.SIAM Journal on Scientific Computing，2000，21(1):24－45.

［13］QU Y，F(xiàn)ISH J.Multifrontal incomplete factorization for indefinite and complex symmetric systems［J］.International journal for numerical methods in engineering，2002，53(6):1433－1459.