1.   引言

矩量法(Method of Moments, MoM)是电磁数值分析方法中的经典算法之一，也是高精度的数值算法之一，在电磁仿真计算中得到了广泛的应用^[1]。频域矩量法需要建立并求解复数稠密矩阵方程，为了保证求解精度，采用基于高斯消元法的LU分解算法直接得到矩阵方程的解，其存储复杂度和计算复杂度分别为${\rm O}({N^2})$和${\rm O}({N^3})$，其中N为矩阵维数。当计算电大尺寸目标的电磁问题时会生成一个未知量十分庞大的复数稠密矩阵，此时，矩阵LU分解的时间占整个矩量法计算过程的90%以上^[2]。因而，巨大的存储需求和过长的矩阵求解时间成为限制矩量法计算电大尺寸电磁问题的主要因素。突破该瓶颈的常用方法之一是利用超级计算机集群进行大规模并行计算^[3,4]。

随着计算机技术的飞速发展，针对高度并行工作负载而设计的众核处理器，依靠其集成的大量计算单元，以更小的功耗，展现出比同期多核处理器更加强大的计算能力，特别适合处理矩阵LU分解这种数据密集型、计算密集型的高并发度任务。此前，矩量法已在国内超级计算机平台中的国产处理器^[5,6]，通用处理器^[7,8]以及异构众核处理器^[9－12]上开展过大量的研究工作。但是，利用的众核处理器平台都是采用主机+设备(Host+Device)的异构模式，即Host端采用多核的CPU(Center Process Unit)控制任务的划分和数据传输，同时执行部分任务的计算。Device端的众核处理器以协处理器或加速卡的形式出现，执行大部分的高并发度任务，从而达到加速的效果^[13]。但是，该种异构计算模式的缺陷在于Host端和Device端的通信都是通过PCI-E接口实现的，通信速率较慢，限制了计算性能的提升。此外，在Host+Device异构超算平台进行算法的开发难度大、周期长，不利于应对复杂多变和迫切的科研需求。最为关键的是，国内大多数异构众核平台其加速卡或协处理器过分依赖国外，算法的性能提升和安全“受制于人”。

为保证科研工作的安全，近年来，国内加快了对超级计算机特别是众核超级计算机的研制。“天河二号”超级计算平台，采用自主研发的FT2000+众核处理器，单核性能强大，是国内首款独立自主研发的同构众核超算平台，其同构计算系统使算法的研发和优化大为简化，在保证科研成果安全的同时，大大加快了科学研究的进度。矩量法在“天河二号”的研发，弥补了国产超算平台中电磁软件的不足，实现了电磁仿真的自主可控和安全可靠。

然而，当矩量法在超算平台进行大规模并行计算时，对于计算复杂、通信复杂的矩阵求解过程，进程间密集而频繁地通信对计算平台造成的通信压力和通信延迟是我们应该正视和亟待解决的问题。虽然，文献[14], [15]提出单向的避免通信矩阵LU分解(Communication Avoiding LU decomposition, CALU)算法，通过改变选主元策略，减少了panel列分解过程中的通信次数和通信量。而且，在通用处理器平台对CALU算法做了大量的测试工作和深入分析，其测试结果表明，CALU算法比公认的最快的商业数学库(Math Kernel Library, MKL)^[16]表现出更快的矩阵求解速度和更好的扩展性。但是，CALU算法增加了计算量，主元依然是在整个panel列中进行选取。

本文基于FT2000+同构众核处理器平台，通过分析矩量法生成的矩阵具有对角占优特性，提出一种新型的LU分解算法，即对角块矩阵选主元LU分解(Block Diagonal Pivoting LU decomposition, BDPLU)算法，该算法结合矩量法特性，改变了传统选主元策略，在未增加panel列分解计算量的前提下，完全消除了选主元过程的通信，进一步加速了矩阵求解过程，并验证了算法的正确性。

2.   算法分析

2.1   并行矩量法

本文矩量法采用的基函数为RWG(Rao-Wilton-Glisson)基函数，是现今广泛使用的一种矩量法基函数，RWG基函数可模拟任意形状物体的表面电、磁流分布，采用伽略金方法可得到矩量法矩阵方程^[4]

