自適應交叉近似壓縮的高階矩量法的并行實現

袁浩波　何力　黨曉杰　王志軍

(1.中北大學,太原 030051;2.西安電子科技大學電子工程學院,西安 710071)

?

自適應交叉近似壓縮的高階矩量法的并行實現

袁浩波1,2何力2黨曉杰2王志軍1

(1.中北大學,太原 030051;2.西安電子科技大學電子工程學院,西安 710071)

摘要高階矩量法在計算電磁學中的應用越來越廣泛,為了進一步提高其計算規模,引入并行的自適應交叉近似壓縮算法(Adaptive Cross Approximation algorithm, ACA).該算法首先采用非均勻有理B樣條建模(Non-Uniform Rational B-Splines, NURBS)的方法進行面片分組;然后利用矩量法中遠區阻抗矩陣的低秩特性進行ACA壓縮;最后采用稀疏近似逆預條件(Sparse Pattern Approximate Inverse preconditioning, SPAI)的共軛梯度法(Conjugate Gradient method, CG) 快速求解矩陣方程.該算法中的ACA壓縮過程和迭代求解過程都特別適合并行計算.數值實驗表明,對于電大尺寸問題,ACA壓縮后的矩陣占用的內存遠遠低于原矩陣,而預條件的共軛梯度法可以很快收斂.此外該算法在大規模并行時的效率較高.

關鍵詞高階矩量法;ACA壓縮算法;共軛梯度法;并行計算

DOI10.13443/j.cjors.2015020701

A parallelized higher order moment method combined with the ACA compressing

YUAN Haobo1,2HE Li2DANG Xiaojie2WANG Zhijun1

(1.SchoolofMechano-ElectronicEngineering,NorthUniversityofChina,Taiyuan030051,China;2.SchoolofElectronicEngineering,XidianUniversity,Xi’an710071,China)

AbstractThe higher order moment method is widely applied in the computational electromagnetics. In order to compute the electrically massive problems, this paper introduces a parallel adaptive cross approximation algorithm(ACA) to accelerate the higher order moment method. At first, the non-uniform rational B-Splines modeling (NURBS) is applied to divide the patches into groups. Then the ACA algorithm is used to compress the impedance matrix in the far zone, which is low in rank. Finally, the conjugate gradient method(CG) combined with the sparse pattern approximate inverse preconditioning(SPAI) is used to solve the matrix equation. Both the ACA compressing and the CG method are suitable for parallel computation. Numerical experiments show that the memory of the compressed matrix is much less than that of the original matrix, and the preconditioned CG method converges very fast. Besides, the massively parallel method often has a high efficiency.

Keywords higher order moment method; ACA compressing; conjugate gradient method; parallel computing

引言

盡管多層快速多極子技術[1]使得傳統低階矩量法可以求解大規模的電磁問題,但是該技術與具體問題的積分核相關,而且并行化難度很高.而自適應交叉近似算法(Adaptive Cross Approximation algorithm, ACA)是一種非常簡單的線性代數算法,與積分核無關,可以很方便地移植到任何矩量法代碼中,特別適合并行計算.ACA算法于2000年由Bebendorf首次提出[2],它將大的矩陣分解為多層塊矩陣,其中低秩的塊矩陣可以通過一個類似LU分解的過程進行壓縮.2005年李金發首次將ACA算法應用于矩量法中[3],他所給出的ACA算法流程可以很簡單地移植到任何新算法中.2008年Astner解決了并行ACA壓縮在低階矩量法中使用的負載均衡問題[4].2009年麻連鳳提出對高階矩量法的矩陣采用一種局部ACA方法[5]進行壓縮,從而提高壓縮效果.2013年吳君輝采用并行核外技術[6]提高ACA壓縮的低階矩量法的計算效率.2014年晏嬰[7]采用并行ACA技術加速時域矩量法,對于三角形面片采用八叉樹分組.但是上述工作[4,6-7]中的并行規模都很小,難以用于實際的電大尺寸問題.

文獻[8]的高階矩量法采用非均勻有理B樣條建模(Non-Uniform Rational B-Splines, NURBS)結合多層高階基函數求解電場積分方程,其優勢是產生的未知數可以比低階矩量法的未知數少一個數量級.在此基礎上,建立并行ACA壓縮的高階矩量法,目的是將ACA壓縮算法移植到高階矩量法中,并通過大規模并行計算使其能夠求解電大尺寸模型的電磁散射問題.

1ACA壓縮的理論

將ACA算法用于矩量法中時,假定有兩組相距較遠的面片.第一組的若干個面片上定義m個基函數,第二組的若干個面片上定義n個基函數,它們之間的互阻抗矩陣為Zm×n.該矩陣可以近似為兩個矩陣的乘積

Zm×n≈Um×rVr×n,

(1)

式中r稱為矩陣Zm×n的有效秩.ACA算法的目標是使得近似矩陣的相對誤差低于某個門限ε,即

‖Z-UV‖≤ε‖Z‖,

(2)

