一種基于空間密鋪的星型Stencil并行算法

曹 杭 袁 良 黃 珊 張云泉 徐勇軍 陸鵬起 張廣婷

1(計算機體系結(jié)構(gòu)國家重點實驗室(中國科學院計算技術研究所) 北京 100190)

2(中國科學院大學 北京 100049)

Stencil計算(模板計算)是科學工程應用中一類常見的嵌套循環(huán)算法，也被稱為“結(jié)構(gòu)化網(wǎng)格計算”，是13個伯克利核心計算模式之一.在動態(tài)規(guī)劃、圖像處理等領域的很多算法中也包含相似的依賴模式.Stencil定義了一種更新結(jié)點所需鄰居結(jié)點的模式.Stencil計算在時間維度上迭代更新規(guī)則的d維網(wǎng)格(也稱為數(shù)據(jù)空間)，數(shù)據(jù)空間沿著時間維度更新，產(chǎn)生的d+1維空間被稱為迭代空間.

與其他的數(shù)值算法不同，例如稠密或稀疏矩陣計算類型的多樣性主要來源于輸入數(shù)據(jù)如稠密或稀疏矩陣的數(shù)據(jù)量,Stencil的計算類型本質(zhì)上取決于問題的類型.Stencil可以從不同的角度進行分類，如格點維度(1維，2維，…)，鄰居數(shù)目，即階數(shù)(3點，5點，…)，形狀(盒型，星型，…)，依賴類型(高斯-賽德爾,雅可比,…)和邊界條件(常量的，周期的，…)等.

Stencil計算具有計算密集度低的特征，d維Stencil的簡單實現(xiàn)由d+1個循環(huán)組成，其中最外層循環(huán)遍歷時間維度，內(nèi)層循環(huán)更新d維空間的所有格點.簡單實現(xiàn)的數(shù)據(jù)重用度較差且典型的帶寬受限.

分塊是提高多重嵌套循環(huán)的數(shù)據(jù)局部性和并行度最有效的優(yōu)化技術之一，目前有大量針對Stencil計算分塊方案的研究.空間分塊算法促進了2維至更高維度Stencil在一個時間步內(nèi)的數(shù)據(jù)重用.為了進一步增強數(shù)據(jù)重用度，往往同時考慮時間維度和空間維度，后文將詳細介紹這些技術.另一類工作采用簡單的超矩形分塊形狀[1-4]并引入冗余計算來解決塊間數(shù)據(jù)依賴.這些研究包括自動調(diào)優(yōu)框架[5-12]、編程模型[13-14]、CPU[2,15-18]、GPU[19-23]和Phi[24-25]上的手動調(diào)優(yōu)實現(xiàn).

不同的Stencil類型也衍生出不同的解決方案.Pluto系統(tǒng)[26]延伸到周期性邊界條件Stencil的處理[27].菱形分塊最初是針對1維Stencil的分塊[28]，然后擴展到更高維度的Stencil[29].高速緩存參數(shù)無關(cache oblivious)方案從最初的算法[30-31]，CORALS(cache oblivious parallelograms)方案[32]，發(fā)展到能處理任意Stencil并達到最大并發(fā)度的Pochoir系統(tǒng)[33].

據(jù)我們所知，目前沒有針對不同Stencil形狀的解決方案，特別是盒型和星型Stencil.星型Stencil只依賴于每個坐標軸上的點，盒型Stencil在此之外還有對角線上的數(shù)據(jù)依賴.不區(qū)分盒型和星型Stencil的根本原因可能是盒型Stencil包含星型Stencil，任何適用于盒型Stencil的算法也能應用到星型Stencil.

本文首先在數(shù)據(jù)空間上提出“自然塊”的概念來確定星型和盒型Stencil之間的區(qū)別.然后針對星型Stencil提出了一個新的2層密鋪方案，在此方案中，自然塊和它的后繼塊可以密鋪整個空間，而他們沿著時間維度擴展后形成對迭代空間的密鋪.密鋪框架類似于文獻[34]中提出的方法.我們將統(tǒng)一闡述基于自然塊概念的2種方案并將已有的方案稱為盒型密鋪而將新提出的方案稱為星型密鋪.

