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

1 快速多极子方法的并行技术 国家973项目高性能科学计算研究 大规模并行计算研究 2 纲要 FMM Data Structures Parallelization 3 纲要 FMM Data Structures Parallelization 4 FMM in Computational Electromagnetics EFIE MFIE CFIE Green函数 5 积分方程的离散 RWG矢量基函数 MOM 离散 Rao-Wilson-Glisson 6 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 7 线性系统： 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,768 Finding:快速矩阵乘向量的算法(NlogN)； 并行实施。 8 Fast Multipole Methods(FMM) Introduced by Rokhlin 需要构建的函数，如Parent(n), ChildrenAll(n), Children(X;n,l), NeighborsAll(n,l), Neighbors(X;n,l) 24 Grouping 每个盒子在l=0,L层的指标n=Number=0,1,2ld-1. 由于E2(n,l)和 E3(n,l)是互补的， 因此对任意的n,l 25 纲要 FMM Data Structures Parallelization 26 Data Structure 构造树 压缩八叉树 建立连接 morton键 排序 27 构造树 离散点的空间层次划分 对应的四叉树及其压缩四叉树 28 确定层数L 根据读入模型的最大几何尺寸确定计算区域D,根据问题的参数确定最小盒子 尺寸d ,树结构的层数为L=log2(D/d) 第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) . 宽度周游 Eliminating empty leaves . 去掉空的叶节点 Compress Octree . 压缩八叉树 . 如果某中间节点仅包含一个子节点，则被其替换 每个节点用两个键值描述: 包含相同数据的最小单元和最大单元 30 构造树八叉树（3） 包含N个叶节点的压缩八叉树节点总数不超过2N-1 因此可以采用数组而不是链表来存储压缩八叉树 MLFMM基于后序周游的压缩八叉树数据结构 从键值最小的叶节点开始周游 每个节点都存储在其子节点之后，且紧挨其子节点存储 节点排序时，使用的是其所表示的最小单元键值 对于兄弟节点，按键值从小到大排序 各节点所表示的单元指的是它所包含的最小单元 后序周游存储方式是实现与分布无关的自动负载平衡并行MLFMM的 关键 分布无关自适应树(Distribution-Independent Adaptive Tree, DIAT) 构造2d维DIAT的复杂度为O(dnlogn) 树 节 点 存 储 31 Morton键 为什么不用指针？ 用指针指向Children的指针可以很方便的建立父子节点之间的关系， 从而构成树结构的拓扑结构。在串行程序，指针可以在全局存储空间 中寻址，效率很高也很准确。但在内存分布式并行环境中，一个计算 节点不能直接访问另一个计算节点上的存储空间，因此用于联系树结 构拓扑结构的指针只能在其所在的计算节点上才有意义，如果要让指 针所指向的树节点能够存储在其他节点上，就必须小心处理指针的变 换关系。以便将指针的指向正确地映射到其他计算节点上的内存空间 。这种转换，使得基于指针的树拓扑方法在分布式并行环境中实现起 来不仅很复杂，而且效率也不高。 Morton键技术是实现并行多层快速多极子的关键技术之一！ 原理：将空间坐标的二进制序列按位交叉，把空间中某一点在x、y、 z方向的坐标/序号映射为一个值，这个值就是morton键值。 32 Morton次序 位于第m层,在该层中编号为n的盒子, 一般采用Morton次序编号为(n, m). 左图为第三层的Morton次序 ， 右图为每一层编号,前三层分别有1, 4, 16个点. 33 Morton编码 小的灰盒子在第3层,本层编 号为11, 于是其Morton编号为 (11,3); (023)4=(11)10 采用四进制编号为(0,2,3); Morton编号(Num,l)， 在2d叉树中可以得到某盒 子对应的 2d进制编号: (N1, N2,Nl)2d 再按照下面的公式计算其Morton编号： 算法 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， xmin=(x1,min,，xd,min), 称x为真实的坐标， 为该点标准化的坐标 目的：实施二进制编码转换，方便查找！ 36 二进制编码 For example d=1: 用十进制表示为 对于 不仅可以表示为 也可以表示为 用二进制表示 为： 对于 37 位交错 在d维单位盒子里, 点坐标 其中第k个分量 也可以写成交错二进制(bit interleaving)的形式: 38 解交错 d维盒子在第l层的Morton键值为 可以解交错(bit deinterleaving)为d个数值代表该盒子的d维坐标 Number=7689310= =1112=710 =1110102=5810 =10012=910 (9, 58, 7) Number3= Number2= Number1= 39 寻找给定点所在的盒子 在d维单位盒子里, 点坐标的交错2d进制形式: 可以利用2d进制移位取整寻找给定点在第l层所在的盒子: 或者写成dl位的平移形式: 查找分为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 3*5=15 (Number,5)=1001010010000102=19010 40 查找盒子中心 单位盒子第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+0.5)=0.078125; x1c(533, 5)=2-5(1+0.5)=0.04875; 因此xc=(0.04875, 0.078125, 0.296875) 算法 41 父盒子编码 计算某盒子的父盒子 父盒子编码与层次无关, 其计算方法: 除法的商取整, 比如11/4=2，于是 例如： 于是#5981的父盒子是#747 除以2d相当于作d次2进制位移就可以了 算法 42 子盒子编码 计算某盒子的子盒子 子盒子是个集合,与所在的层无关: 其计算方法为: For example 43 查找邻居盒子（1） 基于交错二进制的邻居查询算法 Step1.解交错(deinterleaving).十进制(或二进制)编号可转为d维坐标: Step2.每个分坐标的邻居: 于是邻居集为: 其中nk定义为 44 比如编号为#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, 100101, 110000, 110001 十进制编号为:13,24,25, 15,27, 37,48,49. 查找邻居盒子（2） 45 纲要 FMM Data Structures Parallelization 46 并行实施技术(parallel implementation) MLFMM具有很好的并行效率 大多数操作都在本地机上进行 例如，Parant to Children,box to its interaction list 八叉树结构是很大的 需要生成和分布式存储 每个处理器只保持本地的基本树 大多数工作只在本地的基本树上进行 数据需要同步 父节点需要子节点上的值，但这两个节点在 不同的处理器上 荷载平衡问题 但是，还存在： 分布式八叉树 负载平衡 相互作用表列 相邻结点的通信 次相邻点的通信 需要解决 47 并行计算步骤 构造压缩八叉树 近场矩阵计算 建立转移节点列表 远场矩阵计算 聚合 转移 发散 48 树结构的并行划分（1） 49 二维计算区域对应的分布式四叉树 树结构的并行划分（2） 50 构造分布式压缩八叉树(1) 层数L=log2(D/d), D和d为计算区域划分的最大和最小盒子尺寸 将n个基函数在p个处理机上按次序平均分配，再按照基函数的位置生 成包含这些基函数的叶节点 由于不同基函数可包含于同一叶节点，因此这样的叶节点会同时存 储在不同处理机上 RWG基函数定义在各边上，并包含一对完整的三角形 P1 P2 P3 A1 A2 A3 用边的中点代表整个边, 各点都有相应的最底层 非空盒子, 即叶节点.将 叶节点的Morton键值赋 给中点所在的边 51 构造分布式压缩八叉树(2) 采用并行排序算法对所有处理机中的基函数和叶节点排序，使个处理 机包含相同数量的基函数 对每个处理机里的N/p个键值，采用快速排序(quicksort) 全局并行排序采用取样排序(samplesort), 它用到位排序(bitornic sort) 排序时用到的通信为MPI_Allgather和MPI_Alltoall 每个处理机还包含下一处理机的第一个叶节点, 并根据这些叶节点建立 本地压缩八叉树，通过后序周游存储 各处理机将本地压缩树中位于从下一处理机得到的叶节点之后的节点 发送到下一处理机，并按后序周游插到对应的位置 共享叶节点个数不超过L, 每个处理机接收的非本地节点不超过7L 各处理机从下而上构造本地树的复杂度为O( (N/p)log(N/p) 52 近场计算 Morton次序保证兄弟节点的键值是相邻的，但并不是只有兄弟节点 相邻, 因此需要调用位交错和解交错函数去查找邻居节点 近场矩阵Anear只需要考虑最细层(第L层)每个节点及其相邻节点所 包含的基函数相互作用; 必须按照MOM计算, 并在迭代前存储 每个叶节点最多有26个相邻节点。若最细层每个盒子至多包含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 M2L B L2L 外插值 二叉树的例子 55 建立相互作用表列 相互作用表列(interaction list)包含每一层、每个节点的次相邻节点 次相邻节点指它们本身不相邻，但它们的父节点相邻 因此每个节点最多有63-33=189个次相邻节点，远多于相邻节点（26） 需要在表列中注明次相邻点是否位于其它处理机 每层都要建立次相邻表列，但次相邻点可能不是物理意义的同一层 empty M2L 相互作用表列 在迭代前存储, 每次远场作用 的转移项都会 调用该表列 56 聚集 对于每一层的每个盒子, 聚集相当于将子层组(t)中心的平面波移 置到父层组(Pt )的中心(Ct ), 并通过内插值得到大数目的平面波: 父节点(或部分子节点)在其它处理机上的节点构成剩余八叉树, 它至多包含8pL个子节点, 因此数据交换的量级为O(logp+logL) 对压缩八叉树的所有节点(包括剩余节点)应用并行聚合算法,因 此每个处理机的计算量