自適應(yīng)交叉近似壓縮的高階矩量法的并行實(shí)現(xiàn)

袁浩波　何力　黨曉杰　王志軍

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

?

自適應(yīng)交叉近似壓縮的高階矩量法的并行實(shí)現(xiàn)

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

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

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

關(guān)鍵詞高階矩量法;ACA壓縮算法;共軛梯度法;并行計(jì)算

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

引言

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

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

1ACA壓縮的理論

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

Zm×n≈Um×rVr×n,

(1)

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

‖Z-UV‖≤ε‖Z‖,

(2)

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

圖1　用于分組的模型A

圖2　需進(jìn)行電磁計(jì)算的模型B

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

2并行ACA實(shí)現(xiàn)

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

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

圖3　并行ACA算法流程

圖4　并行ACA壓縮時(shí)各進(jìn)程的任務(wù)分配

圖5　進(jìn)程0中矩陣向量乘積運(yùn)算

3計(jì)算實(shí)例

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

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

表1　不同壓縮門限時(shí)采用768核并行計(jì)算導(dǎo)體球RCS

圖6　不同壓縮門限時(shí)計(jì)算的導(dǎo)體球在xoz面的RCS

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

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

圖8　導(dǎo)彈模型

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

表2　不同壓縮門限時(shí)采用980核并行計(jì)算導(dǎo)彈RCS

圖9　不同壓縮門限時(shí)計(jì)算的導(dǎo)彈在xoz面的RCS

圖10　導(dǎo)彈在ACA壓縮時(shí)的并行效率(ε=0.001)

4結(jié)論

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

參考文獻(xiàn)

[1]袁軍, 邱揚(yáng), 劉其中, 等. 自適應(yīng)多層快速多極子算法及其并行算法[J]. 電波科學(xué)學(xué)報(bào), 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] 吳君輝, 曹祥玉, 袁浩波, 等. 一種電大目標(biāo)散射特性的核外并行快速算法[J]. 電波科學(xué)學(xué)報(bào), 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]. 北京: 科學(xué)出版社, 2011:393-396.

袁浩波(1980-),男,湖北人,西安電子科技大學(xué)副教授,博士,研究方向?yàn)殡姶艌?chǎng)數(shù)值計(jì)算.

何力(1989-),男,四川人,碩士研究生,研究方向?yàn)殡姶艌?chǎng)數(shù)值計(jì)算.

黨曉杰(1980-),男,內(nèi)蒙古人,西安電子科技大學(xué)講師,研究方向?yàn)殡姶判虏牧霞夹g(shù).

王志軍(1963-),男,山西人,中北大學(xué)教授,博士生導(dǎo)師,研究方向?yàn)殪`巧彈藥技術(shù)、彈箭毀傷控制技術(shù)、計(jì)算機(jī)仿真與實(shí)驗(yàn)研究等.

作者簡(jiǎn)介

中圖分類號(hào)TN011

文獻(xiàn)標(biāo)志碼A

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

收稿日期：2015-02-07

袁浩波, 何力, 黨曉杰, 等. 自適應(yīng)交叉近似壓縮與高階矩量法的并行實(shí)現(xiàn)[J]. 電波科學(xué)學(xué)報(bào),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

資助項(xiàng)目: 國(guó)家自然科學(xué)基金 (61072018，60901030)； 中國(guó)博士后基金(2014M 561211)； 中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金資助(JB150223，WRYB142105).

聯(lián)系人： 袁浩波 E-mail：useryuanhaobo@163.com