式中的矩陣范數都是F范數.由于矩量法中遠區互阻抗矩陣的有效秩一般滿足r?min(m,n),因此不需要存儲整個分塊陣的m×n個元素,而只要存儲近似矩陣的(m+n)×r個元素,由此降低存儲空間.ACA壓縮算法一般按照文獻[3]的流程實現,是一種簡單的純線性代數算法,用于低階矩量法時壓縮效果很好.

圖1　用于分組的模型A

圖2　需進行電磁計算的模型B

ACA壓縮算法用于高階矩量法中與用于低階矩量法中有不少區別,其中最大的區別在于面片分組方法不同.文獻[5]中使用一種八叉樹分組方法,其缺點是各組包含的面片數目差距很大,導致并行計算時難以達到負載均衡.本文提出采用兩次NURBS建模的方法進行分組.例如對于一個平板,首先用如圖1所示的9個較大的面片建立模型A,然后將A的每個面片剖分為4個面片從而構成如圖2所示的模型B.其中模型B是需要進行電磁計算的模型,而模型A專門用于給模型B的36個面片分組.只要判斷模型B中每個面片的中心點處于模型A中的第幾個面片上就分到第幾個組.如果兩個模型都剖分得比較均勻,那么每組中包含的面片數目就差不多,因而容易達到負載均衡.

2并行ACA實現

如圖3所示,并行ACA算法主要包括五個步驟,其中關鍵是阻抗矩陣的ACA壓縮和共軛梯度法(ConjugateGradientmethod,CG)迭代求解[9]兩個過程,這兩個過程都特別適合并行計算.在ACA壓縮過程中,假定矩量法的未知數有N個,并將這些未知數分成9組,同時假定進程數目為3個.在并行程序中只需將9×9的分塊矩陣再平均分成如圖4所示的3個橫向條帶,每個進程采用ACA壓縮算法依次填充對應的那個條帶中的27個子陣并存儲.

在矩陣方程的CG求解過程中,由主進程負責耗時較少的主流程計算,而由所有進程共同完成核心的矩陣與向量的乘積運算.在計算矩陣與向量的乘積時,進程0只需要計算其本身存儲的條帶上的壓縮矩陣與向量的乘積,計算結果發送給主進程,如圖5所示.顯然此過程只需少量通信,并行效率很高.為了加快迭代收斂速度,采用了稀疏近似逆預條件技術(SparsePatternApproximateInversepreconditioning,SPAI)[10].計算該預條件矩陣的各個列向量就是求解N個獨立的均方問題,這N個均方問題在并行程序中平均分配給所有進程,各進程之間不需要通信.

圖3　并行ACA算法流程

圖4　并行ACA壓縮時各進程的任務分配

圖5　進程0中矩陣向量乘積運算

3計算實例

首先計算一個半徑為1m的導體球面的散射問題.激勵為x方向極化z方向入射的平面波,波長為0.02m.為了計算模型的雙站雷達散射截面(RadarCrossSection,RCS),首先將該模型剖分成24 576個面片,最大電尺寸為0.57個波長,采用3階基函數,一共得到442 368個未知數.然后將導體球模型剖分為1 536個面片用于ACA分組.并行程序在國家超級計算深圳中心的曙光6 000上進行測試,每個計算節點配置4顆AMD6136八核處理器,主頻2.4GHz,內存128GB.編譯環境采用IntelFortran12.1編譯器和openMPI并行庫.

表1對比了768核的并行ACA算法在3種不同壓縮門限時的求解結果.不同壓縮門限時得到的RCS如圖6所示,將其與MIE級數[11]得到的解析結果對比算出均方根誤差,如表1第6列所示.可見,壓縮門限ε越大則壓縮矩陣占用內存越小,但是所得RCS的精度越低.從表1的第4列可見,預條件的CG只需幾十步迭代即可收斂.從表1的第3列和第5列可見,ACA壓縮過程占用了算法的絕大多數時間,因此該過程的并行效率決定了整個并行算法的計算速度.圖7給出了不同核數時ACA壓縮的并行效率.由于串行程序的計算時間太長而無法得到,這里以64核的計算時間(15.3h)作為基準計算并行效率.可見,核數越多則計算效率越低.

表1　不同壓縮門限時采用768核并行計算導體球RCS

圖6　不同壓縮門限時計算的導體球在xoz面的RCS

圖7　導體球在ACA壓縮時的并行效率(ε=0.001)

接著分析如圖8所示的導彈模型,長3.5 m,機翼寬1.8 m,整個模型的表面積為7.2 m2.激勵為x方向極化z方向入射的平面波,波長為0.02 m.為了計算RCS,首先將該模型剖分成131 220個面片,最大電尺寸為0.48個波長,采用2階基函數,一共得到1 048 896個未知數.然后將圖8中模型剖分為 1 620個面片用于ACA分組.

圖8　導彈模型