其中，[Z_mn]表示大小为N×N的阻抗矩阵，[I_n]和[V_m]分别表示电流列向量和激励电压列向量。为了保证计算精度，RWG基函数三角形剖分的网格边长一般在$[\lambda /12,\lambda /12]$范围内^[4]。为了寻求更为高效地求解矩量法大型复数稠密矩阵方程的新方法，以下从矩量法产生的阻抗矩阵特性进行分析和推导。

2.2   矩阵特性分析

矩量法生成的阻抗矩阵[Z_mn]如公式(2)所示^[4]，其中，${r_t}$和${r_s}'$分别为场三角形和源三角形上的采样点，W(t)和W(s)为采样点对应的权值，${{v_m}^ \pm }$和${{v_n}^ \pm }$分别是第m和n条公共边对应的顶点坐标，$G\left( R \right) =$${{\rm{e}}^{ - {\rm{j}}kR}}/4{{π}} R$为自由空间格林函数，$R = \left| {{r_t} - {r_s}'} \right|$。分析式(2)可知，当${r_t} = {r_s}'$时，即场三角形和源三角形重合时，[Z_mn]的值最大，因而可得MoM阻抗矩阵具有对角占优特性。为进一步证明该结论的正确性，图1给出了积分方程离散生成的复数稠密矩阵元素分布。由图1对矩阵元素的描述可以看出，LU分解前，原始矩阵其较大值主要分布在对角线附近。当完成一次panel列分解，panel行更新、trailing更新后，矩阵的元素分布仍然具有对角占优特性，直到整个矩阵分解完成。对于矩量法生成的矩阵，主元主要分布在对角线附近，而CALU算法在整个panel列中进行选主元无疑产生了冗余的通信量和计算量。

图 1  LU分解过程矩阵特性分布

下载: 全尺寸图片 幻灯片

2.3   BDPLU算法

