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

1快速多极子方法的并行技术冯仰德王武迟学斌中科院计算机网络信息中心超级计算中心[e-mail]2007年2月5日国家973项目——高性能科学计算研究大规模并行计算研究2纲要•FMM•DataStructures•Parallelization3纲要•FMM•DataStructures•Parallelization4FMMinComputationalElectromagnetics•EFIE•MFIE•CFIE•Green函数2()(),,()/()()';iSLSLjkkgdStJtErrJIr,r'Jr'|-'|(,')/4|-'|jkrrgrrerr()/2()(),()()()';iSKKgdStJrtnJtnHrJJr'r,r'(1)CFIEEFIEMFIE5积分方程的离散•RWG矢量基函数•MOM离散1()()NiiijJrfr+-()/2,(),;()()/2,(),.iiiiiiiiiiiilAifTlAifTiρrρrrrrfrρrρrrrr()()(1)[()/2()],()()(1)().mnmnnnSiimmSzLKdSvdSfrffrnffrErnHrn=1,1,2,,.NmnnmzjvmNRao-Wilson-Glisson6MethodofMoments•SurfaceisDiscretizedintoPatches(BasisFunctions)•BasisFunctionsInteractthroughtheGreen’sFunction,Grr•GeneratesaDenseMethodofMomentsMatrix11nnnnZIV,,,''ZrrrjiSjifGfds()rf()rtPulse7线性系统：Mx=sM是N×N矩阵，x、s是N矢量●Directsolution(Gausselimination,LUDecomposition,SVD,…)空间复杂度为O(N2)，需要O(N3)次运算；●Iterativemethods,空间复杂度仍为O(N2)，如果K(kN)步收敛，每步需要的矩阵乘向量的运算为O(N2)；—Moore’Law:processorspeeddoublesevery18months—amillionvariable,16generationsofMoore’LawbeforeaO(N2)algorithmwascomparablewithaO(N)algorithm—1GBRAM=10243=1,073,741,824bytes=largestN=32,768●Finding:快速矩阵乘向量的算法(NlogN)；并行实施。8FastMultipoleMethods(FMM)•IntroducedbyRokhlin&Greengardin1987•Calledoneofthe10mostsignificantadvancesincomputingof20thcentury•Speedsupmatix-vectorproducts(sums)ofaparticulartype•以上求和要求O(MN)运算复杂度•对给定的精度，FMM可以获得O(M+N)运算复杂度•可以加速matix-vectorproducts，使O(N2)变为O(NlogN)•加速线性系统求解，如果用迭代方法，k步收敛，每步用矩阵矢量相乘，使计算复杂度由O(N3)变为O(kNlogN)1()(),{}[]{}Njijijjijisxxxs9FMM:Application•MolecularandStellardynamics—computationofforcefieldsanddynamics•Solutionofacousticalscatteringproblems—HelmholtzEquation•ElectromagneticWaveScattering—Maxwell’sEquations•FluidMechanics:Potentialflow,vertexflow—Laplace/PoissonEquations10FMM:Fundament•格林函数的加法定理•jlpl平面波展开||(2)1e(1)(21)()()()||jkllllljkljkdhkrPr+ddrr+djl_—第一类球面Bessel函数hl(2)—第二类球面Hankel函数Pl—Legendre多项式24()()()e()njknnnSajjkdPPdkddrkrk注意到lz时，函数jl(z)和hl(2)(z)幅值大致保持为常数；lz时，函数jl(z)衰减非常快，而hl(2)(z)递增非常快。当dr时，上式在保证精度的情况下截断。则上式可以写为：||(2)1e(1)(21)()()()||jkLllllljkljkdhkrPr+ddrr+dL=kd+cln(kd+)Kd-源点到观察点的最大半径c-是一个依赖希望精度的常数=1最小的相对误差小于0.1=5相对误差小于10-6=10准确到双精度11FastMultipoleBasics直接求解：分解：复杂度：O(MN)复杂度：O(pN+pM){cm}12•FMM形式的矩阵向量乘积1'()'221'()'[()]16NijjijjnmGnmjGmKpipppjjpmGnmjGmzjzjkwTjDA'2''(')(,),ln(),2;e()(')';mpmmpLmmjkpppjjjTkrLkdckdKLdSkrrTkrAI-kkfr()e[()()(1)()]mjkpppipiiidSkrrDI-kkfrkfrn近场部分远场部分电磁场积分方程的多极子展开13FMM{}{}[][][]{}NearyyDTAx1{}[]{}xAx21{}[]{}xTx——聚集过程，将基函数聚集成平面波函数，其结果表示K个平面波12{}[]{}yDx——转移过程，将得到的{x1}平移到另一个盒子的中心，其结果表示该盒子中心的K个平面波——发散过程，将得到的{x2}发散到盒子中的基函数。M2M,M2L,L2L：聚集–转移–发散M:多极子展开L:局部展开14FMM123()1()()2()()3()()(,),()(,),()(,)iiinixEnnixEnnixEnuuuiiiyyxyyxyyx由于E2(n)和E3(n)是互补的，因此对任意的n,23()()23()()()(,)()()innixEnEnuiyyxyy15FMMAlgorithmStep1M2M16Step2M2L17Step3L2LSummation18MLFMMAlgorithm•UpwardPass:mergescatteringmatrices•DownwardPass:constructsplittingandexchangematrices•FinalSummationBasedOn:◆HierarchicalGrouping◆DataStructure19L2LM2MM2MM2LM2LM2LMultilevelMultipoleOperatorsFinestLevelFinest-1LevelL2LUpTreeAcrossTreeDownTree20MLFMMAlgorithm—UpwardPass•Step1:在最细的层盒子求解远场展开系数，xi∈E1(n,l),得到C(n,l)或，这也可以用于xi∈E3(n,l)•Step2:对于l=L-1,…,2,从step1得到值进行递归得到。同样适合xi∈E3(n,l)结果：得到分层组聚集系数(,)1()nly21MLFMMAlgorithm—DownwardPass•Step1:l=2,..L递归进行E4•Step2:任意一个非空组自身加上其邻接组最多可有3d个，其父组及父组的邻接组最多可形成3d×2d，因此次相邻数目最多为p=3d(2d-1)；三维是189，二维是27，一维是3。局部到局部的转移，E3(n,l)和E4(n,l+1)产生E3(n,l+1)结果：可以得到各分层组的转移系数22KeyWords•空间多层组划分•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,...,21}\(,)(,){(((),1),1))}\{(,),((,),)}ldInlnlInlnlNeighbornllInlInlInlChildrenNeighborParentnllnlNeighbornll12312422(,)(,)(,)(,){((,),)}(,)(0,0)\(,)(,){((,),)}(,)((),1)\(,)((),1)ddEnlRnlEnlRnlNeighbornllEnlEEnlnlNeighbornllEnlEParentnlEnlParentnlNeighb表示盒子的空间区域;表示盒子与其相邻的空间区域;表示盒子及其相邻的外部空间区域；表示父盒子及其相邻{(),1),1}orParentnll(的空间区域；1234(,)(,)12(,)(,)(,)(,)34(,)(,)()(,),()(,),()(,),()(,)iiiinlnliixEnlxEnlnlnliixEnlxEnluuuuiiiiyyxyyxyyxyyx由于E2(n,l)和E3(n,l)是互补的，因此对任意的n,l(,)(,)23()()()nlnlyyy25纲要•FMM•DataStructures•Parallelization26DataStructure•构造树—压缩八叉树•建立连接—morton键—排序27构造树•离散点的空间层次划分•对应的四叉树及其压缩四叉树28•确定层数L—根据读入模型的最大几何尺寸确定计算区域D,根据问题的参数确定最小盒子尺寸d,树结构的层数为L=log2(D/d)•第l-1层立方体等分为八个子立方体,形成第l层更小的立方体,l-1是l层的父层,l层是l-1层的子层.•形成相邻组、次相邻组、远场组•第l层的节点数不超过2dl个构造树—八叉树(1)29构造树—八叉树（2）procedureOctree_BuildOctree={empty}Fori=1ton...对所有的点做循环Octree_Insert(i,root)...将点i插入八叉树相应的位置endFor...八叉树中可能有很多空的叶节点,但它们的兄弟节点非空Traversethetree(via,say,breadthfirstsearch)...宽度周游Eliminatingemptyleaves...去掉空的叶节点CompressOctree...压缩八叉树...如果某中间节点仅包含一个子节点，则被其替换每个节点用两个键值描述:包含相同数据的最小单元和最大单元30构造树—八叉树（3）•包含N个叶节点的压缩八叉树节点总数不超过2N-1–因此可以采用数组而不是链表来存储压缩八叉树•MLFMM基于后序周游的压缩八叉树数据结构–从键