在密鋪方案中區(qū)別盒型和星型Stencil有2個主要的優(yōu)點，分別利用不同層次上的數(shù)據(jù)局部性.首先，當給定緩存大小時，相比盒型密鋪，用星型Stencil的自然塊密鋪能得到更大的塊，這樣可以更有效地減少內(nèi)存系統(tǒng)的壓力；其次，我們針對第二后繼塊中的元素提出了一個新穎的“2次更新”技術，可對1個元素連續(xù)進行2次更新，進一步提升了核心數(shù)據(jù)的重用.

1 相關工作

1.1 分塊技術

分塊技術[35-39]是探究多層嵌套循環(huán)數(shù)據(jù)局部性和并行度最有效的轉(zhuǎn)換技術之一.Wonnacott和Strout[40]對比了一些現(xiàn)有分塊方案的可擴展性，下面我們將總結(jié)一些和Stencil計算有關的技術.

1) 重疊分塊(overlapped tiling).高性能領域中經(jīng)常在手動調(diào)優(yōu)的實現(xiàn)中采用超矩形分塊[1-4].為了執(zhí)行多個時間步，往往采用冗余計算[41]來解決分塊之間的依賴，這就是重疊分塊[28,42].Philips和Fatica[22]提出了GPU上手動調(diào)優(yōu)的手寫3.5維分塊實現(xiàn)方案，在2.5維分塊方案[4]上增加了時間維度上的劃分.Demmel等人[43-44]減少了2倍的冗余計算.形狀規(guī)則的超矩形可以獲得更高的并發(fā)度和更多細粒度優(yōu)化.

2) 時間偏移(time skewing).時間偏移分塊[45-47](分塊形狀在2維上是平行四邊形，在3維上是平行六面體，在更高維上是超平行體)消除了冗余計算，但是大多數(shù)方法只用單一形狀的塊填充空間，導致流水線啟動和有限的并發(fā)度[48-49].

3) 菱形分塊(diamond tiling).Bondhugula等人[26]對1維Stencil采用菱形分塊.Bandishti等人[29,50]將此方法擴展到更高空間維度的Stencil.Grosser等人[51]對2維的菱形尖端粗粒度化形成六邊形，在3維上演化成截掉頂端的八面體.

4) 緩存無關分塊(cache oblivious tiling).緩存無關技術旨在內(nèi)存層次結(jié)構(gòu)參數(shù)未知時充分利用數(shù)據(jù)局部性.Frigo和Strumpen提出了第1個串行[30]和并行[31]緩存無關Stencil算法.緩存無關平行四邊形方法[32]同時分離時間空間維度，但是會導致波陣面并行.Pochior[33]實現(xiàn)了多維空間劃分，能夠同時劃分所有空間維度.

5) 分裂分塊(split tiling).分裂分塊方法發(fā)掘每個塊中的獨立子塊，然后將這些子塊發(fā)送給他們的后繼塊，使得剩下的區(qū)域同樣可以并發(fā)執(zhí)行[28,52-54].Grosser等人[55]提出了另一個類似于具有緩存無關范式的多維空間劃分的分裂分塊方法.嵌套分裂分塊(the nested split-tiling)[56]能在所有空間維度遞歸分裂.

6) 混合分塊(hybrid tiling).Strzodka等人[57]提出的CATS(cache accurate time skewing)算法結(jié)合了菱形分塊、平行四邊形分塊和流水線處理.Grosser等人[58]采用混合六邊形和平行四邊形的分塊算法.這些算法取某一維空間和時間維度組成平面，用六邊形分塊和菱形分塊進行劃分，在剩下的空間維度中采用時間偏移分塊，從而分解迭代空間.混合分裂分塊方法[56]結(jié)合了嵌套分裂分塊和時間偏移分塊.

7) 密鋪分塊(tessellating).密鋪方案[34]適用于包括星型和盒型等不同Stencil類型.空間被一系列分塊密鋪，分塊之間可以無冗余并行執(zhí)行.相應地，這些分塊在時間維度上擴展后得到的擴展塊可以在迭代空間中形成密鋪.簡單明了的數(shù)學特性使得這種方法可以與細粒度優(yōu)化方法相結(jié)合.

1.2 優(yōu)化方法