BDPLU算法的前提是矩阵具有对角占优特性。算法原理如图2所示，假设图2(a)中矩阵的大小为N×N，分块矩阵大小为n_b×n_b，采用2维块循环分布方式分布在进程网格P_r×P_c中，当前panel列矩阵的对角分块矩阵B_ij分布在第P_r(i)行进程和第P_c(j)列进程中，主元只从B_ij分块矩阵中进行选取，参与选主元的进程仅为(P_r(i), P_c(j)，选主元过程完全避免了通信。而且行交换只有panel行中的进程参与执行，其他行进程不参与行交换，减少了行交换的通信。图2(b)中的(P_r(i), P_c(j)进程向下广播U_kk矩阵到其他行进程，$\left( {{P_r}_{(i + 1)} \sim {P_r}_{(i + n)},{P_c}_{\left( j \right)}} \right)$进程执行完图2(c)所示的消元与更新操作后，整个panel列分解计算完成。

图 2  BDPLU算法原理图

下载: 全尺寸图片 幻灯片

以矩量法生成的双精度复数稠密矩阵LU分解为例，分析CALU算法和BDPLU算法的通信时间。假设计算机发送1次消息的通信延迟为$\alpha$，发送1个矩阵元素的时间为$\beta$。式(3)和式(4)分别给出了CALU算法和BDPLU算法panel列分解的通信时间T_CALU^[4]和T_BPULU。可以得出，BDPLU算法其panel列分解的通信时间比CALU算法减少了$\ge {1 / 3}$，这对LU分解性能的提升至关重要。

为进一步证明BDPLU算法的正确性和高效性。此处以飞机模型I的散射特性计算为例进行研究^[17,18]。飞机I仿真模型如图3(a)所示，平面波沿机头方向入射，沿Z轴方向极化，频率为500 MHz，矩量法所产生的未知量为58646，矩阵方程分别用BDPLU算法和CALU算法进行求解，飞机I的双站雷达散射截面(Radar-Cross Section, RCS)结果如图3(b)和图3(c)所示，并与商业软件FEKO计算结果完全吻合，证明了BDPLU算法和CALU算法在“天河二号”国产超算平台进行电磁散射计算的正确性。表1给出了BDPLU算法和CALU算法求解矩阵的时间和效率，可以看出采用不同并行规模时，BDPLU算法总能表现出比CALU算法更快的矩阵求解速度和更好的扩展性。

图 3  飞机I仿真模型和双站RCS结果

下载: 全尺寸图片 幻灯片

表 1  CALU算法与BDPLU算法矩阵求解时间对比

FT2000+核数	矩阵求解时间(s)			并行效率(%)
	CALU	BDPLU		CALU	BDPLU
2000	796.54	742.57		100	100
3000	567.78	518.68		93.53	95.44
4000	463.93	421.12		85.85	88.17
5000	386.89	338.24		82.35	87.82
10000	226.57	187.83		70.31	79.07
15000	172.91	139.05		61.42	71.20
20000	133.97	118.38		59.46	62.73
40000	72.83	64.61		54.68	57.47

下载: 导出CSV 

| 显示表格

3.   计算平台简介

本文使用的“天河二号”国产同构众核处理器超算平台，共有12个机柜，每个机柜有1440个FT2000+计算节点，全系统共17280个FT2000+计算节点，节点间采用25 G高速互连网络，全系统的峰值性能为10.17 PFlops (PetaFlops, 10¹⁵次浮点运算/s)。

FT2000+是国内首款自主研发的同构众核处理器，采用16 nm工艺，工作频率为2.3 GHz，每个处理器64核心，双精度峰值为588.8 GFlops，片内集成8个DDR4存储访问通道，Stream带宽初测为89 GB/S, DGEMM效率为92%, LINPACK效率为86%，操作系统为具有完全自主版权的“麒麟”(kylin)高性能操作系统^[19]。

4.   数值算例

本文基于“天河二号”国产众核超算平台，利用并行矩量法，采用BDPLU算法求解矩阵方程，计算飞机模型II的散射特性，验证了并行矩量法可在“天河二号”国产众核超算平台进行大规模并行计算且能够有效解决电大尺寸目标电磁问题，并给出矩阵求解并行效率。

飞机II仿真模型如图4(a)所示，其几何尺寸为18.92 m×14.56 m×5.05 m，平面波入射频率为400 MHz，相应的电尺寸为$25.23\lambda \times 19.41\lambda \times 6.73\lambda$，平面波沿机头方向入射，极化方向为垂直极化，计算其双站RCS，在该频率下矩量法生成一个未知量为190542×190542的复数稠密矩阵。图4(b)、图4(c)给出了xoy面和xoz面的2维双站RCS结果^[17,18]。表2给出了求解该矩阵的时间和并行效率，60万核规模是目前在国产超级计算平台中实现的最大规模的并行矩量法计算。

图 4  飞机II双站RCS结果

下载: 全尺寸图片 幻灯片

表 2  BDPLU算法求解矩阵的加速比和并行效率

FT2000+核数	矩阵求解时间(s)	加速比	并行效率(%)
9600	29183.33	1	100
48000	6336.53	4.61	92.11
96000	3501.73	8.33	83.34
192000	2035.35	14.34	71.69
240000	1764.53	16.54	66.16
336000	1328.13	21.97	62.78
384000	1227.91	23.77	59.42
432000	1133.64	25.74	57.21
480000	1043.45	27.97	55.94
504000	997.94	29.24	55.70
552000	937.82	31.12	54.12
600000	898.81	32.49	51.95

下载: 导出CSV 

| 显示表格

图5给出了该算法的并行性能曲线，以9600核为基准，当并行规模扩展到60万核时，并行效率仍能达到51.95%，展现出了BDPLU算法良好的可扩展性。

图 5  加速比和并行效率

下载: 全尺寸图片 幻灯片

5.   结束语

本文基于“天河二号”国产众核超级计算平台，针对并行矩量法中通信复杂、计算复杂的矩阵求解过程，分析矩量法电场积分方程离散生成的矩阵具有对角占优特性，提出了BDPLU算法，完全消除了LU分解过程中panel列分解选主元的通信开销，加速了复数稠密矩阵的求解，缓解了大规模并行计算时计算集群的通信压力，实现了国产众核超算平台中并行矩量法的大规模并行计算。测试结果表明：BDPLU算法在保证结果精确的同时，比单向通信的CALU算法具有更好的性能。利用BDPLU算法，并行矩量法在“天河二号”国产众核超算平台率先突破了6×10⁵CPU核的并行计算，并展现出了良好的并行性能。本文的研究工作为今后矩量法在完全自主研发的国产超算平台中解决更多具有实际工程意义的复杂电磁问题奠定了基础。