專利名稱：一種基于迭代空間條塊的并行有限差分模版方法
技術領域：
本發(fā)明屬于工程計算領域，涉及一種基于迭代空間條塊的并行有限差分模版方法。
背景技術：
隨著科學計算需要解決的問題漸趨復雜和并行計算機的迅速發(fā)展，并行計算已成為解決大規(guī)模并行科學計算問題必不可少的手段。在計算數(shù)學和計算物理等科學與工程計算領域，很多問題最終都歸結為求解稀疏線性代數(shù)方程組。因此，在并行計算機上高效求解大規(guī)模稀疏線性代數(shù)方程組已成為當前科學計算領域的一項重要任務。由于受舍入誤差、計算機內存和計算復雜度的限制，對大規(guī)模問題，直接求解該類方程組幾乎是不可能的，通常采用有限差分模版計算方法，迭代法的主要思想是通過構造有效的迭代格式，在有限步數(shù)內收斂于方程的精確解。有限差分模版并行方法的實現(xiàn)一直是數(shù)值計算的重要研究對象。傳統(tǒng)的做法是通過使用基于區(qū)域分解的多色排序方法來實現(xiàn)面向集群的并行GS方法，但是當數(shù)據(jù)量增大時，數(shù)據(jù)局部性成下降趨勢，并且在每次迭代計算過程，都需要通信和同步。另一方面，循環(huán)分塊技術可以提高循環(huán)并行度和數(shù)據(jù)局部性優(yōu)化。為此，研究人員對迭代空間分塊做了很多的研究。這些研究主要集中在兩個方面(1)數(shù)據(jù)局部性優(yōu)化研究；(2)提高并行效率的分塊尺寸/形狀研究。通過數(shù)據(jù)依賴向量實現(xiàn)數(shù)據(jù)分塊使其能夠最大化復用，但是分塊技術會引入數(shù)據(jù)塊的空間依賴性，導致其基于空間條塊的流水線執(zhí)行方式不利于并行執(zhí)行。 但模版迭代方法并行化還有4個問題需要進一步的優(yōu)化
數(shù)據(jù)局部性問題。大多數(shù)方法能夠有效的提高迭代內數(shù)據(jù)局部性，但對迭代間的數(shù)據(jù)局部性優(yōu)化效果有限?？蓴U展性問題。傳統(tǒng)的并行化迭代方法在迭代內和迭代間都需要同步操作以維護數(shù)據(jù)依賴關系。處理機進行全局的同步，會增加開銷時間，當處理機的臺數(shù)增多時，全局同步的代價變得更加重要，并且影響方法的可擴展性。通信和同步開銷問題。由于傳統(tǒng)的并行化方法需要在每次迭代過程中通過通信操作得到邊界數(shù)據(jù)，通信開銷制約了并行方法的效率。并且當問題給定時，隨著處理機臺數(shù)的增大，并行純計算時間在減少，而通信時間在不斷增加，這必將影響并行方法的可擴展性。依賴關系問題。迭代空間劃分所引入的條塊依賴關系嚴重阻礙了循環(huán)的并行化。