為了使核心程序達到更高的性能，核心內(nèi)的優(yōu)化技術如計算重用[59-60]、計算重組[61]、向量化[15,62-64]和數(shù)據(jù)布局轉(zhuǎn)換[65-66]等十分關鍵.但是這些方法大多沒有對不同的Stencil形狀進行區(qū)分.本文設計實現(xiàn)了一個新穎的技術來提升星型Stencil核心程序級別上的數(shù)據(jù)重用.這種技術也可以應用到其他的分塊方法上.

2 算法描述

本節(jié)首先介紹一個新的概念，即自然塊，來區(qū)分不同的Stencil形狀特征.然后圍繞自然塊和后繼塊這2個方面，詳細描述新提出的星型Stencil密鋪方案以及已有的盒型密鋪方案[34].最后將密鋪方案應用于高階Stencil和不同的邊界條件.

Gauss-Seidel Stencil需要執(zhí)行逐點45°流水線啟動，使用分塊技術也不能使其完全并發(fā)啟動[28].因此我們僅考慮Jacobi Stencil.此外，星型Stencil只能應用到2維數(shù)據(jù)空間中，我們只討論2維Stencil.在第5節(jié)的實驗中，對于3維星型Stencil將保留1維數(shù)據(jù)空間不進行切分.

密鋪3維迭代空間對應于計算2維Stencil，我們首先將其劃分成多個相同的時間片.在一個時間片開始時所有格點處于同一時間維度，結(jié)束時所有格點都更新了t個時間步(在下面的例子中t=4).不失一般性，我們只討論第1個時間片的密鋪，即密鋪前所有格點的時間維度均為0，密鋪后所有格點的時間維度為t.

Fig. 1 Natural block B0(0)圖1 自然塊B0(0)

2.1 自然塊B0(0)

當所有格點在時間維度上的值相同時，更新某一格點t個時間步所需的最小鄰居點集被稱為時間步為t的自然塊，簡稱自然塊.

圖1(a)和圖1(b)分別展示了2D9P盒型Stencil和2D5P星型Stencil時間步長為4的自然塊.

最大更新方法[34]與自然塊概念相近.具體定義為：在數(shù)據(jù)空間中給定B，沿著時間維度最大更新每個格點，直到不滿足Stencil定義的依賴關系為止.形式上而言，對所有格點b∈B，滿足以下2個條件中任一個時，停止更新：

1)time[b]=t，t是給定的最大時間更新步數(shù)；

