快速多極子方法的并行技術(shù).ppt

1、快速多極子方法的并行技術(shù),國(guó)家973項(xiàng)目高性能科學(xué)計(jì)算研究 大規(guī)模并行計(jì)算研究,綱要,FMM Data Structures Parallelization,綱要,FMM Data Structures Parallelization,FMM in Computational Electromagnetics,EFIE MFIE CFIE Green函數(shù),積分方程的離散,RWG矢量基函數(shù) MOM 離散,Rao-Wilson-Glisson,Method of Moments,Surface is Discretized into Patches (Basis Functions),Basis

2、Functions Interact through the Greens Function,Generates a Dense Method of Moments Matrix,線性系統(tǒng)：,Mx=s M是NN矩陣，x、s是N矢量 Direct solution (Gauss elimination,LU Decomposition,SVD,) 空間復(fù)雜度為O(N2) ，需要O(N3)次運(yùn)算； Iterative methods,空間復(fù)雜度仍為O(N2)，如果K(k largest N =32,768 Finding:快速矩陣乘向量的算法(NlogN)； 并行實(shí)施。,Fast Multipol

3、e Methods(FMM),Introduced by Rokhlin 需要構(gòu)建的函數(shù)，如Parent(n), ChildrenAll(n), Children(X;n,l), NeighborsAll(n,l), Neighbors(X;n,l),Grouping,每個(gè)盒子在l=0,.,L層的指標(biāo)n=Number=0,1,2ld-1.,由于E2(n,l)和E3(n,l)是互補(bǔ)的，因此對(duì)任意的n,l,綱要,FMM Data Structures Parallelization,Data Structure,構(gòu)造樹(shù) 壓縮八叉樹(shù) 建立連接 morton鍵 排序,構(gòu)造樹(shù),離散點(diǎn)的空間層次劃分,對(duì)應(yīng)

4、的四叉樹(shù)及其壓縮四叉樹(shù),確定層數(shù)L 根據(jù)讀入模型的最大幾何尺寸確定計(jì)算區(qū)域D,根據(jù)問(wèn)題的參數(shù)確定最小盒子 尺寸d ,樹(shù)結(jié)構(gòu)的層數(shù)為L(zhǎng)=log2(D/d) 第l-1層立方體等分為八個(gè)子立方體,形成第l層更小的立方體,l-1是l層的父層,l層是l-1層的子層. 形成相鄰組、次相鄰組、遠(yuǎn)場(chǎng)組 第l層的節(jié)點(diǎn)數(shù)不超過(guò)2dl個(gè),構(gòu)造樹(shù)八叉樹(shù)(1),構(gòu)造樹(shù)八叉樹(shù)（2）,procedure Octree_Build Octree = empty For i = 1 to n . 對(duì)所有的點(diǎn)做循環(huán) Octree_Insert(i, root) . 將點(diǎn)i插入八叉樹(shù)相應(yīng)的位置 end For . 八叉樹(shù)中可能有很

5、多空的葉節(jié)點(diǎn), 但它們的兄弟節(jié)點(diǎn)非空 Traverse the tree (via, say, breadth first search) . 寬度周游 Eliminating empty leaves . 去掉空的葉節(jié)點(diǎn) Compress Octree . 壓縮八叉樹(shù) . 如果某中間節(jié)點(diǎn)僅包含一個(gè)子節(jié)點(diǎn)，則被其替換 每個(gè)節(jié)點(diǎn)用兩個(gè)鍵值描述: 包含相同數(shù)據(jù)的最小單元和最大單元,構(gòu)造樹(shù)八叉樹(shù)（3）,包含N個(gè)葉節(jié)點(diǎn)的壓縮八叉樹(shù)節(jié)點(diǎn)總數(shù)不超過(guò)2N-1 因此可以采用數(shù)組而不是鏈表來(lái)存儲(chǔ)壓縮八叉樹(shù) MLFMM基于后序周游的壓縮八叉樹(shù)數(shù)據(jù)結(jié)構(gòu) 從鍵值最小的葉節(jié)點(diǎn)開(kāi)始周游 每個(gè)節(jié)點(diǎn)都存儲(chǔ)在其子節(jié)點(diǎn)之后，且緊

6、挨其子節(jié)點(diǎn)存儲(chǔ) 節(jié)點(diǎn)排序時(shí)，使用的是其所表示的最小單元鍵值 對(duì)于兄弟節(jié)點(diǎn)，按鍵值從小到大排序 各節(jié)點(diǎn)所表示的單元指的是它所包含的最小單元 后序周游存儲(chǔ)方式是實(shí)現(xiàn)與分布無(wú)關(guān)的自動(dòng)負(fù)載平衡并行MLFMM的關(guān)鍵 分布無(wú)關(guān)自適應(yīng)樹(shù)(Distribution-Independent Adaptive Tree, DIAT) 構(gòu)造2d維DIAT的復(fù)雜度為O(dnlogn),樹(shù) 節(jié) 點(diǎn) 存 儲(chǔ),Morton鍵,為什么不用指針？ 用指針指向Children的指針可以很方便的建立父子節(jié)點(diǎn)之間的關(guān)系，從而構(gòu)成樹(shù)結(jié)構(gòu)的拓?fù)浣Y(jié)構(gòu)。在串行程序，指針可以在全局存儲(chǔ)空間中尋址，效率很高也很準(zhǔn)確。但在內(nèi)存分布式并行環(huán)境中，一

