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

1、1快速多极子方法的并行技术快速多极子方法的并行技术冯仰德冯仰德 王武王武 迟学斌迟学斌中科院计算机网络信息中心中科院计算机网络信息中心 超级计算中心超级计算中心 2007年年2月月5日日 国家国家973项目项目高性能科学计算研究高性能科学计算研究大规模并行计算研究大规模并行计算研究2纲要 fmm data structures parallelization3纲要 fmm data structures parallelization4fmm in computational electromagnetics efie mfie cfie green函数函数2( )( ) , , ( )/()

2、( );islsljkkgds tjt e rrjir,rj r| - |( , )/4| - |jk r rg r rer r( )/2( )( ) ,( )( )() ;iskkgds t j rt njt n h rjj rr,r(1)cfieefiemfie 5积分方程的离散积分方程的离散 rwg矢量基函数矢量基函数 mom 离散离散1( )( )niiijj rf r+-( )/2, ( ), ;( )( )/2, ( ), .iiiiiiiiiiiilaiftlaifti r rrrrf r r rrrr( )( )(1) ( )/2( ),( )( )(1)( ).mnmnnns

3、iimmszlkdsvds frff rnffre rn h rn=1 , 1,2,.nmnnmzjvmnrao-wilson-glisson6method 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

4、 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:快速矩阵乘向量的算法(nlogn)； 并行实施。8fast multipole methods(fmm)introduced by rokhlin &greengard in 1987called one of the 10 most significant advances in

5、 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)加速线性系统求解，如果用迭代方法，k步收敛，每步用矩阵矢量相乘，使计算复杂度由o(n3)变为o(knlogn)1( )(), njijijjijis xxxs 9fmm:application molecular and stellar dynamics com

6、putation 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 equations10fmm: fundament 格林函数的加法定理 jlpl平面波展开|(2)1e( 1) (21) ()()()|jkllllljklj kd hkr p r

7、+dd rr+djl_第一类球面bessel函数hl(2) 第二类球面hankel函数pl legendre多项式24()()()e()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 mult

8、ipole basics直接求解：分解：复杂度：o(mn)复杂度：o(pn+pm)cm12 fmm形式的矩阵向量乘积形式的矩阵向量乘积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

9、聚集过程，将基函数聚集成平面波函数，其结果表示k个平面波12 yd x转移过程，将得到的x1 平移到另一个盒子的中心，其结果表示该盒子中心的k个平面波发散过程，将得到的x2发散到盒子中的基函数。m2m , m2l , l2l： 聚集 转移 发散 m: 多极子展开 l: 局部展开14fmm123( )1( )( )2( )( )3( )( )( ,),( )( ,),( )( ,)iiinixe nnixennixenuuuiiiyy xyy xyy x由于e2(n)和e3(n)是互补的，因此对任意的n,23( )( )23( )( )( )( ,)( )( )innixenenuiyy xyy

10、15fmm algorithmstep1 m2m16step2 m2l17step3 l2lsummation18mlfmm algorithm upward pass: merge scattering matrices downward pass: construct splitting and exchange matrices final summationbased on: hierarchical grouping data structure19l2lm2mm2mm2lm2lm2lmultilevel multipole operatorsfinestlevelfinest -

11、1levell2lup treeacross treedown tree20mlfmm algorithmupward pass step1: 在最细的层盒子求解远场展开系数，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，因此次相邻

12、数目最多为p=3d(2d-1)；三维是189，二维是27，一维是3。局部到局部的转移，e3(n,l)和e4(n,l+1)产生e3(n,l+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,.,

13、l层的指标n=number=0,1,2ld-1.12324( , )( , )( , )( , ),( , ), )( , )0,.,21( , )( , )( ),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 nlrnlneighb

14、or nl le nlee nlnlneighbor nl le nle parent n le nlparent n lneighb表 示 盒 子的 空 间 区 域 ;表 示 盒 子与 其 相 邻的 空 间 区 域 ;表 示 盒 子及 其 相 邻的 外 部 空 间 区 域 ；表 示 父 盒 子及 其 相 邻 ( ),1),1 orparent n ll(的 空 间 区 域 ；1234( , )( , )12( , )( , )( , )( , )34( , )( , )( )( ,),( )( ,),( )( ,),( )( ,)iiiin ln liixen lxen ln ln liix

15、en lxen luuuuiiiiyy xyy xyy xyy x由于e2(n,l)和e3(n,l)是互补的，因此对任意的n,l( , )( , )23( )( )( )n ln l yyy25纲要 fmm data structures parallelization26data structure 构造树 压缩八叉树 建立连接 morton键 排序27构造树 离散点的空间层次划分离散点的空间层次划分 对应的四叉树及其压缩四叉树对应的四叉树及其压缩四叉树 28 确定层数l 根据读入模型的最大几何尺寸确定计算区域d,根据问题的参数确定最小盒子 尺寸d ,树结构的层数为l=log2(d/d) 第

16、l-1层立方体等分为八个子立方体,形成第l层更小的立方体,l-1是l层的父层,l层是l-1层的子层. 形成相邻组、次相邻组、远场组 第l层的节点数不超过2dl个 构造树八叉树八叉树(1)29构造树八叉树（八叉树（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) . 宽