2) 存在至少一個b的鄰居b′，b′的時間維度嚴格小于b(time[b′]

但是這2種概念并不完全等價，實際上它們互為補充.一方面，最大更新沒有給出密鋪數(shù)據(jù)空間塊的具體形狀.盒型密鋪[34]采用盒型Stencil的自然塊(超立方體)來密鋪盒型和星型Stencil，盒型Stencil包含星型Stencil，因此可以保證結(jié)果的正確性.在本工作中，我們將采用不同類型Stencil各自的自然塊.

2.2 用B0(0),B1(1)和B2(2)密鋪

盒型Stencil自然塊B0的形狀是正方形，星型Stencil自然塊的形狀是菱形，這2種形狀都能密鋪2維數(shù)據(jù)空間.圖2(a)和圖2(b)(忽略灰色區(qū)域)分別表示盒型Stencil的4個自然塊和星型Stencil的12個自然塊.注意這些自然塊可以并行執(zhí)行.

Fig. 2 Tessellating data space with B0圖2 用B0密鋪數(shù)據(jù)空間

(1)

Fig. 3 Determining the center points of B1圖3 由0確定B1塊的中心點

后繼塊B1,B2等是對自然塊B0進行分裂而后重組格點所產(chǎn)生的.分裂和重組的關鍵在于確定這些B1的中心點，之后每個格點可以簡單地劃分給離它最近的中心點.這種方法保證除邊界上的點外，每個格點只屬于某一個B1，并且B1可以密鋪數(shù)據(jù)空間.

自然塊的一個性質(zhì)是中心點在時間維度更新步數(shù)最多，因此最顯然的處理方式是將更新次數(shù)為0的格點凸集的中心點作為后繼塊的中心點.需要注意的是，后繼塊的中心點總在前驅(qū)塊的邊界上，只用考慮邊界上標為0的格點凸集.

由圖3可知，盒型和星型Stencil中所有標為0的格點分別在正方形和菱形的邊界上.根據(jù)凸集的要求，4條邊被視為4個凸集，因此方框中的4個格點0即為對應4個B1的中心點.

圖4(a)和圖4(b)中展示了由B0到B1密鋪數(shù)據(jù)空間的轉(zhuǎn)變過程，其中虛線代表對B0的劃分，可以看到每個B0被劃分成4個子塊.圖4(c)和圖4(d)展示了用B1密鋪數(shù)據(jù)空間的方法，其中共享同一個中心點的B0子塊合并組成1個B1.經(jīng)過此轉(zhuǎn)變過程，盒型Stencil由4個B0轉(zhuǎn)換為12個B1，星型Stencil由12個B0轉(zhuǎn)換為16個B1.

Fig. 4 Tessellating data space with B1圖4 用B1塊密鋪數(shù)據(jù)空間

(2)

Fig. 5 1圖5 1

Fig. 6 Determining the center points of B2圖6 由0和1確定B2的中心點

通過相似的步驟確定圖6(a)和圖6(b)中所示的2個0點是B2的中心點.從B1到B2的轉(zhuǎn)換過程如圖7(a)和圖7(b)中虛線所示，由B2密鋪數(shù)據(jù)空間的方式如圖7(c)和圖7(d)所示.盒型Stencil的密鋪方式由12個B1轉(zhuǎn)換為9個B2，星型Stencil由16個B1轉(zhuǎn)換為13個B2.

Fig. 7 Tessellating data space with B2圖7 用B2塊密鋪數(shù)據(jù)空間

(3)

Fig. 8 Natural block B2(2)圖8 自然塊B2(2)

2.3 高階Stencil和邊界條件

解決周期性邊界條件問題時，我們可能需要在1維Stencil中對1個塊進行拉伸，或者在高維Stencil中拉伸1組塊.采用拉伸過的六邊形塊來保證在目標維度2邊的2個子塊可以組成1個自然塊.

3 實現(xiàn)方法

基于第2節(jié)描述的數(shù)學框架可以手工實現(xiàn)代碼.本節(jié)將介紹手工實現(xiàn)代碼的具體細節(jié).自動代碼生成工具的實現(xiàn)將作為未來的工作.

3.1 粗粒度化

Fig. 9 Coarsening star tessellation圖9 粗粒度化星型密鋪

如圖9所示，B0和B2塊在數(shù)據(jù)空間中具有同樣的形狀.因此可以將相鄰時間片的B0和B2塊合并，減少一次同步并增加數(shù)據(jù)重用.

此外，從圖9中可以看出，2個B0在某一維度的距離等于B0結(jié)束區(qū)域的大小，這使得B0的開始區(qū)域和B2的結(jié)束區(qū)域完全相同.

3.4 并行化和向量化

在同一階段的所有塊可以并發(fā)執(zhí)行，因此很容易實現(xiàn)并行化密鋪方案.因此，我們只在共享內(nèi)存機器上簡單使用OpenMP編譯指示parallel for.對于分布式內(nèi)存計算機，清晰的密鋪方案使得我們可以實現(xiàn)一個簡單的數(shù)據(jù)/計算分布以及高效的數(shù)據(jù)通信,留做下一步工作.

星型密鋪方案只適用于2維Stencil，而且B1的計算很難向量化，因此我們只對此方案在3維星型Stencil上的計算性能進行評估.因為3維Stencil中每個維度上的大小通常小于1 000，我們選擇不劃分單位跨度維度來利用硬件預取的能力,采用pragma simd對最內(nèi)層循環(huán)進行向量化.

3.5 代碼實現(xiàn)

代碼實現(xiàn)主體包含5層循環(huán)，如算法1所示.最外層循環(huán)控制時間維度上時間片的切換，與分塊相關的時間片大小為tB，t0即為當前時間片的開始時間.內(nèi)部包含2個4層循環(huán)，第1個循環(huán)用來進行B0塊和B2塊的更新，第2個循環(huán)更新B1塊.

第2個循環(huán)的第1層循環(huán)遍歷所有B1塊，在星型密鋪中，B1塊是沿2個不同的對角線方向合并產(chǎn)生的，可以按照對角線方向的不同將B1塊進一步細分為B11塊和B12塊，對應的dy值分別為1和-1，同時根據(jù)B1塊的位置可以確定對應數(shù)據(jù)空間區(qū)域邊界xmin,ymax.第2層循環(huán)遍歷B1所在的時間維度，跨度為tB，并根據(jù)t計算特定時間維度上對角線區(qū)域邊界xmax,ymin.第3層循環(huán)遍歷對角線區(qū)域中每條對角線，確定每條對角線的起始坐標x,y和結(jié)束橫坐標xr.第4層循環(huán)遍歷某一條對角線上的所有格點，每個格點被連續(xù)更新2次.

算法1.2D星型Stencil算法.

輸入：時間步長tB;迭代空間上下界0,T;數(shù)據(jù)空間菱形邊界xmin,xmax,ymin,ymax;數(shù)據(jù)空間(t,x,y);

輸出：更新后的數(shù)據(jù)空間(t,x,y).

① fort0= -tBtoTsteptBdo

② #pragma omp parallel for private(xmin,

xmax,ymin,ymax,t,x,y)