7、個(gè)計(jì)算節(jié)點(diǎn)不能直接訪問(wèn)另一個(gè)計(jì)算節(jié)點(diǎn)上的存儲(chǔ)空間，因此用于聯(lián)系樹(shù)結(jié)構(gòu)拓?fù)浣Y(jié)構(gòu)的指針只能在其所在的計(jì)算節(jié)點(diǎn)上才有意義，如果要讓指針?biāo)赶虻臉?shù)節(jié)點(diǎn)能夠存儲(chǔ)在其他節(jié)點(diǎn)上，就必須小心處理指針的變換關(guān)系。以便將指針的指向正確地映射到其他計(jì)算節(jié)點(diǎn)上的內(nèi)存空間。這種轉(zhuǎn)換，使得基于指針的樹(shù)拓?fù)浞椒ㄔ诜植际讲⑿协h(huán)境中實(shí)現(xiàn)起來(lái)不僅很復(fù)雜，而且效率也不高。 Morton鍵技術(shù)是實(shí)現(xiàn)并行多層快速多極子的關(guān)鍵技術(shù)之一！ 原理：將空間坐標(biāo)的二進(jìn)制序列按位交叉，把空間中某一點(diǎn)在x、y、z方向的坐標(biāo)/序號(hào)映射為一個(gè)值，這個(gè)值就是morton鍵值。,Morton次序,位于第m層,在該層中編號(hào)為n的盒子, 一般采用Morton次

8、序編號(hào)為(n, m). 左圖為第三層的Morton次序 ， 右圖為每一層編號(hào),前三層分別有1, 4, 16個(gè)點(diǎn).,Morton編碼,小的灰盒子在第3層,本層編號(hào)為11, 于是其Morton編號(hào)為(11,3); (023)4=(11)10 采用四進(jìn)制編號(hào)為(0,2,3); Morton編號(hào)(Num,l)， 在2d叉樹(shù)中可以得到某盒子對(duì)應(yīng)的 2d進(jìn)制編號(hào): (N1, N2,Nl)2d,再按照下面的公式計(jì)算其Morton編號(hào)：,算法,空間編碼,尺度轉(zhuǎn)換 二進(jìn)制編碼 d維空間編碼 二進(jìn)制交錯(cuò)及其解交錯(cuò) 涉及到的算法：查找Parent、Children、Neighbor、盒子中心坐標(biāo)、盒子編碼等,尺度轉(zhuǎn)

9、換,2d樹(shù)結(jié)構(gòu)每個(gè)參數(shù)都有自己的參數(shù)Dj映射原始的盒子的大小 x1,min, x1,max x1,min, x1,max ,xj,max-xj,min=Dj,j=1,d 把原始的盒子映射到單位盒子 上 實(shí)際的三維物理空間Dj=D， xmin=(x1,min,，xd,min), 稱x為真實(shí)的坐標(biāo)， 為該點(diǎn)標(biāo)準(zhǔn)化的坐標(biāo) 目的：實(shí)施二進(jìn)制編碼轉(zhuǎn)換，方便查找！,二進(jìn)制編碼,For example d=1: 用十進(jìn)制表示為 對(duì)于 不僅可以表示為 也可以表示為 用二進(jìn)制表示 為： 對(duì)于,位交錯(cuò),在d維單位盒子里, 點(diǎn)坐標(biāo) 其中第k個(gè)分量 也可以寫(xiě)成交錯(cuò)二進(jìn)制(bit interleaving)的形式:,解

