快速多极子方法的并行技术ppt课件

1、1快速多极子方法的并行技术快速多极子方法的并行技术国家国家973工程工程高性能科学计算研讨高性能科学计算研讨大规模并行计算研讨大规模并行计算研讨2纲要 FMM Data Structures Parallelization3纲要 FMM Data Structures Parallelization4FMM in Computational Electromagnetics EFIE MFIE CFIE Green函数函数2( )( ) , , ( )/()( );iSLSLjkkgdS tJt E rrJIr,r J r| - |( , )/4| - |jk r rg r rer r( )/

2、2( )( ) ,( )( )() ;iSKKgdS t J rt nJt n H rJJ rr,r(1)CFIEEFIEMFIE 5积分方程的离散积分方程的离散 RWG矢量基函数矢量基函数 MOM 离散离散1( )( )NiiijJ rf r+-( )/2, ( ), ;( )( )/2, ( ), .iiiiiiiiiiiilAifTlAifTirrrrrf rrrrrr( )( )(1) ( )/2( ),( )( )(1)( ).mnmnnnSiimmSzLKdSvdS frff rnffrE rn H rn=1 , 1,2,.NmnnmzjvmNRao-Wilson-Glisson6

3、Method of Moments Surface is Discretized into Patches (Basis Functions) Basis Functions Interact through the Greens Function,Gr r Generates a Dense Method of Moments Matrix11n n nnZIV , Zr rrjiSj ifGfds( )rf( )rtPulse7线性系统： Mx=s M是NN矩阵，x、s是N矢量Direct solution (Gauss elimination,LU Decomposition,SVD,)

4、 空间复杂度为O(N2) ，需求O(N3)次运算；Iterative methods,空间复杂度仍为O(N2)，假设K(k largest N =32,768Finding:快速矩阵乘向量的算法(NlogN)； 并行实施。8Fast Multipole Methods(FMM)Introduced by Rokhlin &Greengard in 1987Called one of the 10 most significant advances in computing of 20th centurySpeeds up matix-vector products (sums) of

5、a particular type以上求和要求O(MN)运算复杂度对给定的精度，FMM可以获得O(M+N)运算复杂度可以加速matix-vector products ，使O(N2)变为O(NlogN)加速线性系统求解，假设用迭代方法，k步收敛，每步用矩阵矢量相乘，使计算复杂度由O(N3)变为O(kNlogN)1( )(), Njijijjijis xxxs 9FMM:Application Molecular and Stellar dynamics computation of force fields and dynamics Solution of acoustical scatter

6、ing problems Helmholtz Equation Electromagnetic Wave Scattering Maxwells Equations Fluid Mechanics:Potential flow,vertex flow Laplace/Poisson Equations10FMM: Fundament 格林函数的加法定理 jlpl平面波展开|(2)1e( 1) (21) ()()()|jkllllljklj kd hkr P r+dd rr+djl_第一类球面Bessel函数hl(2) 第二类球面Hankel函数Pl Legendre多项式24()()()e()

7、njknnnSajjkd PPdk dd rkrk留意到lz时，函数jl(z)衰减非常快，而hl(2)(z)递增非常快。当dr时，上式在保证精度的情况下截断。那么上式可以写为：|(2)1e( 1) (21) ()()()|jkLllllljklj kd hkr P r+dd rr+dL=kd+cln(kd+ )Kd-源点到察看点的最大半径c-是一个依赖希望精度的常数=1 最小的相对误差小于0.1=5 相对误差小于10-6=10 准确到双精度11Fast Multipole Basics直接求解：分解：复杂度：O(MN)复杂度：O(pN+pM)cm12 FMM方式的矩阵向量乘积方式的矩阵向量乘积

8、1()221()()16NijjijjnmGn mj GmKpipppjjpmGn mj Gmz jz jkwTj DA 2( )( , ), ln(), 2;e()( ) ;mpmmpLmmjkpppjjjT krLkdckdKLdSk rrTkrAI-k kf r ()e ()( )(1)( )mjkpppipiiidSk r rDI-k kf rkf rn近场部分远场部分电磁场积分方程的多极子展开13FMM NearyyD TAx1 xA x21 xTx 聚集过程，将基函数聚集成平面波函数，其结果表示K个平面波12 yD x转移过程，将得到的x1 平移到另一个盒子的中心，其结果表示该盒子

9、中心的K个平面波发散过程，将得到的x2发散到盒子中的基函数。M2M , M2L , L2L： 聚集 转移 发散 M: 多极子展开 L: 部分展开14FMM123( )1( )( )2( )( )3( )( )( ,),( )( ,),( )( ,)iiinixE nnixEnnixEnuuuiiiyy xyy xyy x由于E2(n)和E3(n)是互补的，因此对恣意的n,23( )( )23( )( )( )( ,)( )( )innixEnEnuiyy xyy15FMM AlgorithmStep1 M2M16Step2 M2L17Step3 L2LSummation18MLFMM Algo

10、rithm Upward Pass: merge scattering matrices Downward Pass: construct splitting and exchange matrices Final SummationBased On: Hierarchical Grouping Data Structure19L2LM2MM2MM2LM2LM2LMultilevel Multipole OperatorsFinestLevelFinest - 1LevelL2LUp TreeAcross TreeDown Tree20MLFMM AlgorithmUpward Pass St

11、ep1: 在最细的层盒子求解远场展开系数，xiE1(n,l),得到C(n,l)或 ，这也可以用于xiE3(n,l) Step2:对于l=L-1,2,从step1得到值进展递归得到。同样适宜xiE3(n,l) 结果：得到分层组聚集系数( , )1( )nly21MLFMM AlgorithmDownward Pass Step1:l=2,.L递归进展E4 Step2:恣意一个非空组本身加上其邻接组最多可有3d个，其父组及父组的邻接组最多可构成3d2d，因此次相邻数目最多为p=3d(2d-1)；三维是189，二维是27，一维是3。部分到部分的转移，E3(n,l)和E4(n,l+1)产生E3(n,l

12、+1)结果：可以得到各分层组的转移系数22 Key Words 空间多层组划分 Morton编号 相邻组的作用 远场组的上聚 次相邻组中心的转移 远场组的下推23Grouping 每个盒子在第l层的目的号数为n,那么恣意盒子的目的为(n,l); 需求构建的函数，如Parent(n), ChildrenAll(n), Children(X;n,l), NeighborsAll(n,l), Neighbors(X;n,l)24Grouping每个盒子在l=0,.,L层的目的n=Number=0,1,2ld-1.12324( , )( , )( , )( , ),( , ), )( , )0,.,2