③ forn=0 toB0orB2do

④ fort=max(t0,0) to min(t0+2×tB,

T) do

⑤set(xmin,xmax,ymin,ymax)；

⑥update(t,x,y)；

/*在邊界內(nèi)更新B0和B2*/

⑦ end for

⑧ end for

⑨ #pragma omp parallel for private(xmin,

xmax,ymin,ymax,t,x,y,xr,dy,i)

⑩ forn=0 toB11orB12do

4 實驗結(jié)果

4.1 實驗環(huán)境

Fig. 10 Performance of size 2563圖10 問題規(guī)模為2563的性能

我們的實驗運行在Intel Xeon E5-2670處理器上，處理器時鐘頻率為2.70 GHz，實驗規(guī)模由單核擴展到所有12核.每個核獨享32 KB L1 cache和256 KB L2 cache，12個核共享1個統(tǒng)一的30 MB L3 cache.使用版本號為16.0.1的ICC編譯器編譯程序并使用優(yōu)化標志‘-O3 -openmp’.

將我們的方案與1.1節(jié)提及的2種高并發(fā)方案Pluto和Girih進行比較.盒型密鋪、Pochoir和SDSL(Stencil domain specific language)[56]因其較差的性能結(jié)果(特別是3維Stencil)不作考慮.測試用例包括3種3維星型Stencil：1階(3D7P)、2階(3D13P)和4階(3D25P)的具有非周期性邊界條件的情況.測試的問題規(guī)模為2563×500，3843×800，5123×1 000.

Girih提供一種自動調(diào)優(yōu)組件來選擇不同線程數(shù)時的最優(yōu)分塊大小.對于Pluto和星型密鋪方案，我們將測試所有可能的分塊大小，然后選擇12核時性能最優(yōu)的分塊大小并應用于所有1～12核測試中.

Fig. 11 Performance of size 3843圖11 問題規(guī)模為3843的性能

4.2 結(jié)果分析

圖10～12分別顯示了星型密鋪(深色柱)、Girih(白色柱)和Pluto(淺色柱)在3D7P，3D13P，3D25P星型Stencil上大小為2563,3843,5123時的性能.其中性能GStencil的定義為：以最內(nèi)層循環(huán)計算量看為一個Stencil，單位時間完成給定的Stencil迭代更新計算的次數(shù)為GStencil.由圖10～12中可以明顯看出，Pluto對應的性能條柱性能都比我們的代碼和Girih低，且其結(jié)果不是一條直線，特別是在6核或更多核時.這和文獻[29]中的結(jié)果一致.因此，下面主要分析星型密鋪方案和Girih.

Fig. 12 Performance of size 5123圖12 問題規(guī)模為5123的性能