10、交錯(cuò),d維盒子在第l層的Morton鍵值為 可以解交錯(cuò)(bit deinterleaving)為d個(gè)數(shù)值代表該盒子的d維坐標(biāo),Number=7689310=,=1112=710,=1110102=5810,=10012=910,(9, 58, 7),Number3=,Number2=,Number1=,尋找給定點(diǎn)所在的盒子,在d維單位盒子里, 點(diǎn)坐標(biāo)的交錯(cuò)2d進(jìn)制形式: 可以利用2d進(jìn)制移位取整尋找給定點(diǎn)在第l層所在的盒子: 或者寫(xiě)成dl位的平移形式:,查找分為5層八叉樹(shù)結(jié)構(gòu)的點(diǎn)(0.7681, 0.0459,0.3912)在第3層盒子的號(hào) (1)二進(jìn)制轉(zhuǎn)化(0.11000,0.00001,0

11、.01100)2 (2）形成單個(gè)串0.1001010010000102 (3）3*3=9 (Number,3)=1001010012=297 3*5=15 (Number,5)=1001010010000102=19010,查找盒子中心,單位盒子第l層編號(hào)為Num的子盒子中心的二進(jìn)制d維坐標(biāo)形式: 或者不依賴計(jì)數(shù)體系,寫(xiě)為: 例如: 八叉樹(shù)第5層10進(jìn)制編號(hào)為533的盒子，計(jì)算過(guò)程為: (1)將Morton編號(hào)轉(zhuǎn)為2進(jìn)制: 53310=1, 000, 010, 1012 (2)通過(guò)解交錯(cuò)得到該盒子的三維坐標(biāo): Num3=10012=910 , Num2=0102=210 , Num1=0012

12、=110 (3)計(jì)算盒子中心坐標(biāo): 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),算法,父盒子編碼,計(jì)算某盒子的父盒子 父盒子編碼與層次無(wú)關(guān), 其計(jì)算方法: 除法的商取整, 比如11/4=2，于是 例如： 于是#5981的父盒子是#747 除以2d相當(dāng)于作d次2進(jìn)制位移就可以了,算法,子盒子編碼,計(jì)算某盒子的子盒子 子盒子是個(gè)集合,與所在的層無(wú)關(guān): 其計(jì)算方法為: For ex

13、ample,查找鄰居盒子（1）,基于交錯(cuò)二進(jìn)制的鄰居查詢算法 Step1.解交錯(cuò)(deinterleaving).十進(jìn)制(或二進(jìn)制)編號(hào)可轉(zhuǎn)為d維坐標(biāo): Step2.每個(gè)分坐標(biāo)的鄰居: 于是鄰居集為: 其中nk定義為,比如編號(hào)為#26的 2維盒子,其鄰居 查詢過(guò)程如下: 2610=(01 10 10)2 解交錯(cuò)形式為 (011,100) 2=(3,4) 10 其相鄰盒子為 (2,3),(2,4), (2,5) ; (3,3),(3,5); (4,3),(4,4),(4,5); 8個(gè)點(diǎn)的二進(jìn)制: (010,011) 2 , (010,100) 2 , (100,101) 2,相鄰盒子編號(hào)的交錯(cuò)

14、(interleaving)二進(jìn)制形式: 001101, 011000, 011001, 001111, 011011, 100101, 110000, 110001 十進(jìn)制編號(hào)為:13,24,25, 15,27, 37,48,49.,查找鄰居盒子（2）,綱要,FMM Data Structures Parallelization,并行實(shí)施技術(shù)(parallel implementation),MLFMM具有很好的并行效率 大多數(shù)操作都在本地機(jī)上進(jìn)行 例如，Parant to Children,box to its interaction list ,八叉樹(shù)結(jié)構(gòu)是很大的 需要生成和分布式存儲(chǔ)

15、每個(gè)處理器只保持本地的基本樹(shù) 大多數(shù)工作只在本地的基本樹(shù)上進(jìn)行 數(shù)據(jù)需要同步 父節(jié)點(diǎn)需要子節(jié)點(diǎn)上的值，但這兩個(gè)節(jié)點(diǎn)在不同的處理器上 荷載平衡問(wèn)題,但是，還存在：,分布式八叉樹(shù) 負(fù)載平衡 相互作用表列 相鄰結(jié)點(diǎn)的通信 次相鄰點(diǎn)的通信,需要解決,并行計(jì)算步驟,構(gòu)造壓縮八叉樹(shù) 近場(chǎng)矩陣計(jì)算 建立轉(zhuǎn)移節(jié)點(diǎn)列表 遠(yuǎn)場(chǎng)矩陣計(jì)算 聚合 轉(zhuǎn)移 發(fā)散,樹(shù)結(jié)構(gòu)的并行劃分（1）,二維計(jì)算區(qū)域?qū)?yīng)的分布式四叉樹(shù),樹(shù)結(jié)構(gòu)的并行劃分（2）,構(gòu)造分布式壓縮八叉樹(shù)(1),層數(shù)L=log2(D/d), D和d為計(jì)算區(qū)域劃分的最大和最小盒子尺寸 將n個(gè)基函數(shù)在p個(gè)處理機(jī)上按次序平均分配，再按照基函數(shù)的位置生成包含這些基函數(shù)的葉節(jié)

