自适应交叉近似压缩的高阶矩量法的并行实现

袁浩波　何力　党晓杰　王志军

(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