《快速多极子方法的并行技术.ppt》由会员分享，可在线阅读，更多相关《快速多极子方法的并行技术.ppt（60页珍藏版）》请在三一办公上搜索。
1、快速多极子方法的并行技术,冯仰德 王武 迟学斌中科院计算机网络信息中心 超级计算中心 2007年2月5日,国家973项目高性能科学计算研究大规模并行计算研究,纲要,FMMData StructuresParallelization,纲要,FMMData StructuresParallelization,FMM in Computational Electromagnetics,EFIEMFIECFIEGreen函数,积分方程的离散,RWG矢量基函数MOM 离散,Rao-Wilson-Glisson,Method of Moments,Surface is Discretized into P
2、atches(Basis Functions),Basis Functions Interact through the Greens Function,Generates a Dense Method of Moments Matrix,线性系统：,Mx=s M是NN矩阵，x、s是N矢量Direct solution(Gauss elimination,LU Decomposition,SVD,)空间复杂度为O(N2)，需要O(N3)次运算；Iterative methods,空间复杂度仍为O(N2)，如果K(k largest N=32,768Finding:快速矩阵乘向量的算法(Nlog
3、N)；并行实施。,Fast 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 a particular type以上求和要求O(MN)运算复杂度对给定的精度，FMM可以获得O(M+N)运算复杂度可以加速matix-vector products，使O(N2)变为O(NlogN)加速线性系统求解，
4、如果用迭代方法，k步收敛，每步用矩阵矢量相乘，使计算复杂度由O(N3)变为O(kNlogN),FMM:Application,Molecular and Stellar dynamics computation of force fields and dynamicsSolution of acoustical scattering problems Helmholtz EquationElectromagnetic Wave Scattering Maxwells EquationsFluid Mechanics:Potential flow,vertex flow Laplace/Pois
5、son Equations,FMM:Fundament,格林函数的加法定理jlpl平面波展开,jl_第一类球面Bessel函数hl(2)第二类球面Hankel函数Pl Legendre多项式,注意到lz时，函数jl(z)衰减非常快，而hl(2)(z)递增非常快。当dr时，上式在保证精度的情况下截断。则上式可以写为：,Kd-源点到观察点的最大半径c-是一个依赖希望精度的常数=1 最小的相对误差小于0.1=5 相对误差小于10-6=10 准确到双精度,Fast Multipole Basics,直接求解：,分解：,复杂度：O(MN),复杂度：O(pN+pM),cm,FMM形式的矩阵向量乘积,近场部
6、分,远场部分,电磁场积分方程的多极子展开,FMM,聚集过程，将基函数聚集成平面波函数，其结果表示K个平面波,转移过程，将得到的x1 平移到另一个盒子的中心，其结果表示该盒子中心的K个平面波,发散过程，将得到的x2发散到盒子中的基函数。,M2M,M2L,L2L：聚集 转移 发散 M:多极子展开 L:局部展开,FMM,由于E2(n)和E3(n)是互补的，因此对任意的n,FMM Algorithm,Step1 M2M,Step2 M2L,Step3 L2L,Summation,MLFMM Algorithm,Upward Pass:merge scattering matricesDownward
7、Pass:construct splitting and exchange matricesFinal Summation,Based On:Hierarchical Grouping Data Structure,L2L,M2M,M2M,M2L,M2L,M2L,Multilevel Multipole Operators,FinestLevel,Finest-1Level,L2L,Up Tree,Across Tree,Down Tree,MLFMM AlgorithmUpward Pass,Step1:在最细的层盒子求解远场展开系数，xiE1(n,l),得到C(n,l)或，这也可以用于xi
8、E3(n,l)Step2:对于l=L-1,2,从step1得到值进行递归得到。同样适合xiE3(n,l)结果：得到分层组聚集系数,MLFMM AlgorithmDownward Pass,Step1:l=2,.L递归进行E4Step2:,任意一个非空组自身加上其邻接组最多可有3d个，其父组及父组的邻接组最多可形成3d2d，因此次相邻数目最多为p=3d(2d-1)；三维是189，二维是27，一维是3。,局部到局部的转移，E3(n,l)和E4(n,l+1)产生E3(n,l+1),结果：可以得到各分层组的转移系数,Key Words,空间多层组划分Morton编号相邻组的作用远场组的上聚次相邻组中心
9、的转移远场组的下推,Grouping,每个盒子在第l层的指标号数为n,那么任意盒子的指标为(n,l);需要构建的函数，如Parent(n),ChildrenAll(n),Children(X;n,l),NeighborsAll(n,l),Neighbors(X;n,l),Grouping,每个盒子在l=0,.,L层的指标n=Number=0,1,2ld-1.,由于E2(n,l)和E3(n,l)是互补的，因此对任意的n,l,纲要,FMMData StructuresParallelization,Data Structure,构造树 压缩八叉树建立连接 morton键 排序,构造树,离散点的空间
10、层次划分,对应的四叉树及其压缩四叉树,确定层数L 根据读入模型的最大几何尺寸确定计算区域D,根据问题的参数确定最小盒子 尺寸d,树结构的层数为L=log2(D/d)第l-1层立方体等分为八个子立方体,形成第l层更小的立方体,l-1是l层的父层,l层是l-1层的子层.形成相邻组、次相邻组、远场组第l层的节点数不超过2dl个,构造树八叉树(1),构造树八叉树（2）,procedure Octree_Build Octree=empty For i=1 to n.对所有的点做循环Octree_Insert(i,root).将点i插入八叉树相应的位置 end For.八叉树中可能有很多空的叶节点,但它
11、们的兄弟节点非空Traverse the tree(via,say,breadth first search).宽度周游Eliminating empty leaves.去掉空的叶节点Compress Octree.压缩八叉树.如果某中间节点仅包含一个子节点，则被其替换每个节点用两个键值描述:包含相同数据的最小单元和最大单元,构造树八叉树（3）,包含N个叶节点的压缩八叉树节点总数不超过2N-1因此可以采用数组而不是链表来存储压缩八叉树MLFMM基于后序周游的压缩八叉树数据结构从键值最小的叶节点开始周游每个节点都存储在其子节点之后，且紧挨其子节点存储节点排序时，使用的是其所表示的最小单元键值对于
12、兄弟节点，按键值从小到大排序各节点所表示的单元指的是它所包含的最小单元后序周游存储方式是实现与分布无关的自动负载平衡并行MLFMM的关键分布无关自适应树(Distribution-Independent Adaptive Tree,DIAT)构造2d维DIAT的复杂度为O(dnlogn),树节点存储,Morton键,为什么不用指针？用指针指向Children的指针可以很方便的建立父子节点之间的关系，从而构成树结构的拓扑结构。在串行程序，指针可以在全局存储空间中寻址，效率很高也很准确。但在内存分布式并行环境中，一个计算节点不能直接访问另一个计算节点上的存储空间，因此用于联系树结构拓扑结构的指针只
13、能在其所在的计算节点上才有意义，如果要让指针所指向的树节点能够存储在其他节点上，就必须小心处理指针的变换关系。以便将指针的指向正确地映射到其他计算节点上的内存空间。这种转换，使得基于指针的树拓扑方法在分布式并行环境中实现起来不仅很复杂，而且效率也不高。Morton键技术是实现并行多层快速多极子的关键技术之一！原理：将空间坐标的二进制序列按位交叉，把空间中某一点在x、y、z方向的坐标/序号映射为一个值，这个值就是morton键值。,Morton次序,位于第m层,在该层中编号为n的盒子,一般采用Morton次序编号为(n,m).左图为第三层的Morton次序，右图为每一层编号,前三层分别有1,4,
14、16个点.,Morton编码,小的灰盒子在第3层,本层编号为11,于是其Morton编号为(11,3);(023)4=(11)10采用四进制编号为(0,2,3);Morton编号(Num,l)，在2d叉树中可以得到某盒子对应的 2d进制编号:(N1,N2,Nl)2d,再按照下面的公式计算其Morton编号：,算法,空间编码,尺度转换二进制编码d维空间编码二进制交错及其解交错涉及到的算法：查找Parent、Children、Neighbor、盒子中心坐标、盒子编码等,尺度转换,2d树结构每个参数都有自己的参数Dj映射原始的盒子的大小x1,min,x1,max x1,min,x1,max,xj,m
15、ax-xj,min=Dj,j=1,d把原始的盒子映射到单位盒子 上实际的三维物理空间Dj=D，xmin=(x1,min,，xd,min),称x为真实的坐标，为该点标准化的坐标目的：实施二进制编码转换，方便查找！,二进制编码,For example d=1:用十进制表示为 对于 不仅可以表示为 也可以表示为 用二进制表示 为：对于,位交错,在d维单位盒子里,点坐标其中第k个分量 也可以写成交错二进制(bit interleaving)的形式:,解交错,d维盒子在第l层的Morton键值为 可以解交错(bit deinterleaving)为d个数值代表该盒子的d维坐标,Number=768931
16、0=,=1112=710,=1110102=5810,=10012=910,(9,58,7),Number3=,Number2=,Number1=,寻找给定点所在的盒子,在d维单位盒子里,点坐标的交错2d进制形式:可以利用2d进制移位取整寻找给定点在第l层所在的盒子:或者写成dl位的平移形式:,查找分为5层八叉树结构的点(0.7681,0.0459,0.3912)在第3层盒子的号(1)二进制转化(0.11000,0.00001,0.01100)2(2）形成单个串2(3）3*3=9(Number,3)=1001010012=297 3*5=152=19010,查找盒子中心,单位盒子第l层编号为N
17、um的子盒子中心的二进制d维坐标形式:或者不依赖计数体系,写为:例如:八叉树第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),算法,
18、父盒子编码,计算某盒子的父盒子父盒子编码与层次无关,其计算方法:除法的商取整,比如11/4=2，于是例如：于是#5981的父盒子是#747 除以2d相当于作d次2进制位移就可以了,算法,子盒子编码,计算某盒子的子盒子子盒子是个集合,与所在的层无关:其计算方法为:For example,查找邻居盒子（1）,基于交错二进制的邻居查询算法Step1.解交错(deinterleaving).十进制(或二进制)编号可转为d维坐标:Step2.每个分坐标的邻居:于是邻居集为:其中nk定义为,比如编号为#26的2维盒子,其邻居查询过程如下:2610=(01 10 10)2解交错形式为(011,100)2=(
19、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,100101,110000,110001十进制编号为:13,24,25,15,27,37,48,49.,查找邻居盒子（2）,纲要,FMMData StructuresParallelization,并行实施技术(parallel implementation),MLFM
20、M具有很好的并行效率大多数操作都在本地机上进行例如，Parant to Children,box to its interaction list,八叉树结构是很大的 需要生成和分布式存储 每个处理器只保持本地的基本树 大多数工作只在本地的基本树上进行数据需要同步 父节点需要子节点上的值，但这两个节点在不同的处理器上荷载平衡问题,但是，还存在：,分布式八叉树负载平衡相互作用表列相邻结点的通信次相邻点的通信,需要解决,并行计算步骤,构造压缩八叉树近场矩阵计算建立转移节点列表远场矩阵计算聚合转移发散,树结构的并行划分（1）,二维计算区域对应的分布式四叉树,树结构的并行划分（2）,构造分布式压缩八叉树
21、(1),层数L=log2(D/d),D和d为计算区域划分的最大和最小盒子尺寸将n个基函数在p个处理机上按次序平均分配，再按照基函数的位置生成包含这些基函数的叶节点由于不同基函数可包含于同一叶节点，因此这样的叶节点会同时存储在不同处理机上RWG基函数定义在各边上，并包含一对完整的三角形,P1,P2,P3,A1,A2,A3,用边的中点代表整个边,各点都有相应的最底层非空盒子,即叶节点.将叶节点的Morton键值赋给中点所在的边,构造分布式压缩八叉树(2),采用并行排序算法对所有处理机中的基函数和叶节点排序，使个处理机包含相同数量的基函数对每个处理机里的N/p个键值，采用快速排序(quicksort
22、)全局并行排序采用取样排序(samplesort),它用到位排序(bitornic sort)排序时用到的通信为MPI_Allgather和MPI_Alltoall每个处理机还包含下一处理机的第一个叶节点,并根据这些叶节点建立本地压缩八叉树，通过后序周游存储各处理机将本地压缩树中位于从下一处理机得到的叶节点之后的节点发送到下一处理机，并按后序周游插到对应的位置共享叶节点个数不超过L,每个处理机接收的非本地节点不超过7L各处理机从下而上构造本地树的复杂度为O(N/p)log(N/p),近场计算,Morton次序保证兄弟节点的键值是相邻的，但并不是只有兄弟节点相邻,因此需要调用位交错和解交错函数去
23、查找邻居节点近场矩阵Anear只需要考虑最细层(第L层)每个节点及其相邻节点所包含的基函数相互作用;必须按照MOM计算,并在迭代前存储每个叶节点最多有26个相邻节点。若最细层每个盒子至多包含c个基函数,则每个叶节点上的计算量为27c2如果同一棵子树上的相邻节点位于同一处理机,则无需通信 如果相邻节点位于不同的处理机上,则使用MPI_Allgather获得每一处理机的第一个和最后一个叶节点的键值.,并存储在长为2p的数组中,通过该数组的检索得到相邻叶节点调用MPI_Alltoall实现数据(叶结点及其包含的基函数)的分发;整个近场计算复杂度为O(N/p)log(N/p),远场计算,局部的不变项(
24、如最细层的D和A)只需要分配到它对应的子树所在处理机上，全局不变项(如常数)则分配到所有处理机上；由于下一层的计算依赖于上一层的信息，因此完成每一层的计算时，各个处理机都需要进行一次同步存储时采用内存循环策略，它依赖于数据的相关性。聚集项S在层层上聚时分配内存，每当层层下推时某层的发散项B计算完毕，就将该层的S的内存释放掉；发散项B在层层下推时分配内存,每当处理机上某层的所有的B计算完毕，就将其父层的B释放掉归并各处理机得到的结果,即为远场矩阵向量乘通过归约得到整个系数矩阵与向量的乘积通过并行的迭代法得到计算结果(等效电流),每次迭代的矩阵向量乘法计算完成时，需要进行一次同步完成计算结果的后处
25、理，比如RCS的计算。,树结构代码,上聚A M2M内插值下推C M2LB L2L外插值二叉树的例子,建立相互作用表列,相互作用表列(interaction list)包含每一层、每个节点的次相邻节点次相邻节点指它们本身不相邻，但它们的父节点相邻因此每个节点最多有63-33=189个次相邻节点，远多于相邻节点（26）需要在表列中注明次相邻点是否位于其它处理机每层都要建立次相邻表列，但次相邻点可能不是物理意义的同一层,empty,M2L,相互作用表列在迭代前存储,每次远场作用的转移项都会调用该表列,聚集,对于每一层的每个盒子,聚集相当于将子层组(t)中心的平面波移置到父层组(Pt)的中心(Ct),
26、并通过内插值得到大数目的平面波:父节点(或部分子节点)在其它处理机上的节点构成剩余八叉树,它至多包含8pL个子节点,因此数据交换的量级为O(logp+logL)对压缩八叉树的所有节点(包括剩余节点)应用并行聚合算法,因此每个处理机的计算量为O(N/p+L)后序遍历保证父节点在子节点之后,且紧挨子节点存储,转移,对于每一层的盒子,其远场作用只需要考虑次相邻组中心间的作用,即转移.假设组(Pt)和(Ps)的中心分别为(Ct,Cs),则:在第一次迭代时,若相互作用表列中的次相邻组中心不在同一处理机,则通过MPI_Alltoall发送到各个相应的处理机按照上面的公式计算时,通过MPI_Alltoall
27、接收次相邻组的平面波对于压缩八叉树,次相邻组的盒子大小可能不一样,即插值取样点数可能不同,因此需要将取样点少的平面波转换为取样点多的,发散,发散过程是聚集过程的逆操作,对于每一层的每个盒子,发散相当于将父层组(Pr)中心的平面波移置到子层组(r)的中心(Ct),并通过外插值得到小数目的平面波,并结合求和得到远场作用:父节点(或部分子节点)在其它处理机上的节点构成剩余八叉树,它至多包含8pL个子节点,因此数据交换的量级为O(logp+logL)对压缩八叉树的所有节点(包括剩余节点)应用并行发散算法,因此每个处理机的计算量为O(N/p+L)由于后序遍历能保证父节点在子节点之后,且紧挨子节点存储因此采用反向遍历后序存储的树节点,计算流程,网格剖分与几何信息读入,最细层的近场信息,聚集（内插值）,近场矩阵计算,近场作用,转移(次相邻组),矩量法离散,迭代法求解,向量运算,矩阵-向量乘积,BLAS,LAPACK,构造分布式八叉树,发散（外插值）,各层的远场信息,远场作用,预条件子,线性方程组,电磁场积分方程,THANKS,
