2. 英特尔(中国)有限公司, 北京 100013
2. Intel(China) Limited, Beijing 100013, China
随着计算机技术的飞速发展, 针对高度并行工作负载而设计的一种新型处理器, 即众核处理器, 依靠其集成的大量计算单元, 展现出强大的数据处理能力, 在高性能计算领域发挥着越来越重要的作用.与传统通用多核处理器相比, 众核处理器具有成本低、功耗小等显著优势.同时其优势还在于集成的核心数量超出同期通用多核处理器一个数量级, 在处理计算密集型、数据密集型任务时, 例如求解矩量法生成的复数稠密矩阵, 能够展现出强大的优势.
近年来, 国内已有若干利用众核处理器加速矩量法的文献[1-3], 并且取得了较好的成果.但是, 这些文献中利用的众核处理器都是作为协处理器出现的, 如CPU/MIC异构并行矩量法的研究, 由CPU控制任务分配和数据传输, 同时也负责部分计算任务, 将高并发度的代码段在MIC(Many Integrated Core)协处理器执行, 达到加速并行矩量法的目的[3].但是, CPU与MIC协处理之间的通信都是通过PCI-E(Peripheral Component Interconnect Express)接口实现的, 通信速率较慢, 限制了计算性能的提升; 尽管MIC协处理器具有独立的内存和存储, 但容量一般都较小, 因此异构协同计算能解决的实际问题受到协处理器存储容量的限制[4].针对上述问题, Intel于2016年推出了全新一代至强融核处理器, 代号为"骑士登陆"(KNL), 通过消除对PCI-E总线的依赖, 提供更高的可扩展性, 应对更广泛的工作负载和配置.所以, 并行矩量法在KNL众核处理器平台上的计算和优化具有很高的研究价值.
KNL众核处理器是Intel首款专门针对高度并行工作负载而设计的可独立自启动的主处理器:能效型集成架构提供了比类似平台高得多的单元计算能力, 有效减少了计算成本; 内存和结构的集成最大限度地提升了内存容量, 并首次实现了内存与高速互连技术的集成, 为大规模并行和矢量化服务提供了有效平台.
为了充分利用KNL的众核计算资源以及发挥其超宽的矢量宽度优势, 本文基于KNL众核处理器平台对并行矩量法程序开展了优化工作[5-7]; 通过与商业软件计算结果的对比, 验证了优化算法的正确性; 利用KNL计算集群, 仿真了未知量超过20万的飞机模型的散射特性, 程序优化后的性能有较大提升.值得指出的是, 并行矩量法在KNL众核处理器平台上的计算和优化此前尚未见有公开的文献报道.
1 并行矩量法 1.1 矩阵填充分析本文中矩量法采用的基函数为RWG(Rao-Wilton-Glisson)基函数, 是现今广泛使用的一种矩量法基函数, 它被定义在具有公共边的两个相邻三角形上, 可模拟任意形状物体的表面电、磁流分布, 采用伽辽金方法可得到矩量法矩阵方程[8]
| $ \begin{align} ZI=V, \end{align} $ | (1) |
其中, $Z$表示大小为$N\times N$的阻抗矩阵, $I$和$V$均是$N\times 1$的列向量, 二者分别表示电流列向量和激励电压列向量. $Z$中矩阵元素的内外两层面积分可由二维高斯积分来近似得出
| $ \begin{align} Z_{mn}^{\pm \pm } &=\frac{jk\eta }{4\pi }\frac{l_m l_n }{4A_m A_n } \iint{{{\Big[{( {r_t-v^\pm_m } )( {r'_s-v^\pm_n }) -\frac{4}{k^2}} \Big]} \frac{{\rm e}^{-jkR}}{R}{\rm d}s'{\rm d}s \notag}} \\ &=\frac{jk\eta }{4\pi }l_m l_n \Big[{\sum\limits_{t=1}^7 {\sum\limits_{s=1}^7 {W(t)} } W(s)(r_t r'_s-\frac{4}{k^2})\frac{{\rm e}^{-jkR}} {R}-v^\pm_n \sum\limits_{t=1}^7 {\sum\limits_{s=1}^7 {W(t)} } W(s)r_t } \frac{{\rm e}^{-jkR}}{R}\notag\\ &-v^\pm_m\sum\limits_{t=1}^7 {\sum\limits_{s=1}^7 {W(t)} } W(s)r'_s\frac{{\rm e} ^{-jkR}}{R}+v^\pm_m v^\pm_n \sum\limits_{t=1}^7 {\sum\limits_{s=1}^7 {W(t)} } W(s)\frac{{\rm e}^{-jkR}}{R} \Big]. \end{align} $ | (2) |
式(2)中高斯采样点数为7个, $r_t $和${r}'_s $分别为场三角形和源三角形上的采样点, $W(t)$和$W(s)$为采样点对应的权值, $v^\pm_m $和$v^\pm_n $分别是第$m$和$n$条公共边对应的顶点坐标, $R=|r_t-r'_s|$.对于并行矩量法, 由于MPI分布式内存并行策略的引入会使进程间产生冗余积分.
1.2 冗余计算在矩阵填充过程中, 相同的积分计算被不同进程反复多次执行, 称为冗余积分.如图 1所示, 假设为公共边编号与剖分后三角形网格的对应关系, 公共边个数为$N$. 图 2是将$N\times N$个矩阵元素按照二维块循环的方式分布到二维MPI虚拟拓扑网格为$3\times 2$、分块大小为$2\times 2$的矩阵分布情况, 框中的数字代表元素的索引.如果只关心(7, 2)和(7, 3)两个矩阵元素, 则会用到图 1中三角形$KI$的二重积分, 图 2中(7, 2)分配到了进程$P_{00}$上, (7, 3)分配到了$P_{01}$进程上, 这两进程都要计算这一积分, 造成冗余计算.
|
图 1 三角形和公共边编号对应关系 Fig.1 Correspondence between triangle and common edge number |
|
图 2 二维块循环分布 Fig.2 2D block-cyclic distribution |
根据上述分析可知, 矩阵填充过程中的MPI进程越多, 进程间产生的冗余积分计算量越大, 影响了矩阵填充的效率.
1.3 矩阵求解采用LU分解算法对复数稠密矩阵进行求解, LU分解过程包含panel列分解、panel行更新和trailing更新这3个部分[8-9].假设矩阵大小为$N\times N$、分块矩阵的大小为$n_b \times n_b $、进程网格选择为$P_r \times P_c $, 计算机发送一次消息的通信延迟为$\alpha $, 发送一个矩阵元素的时间为$\beta $, 则一次性发送$L$个矩阵元素所需要的通信时间$(T)$为
| $ \begin{align} T=\alpha +\beta L. \end{align} $ | (3) |
以双精度复数矩阵LU分解为例, 参考文献[8], 得出panel列分解、panel行更新、trailing更新的通信时间分别为
| $ T_{\rm panel, Column, comm} =2\alpha n_b \log _2 P_r +2\beta n_b ^2\log _2 P_r, $ | (4a) |
| $ T_{\rm panel, Row, comm} =( {P_c -1} )\Big( {\alpha +\frac{n_{_b }^2 }{2}\beta} \Big), $ | (4b) |
| $ T_{\rm trailing, comm} =( {P_c -1})\Big( {\alpha +\frac{(N-kn_b )n_b }{P_r }\beta }\Big)+( {\log _2 P_r -1} )\Big( {\alpha +\frac{(N-kn_b )n_b }{P_c }\beta }\Big). $ | (4c) |
由式(4a)-式(4c)可以看出, 进程网格$P_r \times P_c $、分块矩阵的大小$n_b \times n_b $、矩阵维数$N$决定了矩阵求解过程的MPI通信开销.如果能够降低LU分解过程的通信时间, 则会大幅提升矩阵求解的效率[10].
2 计算平台简介本文使用的KNL众核处理器编号为7210(简称KNL7210), 采用Slivermont的微架构和光刻14nm工艺, 处理器基本频率为1.3GHz, 拥有64核, 每核支持四线程, 即总计256个线程.每个处理器有32个核片, 每个核片由两个核组成, 并共享1MB二级缓存, 每个核心内搭载两个512位宽的VPU(Vector Processing Units), 可以同时处理8个双精度浮点或者16个单精度浮点运算, 超宽的矢量宽度, 提高了高度并行计算性能的标准.
单个KNL7210能提供超过3TFlops的双精度浮点处理性能或大于6TFlops的单精度浮点处理性能.最大内存达到384GB, 内存类型为DDR4-2133, 最大内存带宽为102GB/s, 同时集成16GB高带宽内存-MCDRAM(Multi-Channel DRAM), 可为内存访问密集型的工作负载提供高达500GB/s的可持续高内存带宽[11].
3 优化策略分析 3.1 向量化优化分析为了充分发挥KNL7210超宽的矢量宽度优势, 进一步提高循环结构的执行效率, 对公式(2)中的高斯数值积分过程进行向量化优化[9, 12], 加速矩阵填充过程. 图 3给出了向量化过程的伪代码.
|
图 3 矩阵填充向量化优化伪代码 Fig.3 Pseudo-code of vector optimization for the filling matrix |
并行矩量法向量化优化前后, 利用单个KNL7210, 计算未知量为58652的飞机模型的散射特性, 性能的提升情况如表 1所示.
|
表 1 矩阵填充向量化优化结果 Tab.1 Results of vector optimization for the filling matrix |
由表 1可以得出, 向量化后矩阵填充效率提升11.09%.所以, 在KNL众核处理器平台进行向量化优化对提升程序效率至关重要.
3.2 MPI+OpenMP混合编程根据前文分析的矩阵填充算法以及存在的问题, 综合考虑实现数据的局部性、提高Cache命中率[13-15]、降低冗余积分计算量以及充分利用KNL7210的计算资源等条件, 在MPI进程内部利用OpenMP共享内存的编程方式开启超线程, 获取最优性能[16].
图 4给出了矩阵填充过程优化后的伪代码.
|
图 4 矩阵填充MPI+OpenMP优化伪代码 Fig.4 Pseudo-code of optimization using MPI+OpenMP for the filling matrix |
对于计算密集型、数据密集型的矩阵求解部分, CPU利用率高, 求解速度与处理器的物理核数有关[17].在不改变LU分解算法的基础上, 降低分解过程的通信时间, 将会大幅提升矩量法求解电磁问题的效率.分析式(4a)-(4c)可以得出panel列分解和panel行更新过程中通信时间和进程网格$P_r \times P_c $成正比的结论. trailing更新过程中同时存在行向和列向通信, 为方便讨论, 将根据式(4a)-式(4c)分析trailing更新过程的通信时间换算为分析如下两式.
| $ {\rm coefficient1}=\frac{P_c -1}{P_r }+\frac{\log _2 P_r -1}{P_c }, $ | (5a) |
| $ {\rm coefficient2}=P_c +\log _2 P_r -2. $ | (5b) |
由式(5a)和式(5b)可以看出, coefficient1的第二项${({\log _2 P_r -1})}/{P_c }$值小于1, 考虑第一项${(P_c -1)}/{P_r }$和$P_r \le P_c $条件[8], 得出当$P_r $和$P_c $较小时, coefficient1的值较小; 同样, coefficient2的值也是随着$P_r $和$P_c $的减小而减小.通过分析可以得出, 矩阵求解过程的通信时间随着$P_r $和$P_c $的变小而变短, 所以求解过程可通过引入线程来减少MPI进程数, 达到减少LU分解通信开销的目的.
4 数值算例首先在单个KNL7210上分析飞机模型Ⅰ的散射特性, 与商业软件FEKO的计算结果对比, 验证程序优化的正确性[18]; 其次, 在KNL7210集群(8个KNL7210组成)上分析未知量超过20万的飞机模型Ⅱ的散射特性, 研究程序优化后在集群上性能的提升.
4.1 飞机模型Ⅰ的散射特性飞机Ⅰ电磁仿真模型如图 5(a)所示, 入射平面波频率为500MHz(沿机头方向), 极化方向为垂直极化, 计算其双站RCS(Radar-Cross Section)[19].该模型被剖分为39102个三角形, 公共边个数为58652, 故阻抗矩阵大小为58652$\times $58652.计算得到飞机的3D双站RCS结果如图 5(b)所示. 图 6给出了飞机的$xoy$面和$xoz$面的2D双站RCS结果, 可见程序优化后的计算结果与商业软件FEKO的仿真结果吻合, 进而证明了程序在KNL7210众核处理器上优化的正确性.
|
图 5 飞机Ⅰ仿真模型及3D双站RCS结果 Fig.5 Model of aircraftⅠand 3D bistatic RCS |
|
图 6 飞机Ⅰ的2D双站RCS结果 Fig.6 2D bistatic RCS of aircraft Ⅰ |
在该算例中, 表 2给出了并行矩量法优化前后其矩阵填充过程中的总积分次数和进程间产生的冗余积分情况.由表 2的数据可得, 程序优化前冗余比例为70.18%, 明显高于优化后冗余积分所占比例, 这会大大影响矩阵填充的速度. 表 3列出了并行矩量法优化前后各部分的计算时间情况.由表 3可以看出, 采用4个MPI进程、矩阵填充和矩阵求解分别开启64和16个OpenMP线程时, 仿真该模型的效率最佳, 比优化前加速了2.62倍, 其中矩阵填充加速了5.81倍.由于该模型的计算规模较小, 矩阵求解过程的通信时间较短, 所以, 仿真该模型时矩阵求解过程的加速效果不明显.
|
|
表 2 矩阵填充过程积分次数测试 Tab.2 Integral numbers of the filling matrix |
|
|
表 3 运行时间对比 Tab.3 Comparison of running time |
飞机模型Ⅱ的仿真模型如图 7(a)所示, 平面入射波频率为450MHz(沿机头方向), 极化方向为垂直极化, 计算其双站RCS.该模型被剖分为151452个三角形, 公共边个数为227178, 故阻抗矩阵大小为227178$\times $227178.计算得到飞机的3D双站RCS如图 7(b)所示.图 8给出了飞机的$xoy$面和$xoz$面的2D双站RCS结果.
|
图 7 飞机Ⅱ仿真模型及3D双站RCS结果 Fig.7 Model of aircraftⅡand 3D bistatic RCS |
|
图 8 飞机Ⅱ的2D双站RCS结果 Fig.8 2D bistatic RCS of aircraftⅡ |
表 4给出了该算例矩阵填充过程中的总积分次数和进程间产生的冗余积分情况.由表 4的数据可得, 程序优化前冗余比例为73.77%, 均高于优化后冗余积分所占的比例.矩阵求解部分开启OpenMP线程后, 减少了MPI进程数; 根据公式(7)得出减少了通信时间, 加速了矩阵求解过程. 表 5列出了并行矩量法优化前后各部分的计算时间情况, 其中矩阵填充和矩阵求解每个节点使用4个MPI进程、分别开启64和16个OpenMP线程, 可以获得最优性能.该算例中进程和线程的最优分配策略和上例中的计算资源分配策略一致.所以, 本文中的计算资源分配策略对今后矩量法在KNL系列上的优化具有重要的借鉴和指导意义.由表 5可得, 优化后矩阵填充速度加速13.27倍, 矩阵求解速度加速1.26倍, 总的仿真计算速度加速了3.62倍.程序优化后在KNL7210集群上的运行性能得到了明显提升.
|
|
表 4 矩阵填充过程积分次数测试 Tab.4 Integral numbers of the filling matrix |
|
|
表 5 运行时间对比 Tab.5 Comparison of running time |
本文基于KNL7210众核处理器计算平台, 对并行矩量法进行了优化, 通过引入OpenMP共享内存的并行编程策略, 显著减少了矩阵填充过程的冗余积分计算量, 降低了矩阵求解过程的通信时间, 向量化优化进一步提高了高斯数值积分过程中循环结构的执行效率.通过对飞机模型Ⅰ和飞机模型Ⅱ散射特性的仿真分析, 验证了优化算法的正确性, 并分析了并行矩量法优化后其矩阵填充和求解的加速情况.测试结果表明, 基于KNL7210众核处理器平台, 优化后的并行矩量法各部分计算均有较好的性能提升, 对后续在KNL系列平台开展矩量法研究和优化, 解决更多具有实际意义的电磁问题奠定了基础.
| [1] |
CHEN Y, ZHANG G, LIN Z, et al. Solution of EM problems using hybrid parallel MIC/CPU implementation of higher-order MoM[C]//IEEE, International Symposium on Microwave, Antenna, Propagation, and Emc Technologies. IEEE, 2016: 789-791.
|
| [2] |
张光辉. CPU/MIC异构平台中矩量法与时域有限差分法的研究[D].西安: 西安电子科技大学, 2015.
|
| [3] |
左胜, 陈岩, 张玉, 等. 一种可扩展异构并行核外高阶矩量法[J]. 西安电子科技大学学报(自然科学版), 2017, 44(1): 146-151. DOI:10.3969/j.issn.1001-2400.2017.01.026 |
| [4] |
赖明澈.数据并行协处理器体系结构的研究与实现[D].长沙: 国防科学技术大学, 2005.
|
| [5] |
HARRINGTON R F, HARRINGTON J L. Field Computation by Moment Methods[M]. NewYork: Oxford University Press, 1996.
|
| [6] |
ZHANG Y, SARKAR T K. Parallel Solution of Integral Equation Based EM Problems in the Frequency Domain[M]. Hoboken, NJ: Wiley-IEEE Press, 2009.
|
| [7] |
RAO S M, WILTON D R, GLISSON A W. Electromagnetic scattering by surfaces of arbitrary shape[J]. IEEE Transactions on Antennas & Propagation, 1982, 30(3): 409-418. |
| [8] |
张玉, 赵勋旺, 陈岩, 等. 计算电磁学中的超大规模并行矩量法[M]. 西安: 西安电子科技大学出版社, 2016.
|
| [9] |
RANA V S, LIN M, CHAPMAN B. A scalable task parallelism approach for LU decomposition with multicore CPUs[C]//Proceedings of the 2nd Internationsl Workshop on Extreme Scale Programming Models and Middleware. Piscataway, NJ, USA: IEEE Press, 2016: 17-23.
|
| [10] |
ZHANG G, CHEN Y, ZHANG Y, et al. MIC accelerated LU decomposition for method of moments[C]//IEEE International Symposium on Antennas and Propagation & Usnc/ursi National Radio Science Meeting. IEEE, 2015: 756-757.
|
| [11] |
JEFFERS J, REINDERS J. Intel Xeon Phi协处理器高性能编程指南[M].陈健, 李慧, 杨昆, 等, 译.北京: 人民邮电出版社, 2014.
|
| [12] |
高伟, 赵荣彩, 韩林, 等. SIMD自动向量化编译优化概述[J]. 软件学报, 2015, 26(6): 1265-1284. |
| [13] |
周领良, 朱延超, 刘轶, 等.基于Cache命中率校准的并行程序性能预测[C]//2014全国高性能计算学术年会论文集.中国计算机学会, 2015: 814-817.
|
| [14] |
艾维丽. 浅析Cache命中率与块的大小之间的关系[J]. 价值工程, 2011, 32: 153. |
| [15] |
叶凝, 应忍冬, 朱新忠, 等. 众核处理器系统可靠性优化方案[J]. 计算机与现代化, 2013, 218(10): 143-148. DOI:10.3969/j.issn.1006-2475.2013.10.036 |
| [16] |
MIWA M, NAKASHIMA K. Progression of MPI Non-blocking Collective Operations Using Hyper-Threading[C]//201523rd Euromicro International Conference on Parallel, Distributed and Network-Based Processing (PDP). IEEE, 2015: 163-171.
|
| [17] |
QUN N H, KHALIB Z I A, WARIP M N, et al. Hyper-threading technology: Not a good choice for speeding up CPU-bound code[C]//International Conference on Electronic Design. IEEE, 2017: 578-581.
|
| [18] |
RAJESH N, MALATHI K, RAJU S, et al. Design of vivaldi antenna with wideband radar cross section reduction[J]. IEEE Transactions on Antennas and Propagation, 2017, 65(4): 2102-2105. DOI:10.1109/TAP.2017.2670566 |
| [19] |
HU C F, LI N J, CHEN W J, et al. High-precision RCS measurement of aircraft's weak scattering source[J]. Chinese Journal of Aeronautics, 2016, 29(3): 772-778. DOI:10.1016/j.cja.2016.03.003 |