發(fā)明內容
本發(fā)明的目的在于鑒于現(xiàn)有方法存在的上述問題，提出基于迭代空間條塊的并行有限模版方法。在不增加數(shù)據(jù)通信量的情況下，通過負超平面對迭代空間進行劃分，減少通信啟動開銷和同步時間，提高數(shù)據(jù)局部性。并且通過對空間迭代塊的重新排序，改變塊間數(shù)據(jù)依賴性，有效地實現(xiàn)迭代方法的并行化。本發(fā)明解決技術問題所采用的技術方案是
本發(fā)明首先通過對迭代空間進行時間軸方向劃分成網(wǎng)格條塊，實現(xiàn)對同一網(wǎng)格塊進行遞歸式多次迭代步更新，從而在不改變串行模版迭代方法性質的同時，提高條塊內數(shù)據(jù)局部性。然后通過對網(wǎng)格條塊重排序，實現(xiàn)交錯條塊有限差分模版方法的并行化。步驟1根據(jù)處理器拓撲進行空間計算區(qū)域劃分
依據(jù)處理器數(shù)目P及處理器的拓撲結構，將空間計算區(qū)域劃分為Pl XP2個子空間 domain，其中 Pl XP2=P 且 P1=P2 ； 步驟2指定迭代方向
交錯條塊迭代方法分為奇數(shù)k次迭代和偶數(shù)k次迭代，設定奇數(shù)k次迭代與偶數(shù)k次迭代執(zhí)行方向相反；
步驟3對子空間進行網(wǎng)格條塊劃分
通過時間軸將各子空間劃分為多個網(wǎng)格條塊，同串行迭代方法，并且奇數(shù)k次和偶數(shù)k 次劃分方向相反；
步驟4對網(wǎng)格條塊重新排序
根據(jù)處理器數(shù)目P按照CydiC(P)排列方法對所有的網(wǎng)格條塊進行重排序； 步驟5以網(wǎng)格條塊為單位執(zhí)行奇數(shù)k次迭代數(shù)據(jù)更新
執(zhí)行奇數(shù)k次迭代更新，更新順序參照重排序之后的網(wǎng)格條塊序；當更新奇數(shù)k次中的發(fā)送條塊和混合條塊后，將邊界數(shù)據(jù)發(fā)送給相應的接收條塊和混合條塊，而接收數(shù)據(jù)的條塊在接收完數(shù)據(jù)后再進行迭代更新；
步驟6以網(wǎng)格條塊為單位執(zhí)行偶數(shù)k次迭代數(shù)據(jù)更新
執(zhí)行偶數(shù)k次迭代更新，更新順序參照重排序之后的網(wǎng)格條塊序；當更新偶數(shù)k次中的發(fā)送條塊和混合條塊后，將邊界數(shù)據(jù)發(fā)送給相應的接收條塊和混合條塊，而接收數(shù)據(jù)的條塊在接收完數(shù)據(jù)后再進行迭代更新。本發(fā)明具有的有益效果是
1、本發(fā)明方法通過網(wǎng)格條塊重排序減少了 cache缺失率、通信啟動和同步次數(shù)。2、本發(fā)明方法比傳統(tǒng)的區(qū)域分解方法和紅黑排序并行方法具有更好的數(shù)據(jù)局部性，并行效率和可擴展性。