17、度周游eliminating empty leaves . 去掉空的叶节点compress octree . 压缩八叉树 . 如果某中间节点仅包含一个子节点，则被其替换每个节点用两个键值描述: 包含相同数据的最小单元和最大单元30构造树构造树八叉树（八叉树（3）包含n个叶节点的压缩八叉树节点总数不超过2n-1 因此可以采用数组而不是链表来存储压缩八叉树mlfmm基于后序周游的压缩八叉树数据结构 从键值最小的叶节点开始周游 每个节点都存储在其子节点之后，且紧挨其子节点存储节点排序时，使用的是其所表示的最小单元键值 对于兄弟节点，按键值从小到大排序 各节点所表示的单元指的是它所包含的最小单元后序周

18、游存储方式是实现与分布无关的自动负载平衡并行mlfmm的关键 分布无关自适应树(distribution-independent adaptive tree, diat) 构造2d维diat的复杂度为o(dnlogn)树树节节点点存存储储31morton键为什么不用指针？ 用指针指向children的指针可以很方便的建立父子节点之间的关系，从而构成树结构的拓扑结构。在串行程序，指针可以在全局存储空间中寻址，效率很高也很准确。但在内存分布式并行环境中，一个计算节点不能直接访问另一个计算节点上的存储空间，因此用于联系树结构拓扑结构的指针只能在其所在的计算节点上才有意义，如果要让指针所指向的树节点能

19、够存储在其他节点上，就必须小心处理指针的变换关系。以便将指针的指向正确地映射到其他计算节点上的内存空间。这种转换，使得基于指针的树拓扑方法在分布式并行环境中实现起来不仅很复杂，而且效率也不高。morton键技术是实现并行多层快速多极子的关键技术之一！ 原理：将空间坐标的二进制序列按位交叉，把空间中某一点在x、y、z方向的坐标/序号映射为一个值，这个值就是morton键值。32morton次序次序位于第m层,在该层中编号为n的盒子,一般采用morton次序编号为(n, m).左图为第三层的morton次序 ，右图为每一层编号,前三层分别有1, 4, 16个点. 33morton编码小的灰盒子在第

20、3层,本层编号为11, 于是其morton编号为(11,3); (023)4=(11)10采用四进制编号为(0,2,3); morton编号(num,l)，在2d叉树中可以得到某盒子对应的 2d进制编号:(n1, n2,nl)2d再按照下面的公式计算其morton编号：111012112()(.)(2 ). (2 ).ddldlllnumn nnnumnnn算法34空间编码 尺度转换 二进制编码 d维空间编码 二进制交错及其解交错 涉及到的算法：查找parent、children、neighbor、盒子中心坐标、盒子编码等35尺度转换 2d树结构每个参数都有自己的参数dj映射原始的盒子的大小x

21、1,min, x1,max x1,min, x1,max ,xj,max-xj,min=dj,j=1,d把原始的盒子映射到单位盒子 上 实际的三维物理空间dj=d， xmin=(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

22、101(0.999999.)x 0,1x 1 2 32(0.) ,0,1;1,2,.jxbbbbj1x 21(0.11111.)x 37位交错位交错 在d维单位盒子里,点坐标其中第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层的mort

23、on键值为 可以解交错解交错(bit deinterleaving)为d个数值代表该盒子的d维坐标11 21112222122(.|.|.|.)ddlldlnumb bbb bbb bbnumber=7689310=1112=710=1110102=5810=10012=910(9, 58, 7)11112122212222122(.) ;(.) ;.(.) ;llddddlnumb bbnumb bbnumb bbnumber3=number2=number1=39寻找给定点所在的盒子寻找给定点所在的盒子在d维单位盒子里, 点坐标的交错交错2d进制进制形式:可以利用2d进制进制移位取整移位取

24、整寻找给定点在第l层所在的盒子:或者写成dl位的平移形式:11 21112222122(.|.|.|.) , 0,1ddljdlkjnumb bbb bbb bbb121 21222(0.)(,.,) , )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(3）3*3=9 (number,3)=1001010012=297

25、3*5=15 (number,5)=1001010010000102=1901040查找盒子中心查找盒子中心单位盒子第l层编号为num的子盒子中心的二进制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

26、+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,.,kckkklxnum 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 n

27、umnnn(, )(),1),parentnum lparent num l()/2 dparent numnum104410(11)(23)(2)(2).parentparent3103,(5981)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.

28、ddchildren numnumjj104444410101010(11 )(23 )230 ,231 ,232 ,233 44 ,45 ,46 ,47 childrenchildren1011 2111()(.) ()(2 ).(2 ),lldldllnumnn nchildren numnnn43查找邻居盒子（查找邻居盒子（1） 基于交错二进制的邻居查询算法基于交错二进制的邻居查询算法step1.解交错解交错(deinterleaving).十进制(或二进制)编号可转为d维坐标: step2.每个分坐标的邻居每个分坐标的邻居: 于是邻居集为: 其中nk定义为1,.,dnumnumnum1

