LecNote-14-并行算法設計(二)

第十四講并行算法設計(二)分治并行(partitioninganddivide-and-conquerstrategies)一維FFT問題多體問題N-Body的Barnes-Hut算法流水并行(pipelinedcomputation)Gauss-Seidel迭代法求解線性方程組一維FFT在前面的二維FFT算法中，我們沒有考慮矩陣的每一行(列)是如何實現，只是假定已經有了一個串行的一維FFT/DFT函數，并沒有考慮每一行(列)內部是否有并行性但是，如果計算任務是對一個一維的向量進行DFT，且n非常大，比如n=1M/G并行性：每個bk都是可以并行計算的局部性：每個bk的計算都需要整個的向量(a0,a1,…,an-1)從上一講中，我們知道，為了提高并行的性能，必須對(a0,a1,…,an-1)進行劃分提高程序對問題規(guī)模的scalability：最大可解問題的規(guī)模可以隨處理器數量增加提高程序中的數據訪問效率：通常，涉及的數據規(guī)模越小，訪問的命中率就越高通信問題？bk=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)(0,2,4,6,8,10,12,14)(0,4,8,12)(0,8)wk(4,12)wk(2,6,10,14)(2,10)wk(6,14)wk(1,3,5,7,9,11,13,15)(1,5,9,13)(1,9)wk(5,13)wk(3,7,11,15)(3,11)wk(7,15)a0a1a2a3a4a5a6a7a8a9a10a11a12a13a14a15b0b1b2b3b4b5b6b7b8b9b10b11b12b13b14b15a0a1a2a3a4a5a6a7a8a9a10a11a12a13a14a15b0b1b2b3b4b5b6b7b8b9b10b11b12b13b14b15bk=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14)(0,2,4,6,8,10,12,14)(0,4,8,12)(0,8)wk(4,12)wk(2,6,10,14)(2,10)wk(6,14)wk(1,3,5,7,9,11,13)(1,5,9,13)(1,9)wk(5,13)wk(3,7,11)(3,11)wk(7)a0a1a2a3a4a5a6a7a8a9a10a11a12a13a14b0b1b2b3b4b5b6b7b8b9b10b11b12b13b14一維FFT并行特征如果在存在關系的兩個數據之間連接一條邊，所有數據之間形成一個全連同的圖每個處理器要執(zhí)行l(wèi)og(p)次通信，其中p是處理器的數量每次通信交換數據量為n/p每一次通信時，都要對處理器進行重新分組每個處理器上的循環(huán)空間nlog(n)p每個循環(huán)的計算量：一維并行FFT的scalability每個處理器上的負載n/p個bk計算每個處理器上的數據量于p成反比通信的次數log(p)消息大小n/p分治并行partitioninganddivide-and-conquerstrategies劃分策略(partitioning)：把一個問題分解成若干個組成部分通常，需要對每個組成部分的結果進行合成(combine)后，才能夠獲得整個問題的結果例如：理想并行采用的是一種劃分策略，這種策略下，在把各個部分的結果合成為整個問題的結果時，幾乎不需要什么額外的運算。劃分策略的分類數據劃分(datapartitioning)/域分解(domaindecomposition)：把對計算的數據分解成一組數據子集，在不同的子集上并行的執(zhí)行處理。要求每個數據子集上執(zhí)行的運算沒有依賴關系不同的數據子集可以采用數據復制(datareplication)的策略，使它們有數據“重疊”例如在二維的FFT計算中，每一個super-step上，都執(zhí)行的是數據分解，把整個二維陣列按照“行”/“列”劃分，在每一“行”/“列”上并發(fā)執(zhí)行一維FFT功能分解(functionaldecomposition)：把計算任務劃分成一組獨立的功能模塊，并發(fā)執(zhí)行每一個功能模塊。不同的功能模塊所需要的初始數據可以有“重疊”例如Jacobi迭代(求解方程組AX=B)，每個X[i](t)的計算作為一個功能模塊，每個功能模塊出了涉及A的第i行、B[i]為。，都需要X(t-1)(“重疊”的數據部分

)分治策略(divide-and-conquer)：把一個復雜的問題分解成一組子問題，每個子問題與原問題的形式相同，但規(guī)模比原問題小，而且對子問題可以采用這一策略進一步細分這是另一種特殊情形的partitioning策略：子問題與原問題除了問題規(guī)模外，完全相同例如：一維FFT(A)，其中A是一個長度為n的向量。我們采用的并行策略是令FFT(A)=FFT([FFT(A0dd),FFT(Aeven)])把FFT(A)劃分為三個一維FFT子問題：FFT(A0dd)、FFT(Aeven)和FFT([FFT(A0dd),FFT(Aeven)]FFT(A0dd)執(zhí)行對A的奇數位置的元素組成的向量的一維FFT，問題規(guī)模是原問題規(guī)模的一半FFT(Aeven)執(zhí)行對A的偶數位置的元素組成的向量的一維FFT，問題規(guī)模是原問題規(guī)模的一半FFT([FFT(A0dd),FFT(Aeven)])執(zhí)行對一個由兩個元素所組成向量的一維FFTM-ary的分治策略把一個復雜的問題分解成M個子問題，M2，每個子問題與原問題的形式相同當子問題的規(guī)模還是太大時，進一步把每個子問題分解成M個粒度更細的子問題如此遞歸，直到每個子問題的規(guī)模合適為止一個4-ary分治的例子用分治策略解決N_Body問題有N個粒子(body)

在天體物理學里代表星體，例如地球、太陽、火星等在分子動力學里代表構成一個分子的各個原子……每個粒子有一定的狀態(tài)：位置、速度、加速度、能量、溫度等，這些狀態(tài)是隨時間變化的影響一個粒子狀態(tài)變化的因素是粒子的初始狀態(tài)其它粒子對它的作用力任意兩個粒子之間都有牛頓萬有引力，與兩個粒子的距離有關對不同物理問題，粒子之間還存在其它與雙方狀態(tài)有關的作用力。例如在兩個原子間還存在勢能力(由原子的自轉速度、原子間的距離等有關)。問題：對一個N_Body系統，給定其中各個粒子的初始狀態(tài)最終狀態(tài)是什么樣的，即其中每個粒子的狀態(tài)不再變化時，各個粒子的狀態(tài)是什么樣的？從初態(tài)到終態(tài)的變化過程中，變化的軌跡是什么樣的？即對于我們關心的某(幾)個狀態(tài)量(例如電磁場勢能的分布)，隨時間變化的規(guī)律是什么樣的？參考《DesigningandBuildingParallelPrograms》第1.4.2節(jié)把N個粒子按照“block”方式，劃分成M組，M是處理器的數量每組作為一個粒子的子集，分配給一個處理器在每個處理器上，除了存儲本地粒子的子集local_body_set外，開辟一個buffer，其容量大小是能夠存儲一個粒子的子集。在每個處理器上，用它的local_body_set對buffer進行初始化在每個時間步上，執(zhí)行下列循環(huán)從1到M，執(zhí)行循環(huán)對local_body_set中的每個粒子，計算buffer中各粒子對它的作用力把buffer中的內容發(fā)送給“左”鄰居(0#處理器的“左”鄰居是M-1號處理器)從“右”鄰居接收消息，把消息數據存儲在buffer中(M-1#處理器的“右”鄰居是0#處理器)根據得到的作用力，更新local_body_set中的各個粒子的狀態(tài)以100各粒子、4個處理器為例每個時間步上粒子0~24粒子0~24粒子25~49粒子25~49粒子50~74粒子50~74粒子75~99粒子75~99粒子0~24粒子25~49粒子25~49粒子50~74粒子50~74粒子75~99粒子75~99粒子0~24粒子0~24粒子50~74粒子25~49粒子75~99粒子50~74粒子0~24粒子75~99粒子25~49粒子0~24粒子75~99粒子25~49粒子0~24粒子50~74粒子25~49粒子75~99粒子50~740#處理器0#處理器0#處理器0#處理器N_Body問題的一種并行算法：性能分析對問題規(guī)模的scalability數據的scalability：計算數據以“block”方式劃分在各個處理器上，數據存儲的能力與處理器數量成線性關系速度的scalability：每個處理器上的主要運算開銷是計算local_body_set各個粒子所受作用力，每個時間步上的復雜度為O(N2M)=O(N2)，其中N是粒子數量、M是處理器數量通信復雜度：在每個時間步上通信啟動O(M)每次通信的數據傳遞開銷O(NM)總的數據傳遞開銷O(N)N是粒子數量、M是處理器數量因此，采用這種算法解決N_Body問題：隨著問題規(guī)模的上升，無論是否增加處理器的數量，計算的性能都按照O(N2)的比例下降在大的分子模擬或天體物理學計算中，N的數量會很大，104或者更多怎樣才能夠提高實際應用問題的計算性能？N_Body問題一種優(yōu)化并行算法Barnes-Hut算法：采用分治策略、進行近似計算基本思想：對任何一個粒子A，當一群粒子與A的距離“足夠遠”，這群粒子對A的影響可被一個“個體”來(近似)代替。在牛頓力學中，任何粒子B對A的作用與r2成反比，r是A與B間的距離對于一組粒子B1、B2、…、Bk，如果它們在空間位置上位于一個邊長為d的立方體范圍內，該立方體中心與A的距離為r。只要r“足夠”大、d“足夠”小，那么在計算B1、B2、…、Bk對A的影響時，可以把它們看作單個“個體”B來，用B對A的影響近似代替B1、B2、…、Bk對A的影響關鍵問題是：什么叫“足夠遠”、如何確定“一群”、用哪“一個體”來取代？N_Body問題一種優(yōu)化并行算法：實現思路建樹：確定一個包含所涉及空間的立方體；用一棵8叉樹來表示對該空間的如下劃分，把空間中的粒子分布的區(qū)域性和層次性同時表現出來。該立方體為根，如果其中只有一個粒子，停止；否則按照自然的方式，將它劃分為8個相等的子立方體，如果其中有粒子，則用一個子節(jié)點代表。分別以子節(jié)點為根，重復上述過程。結果：每個葉節(jié)點代表一個粒子，每個粒子由唯一一個葉節(jié)點代表。內節(jié)點代表一個空間單元?？梢钥紤]為節(jié)點實現為一個數據結構節(jié)點所代表立方體的質心的屬性，包括空間位移、質量、速度、加速度等節(jié)點所代表立方體的邊長d如果粒子分布在三維空間中，得到一棵8叉樹如果粒子分布在平面空間中，得到一棵4叉樹N_Body問題一種優(yōu)化并行算法：實現思路計算一個粒子(葉節(jié)點)所受的力：從根開始做樹的遍歷，如果一個空間單元的質心距離本粒子足夠遠，則用在該質心的一個等價體來計算，不再考察空間單元下面的子樹。足夠遠的標準：設一個立方體空間的邊長為d，本粒子到它的質心的距離為r若則可以用該立方體包含星體的總質量和質心來計算，不需要再考慮下面的個別星體。其中，選在0.5—1.2之間，它越小，意味著近似的精度越高。這個式子也表達了“距離越遠，能被近似的空間越大”的含義。

注意不同粒子的計算量不同對一個粒子而言，在遍歷樹的時候，“鄰近”的分枝會遍歷較深。由于樹的平均高度為logn，計算一個星體的受力平均也就是logn。于是整個一個時間步的計算就是大約nlogn。N_Body問題一種優(yōu)化并行算法：并行性分析采用Barnes-Hut算法解決N_Body問題時，在每個時間步上，都要依次執(zhí)行下列四個super-step建樹計算數上各個內節(jié)點的參數(質心、質量等)遍歷樹，計算各個粒子所受到的力更新粒子的屬性

實際上，每個super-step都是可以并行的建樹：不同的子空間可由不同的處理器負責向下分解，獨立的操作計算數上各個內節(jié)點的參數：每個內節(jié)點的參數計算都完全獨立遍歷樹，計算各個粒子所受到的力：以粒子為單位，完全獨立的更新粒子的屬性：以粒子為單位，完全獨立的從各super-step的時間開銷比例來看，計算各個粒子所受到的力是主要的時間開銷N_Body問題一種優(yōu)化并行算法：實現的難點涉及兩個數據結構粒子，包括每個粒子的空間位置、質量、速度、加速度等屬性通過建樹過程得到的8(或4)叉樹：在每個時間步上，得到的樹不同如何劃分這兩個數據結構、如何分配計算任務建樹過程中，自然的想法是：相鄰的粒子位于同一個處理器上。但粒子是運動的，每個時間步上，粒子間的鄰接關系會改變計算任何一個粒子所受的力時，都需要從數的根開始遍歷，如何劃分樹？每個處理器上都保存一棵完整的樹嗎？計算各個粒子所受到的力時，自然的想法是每個處理器負責一組粒子的作用力計算，但是每個粒子作用力的計算復雜性不同在不同的時間步上，同一個粒子作用力的計算復雜性不同更新粒子的屬性時，最好的分配辦法是：把粒子以“block”方式均勻劃分給各個處理器請考慮考慮數據劃分的最簡單辦法是：每個處理器都保留這兩個完整的數據結構但是，粒子的數量越大，樹中節(jié)點的數量也越多，這樣損傷了對問題規(guī)模的scalability比較二維FFT、一維FFT和FOX算法三者都可以用BSP模型刻畫，而且super-step數量與被計算的數據值無關二維FFT可以表示成兩個super-step一維FFT可以劃分成log(p)個super-stepFOX算法中super-step取決于并行程序的進程拓撲結構Gauss-Seidel迭代法一種線性方程組的數值解法，比Jacobi迭代法的收斂速度快什么是Gauss-Seidel迭代法：AX=B表示一個方程組A是mn的矩陣，且

A(i,i)0X是長度n的未知向量B是長度為m的已知向量令X(0)=(c0,c1,…,cm-1)當max(|X[i](t)-X[i](t-1)|)<時，X=X(t)Gauss-Seidel迭代法的特征Gauss-Seidel迭代法的并行性分析與Jacobi迭代法不同的是，在每個時間步上，X[0](t)、X[1](t)、…、X[m-1](t)的計算必須是順序執(zhí)行的下列計算可以并行執(zhí)行X[i](t)=(b[i]-left(i)(t)-right(i)(t-1))a[i,i]對于j>i，left(j)(t)=left(j)(t)+a[j,i-1]X[i-1](t)對于j<i，right(j)(t)=right(j)(t)+a[j,i-1]X[i-1](t)并行算法進程的拓撲結構是一個線性序列A按照(block,*)方式劃分，left、right、B按照“block”方式劃分X(t)按照“block”方式劃分通信：每個時間步上，每個處理器執(zhí)行p次廣播MPI_Bcast其中一次是把局部的X(t)片段廣播給其他進程p次用于接收其它處理器廣播的X(t)片段流水并行(pipelinedcomputation)什么是流水并行(pipelinedcomputation)把一個問題分解成一組功能相同的子任務task，每個子任務處理的數據不同，這些子任務需要依次執(zhí)行把每個子任務的執(zhí)行劃分成多個階段stage每個進程負責實現一個stageP0stage0P1stage1P2stage2P3stage3P4stage4data0data1data2data3data4data5data6data0data1data2data3data4data5data0data1data2data3data4data0data1data2data3data0data1data2時間Gauss-Seidel迭代的流水并行什么是Gauss-Seidel迭代法：AX=B表示一個方程組A是mn的矩陣，且

A(i,i)0X是長度n的未知向量B是長度為m的已知向量令X(0)=(c0,c1,…,cm-1)當max(|X[i](t)-X[i](t-1)|)<時，X=X(t)X[0](t)、X[1](t)、…、X[m-1](t)的計算順序執(zhí)行下列計算并行執(zhí)行X[i](t)=(b[i]-left(i)(t)-right(i)(t-1))a[i,i]對于j>i，left(j)(t)=left(j)(t)+a[j,i-1]X[i-1](t)對于j<i，right(j)(t)=right(j)(t)+a[j,i-1]X[i-1](t)并行算法進程的拓撲結構是一個線性序列A按照(block,*)方式劃分，left、right、B按照“block”方式劃分X(t)按照“block”方式劃分第一步第二步第三步第四步第五步第六步第七步第八步規(guī)則計算與非規(guī)則計算對并行計算而言，我們可以這樣理解(注意：我沒有試圖給一個精確的定義)規(guī)則計算：訪問異地數據時，不需要使用間接地址劃分子任務的數據時，不需要使用數據劃分的索引非規(guī)則計算：不是規(guī)則計算的其他計算非規(guī)則計算的例子：線性稀疏方程求解一個大的稀疏方程，為了提高性能，常使用壓縮存儲的辦法避免“0”元素的乘法：減少Jacobi迭代法或Gauss-Seidel迭代法中乘法運算量避免“0”元素占用的存儲空間：提高數據訪問效率例如用一個一維的數組A存儲系數矩陣中的非“0”值，用另一個與A相同長度的數組B作為地址索引，B[i]記錄A[i]在系數矩陣中的行號和列號采用前面的Jacobi迭代法或Gauss-Seidel迭代法并行算法計算，在劃分方程組的系數矩陣時，需要根據B的值劃分A的數據，B的作用是訪問A的索引為什么會出現非規(guī)則計算的問題原因一：問題空間本身具有不規(guī)則性，通常表現為：有一組數據并行的子任務它們分別對不同的數據集合執(zhí)行相同的運算但每個數據集合的規(guī)模不同例如：非連續(xù)彈簧體問題(multi-elasticbodysystem)每個彈簧體被離散化成一組節(jié)點，由于彈簧體的大小不同，節(jié)點的數量不同根據各個彈簧體邊界上的受力，計算彈簧體的各個節(jié)點的移動情況12345678原因二：為了性能優(yōu)化例如：稀疏線性方程的求解性能優(yōu)化時常需要用到領域知識例如：對于分塊矩陣形式的稀疏方程，我們可以把它劃分成若干個非稀疏矩陣的小方程對于一般的稀疏矩陣，則必須采用索引的方法非規(guī)則計算的表現用兩個不同的數組(數據文件)描述原始數據索引存儲原始數據的數組/文件：一組定長的記錄稀疏矩陣的一個非零元素：同一行/列的非零元素連續(xù)存放多彈簧體中，彈簧體上一個“剛點”：同一彈簧體的“剛點”連續(xù)存放存儲索引的數組/文件：一組定長的記錄稀疏矩陣：每一行/列用一個記錄，這一行在原始數據(的數組/文件)中起始記錄號有多少個記錄多彈簧體：每個彈簧體用一個記錄，這個彈簧體原始數據(的數組/文件)中起始記錄號有多少個記錄