圖1為迭代空間交錯條塊串行模版方法，⑴為從1到k次的正向執(zhí)行，(2)為從 k+Ι到業(yè)次的反向執(zhí)行。其中不同線條表示不同時間的迭代范圍。圖2為本發(fā)明方法的系統(tǒng)框圖。圖3為三維迭代空間的兩類條塊劃分方法，(1)為對稱遍歷條塊串行執(zhí)行的前半部分對稱遍歷，其中迭代步是從第1到k次迭代，(2)為對稱遍歷條塊串行執(zhí)行的后半部分對稱遍歷，其中迭代步是從第k+l到業(yè)次迭代，(3)為對稱遍歷條塊重排后的前半部分對稱遍歷，其中迭代步是從第1到k次迭代，(4)為對稱遍歷條塊重排后后半部分對稱遍歷其中迭代步是從第k+l到業(yè)次迭代。圖4為條塊超平面劃分俯視圖。圖5為圖3所示兩類條塊劃分方法的條塊依賴圖，(1)為圖3(1)和圖3(2)中劃分方法的條塊依賴圖，( 為圖3C3)和圖3(4)中劃分方法的條塊依賴圖。其中數(shù)字前面的標號表示不同的迭代方向f表示奇數(shù)次迭代，b表示為偶數(shù)次迭代。
圖6為四個子區(qū)域的并行化，(1)為對稱遍歷過程的正向部分，(2)為對稱遍歷過程反向部分。
具體實施例方式下面結合附圖和實施例，對本發(fā)明的具體實施方案作進一步詳細描述 傳統(tǒng)的迭代方法執(zhí)行序過程如下
在一次迭代內部依據(jù)網(wǎng)格點的順序，依次對所有網(wǎng)格點進行迭代更新操作 for ( t=l; t<=T; t++ ) for ( i=l; i<N-l; i++ ) for ( j=l; j<N-l; j++ )
A[t+1] [i] [j] = (A[t+l] [i-1] [j]+A[t] [i+1] [j]+A[t] [i] [j] + A[t+l][i][j-l] + A[t][i][j+l])/5
其中T表示迭代次數(shù)，NXN表示矩陣的大小。在上述方法中，數(shù)據(jù)復用的方式共有以下四種
行優(yōu)先存儲格式下，在j循環(huán)中復用A[i] [j]，A[i] [j+1]…… 在j循環(huán)中，更新A[i] [j-1]、A[i] [j]和A[i] [j+1]過程中復用A[i] [j]。在i 循環(huán)中，更新 A[i_l] [j]、A[i] [j]和 A[i+1] [j]過程中復用 A[i] [j]。在時間步中，復用A[i][j]。由于數(shù)據(jù)存取位置相近，因此出現(xiàn)第一類和第二類數(shù)據(jù)復用。當cache的容量C 大于2N (N為一維邊界長度)時，A[i+l][j]與A[i-l][j]在行優(yōu)先存儲格式下相鄰2N個數(shù)據(jù)，因此會出現(xiàn)第三類數(shù)據(jù)復用。當N2小于C時，所有數(shù)據(jù)都會一次性讀取到cache中， 因此會出現(xiàn)第四類數(shù)據(jù)復用即提高了數(shù)據(jù)空間局部性。但這種“串行”執(zhí)行序不僅導致其數(shù)據(jù)依賴關系阻礙了迭代方法并行化執(zhí)行，而且當數(shù)據(jù)量大時，其數(shù)據(jù)局部性較差。其原因在于當數(shù)組大于cache容量時，本次更新的數(shù)據(jù)在下一次更新之前已經(jīng)被寫回內存。而且當數(shù)據(jù)量大時，傳統(tǒng)迭代方法中多次迭代會導致數(shù)據(jù)cache的容量缺失。此外，多次迭代使數(shù)據(jù)的地址轉換信息會周期性地在TLB中進行存取，也是影響其性能的一個重要因素。當矩陣規(guī)模增加時，大量的TLB容量缺失會嚴重影響程序的性能。因此，本發(fā)明采取迭代空間條塊串行迭代執(zhí)行序，其核心思想是改變以往傳統(tǒng)的以迭代次序為執(zhí)行序的特點，將迭代空間分塊引入執(zhí)行序中，如圖1所示。本發(fā)明還通過網(wǎng)格條塊重排序，來實現(xiàn)交錯條塊有限差分模版方法的并行化。其具體步驟描述如圖2所示
步驟1根據(jù)處理器拓撲進行空間計算區(qū)域劃分
依據(jù)處理器數(shù)目P及處理器的拓撲結構，將空間計算區(qū)域劃分為Pl XP2個子空間 i/offlai/ ，其中 Ρ1ΧΡ2=Ρ 且 Ρ1=Ρ2，如圖 3(3)禾Π圖 3 )所示； 步驟2指定迭代方向
交錯條塊迭代方法分為奇數(shù)k次迭代和偶數(shù)k次迭代，設定奇數(shù)k次迭代與偶數(shù)k次迭代執(zhí)行方向相反；例如，使用LB-RT順序(左下到右上)執(zhí)行奇數(shù)k次迭代，使用RT-LB 順序(右上到左下)執(zhí)行偶數(shù)k次迭代。步驟3對子空間進行網(wǎng)格條塊劃分將各子空間劃分為多個網(wǎng)格條塊，同串行迭代方法，并且奇數(shù)k次和偶數(shù)k次劃分方向相反；
以二維Possion方程為例，迭代空間交錯條塊串行迭代方法執(zhí)行過程如下 1)網(wǎng)格條塊劃分
m維網(wǎng)格點x(il，…，im)組成m維數(shù)據(jù)空間data_space (m)。 n_l維數(shù)據(jù)空間data_Space (n-1)和迭代維T的組合可以被看作η維的迭代空間iter_ space (I1，…，Ilri，Τ)，其中的每個點都可由一個ηΧ 1維的列向量來表示，即
F = Ci1,...,^f，其中.,U從左至右分別代表網(wǎng)格點在數(shù)據(jù)空間中的I1,…，Ilri維的維度
坐標。數(shù)據(jù)空間中網(wǎng)格點Wi1，…，im)在k次迭代的值，在迭代空間中可表示為!^(‘… ，‘)。例如，三維迭代空間iter_SpaCe(I，J，T)由二維網(wǎng)格點x(i，j)和時間維T組成。 uk(i, j)表示網(wǎng)格點x(i，j)在k次迭代的值。用區(qū)域分解方法將迭代空間iter_SpaCe(I，J, Τ)在T=O處進行數(shù)據(jù)劃分。如圖 1所示，定義Λ為子空間Sub_iter_Space(p，q，0)的行數(shù)，定義及為子空間sub_iter_ space (ρ, q, 0)的列數(shù)，劃分后使得每個子空間sub_iter_space (p，q, 0)中的網(wǎng)格點數(shù)為R= I1 χ Λ，A與及滿足式(1)
I1 >k Π I2 >k (1) 其中，k為單向迭代次數(shù)。2)在網(wǎng)格條塊的基礎上，沿時間軸對迭代空間進行劃分
劃分方法采用時滯技術，對每層迭代的子空間修正邊界，圖1中虛邊界表示修正后的邊界。定義有向圖G(V，E)存放相鄰網(wǎng)格條塊的關系.若網(wǎng)格條塊Vi與網(wǎng)格條塊Vj邊界相
連且Vi<\，則<Vi，、> e Ε。定義Pmfe (Vi, Vj, k)為在第k次迭代中屬于條塊Vi但與條塊Vj
相鄰的邊界數(shù)據(jù)。算法描述如下 /*條塊邊界修正算法*/ for (n=0;n<N;n+=2)
for( k= n*K+l ； k<=(n+l)*K; k++) { //奇數(shù) K 次迭代修正邊界 Foreach Ky1, v」> E
sub_iter_space (Vi, k+1) = sub_iter_space (Vi, k) _ Pnode (Vi, Vj, k)
sub_iter_space (v」，k+1) = sub_iter_space (v」，k) + Pnode (Vi, Vj, k)
}
for (k=(n+l)^K+l; k〈=(n+l)*K; k++) { //偶數(shù) K 次迭代修正邊界 Foreach〈ν” ν』> E
sub_iter_space (Vjj k+1) = sub_iter_space (Vij k) - Pnode (Vjj Vij k) sub_iter_space (Vij k+1) = sub_iter_space (Vij k) + Pnode (Vjj Vij k)
ι
3)按空間網(wǎng)格條塊順序執(zhí)行模版迭代方法
網(wǎng)格條塊生成后，按條塊順序執(zhí)行迭代方法。以條塊內部網(wǎng)格點層作為內部執(zhí)行序，條塊內部迭代次數(shù)作為中間執(zhí)行序，條塊序作為外部執(zhí)行序，更新每個網(wǎng)格點的值。執(zhí)行過程如圖1所示，a. 1-a. 4為條塊正向執(zhí)行順序，b. l_b. 4為條塊反向執(zhí)行順序，其中數(shù)據(jù)塊中3種不同的線條框分別代表連續(xù)的3次迭代計算。步驟4對網(wǎng)格條塊重新排序
根據(jù)處理器個數(shù)P按照CyCliC(P)排列對所有的網(wǎng)格條塊進行重排序。傳統(tǒng)的多面體模型通過迭代空間多面體矩陣B，依賴矩陣D，超平面矩陣H，條塊依賴關系矩陣S描述迭代空間的數(shù)據(jù)劃分及依賴關系。但傳統(tǒng)多面體模型中引入超平面執(zhí)行順序會降低迭代空間并行化的效率，因此需要改變超平面方向消除引入的數(shù)據(jù)依賴關系， 提高迭代計算的并行性。在標準條塊迭代方法中超平面由一組法向量組成。法向量定義了條塊在同一超平面方向的兩面，其中Iii (乃定義為I方向第/條塊區(qū)域tile(A I)在I方向超平面的后平面，^(7+1)定義為第/條塊區(qū)域的I方向超平面的前平面。在條塊中前平面不會引入I方向其他條塊的依賴，而后平面中的數(shù)據(jù)在計算時需要I方向相鄰條塊的數(shù)據(jù)，因此引入了 I方向數(shù)據(jù)依賴。例如V <Ρ'·0 , V fik(q,rj€ suhdomain(m)
，且ftKiM)與iAfeJ)相鄰，在計算時，會使用fife(ftZ)的邊界數(shù)據(jù)。在邊界條塊中通過建立I方向的負超平面可以減少條塊之間在I方向的依賴關系。將后平面、(乃替換為負超平面、、滿足
權利要求
1.一種基于迭代空間條塊的并行有限差分模版方法，其特征在于該方法包括如下步驟步驟1根據(jù)處理器拓撲進行空間計算區(qū)域劃分，具體是依據(jù)處理器數(shù)目P及處理器的拓撲結構，將空間計算區(qū)域劃分為P1XP2個子空間SizliZofflai/ ，其中P1XP2=P且P1=P2 ； 步驟2指定迭代方向，具體是交錯條塊迭代方法分為奇數(shù)k次迭代和偶數(shù)k次迭代， 設定奇數(shù)k次迭代與偶數(shù)k次迭代執(zhí)行方向相反；步驟3對子空間進行網(wǎng)格條塊劃分，具體是將各子空間劃分為多個網(wǎng)格條塊，并且奇數(shù)k次和偶數(shù)k次劃分方向相反；步驟4對網(wǎng)格條塊重新排序，具體是根據(jù)處理器數(shù)目P按照cyclic (P)排列方法對所有的網(wǎng)格條塊進行重排序；步驟5以網(wǎng)格條塊為單位執(zhí)行奇數(shù)k次迭代數(shù)據(jù)更新，具體是執(zhí)行奇數(shù)k次迭代更新，更新順序參照重排序之后的網(wǎng)格條塊序；當更新奇數(shù)k次中的發(fā)送條塊和混合條塊后， 將邊界數(shù)據(jù)發(fā)送給相應的接收條塊和混合條塊，而接收數(shù)據(jù)的條塊在接收完數(shù)據(jù)后再進行迭代更新；步驟6以網(wǎng)格條塊為單位執(zhí)行偶數(shù)k次迭代數(shù)據(jù)更新，具體是執(zhí)行偶數(shù)k次迭代更新，更新順序參照重排序之后的網(wǎng)格條塊序；當更新偶數(shù)k次中的發(fā)送條塊和混合條塊后， 將邊界數(shù)據(jù)發(fā)送給相應的接收條塊和混合條塊，而接收數(shù)據(jù)的條塊在接收完數(shù)據(jù)后再進行迭代更新。
全文摘要
本發(fā)明涉及一種基于迭代空間條塊的并行有限差分模版方法。傳統(tǒng)的并行化迭代方法在迭代內和迭代間都需要同步操作以維護數(shù)據(jù)依賴關系。本發(fā)明首先通過對迭代空間進行時間軸方向劃分成網(wǎng)格條塊，實現(xiàn)對同一網(wǎng)格塊進行遞歸式多次迭代步更新，從而在不改變串行模版迭代方法性質的同時，提高條塊內數(shù)據(jù)局部性。然后通過對網(wǎng)格條塊重排序，實現(xiàn)交錯條塊有限差分模版方法的并行化。本發(fā)明方法比傳統(tǒng)的區(qū)域分解方法和紅黑排序并行方法具有更好的數(shù)據(jù)局部性，并行效率和可擴展性。
文檔編號G06F17/00GK102200962SQ20111005128
公開日2011年9月28日 申請日期2011年7月25日 優(yōu)先權日2011年7月25日
發(fā)明者萬健, 任永堅, 張偉, 張紀林, 徐向華, 蔣從鋒 申請人:杭州電子科技大學