29、,1,1, .,.kkkkn u mn u mn u mn u mkd1()( ,.,)dneighbor numnn, 0,21, 0;, 21.dkkkkkkkkkdkkkknnumnumnumif 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)

30、 2 , (100,101) 2相邻盒子编号的交错交错 (interleaving)二进制形式二进制形式: 001101, 011000, 011001, 001111, 011011, 100101, 110000, 110001十进制编号为:13,24,25, 15,27, 37,48,49.查找邻居盒子（查找邻居盒子（2）45纲要 fmm data structures parallelization46并行实施技术并行实施技术(parallel implementation) mlfmm具有很好的并行效率 大多数操作都在本地机上进行例如，parant to children,box t

31、o its interaction list 八叉树结构是很大的 需要生成和分布式存储 每个处理器只保持本地的基本树 大多数工作只在本地的基本树上进行数据需要同步 父节点需要子节点上的值，但这两个节点在不同的处理器上荷载平衡问题但是，还存在： 分布式八叉树 负载平衡 相互作用表列 相邻结点的通信 次相邻点的通信需要解决47并行计算步骤 构造压缩八叉树 近场矩阵计算 建立转移节点列表 远场矩阵计算 聚合 转移 发散48 树结构的并行划分（树结构的并行划分（1）49二维计算区域对应的分布式四叉树二维计算区域对应的分布式四叉树 树结构的并行划分（树结构的并行划分（2）50构造分布式压缩八叉构造分布式

32、压缩八叉树树(1) 层数l=log2(d/d), d和d为计算区域划分的最大和最小盒子尺寸将n个基函数在p个处理机上按次序平均分配，再按照基函数的位置生成包含这些基函数的叶节点 由于不同基函数可包含于同一叶节点，因此这样的叶节点会同时存储在不同处理机上 rwg基函数定义在各边上，并包含一对完整的三角形p1p2p3a1a2a3用边的中点代表整个边,各点都有相应的最底层非空盒子, 即叶节点.将叶节点的morton键值赋给中点所在的边51构造分布式压缩构造分布式压缩八叉树八叉树(2)采用并行排序算法对所有处理机中的基函数和叶节点排序，使个处理机包含相同数量的基函数 对每个处理机里的n/p个键值，采用

33、快速排序(quicksort) 全局并行排序采用取样排序(samplesort), 它用到位排序(bitornic sort) 排序时用到的通信为mpi_allgather和mpi_alltoall每个处理机还包含下一处理机的第一个叶节点, 并根据这些叶节点建立本地压缩八叉树，通过后序周游存储各处理机将本地压缩树中位于从下一处理机得到的叶节点之后的节点发送到下一处理机，并按后序周游插到对应的位置共享叶节点个数不超过l, 每个处理机接收的非本地节点不超过7l各处理机从下而上构造本地树的复杂度为o( (n/p)log(n/p)52近场计算近场计算morton次序保证兄弟节点的键值是相邻的，但并不是

34、只有兄弟节点相邻, 因此需要调用位交错和解交错函数去查找邻居节点近场矩阵anear只需要考虑最细层(第l层)每个节点及其相邻节点所包含的基函数相互作用; 必须按照mom计算, 并在迭代前存储每个叶节点最多有26个相邻节点。若最细层每个盒子至多包含c个基函数, 则每个叶节点上的计算量为27c2 如果同一棵子树上的相邻节点位于同一处理机, 则无需通信 如果相邻节点位于不同的处理机上, 则 使用mpi_allgather获得每一处理机的第一个和最后一个叶节点的键值.,并存储在长为2p的数组中, 通过该数组的检索得到相邻叶节点 调用mpi_alltoall实现数据(叶结点及其包含的基函数)的分发; 整

35、个近场计算复杂度为o( (n/p)log(n/p)53远场计算远场计算 局部的不变项(如最细层的d和a)只需要分配到它对应的子树所在处理机上，全局不变项(如常数)则分配到所有处理机上；由于下一层的计算依赖于上一层的信息，因此完成每一层的计算时，各个处理机都需要进行一次同步 存储时采用内存循环策略，它依赖于数据的相关性。聚集项s在层层上聚时分配内存，每当层层下推时某层的发散项b计算完毕，就将该层的s的内存释放掉；发散项b在层层下推时分配内存, 每当处理机上某层的所有的b计算完毕，就将其父层的b释放掉 归并各处理机得到的结果, 即为远场矩阵向量乘 通过归约得到整个系数矩阵与向量的乘积 通过并行的迭代法得到计算结果(等效电流), 每次迭代的矩阵向量乘法计算完成时，需要进行一次同步 完成计算结果的后处理，比如rcs的计算。54树结构代码树结构代码 上聚上聚 a m2m 内插值内插值 下