13、1( , )( , )( ),1),1)( , ),( , ), )ldI n ln lIn ln lNeighbor n l lI n lIn lIn lChildren Neighbor Parent n lln lNeighbor n l l12312422( , )( , )( , )( , )( , ), )( , )(0,0)( , )( , )( , ), )( , )( ),1)( , )( ),1)ddE nlRnlE nlRnlNeighbor nl lE nlEE nlnlNeighbor nl lE nlE Parent n lE nlParent n lNeighb表

14、 示 盒 子的 空 间 区 域 ;表 示 盒 子与 其 相 邻的 空 间 区 域 ;表 示 盒 子及 其 相 邻的 外 部 空 间 区 域 ；表 示 父 盒 子及 其 相 邻 ( ),1),1 orParent n ll(的 空 间 区 域 ；1234( , )( , )12( , )( , )( , )( , )34( , )( , )( )( ,),( )( ,),( )( ,),( )( ,)iiiin ln liixEn lxEn ln ln liixEn lxEn luuuuiiiiyy xyy xyy xyy x由于E2(n,l)和E3(n,l)是互补的，因此对恣意的n,l( ,

15、)( , )23( )( )( )n ln l yyy25纲要 FMM Data Structures Parallelization26Data Structure 构造树 紧缩八叉树 建立衔接 morton键 排序27构造树 离散点的空间层次划分离散点的空间层次划分 对应的四叉树及其紧缩四叉树对应的四叉树及其紧缩四叉树 28 确定层数L 根据读入模型的最大几何尺寸确定计算区域D,根据问题的参数确定最小盒子 尺寸d ,树构造的层数为L=log2(D/d) 第l-1层立方体等分为八个子立方体,构成第l层更小的立方体,l-1是l层的父层,l层是l-1层的子层. 构成相邻组、次相邻组、远场组 第l

16、层的节点数不超越2dl个 构造树八叉树(1)29构造树八叉树2procedure Octree_Build Octree = empty For i = 1 to n . 对一切的点做循环Octree_Insert(i, root) . 将点i插入八叉树相应的位置 end For . 八叉树中能够有很多空的叶节点, 但它们的兄弟节点非空Traverse the tree (via, say, breadth first search) . 宽度周游Eliminating empty leaves . 去掉空的叶节点Compress Octree . 紧缩八叉树 . 假设某中间节点仅包含一个子节

17、点，那么被其交换每个节点用两个键值描画: 包含一样数据的最小单元和最大单元30构造树构造树八叉树八叉树3包含N个叶节点的紧缩八叉树节点总数不超越2N-1因此可以采用数组而不是链表来存储紧缩八叉树MLFMM基于后序周游的紧缩八叉树数据构造从键值最小的叶节点开场周游每个节点都存储在其子节点之后，且紧挨其子节点存储节点排序时，运用的是其所表示的最小单元键值对于兄弟节点，按键值从小到大排序各节点所表示的单元指的是它所包含的最小单元后序周游存储方式是实现与分布无关的自动负载平衡并行MLFMM的关键分布无关自顺应树(Distribution-Independent Adaptive Tree, DIAT)

18、构造2d维DIAT的复杂度为O(dnlogn)树树节节点点存存储储31Morton键为什么不用指针？ 用指针指向Children的指针可以很方便的建立父子节点之间的关系，从而构成树构造的拓扑构造。在串行程序，指针可以在全局存储空间中寻址，效率很高也很准确。但在内存分布式并行环境中，一个计算节点不能直接访问另一个计算节点上的存储空间，因此用于联络树构造拓扑构造的指针只能在其所在的计算节点上才有意义，假设要让指针所指向的树节点可以存储在其他节点上，就必需小心处置指针的变换关系。以便将指针的指向正确地映射到其他计算节点上的内存空间。这种转换，使得基于指针的树拓扑方法在分布式并行环境中实现起来不仅很复

19、杂，而且效率也不高。Morton键技术是实现并行多层快速多极子的关键技术之一！ 原理：将空间坐标的二进制序列按位交叉，把空间中某一点在x、y、z方向的坐标/序号映射为一个值，这个值就是morton键值。32Morton次序次序位于第m层,在该层中编号为n的盒子,普通采用Morton次序编号为(n, m).左图为第三层的Morton次序 ，右图为每一层编号,前三层分别有1, 4, 16个点. 33Morton编码小的灰盒子在第3层,本层编号为11, 于是其Morton编号为(11,3); (023)4=(11)10采用四进制编号为(0,2,3); Morton编号(Num,l)，在2d叉树中可以

20、得到某盒子对应的 2d进制编号:(N1, N2,Nl)2d再按照下面的公式计算其Morton编号：111012112()(.)(2 ). (2 ).dd ldlllNumN NNNumNNN算法34空间编码 尺度转换 二进制编码 d维空间编码 二进制交错及其解交错 涉及到的算法：查找Parent、Children、Neighbor、盒子中心坐标、盒子编码等35尺度转换 2d树构造每个参数都有本人的参数Dj映射原始的盒子的大小x1,min, x1,max x1,min, x1,max ,xj,max-xj,min=Dj,j=1,d 把原始的盒子映射到单位盒子 上 实践的三维物理空间Dj=D， x

21、min=(x1,min,，xd,min), 称x为真实的坐标， 为该点规范化的坐标 目的：实施二进制编码转换，方便查找！0,1 . 0,1x ,min1,1,. ,(,.)jjdjxxxjdxxDxminxxxDx36二进制编码For example d=1: 用十进制表示为 对于 不仅可以表示为 也可以表示为 用二进制表示 为： 对于 0,1x 12310(0.) ,0,.,9,1,2,.jxa a aaj1x 1.000.x 101(0.999999.)x 0,1x 1 2 32(0.) ,0,1;1,2,.jxbbbbj1x 21(0.11111.)x 37位交错位交错 在d维单位盒子里

22、, 点坐标 其中第k个分量 也可以写成交错二进制(bit interleaving)的方式:1(,.,),dxxx122(0.) ,0,1;1,2,.,1,., .kkkkjxb bbjkd11 21112222122(0.|.|.|.|.) , 0,1ddjjdjkjxb bbb bbb bbb121222(0.) ,(.) ,0,.,21ddjjjjdjjxN NNNb bbN38解交错解交错 d维盒子在第l层的Morton键值为 可以解交错(bit deinterleaving)为d个数值代表该盒子的d维坐标11 21112222122(.|.|.|.)ddlldlNumb bbb bb

23、b bbNumber=7689310=1112=710=1110102=5810=10012=910(9, 58, 7)11112122212222122(.) ;(.) ;.(.) ;llddddlNumb bbNumb bbNumb bbNumber3=Number2=Number1=39寻觅给定点所在的盒子寻觅给定点所在的盒子在d维单位盒子里, 点坐标的交错2d进制方式:可以利用2d进制移位取整寻觅给定点在第l层所在的盒子:或者写成dl位的平移方式:11 21112222122(.|.|.|.) , 0,1ddljdlkjNumb bbb bbb bbb121 21222(0.)(,.,

24、) , )ddllxN NN c cbox N NNl121 22(, )2(. . .) ddllNum lxN NNbb查找分为5层八叉树构造的点(0.7681, 0.0459,0.3912)在第3层盒子的号(1)二进制转化(0.11000,0.00001,0.01100)2(2构成单个串0.1001010010000102(33*3=9 (Number,3)=1001010012=297 3*5=15 (Number,5)=1001010010000102=1901040查找盒子中心查找盒子中心单位盒子第l层编号为Num的子盒子中心的二进制d维坐标方式: 或者不依赖计数体系,写为: 例如

25、: 八叉树第5层10进制编号为533的盒子，计算过程为:(1)将Morton编号转为2进制: 53310=1, 000, 010, 1012(2)经过解交错得到该盒子的三维坐标:Num3=10012=910 , Num2=0102=210 , Num1=0012=110(3)计算盒子中心坐标:x3c(533, 5)=2-5(9+0.5)=0.296875; x2c(533, 5)=2-5(2+0.5)=0.078125; x1c(533, 5)=2-5(1+0.5)=0.04875; 因此xc=(0.04875, 0.078125, 0.296875)12(, )(0.1), 1,.,kckk

26、klxNum lb bbkd(, )2(1/2), 1,.,lkckxNum lNumkd算法41父盒子编码父盒子编码12(,.,)lP a ren t NNN 计算某盒子的父盒子 父盒子编码与层次无关, 其计算方法: 除法的商取整, 比如11/4=2，于是 例如： 于是#5981的父盒子是#747 除以2d相当于作d次2进制位移就可以了1111()(2 ).(2 ).dldllParent NumNNN(, )(),1),ParentNum lParent Num l()/2 dParent NumNum104410(11)(23)(2)(2).ParentParent3103,(5981)

27、5981/2 747.dParent算法42子盒子编码子盒子编码12(,.,)lChildren N NN 计算某盒子的子盒子 子盒子是个集合,与所在的层无关: 其计算方法为: For example121211(,.,)(,.,),0,.,21.dllllChildren N NNN NN NN(, )(),1),Children Num lChildren Num l()2,0,1,.,21.ddChildren NumNumjj104444410101010(11 )(23 )230 ,231 ,232 ,233 44 ,45 ,46 ,47 ChildrenChildren1011 2

28、111()(.) ()(2 ).(2 ),lldldllNumNN NChildren NumNNN43查找邻居盒子查找邻居盒子1 基于交错二进制的邻居查询算法基于交错二进制的邻居查询算法 Step1.解交错解交错(deinterleaving).十进制十进制(或二进制或二进制)编号可转为编号可转为d维坐标维坐标: Step2.每个分坐标的邻居每个分坐标的邻居: 于是邻居集为于是邻居集为: 其中其中nk定义为定义为1,.,dNumNumNum1,1,1, .,.kkkkN u mN u mN u mN u mkd1()( ,.,)dNeighbor Numnn, 0,21, 0;, 21.dk

29、kkkkkkkkdkkkknNumNumNumif NumnNumNumif NumnNumNumif Num44比如编号为#26的2维盒子,其邻居查询过程如下: 2610=(01 10 10)2解交错方式为(011,100) 2=(3,4) 10其相邻盒子为(2,3),(2,4), (2,5) ; (3,3),(3,5); (4,3),(4,4),(4,5); 8个点的二进制:(010,011) 2 ,(010,100) 2 , (100,101) 2相邻盒子编号的交错 (interleaving)二进制方式: 001101, 011000, 011001, 001111, 011011,

30、100101, 110000, 110001十进制编号为:13,24,25, 15,27, 37,48,49.查找邻居盒子查找邻居盒子245纲要 FMM Data Structures Parallelization46并行实施技术并行实施技术(parallel implementation) MLFMM具有很好的并行效率 大多数操作都在本地机上进展 例如，Parant to Children,box to its interaction list 八叉树构造是很大的 需求生成和分布式存储 每个处置器只坚持本地的根本树 大多数任务只在本地的根本树上进展数据需求同步 父节点需求子节点上的值，但这

31、两个节点在不同的处置器上荷载平衡问题但是，还存在： 分布式八叉树 负载平衡 相互作用表列 相邻结点的通讯 次相邻点的通讯需求处理47并行计算步骤 构造紧缩八叉树 近场矩阵计算 建立转移节点列表 远场矩阵计算 聚合 转移 发散48 树构造的并行划分149二维计算区域对应的分布式四叉树二维计算区域对应的分布式四叉树 树构造的并行划分250构造分布式紧缩八叉树构造分布式紧缩八叉树(1) 层数L=log2(D/d), D和d为计算区域划分的最大和最小盒子尺寸 将n个基函数在p个处置机上按次序平均分配，再按照基函数的位置生成包含这些基函数的叶节点 由于不同基函数可包含于同一叶节点，因此这样的叶节点会同时

32、存储在不同处置机上 RWG基函数定义在各边上，并包含一对完好的三角形P1P2P3A1A2A3用边的中点代表整个边,各点都有相应的最底层非空盒子, 即叶节点.将叶节点的Morton键值赋给中点所在的边51构造分布式紧缩八叉树构造分布式紧缩八叉树(2)采用并行排序算法对一切处置机中的基函数和叶节点排序，使个处置机包含一样数量的基函数对每个处置机里的N/p个键值，采用快速排序(quicksort)全局并行排序采用取样排序(samplesort), 它用到位排序(bitornic sort)排序时用到的通讯为MPI_Allgather和MPI_Alltoall每个处置机还包含下一处置机的第一个叶节点,

33、 并根据这些叶节点建立本地紧缩八叉树，经过后序周游存储各处置机将本地紧缩树中位于从下一处置机得到的叶节点之后的节点发送到下一处置机，并按后序周游插到对应的位置共享叶节点个数不超越L, 每个处置机接纳的非本地节点不超越7L各处置机从下而上构造本地树的复杂度为O( (N/p)log(N/p)52近场计算近场计算Morton次序保证兄弟节点的键值是相邻的，但并不是只需兄弟节点相邻, 因此需求调用位交错和解交错函数去查找邻居节点近场矩阵Anear只需求思索最细层(第L层)每个节点及其相邻节点所包含的基函数相互作用; 必需按照MOM计算, 并在迭代前存储每个叶节点最多有26个相邻节点。假设最细层每个盒子

34、至多包含c个基函数, 那么每个叶节点上的计算量为27c2假好像一棵子树上的相邻节点位于同一处置机, 那么无需通讯 假设相邻节点位于不同的处置机上, 那么运用MPI_Allgather获得每一处置机的第一个和最后一个叶节点的键值.,并存储在长为2p的数组中, 经过该数组的检索得到相邻叶节点调用MPI_Alltoall实现数据(叶结点及其包含的基函数)的分发; 整个近场计算复杂度为O( (N/p)log(N/p)53远场计算远场计算 部分的不变项(如最细层的D和A)只需求分配到它对应的子树所在处置机上，全局不变项(如常数)那么分配到一切处置机上；由于下一层的计算依赖于上一层的信息，因此完成每一层的计算时，各个处置机都需求进展一次同步 存储时采用内存循环战略，它依赖于数据的相关性。聚集项S在层层上聚时分配内存，每当层层下推时某层的发散项B计算终了，就将该层的S的内存释放掉；发散项B在层层下推时分配内存, 每当处置机上某层的一切的B计算终了，就将其父层的B释放掉 归并各处置机得到的结果, 即为远场矩阵向量乘 经过归约得到整个系数矩阵与向量的乘积 经过并行的迭代法得到计算结果(等效电流), 每次迭代的矩阵向量乘法计算完成时，需求进展一次同步 完成计算结果的后处置，比如RCS的计算。54树构造代码树构造代码 上聚上聚 A M2M 内插值内插值 下推下推 C M2