快速多极子方法并行的技术.pptx

,#,快速多极子方法的并行技术,冯仰德 王武 迟学斌,中科院计算机网络信息中心 超级计算中心,ydfeng,2007,年,2,月,5,日,国家,973,项目,高性能科学计算研究,大规模并行计算研究,纲要,FMM,Data Structures,Parallelization,纲要,FMM,Data Structures,Parallelization,FMM,in Computational Electromagnetics,EFIE,MFIE,CFIE,Green,函数,积分方程的离散,RWG,矢量基函数,MOM,离散,Rao-Wilson-Glisson,Method of Moments,Surface is Discretized into Patches(Basis Functions),Basis Functions Interact through the Greens Function,Generates a Dense Method of Moments Matrix,Pulse,线性系统：,Mx=s,M,是,NN,矩阵，,x,、,s,是,N,矢量,Direct solution(Gauss elimination,LU Decomposition,SVD,),空间复杂度为,O(N,2,),，需要,O(N,3,),次运算；,Iterative methods,空间复杂度仍为,O(N,2,),，,如果,K(k largest N=32,768,Finding:,快速矩阵乘向量的算法,(,NlogN,),；,并行实施。,Fast Multipole Methods(FMM),Introduced by Rokhlin&Greengard in 1987,Called one of the 10 most significant advances in computing of 20,th,century,Speeds up matix-vector products(sums)of a particular type,以上求和要求,O(MN),运算,复杂度,对给定的精度，,FMM,可以获得,O(M+N),运算复杂度,可以加速,matix-vector products,，使,O(N,2,),变为,O(NlogN),加速线性系统求解，如果用迭代方法，,k,步收敛，每步用矩阵矢量相乘，使计算复杂度由,O(N,3,),变为,O(kNlogN),FMM:Application,Molecular and Stellar dynamics,computation of force fields and dynamics,Solution of acoustical scattering problems,Helmholtz Equation,Electromagnetic Wave Scattering,Maxwells Equations,Fluid Mechanics:Potential flow,vertex flow,Laplace/Poisson Equations,FMM:Fundament,格林函数的加法定理,j,l,p,l,平面波展开,j,l,_,第一类球面,Bessel,函数,h,l,(2),第二类球面,Hankel,函数,P,l,Legendre,多项式,注意到,lz,时，函数,j,l,(z),衰减非常快，而,h,l,(2),(z),递增非常快。当,dr,时，上式在保证精度的情况下截断。则上式可以写为：,Kd-,源点到观察点的最大半径,c-,是一个依赖希望精度的常数,=1,最小的相对误差小于,0.1,=5,相对误差小于,10,-6,=10,准确到双精度,Fast Multipole Basics,直接求解：,分解：,复杂度：,O(MN),复杂度：,O(pN+pM),c,m,FMM,形式的矩阵向量乘积,近场部分,远场部分,电磁场积分方程的多极子展开,FMM,聚集,过程，将基函数聚集成平面波函数，其结果表示,K,个平面波,转移,过程，将得到的,x,1,平移到另一个盒子的中心，其结果表示该盒子中心的,K,个平面波,发散,过程，将得到的,x,2,发散到盒子中的基函数。,M2M,M2L,L2L,：聚集,转移,发散,M:,多极子展开,L:,局部展开,FMM,由于,E,2,(n),和,E,3,(n),是互补的，因此对任意的,n,FMM Algorithm,Step1,M2M,Step2,M2L,Step3,L2L,Summation,MLFMM Algorithm,Upward Pass:,merge scattering matrices,Downward Pass:,construct splitting and exchange matrices,Final Summation,Based On:,Hierarchical,Grouping,Data Structure,L2L,M2M,M2M,M2L,M2L,M2L,Multilevel Multipole Operators,Finest,Level,Finest-1,Level,L2L,Up Tree,Across Tree,Down Tree,MLFMM Algorithm,Upward Pass,Step1:,在最细的层盒子求解远场展开系数，,x,i,E,1,(n,l),得到,C,(n,l),或 ，这也可以用于,x,i,E,3,(n,l),Step2:,对于,l=L-1,2,从,step1,得到值进行递归得到。同样适合,x,i,E,3,(n,l),结果,：得到分层组聚集系数,MLFMM Algorithm,Downward Pass,Step1:,l=2,.L,递归进行,E,4,Step2:,任意一个非空组自身加上其邻接组最多可有,3,d,个，其父组及父组的邻接组最多可形成,3,d,2,d,，因此次相邻数目最多为,p=3,d,(2,d,-1),；三维是,189,，二维是,27,，一维是,3,。,局部到局部的转移，,E,3,(n,l),和,E,4,(n,l+1),产生,E,3,(n,l+1),结果：可以得到各分层组的转移系数,Key Words,空间多层组划分,Morton,编号,相邻组的作用,远场组的上聚,次相邻组中心,的转移,远场组的下推,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,2,ld,-1.,由于,E,2,(n,l),和,E,3,(n,l),是互补的，因此对任意的,n,l,纲要,FMM,Data Structures,Parallelization,Data Structure,构造树,压缩八叉树,建立连接,morton,键,排序,构造树,离散点的空间层次划分,对应的四叉树及其压缩四叉树,确定层数,L,根据读入模型的最大几何尺寸确定计算区域,D,根据问题的参数确定,最小盒子 尺寸,d,树结构的层数为,L=log,2,(D/d),第,l-1,层立方体等分为八个子立方体,形成第,l,层更小的立方体,l-1,是,l,层的父层,l,层是,l-1,层的子层,.,形成相邻组、次相邻组、远场组,第,l,层的节点数不超过,2,dl,个,构造树,八叉树,(1),构造树,八叉树（,2,）,procedure 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 .,压缩八叉树,.,如果某中间节点仅包含一个子节点，则被其替换,每个节点用两个键值描述,:,包含相同数据的最小单元和最大单元,构造树,八叉树（,3,）,包含,N,个叶节点的压缩八叉树节点总数不超过,2N-1,因此可以采用数组而不是链表来存储压缩八叉树,MLFMM,基于后序周游的压缩八叉树数据结构,从键值最小的叶节点开始周游,每个节点都存储在其子节点之后，且紧挨其子节点存储,节点排序时，使用的是其所表示的最小单元键值,对于兄弟节点，按键值从小到大排序,各节点所表示的单元指的是它所包含的最小单元,后序周游存储方式是实现与分布无关的自动负载平衡并行,MLFMM,的关键,分布无关自适应树,(Distribution-Independent Adaptive Tree,DIAT),构造,2,d,维,DIAT,的复杂度为,O(dnlogn),树,节,点,存,储,Morton,键,为什么不用指针？,用指针指向,Children,的指针可以很方便的建立父子节点之间的关系，从而构成树结构的拓扑结构。在串行程序，指针可以在全局存储空间中寻址，效率很高也很准确。但在内存分布式并行环境中，一个计算节点不能直接访问另一个计算节点上的存储空间，因此用于联系树结构拓扑结构的指针只能在其所在的计算节点上才有意义，如果要让指针所指向的树节点能够存储在其他节点上，就必须小心处理指针的变换关系。以便将指针的指向正确地映射到其他计算节点上的内存空间。这种转换，使得基于指针的树拓扑方法在分布式并行环境中实现起来不仅很复杂，而且效率也不高。,Morton,键技术是实现并行多层快速多极子的关键技术之一！,原理：,将空间坐标的二进制序列按位交叉，把空间中某一点在,x,、,y,、,z,方向的坐标,/,序号映射为一个值，这个值就是,morton,键值。,Morton,次序,位于第,m,层,在该层中编号为,n,的盒子,一般采用,Morton,次序编号为,(n,m).,左图为第三层的,Morton,次序，,右图为每一层编号,前三层分别有,1,4,16,个点,.,Morton,编码,小的灰盒子在第,3,层,本层编号为,11,于是其,Morton,编号为,(11,3);,(023),4,=(11),10,采用四进制编号为,(0,2,3);,Morton,编号,(Num,l),，,在,2,d,叉树中可以得到某盒子对应的,2,d,进制编号,:,(N,1,N,2,N,l,),2,d,再按照下面的公式计算其,Morton,编号：,算法,空间编码,尺度转换,二进制编码,d,维空间编码,二进制交错及其解交错,涉及到的算法：查找,Parent,、,Children,、,Neighbor,、盒子中心坐标、盒子编码等,尺度转换,2,d,树结构每个参数都有自己的参数,D,j,映射原始的盒子的大小,x,1,min,x,1,max,x,1,min,x,1,max,x,j,max,-x,j,min,=D,j,j=1,d,把原始的盒子映射到单位盒子 上,实际的三维物理空间,D,j,=D,，,x,min,=(x,1,min,，,x,d,min,),称,x,为真实的坐标，为该点标准化的坐标,目的：,实施二进制编码转换，方便查找！,二进制编码,For example,d=1:,用十进制表示为,对于 不仅可以表示为 也可以表示为,用二进制表示 为：,对于,位交错,在,d,维单位盒子里,点坐标,其中第,k,个分量,也可以写成,交错二进制,(bit interleaving),的形式,:,解交错,d,维盒子在第,l,层的,Morton,键值为,可以,解交错,(,bit deinterleaving,),为,d,个数值代表该盒子的,d,维坐标,Number=76893,10,=,=111,2,=7,10,=111010,2,=58,10,=1001,2,=9,10,(9,58,7),Number,3,=,Number,2,=,Number,1,=,寻找给定点所在的盒子,在,d,维单位盒子里,点坐标的,交错,2,d,进制,形式,:,可以利用,2,d,进制,移位取整,寻找,给定点在第,l,层所在的,盒子,:,或者写成,dl,位的平移形式,:,查找分为,5,层八叉树结构的点,(0.7681,0.0459,0.3912),在第,3,层盒子的号,(1),二进制转化,(,0.11000,0.00001,0.01100),2,(2,）形成单个串,0.100101001000010,2,(3,）,3*3=9 (Number,3)=100101001,2,=297,3*5=15,(Number,5)=100101001000010,2,=19010,查找盒子中心,单位盒子第,l,层编号为,Num,的子盒子中心的二进制,d,维坐标形式,:,或者不依赖计数体系,写为,:,例如,:,八叉树第,5,层,10,进制编号为,533,的盒子，计算过程为,:,(1),将,Morton,编号转为,2,进制,:533,10,=1,000,010,101,2,(2),通过解交错得到该盒子的三维坐标,:,Num,3,=1001,2,=9,10,Num,2,=010,2,=2,10,Num,1,=001,2,=1,10,(3),计算盒子中心坐标,:,x,3c,(533,5)=2,-5,(9+0.5)=0.296875;x,2c,(533,5)=2,-5,(2+0.5)=0.078125;,x,1c,(533,5)=2,-5,(1+0.5)=0.04875;,因此,x,c,=(0.04875,0.078125,0.296875),算法,父盒子编码,计算某盒子的父盒子,父盒子编码与层次无关,其计算方法,:,除法的商取整,比如,11/4=2,，于是,例如：,于是,#5981,的父盒子是,#747,除以,2,d,相当于作,d,次,2,进制位移就可以了,算法,子盒子编码,计算某盒子的子盒子,子盒子是个集合,与所在的层无关,:,其计算方法为,:,For example,查找邻居盒子（,1,）,基于交错二进制的邻居查询算法,Step1.,解交错,(deinterleaving).,十进制,(,或二进制,),编号可转为,d,维坐标,:,Step2.,每个分坐标的邻居,:,于是邻居集为,:,其中,n,k,定义为,比如编号为,#26,的,2,维盒子,其,邻居,查询过程,如下,:,26,10,=(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,100101,110000,110001,十进制编号为,:13,24,25,15,27,37,48,49.,查找邻居盒子（,2,）,纲要,FMM,Data Structures,Parallelization,并行实施技术,(parallel implementation),MLFMM,具有很好的并行效率,大多数操作都在本地机上进行,例如，,Parant to Children,box to its interaction list,八叉树结构是很大的,需要生成和分布式存储,每个处理器只保持本地的基本树,大多数工作只在本地的基本树上进行,数据需要同步,父节点需要子节点上的值，但这两个节点在不同的处理器上,荷载平衡问题,但是，还存在：,分布式八叉树,负载平衡,相互作用表列,相邻结点的通信,次相邻点的通信,需要解决,并行计算步骤,构造压缩八叉树,近场矩阵计算,建立转移节点列表,远场矩阵计算,聚合,转移,发散,树结构的并行划分（,1,）,二维计算区域对应的分布式四叉树,树结构的并行划分（,2,）,构造分布式压缩八叉,树,(1),层,数,L=log,2,(D/d),D,和,d,为计算区域划分的最大和最小盒子尺寸,将,n,个基函数在,p,个处理机上按次序平均分配，再按照基函数的位置生成包含这些基函数的叶节点,由于不同基函数可包含于同一叶节点，因此这样的叶节点会同时存储在不同处理机上,RWG,基函数定义在各边上，并包含一对完整的三角形,P,1,P,2,P,3,A,1,A,2,A,3,用边的中点,代表整个边,各点都有相应的最底层,非空盒子,即叶节点,.,将,叶节点的,Morton,键值赋给中点所在的边,构造分布式压缩,八叉树,(2),采用并行排序算法对所有处理机中的基函数和叶节点排序，使个处理机包含相同数量的基函数,对每个处理机里的,N/p,个键值，采用快速排序,(,quicksort,),全局并行排序采用取样排序,(,samplesort,),它用到位排序,(,bitornic sort,),排序时用到的通信为,MPI_Allgather,和,MPI_Alltoall,每个处理机还包含下一处理机的第一个叶节点,并根据这些叶节点建立本地压缩八叉树，通过后序周游存储,各处理机将本地压缩树中位于从下一处理机得到的叶节点之后的节点发送到下一处理机，并按后序周游插到对应的位置,共享叶节点个数不超过,L,每个处理机接收的非本地节点不超过,7L,各处理机从下而上构造本地树的复杂度为,O(N/p)log(N/p),近场计算,Morton,次序保证兄弟节点的键值是相邻的，但并不是只有兄弟节点相邻,因此需要调用位交错和解交错函数去查找邻居节点,近场矩阵,A,near,只需要考虑最细层,(,第,L,层,),每个节点及其相邻节点所包含的基函数相互作用,;,必须按照,MOM,计算,并在迭代前存储,每个叶节点最多有,26,个相邻节点。若最细层每个盒子至多包含,c,个基函数,则每个叶节点上的计算量为,27c,2,如果同一棵子树上的相邻节点位于同一处理机,则无需通信,如果相邻节点位于不同的处理机上,则,使用,MPI_Allgather,获得每一处理机的第一个和最后一个叶节点的键值,.,并存储在长为,2p,的数组中,通过该数组的检索得到相邻叶节点,调用,MPI_Alltoall,实现数据,(,叶结点及其包含的基函数,),的分发,;,整个近场计算复杂度为,O(N/p)log(N/p),远场计算,局部的不变项,(,如最细层的,D,和,A),只需要分配到它对应的子树所在处理机上，全局不变项,(,如常数,),则分配到所有处理机上；由于下一层的计算依赖于上一层的信息，因此完成每一层的计算时，各个处理机都需要进行一次同步,存储时,采用内存循环策略,，它依赖于数据的相关性。聚集项,S,在层层上聚时分配内存，每当层层下推时某层的发散项,B,计算完毕，就将该层的,S,的内存释放掉；发散项,B,在层层下推时分配内存,每当处理机上某层的所有的,B,计算完毕，就将其父层的,B,释放掉,归并各处理机得到的结果,即为远场矩阵向量乘,通过归约得到整个系数矩阵与向量的乘积,通过并行的迭代法得到计算结果,(,等效电流,),每次迭代的矩阵向量乘法计算完成时，需要进行一次同步,完成计算结果的后处理，比如,RCS,的计算。,树结构代码,上聚,A M2M,内插值,下推,C M2L,B L2L,外插值,二叉树的例子,建立相互作用表列,相互作用,表列,(interaction list),包含每一层、每个节点的次相邻节点,次相邻节点指它们本身不相邻，但它们的父节点相邻,因此每个节点最多有,6,3,-3,3,=189,个次相邻节点，远多于,相邻节点（,26,）,需要在表列中注明次相邻点是否位于其它处理机,每层都要建立次相邻表列，但次相邻点可能不是物理意义的同一层,empty,M2L,相互作用表列,在迭代前,存储,每次远场作用,的转移项都会,调用该表列,聚集,对于每一层的每个盒子,聚集相当于将子层组,(t),中心的平面波移置到父层组,(P,t,),的中心,(C,t,),并通过内插值得到大数目的平面波,:,父节点,(,或部分子节点,),在其它处理机上的节点构成,剩余八叉树,它至多包含,8pL,个子节点,因此数据交换的量级为,O(logp+logL),对压缩八叉树的所有节点,(,包括剩余节点,),应用并行聚合算法,因此每个处理机的计算量为,O(N/p+L),后序遍历保证父节点在子节点之后,且紧挨子节点存储,转移,对于每一层的盒子,其远场作用只需要考虑次相邻组中心间的作用,即转移,.,假设组,(P,t,),和,(P,s,),的中心分别为,(C,t,C,s,),则,:,在第一次迭代时,若相互作用表列中的次相邻组中心不在同一处理机,则通过,MPI_Alltoall,发送到各个相应的处理机,按照上面的公式计算时,通过,MPI_Alltoall,接收次相邻组的平面波,对于压缩八叉树,次相邻组的盒子大小可能不一样,即插值取样点数可能不同,因此需要将取样点少的平面波转换为取样点多的,发散,发散过程是聚集过程的逆操作,对于每一层的每个盒子,发散相当于将父层组,(P,r,),中心的平面波移置到子层组,(r),的中心,(C,t,),并通过外插值得到小数目的平面波,并结合求和得到远场作用,:,父节点,(,或部分子节点,),在其它处理机上的节点构成剩余八叉树,它至多包含,8pL,个子节点,因此数据交换的量级为,O(logp+logL),对压缩八叉树的所有节点,(,包括剩余节点,),应用并行发散算法,因此每个处理机的计算量为,O(N/p+L),由于后序遍历能保证父节点在子节点之后,且紧挨子节点存储因此采用反向遍历后序存储的树节点,计算流程,网格剖分与几何信息读入,最细层的近场信息,聚集（内插值）,近场矩阵计算,近场作用,转移,(,次相邻组,),矩量法离散,迭代法求解,向量运算,矩阵,-,向量乘积,BLAS,LAPACK,构造分布式八叉树,发散（外插值）,各层的远场信息,远场作用,预条件子,线性方程组,电磁场积分方程,THANKS,