在問題規(guī)模大小為2563的3D7P星型Stencil中，星型密鋪方案性能比Girih高，相比于Girih，星型密鋪方案在2核時最高加速1.24倍，平均加速1.19倍.對于圖10中更高階的Stencil類型，我們的方案仍保持較高的性能，但隨著階數(shù)增加，這3種方法對應的性能柱變得更加接近，因為對于給定的自然塊大小，更高階的Stencil會產(chǎn)生更小的時間維度分塊，造成更大的內(nèi)存數(shù)據(jù)傳輸，后面將對此進行定量分析.

類似地，問題規(guī)模變大時3條性能曲線也會相互接近，也是因為Pluto和我們的密鋪方案不對單位跨度的維度進行分塊，導致時間維度分塊大小受限.我們的方案在核數(shù)較少時仍取得了性能提升，在8～12核時雖然性能提升速度比Girih慢，但仍達到了可相比的性能.例如，在大小為3843的3D13P星型Stencil的實驗結(jié)果中，星型密鋪和Girih在12核上的性能分別為4.2,4.4 GStencil.

Fig. 13 Memory transfers of size 2563圖13 問題規(guī)模為2563的數(shù)據(jù)傳輸量

圖13～15中顯示了對應的最后一級cache和內(nèi)存之間數(shù)據(jù)傳輸量.Girih利用整個L3 cache來最大限度開發(fā)數(shù)據(jù)局部性，因此具有最小的數(shù)據(jù)傳輸量.Girih在9個測試中的數(shù)據(jù)傳輸量都不規(guī)則變化，與密鋪或Pluto的變化曲線不同，因為自動調(diào)優(yōu)方案在不同核數(shù)時確定的分塊大小不一樣，例如，大小為2563的3D7P Stencil的Girih測量值在4核上是10.2 GBps，6核上是19.8 GBps，相應的菱形分塊寬度最優(yōu)值分別為32和16.Girih事實上是菱形分塊和波前并行的混合方法，并采用細粒度并行，所有線程同時處理一個塊，因此我們不對它進行理論分析.

Fig. 14 Memory transfers of size 3843圖14 問題規(guī)模為3843的數(shù)據(jù)傳輸量

Fig. 15 Memory transfers of size 5123圖15 問題規(guī)模為5123的數(shù)據(jù)傳輸量

對于3D25P Stencil，星型密鋪方案在3個問題規(guī)模上的最優(yōu)分塊大小分別是36×4，28×3，28×3，對應自然塊大小分別為3.9 MB，4 MB，5.3 MB.因此30 MB的L3 cache將分別在8核，8核和6核時被完全使用，對應圖13(c)～15(c)中星型密鋪方案的數(shù)據(jù)傳輸量開始增長的位置.盡管核數(shù)更大時超出了L3 cache大小，但由于高階依賴，數(shù)據(jù)塊將快速收縮.此外，用時間維度塊大小為2的分塊密鋪方案等價于原始算法，因此它們產(chǎn)生同樣的緩存-內(nèi)存?zhèn)鬏斄考碩×N2.Pluto在3D25P Stencil上也有類似地變化趨勢和原因.

對于另外2個低階Stencil，L3 cache能保留所有12個核的自然塊.例如大小為2563的3D7P Stencil中自然塊占用內(nèi)存為2.5 MB，12核將占用所有30 MB L3 cache.因此，星型密鋪的帶寬變化接近于水平.但是性能在核數(shù)變多時增長變緩，這可能是由于有限的L3帶寬.

4 總 結(jié)

本文提出了一個新的并行框架，能夠高并發(fā)執(zhí)行星型Stencil，并以“自然塊”的新概念來標識不同形狀Stencil的特征.自然塊及其后繼塊能夠密鋪數(shù)據(jù)空間，并且它們沿著時間維度擴展后能形成迭代空間的密鋪.我們的方案能夠?qū)С龊啙嵉牟⑿锌蚣?，揭示一個格點的坐標和它更新時間步之間的關系.此外，專為星型Stencil研發(fā)了一個新穎的實現(xiàn)方式即“2次更新”，能夠在第1個后繼塊的執(zhí)行過程中對每個元素更新2次，從而改進核內(nèi)數(shù)據(jù)重用模式.實驗結(jié)果證明了此方法的有效性.

未來我們將著重于自動調(diào)優(yōu)方法來高效尋找最優(yōu)塊大小,基于提出的框架設計一個工具來自動生成Stencil代碼.