表2對比了980核的并行ACA算法在3種不同壓縮門限時的求解結果.不同壓縮門限時得到的RCS如圖9所示.從表2的第2列可見,壓縮門限0.003時占用的內存只有原始矩陣的3.04%.壓縮門限0.000 1時占用的內存幾乎是壓縮門限0.003時的兩倍,但是前者壓縮矩陣時引入的誤差較小,使得其RCS精度比后者高.從表2的第4列可見,預條件的CG方法能夠在大約290步迭代后收斂.圖10以270核的計算時間(16.7 h)作為基準計算ACA壓縮過程的并行效率.圖中324核的計算效率為102%,表明此時并行效率比270核的并行效率高.這主要是由于324核時各進程的任務分配比270核的任務分配更加均衡.隨著核數進一步增加,并行效率逐漸降低,但仍超過90%.

表2　不同壓縮門限時采用980核并行計算導彈RCS

圖9　不同壓縮門限時計算的導彈在xoz面的RCS

圖10　導彈在ACA壓縮時的并行效率(ε=0.001)

4結論

并行ACA算法結合高階矩量法可以求解電大尺寸問題的RCS.該算法在ACA壓縮過程、SPAI預條件矩陣填充,以及CG迭代求解過程中,各個進程之間都不需要或者僅僅需要極少的通信,因此并行效率很高.該算法可以準確求解電大尺寸問題的RCS,具有良好的工程應用前景.為了求解更大規模的問題,可以進一步在高階矩量法中采用并行的多層ACA壓縮算法.

參考文獻

[1]袁軍, 邱揚, 劉其中, 等. 自適應多層快速多極子算法及其并行算法[J]. 電波科學學報, 2008, 23(3): 454-459.

YUAN J, QIU Y, LIU Q Z, et al.Adaptive multilevel fast multipole algorithm and its parallel algorithm [J]. Chinese journal of radio science, 2008, 23(3): 454-459.(in Chinese)

[2] BEBENDORF M. Approximation of boundary element matrices[J]. Numerische mathematik, 2000, 86(4): 565-589.

[3] ZHAO K Z, VOUVAKIS M, LEE J F. The Adaptive cross approximation algorithm for accelerated method of moment computations of EMC problems[J]. IEEE transactions on electromagnetic compatiability, 2005, 47(4): 763-773.

[4] ASTNER M, BRUNS H D, SINGER H. Simple load balancing in binary-tree based parallel multilevel low-rank compression techniques[C]//IEEE International Symposium on Electromagnetic Compatibility. Detroit, August 18-22, 2008.

[5] MA L F, NIE Z P, HU J, et al. Fast direct solution of high-order MoM accelerated by local AC[C]//Asia Pacific Microwave Conference. Singapore, December 7-10, 2009.

[6] 吳君輝, 曹祥玉, 袁浩波, 等. 一種電大目標散射特性的核外并行快速算法[J]. 電波科學學報, 2013, 28(6):1178-1182.

WU J H, CAO X Y, YUAN H B, et al. A parallel out-of-core fast algorithm for scattering characteristic of electrically large target[J]. Chinese journal of radio science, 2013, 28(6): 1178-1182.(in Chinese)

[7] YAN Y, ZHAO X W, LIANG C H, et al. Parallel adaptive cross approximation for accelerating time-domain method of moments[C]//IEEE International Wireless Symposium. Xi’an, March 24-26, 2014.

[8] YUAN H B, WANG N, LIANG C H. Combining the higher order method of moments [J]. IEEE transactions on antennas and propagation, 2009, 57(11): 3558-3563.

[9] SAAD Y. Iterative methods for sparse linear systaems[M]. Boston: PWS Publishing, 1996: 236-237.

[10]ALLEON G, BENZI M, GIRAUD L. Sparse approximate inverse preconditioning for dense linear systems arising in computational electromagnetics[J]. Numerical algorithms,1997, 16:1-15.

[11]葛德彪, 魏兵. 電磁波理論 [M]. 北京: 科學出版社, 2011:393-396.

袁浩波(1980-),男,湖北人,西安電子科技大學副教授,博士,研究方向為電磁場數值計算.

何力(1989-),男,四川人,碩士研究生,研究方向為電磁場數值計算.

黨曉杰(1980-),男,內蒙古人,西安電子科技大學講師,研究方向為電磁新材料技術.

王志軍(1963-),男,山西人,中北大學教授,博士生導師,研究方向為靈巧彈藥技術、彈箭毀傷控制技術、計算機仿真與實驗研究等.

作者簡介

中圖分類號TN011

文獻標志碼A

文章編號1005-0388(2016)01-0138-05

收稿日期：2015-02-07

袁浩波, 何力, 黨曉杰, 等. 自適應交叉近似壓縮與高階矩量法的并行實現[J]. 電波科學學報,2016,31(1):138-142. DOI: 10.13443/j.cjors.2015020701

YUAN H B, HE L, DANG X J, et al. A parallelized higher order moment method combined with the ACA compressing [J]. Chinese journal of radio science,2016,31(1):138-142. (in Chinese). DOI: 10.13443/j.cjors.2015020701

資助項目: 國家自然科學基金 (61072018，60901030)； 中國博士后基金(2014M 561211)； 中央高校基本科研業務費專項資金資助(JB150223，WRYB142105).

聯系人： 袁浩波 E-mail：useryuanhaobo@163.com