16、點(diǎn) 由于不同基函數(shù)可包含于同一葉節(jié)點(diǎn)，因此這樣的葉節(jié)點(diǎn)會(huì)同時(shí)存儲(chǔ)在不同處理機(jī)上 RWG基函數(shù)定義在各邊上，并包含一對(duì)完整的三角形,P1,P2,P3,A1,A2,A3,用邊的中點(diǎn)代表整個(gè)邊, 各點(diǎn)都有相應(yīng)的最底層 非空盒子, 即葉節(jié)點(diǎn).將 葉節(jié)點(diǎn)的Morton鍵值賦給中點(diǎn)所在的邊,構(gòu)造分布式壓縮八叉樹(shù)(2),采用并行排序算法對(duì)所有處理機(jī)中的基函數(shù)和葉節(jié)點(diǎn)排序，使個(gè)處理機(jī)包含相同數(shù)量的基函數(shù) 對(duì)每個(gè)處理機(jī)里的N/p個(gè)鍵值，采用快速排序(quicksort) 全局并行排序采用取樣排序(samplesort), 它用到位排序(bitornic sort) 排序時(shí)用到的通信為MPI_Allgather和

17、MPI_Alltoall 每個(gè)處理機(jī)還包含下一處理機(jī)的第一個(gè)葉節(jié)點(diǎn), 并根據(jù)這些葉節(jié)點(diǎn)建立本地壓縮八叉樹(shù)，通過(guò)后序周游存儲(chǔ) 各處理機(jī)將本地壓縮樹(shù)中位于從下一處理機(jī)得到的葉節(jié)點(diǎn)之后的節(jié)點(diǎn)發(fā)送到下一處理機(jī)，并按后序周游插到對(duì)應(yīng)的位置 共享葉節(jié)點(diǎn)個(gè)數(shù)不超過(guò)L, 每個(gè)處理機(jī)接收的非本地節(jié)點(diǎn)不超過(guò)7L 各處理機(jī)從下而上構(gòu)造本地樹(shù)的復(fù)雜度為O( (N/p)log(N/p),近場(chǎng)計(jì)算,Morton次序保證兄弟節(jié)點(diǎn)的鍵值是相鄰的，但并不是只有兄弟節(jié)點(diǎn)相鄰, 因此需要調(diào)用位交錯(cuò)和解交錯(cuò)函數(shù)去查找鄰居節(jié)點(diǎn) 近場(chǎng)矩陣Anear只需要考慮最細(xì)層(第L層)每個(gè)節(jié)點(diǎn)及其相鄰節(jié)點(diǎn)所包含的基函數(shù)相互作用; 必須按照MOM計(jì)算

18、, 并在迭代前存儲(chǔ) 每個(gè)葉節(jié)點(diǎn)最多有26個(gè)相鄰節(jié)點(diǎn)。若最細(xì)層每個(gè)盒子至多包含c個(gè)基函數(shù), 則每個(gè)葉節(jié)點(diǎn)上的計(jì)算量為27c2 如果同一棵子樹(shù)上的相鄰節(jié)點(diǎn)位于同一處理機(jī), 則無(wú)需通信 如果相鄰節(jié)點(diǎn)位于不同的處理機(jī)上, 則 使用MPI_Allgather獲得每一處理機(jī)的第一個(gè)和最后一個(gè)葉節(jié)點(diǎn)的鍵值.,并存儲(chǔ)在長(zhǎng)為2p的數(shù)組中, 通過(guò)該數(shù)組的檢索得到相鄰葉節(jié)點(diǎn) 調(diào)用MPI_Alltoall實(shí)現(xiàn)數(shù)據(jù)(葉結(jié)點(diǎn)及其包含的基函數(shù))的分發(fā); 整個(gè)近場(chǎng)計(jì)算復(fù)雜度為O( (N/p)log(N/p),遠(yuǎn)場(chǎng)計(jì)算,局部的不變項(xiàng)(如最細(xì)層的D和A)只需要分配到它對(duì)應(yīng)的子樹(shù)所在處理機(jī)上，全局不變項(xiàng)(如常數(shù))則分配到所有處理機(jī)

19、上；由于下一層的計(jì)算依賴于上一層的信息，因此完成每一層的計(jì)算時(shí)，各個(gè)處理機(jī)都需要進(jìn)行一次同步 存儲(chǔ)時(shí)采用內(nèi)存循環(huán)策略，它依賴于數(shù)據(jù)的相關(guān)性。聚集項(xiàng)S在層層上聚時(shí)分配內(nèi)存，每當(dāng)層層下推時(shí)某層的發(fā)散項(xiàng)B計(jì)算完畢，就將該層的S的內(nèi)存釋放掉；發(fā)散項(xiàng)B在層層下推時(shí)分配內(nèi)存, 每當(dāng)處理機(jī)上某層的所有的B計(jì)算完畢，就將其父層的B釋放掉 歸并各處理機(jī)得到的結(jié)果, 即為遠(yuǎn)場(chǎng)矩陣向量乘 通過(guò)歸約得到整個(gè)系數(shù)矩陣與向量的乘積 通過(guò)并行的迭代法得到計(jì)算結(jié)果(等效電流), 每次迭代的矩陣向量乘法計(jì)算完成時(shí)，需要進(jìn)行一次同步 完成計(jì)算結(jié)果的后處理，比如RCS的計(jì)算。,樹(shù)結(jié)構(gòu)代碼,上聚 A M2M 內(nèi)插值 下推 C M2L

20、 B L2L 外插值 二叉樹(shù)的例子,建立相互作用表列,相互作用表列(interaction list)包含每一層、每個(gè)節(jié)點(diǎn)的次相鄰節(jié)點(diǎn) 次相鄰節(jié)點(diǎn)指它們本身不相鄰，但它們的父節(jié)點(diǎn)相鄰 因此每個(gè)節(jié)點(diǎn)最多有63-33=189個(gè)次相鄰節(jié)點(diǎn)，遠(yuǎn)多于相鄰節(jié)點(diǎn)（26） 需要在表列中注明次相鄰點(diǎn)是否位于其它處理機(jī) 每層都要建立次相鄰表列，但次相鄰點(diǎn)可能不是物理意義的同一層,empty,M2L,相互作用表列 在迭代前存儲(chǔ), 每次遠(yuǎn)場(chǎng)作用 的轉(zhuǎn)移項(xiàng)都會(huì) 調(diào)用該表列,聚集,對(duì)于每一層的每個(gè)盒子, 聚集相當(dāng)于將子層組(t)中心的平面波移置到父層組(Pt )的中心(Ct ), 并通過(guò)內(nèi)插值得到大數(shù)目的平面波: 父節(jié)點(diǎn)(或部分子節(jié)點(diǎn))在其它處理機(jī)上的節(jié)點(diǎn)構(gòu)成剩余八叉樹(shù),它至多包含8pL個(gè)子節(jié)點(diǎn), 因此數(shù)據(jù)交換的量級(jí)為O(logp+logL) 對(duì)壓縮八叉樹(shù)的所有節(jié)點(diǎn)(包括剩余節(jié)點(diǎn))應(yīng)用并行聚合算法,因此每個(gè)處理機(jī)的計(jì)算