郭 鵬，袁 良，張?jiān)迫?，黃 珊,2

1.中國(guó)科學(xué)院 計(jì)算技術(shù)研究所 計(jì)算機(jī)體系結(jié)構(gòu)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100190

2.中國(guó)科學(xué)院大學(xué) 計(jì)算機(jī)與控制學(xué)院，北京 100049

1 引言

處理器計(jì)算性能和存儲(chǔ)器速度之間不斷擴(kuò)大的差距是影響軟件性能的關(guān)鍵瓶頸之一。為了解決這一問(wèn)題，在現(xiàn)代CPU和GPU上，硬件架構(gòu)設(shè)計(jì)人員采用越來(lái)越多的片上資源來(lái)形成硬件或軟件管理的高速緩存，應(yīng)用程序開(kāi)發(fā)人員則利用算法的局部性來(lái)最大化緩存內(nèi)的數(shù)據(jù)重用。

Stencil計(jì)算是科學(xué)和工程應(yīng)用中常見(jiàn)的循環(huán)模式。作為十三種Berkeley motifs之一，Stencil計(jì)算也被稱為“結(jié)構(gòu)化網(wǎng)格計(jì)算”。其他領(lǐng)域的許多算法，如動(dòng)態(tài)規(guī)劃和圖像處理也涉及相似的計(jì)算模式。Stencil是依照預(yù)定義的模式用鄰居節(jié)點(diǎn)更新所有節(jié)點(diǎn)。Stencil計(jì)算是d維空間格點(diǎn)隨時(shí)間迭代遍歷所有格點(diǎn)，這些d維空間格點(diǎn)也被稱為數(shù)據(jù)空間。數(shù)據(jù)空間沿時(shí)間維度更新，形成(d+1)維空間，稱之為迭代空間。Stencil存在多種分類方法，例如格點(diǎn)維度（1D,2D,…），鄰居節(jié)點(diǎn)數(shù)（3-point,5-point,…），形狀（星型，盒型），依賴類型（Gauss-Seidel,Jacobi）和邊界條件（常數(shù)，周期，…）。

Stencil計(jì)算的特征在于其計(jì)算密集度較低。d維Stencil的簡(jiǎn)單實(shí)現(xiàn)由(d+1)個(gè)循環(huán)組成，其中最外層的循環(huán)遍歷時(shí)間維，d層循環(huán)遍歷d維空間中的所有格點(diǎn)。這種方法具有較低的數(shù)據(jù)重用性，因此為帶寬受限。

分塊方法（tiling）是提高多重循環(huán)的數(shù)據(jù)局部性和并行性的強(qiáng)大轉(zhuǎn)換技術(shù)之一?？臻g分塊算法以特定順序遍歷數(shù)據(jù)空間，使得緩存內(nèi)的格點(diǎn)在被替換出cache前可被用于更新其所有鄰居節(jié)點(diǎn)，在2D和更高維度Stencil的單時(shí)間步內(nèi)改進(jìn)數(shù)據(jù)重用。但是空間分塊所能利用的局部性受限于Stencil的鄰居規(guī)模，為了增加數(shù)據(jù)重用，通常同時(shí)考慮時(shí)間維度與空間維度。大部分現(xiàn)有技術(shù)都適用于d維數(shù)據(jù)空間Stencil的(d+1)維迭代空間。

作為代表性的循環(huán)類型，Stencil的分塊和代碼的自動(dòng)生成已經(jīng)被詳細(xì)研究過(guò)了。例如廣泛使用的多面體模型是一個(gè)具有代表性的方法，其處理的循環(huán)嵌套上下界與數(shù)組下標(biāo)為周圍循環(huán)索引和預(yù)定義常數(shù)參數(shù)的函數(shù)。大多數(shù)分塊方法會(huì)使用固定的分塊形狀和大小來(lái)進(jìn)行分塊。例如，時(shí)間傾斜分塊（time skewed tiling）方法（在2D中為平行四邊形，在3D中為平行六面體，在高維中為超平行體）[1-2]消減了分塊之間的依賴關(guān)系。然而，大多數(shù)方法只用一個(gè)分塊形狀填充空間，經(jīng)常導(dǎo)致流水線啟動(dòng)并提供有限的并發(fā)性。

存在多種方法緩解流水線啟動(dòng)產(chǎn)生的開(kāi)銷[3-9]。也有一些允許將空間維度切割為任意數(shù)量的研究[10-11]。這些編譯器生成的代碼通常在編譯時(shí)需要固定的分塊大小。由于運(yùn)行時(shí)參數(shù)化可實(shí)現(xiàn)高效自動(dòng)調(diào)優(yōu)和性能可移植性，因此開(kāi)發(fā)了一些技術(shù)來(lái)允許參數(shù)化，比如矩形[12]和2D diamond[13]。然而，這些方法仍然導(dǎo)致高循環(huán)開(kāi)銷，很少允許精細(xì)優(yōu)化。更多細(xì)節(jié)將在第2章中詳細(xì)討論。

自動(dòng)生成的或手工編寫(xiě)的支持可選優(yōu)化和經(jīng)驗(yàn)自動(dòng)調(diào)優(yōu)參數(shù)的代碼已經(jīng)被證明是實(shí)現(xiàn)許多科學(xué)函數(shù)高性能的有效方法，如ATLAS（automatically tuned linear algebra software）[14]。對(duì)于Stencil計(jì)算，也存在自動(dòng)調(diào)整框架[15-18]和編程模型[19-20]。但是這些方法大多使用單一的超矩形分塊形狀[21-22]，并使用冗余計(jì)算來(lái)解決分塊間的依賴，同時(shí)編譯器通常無(wú)法充分利用硬件資源。

像其他數(shù)學(xué)函數(shù)（例如矩陣乘法和FFT（fast Fourier transformation））一樣，Stencil的定義是清晰簡(jiǎn)單的。因此，對(duì)于不同尺寸、形狀、順序和邊界條件的Stencil，應(yīng)該有一個(gè)簡(jiǎn)潔的并行化框架。然而，高性能社區(qū)通常采用具有冗余計(jì)算的超矩形分塊，而編譯器社區(qū)通常側(cè)重于傳統(tǒng)方法，例如多面體模型，其中大多數(shù)會(huì)產(chǎn)生不可讀和不可修改的代碼。借鑒OpenBLAS[23-25]的思路，提出一個(gè)簡(jiǎn)單而清晰的具備輕量級(jí)循環(huán)條件的并行框架，并允許精細(xì)的內(nèi)核優(yōu)化。

為了實(shí)現(xiàn)這一目標(biāo)，本文為Stencil計(jì)算設(shè)計(jì)了一種新的兩層密鋪算法。與以前直接對(duì)整個(gè)迭代空間進(jìn)行分塊的技術(shù)不同，首先使用不同塊類型密鋪數(shù)據(jù)空間，然后沿時(shí)間維度對(duì)塊進(jìn)行擴(kuò)展密鋪整個(gè)迭代空間。本文貢獻(xiàn)如下：（1）針對(duì)Stencil計(jì)算的一種新的兩層密鋪算法，使得數(shù)據(jù)空間的分塊可以并行執(zhí)行；（2）本文密鋪算法的實(shí)驗(yàn)分析與其他算法的比較。

本文的其余部分安排如下：第2章概述相關(guān)工作，并討論其缺陷。提出的密鋪算法將在第3章描述。第4章介紹實(shí)現(xiàn)和優(yōu)化。第5章介紹性能結(jié)果。第6章為本文的總結(jié)。

2 相關(guān)工作

2.1 分塊方法

分塊方法[26-28]是提高多重循環(huán)嵌套的數(shù)據(jù)局部性和并行性的最有效的轉(zhuǎn)換技術(shù)之一。Wonnacott和Strout對(duì)現(xiàn)有分塊算法的可擴(kuò)展性進(jìn)行了對(duì)比[29]。接下來(lái)本文對(duì)Stencil計(jì)算的一些相關(guān)工作進(jìn)行總結(jié)。

overlapped tiling：hyper-rectangle tiling[21-22]通常用于涉及高性能的手動(dòng)調(diào)優(yōu)。通過(guò)計(jì)算多個(gè)時(shí)間步，分塊之間的依賴將會(huì)被冗余計(jì)算抵消[30]，這種方法也被稱為 overlapped tiling[8]。Phillips和 Fatica[31]提供了一份GPUs上手寫(xiě)的3.5D分塊手動(dòng)調(diào)優(yōu)代碼實(shí)現(xiàn)。這種3.5D分塊方法通過(guò)時(shí)間平鋪提升了空間受限的2.5D切分算法。雖然hyper-rectangle的規(guī)則性使得高并發(fā)和更精細(xì)的優(yōu)化成為可能，但是有時(shí)候冗余操作過(guò)多反而掩蓋了性能提升。因此，本文將專注于零冗余算法。

time skewed tiling：該方法（所產(chǎn)生的塊，在2D時(shí)為平行四邊形，在3D時(shí)為平行六面體，在高維時(shí)為超平行體）[1-2]可以消除冗余計(jì)算，但是大部分的方法只使用一種分塊形狀來(lái)填充整個(gè)空間，并總是導(dǎo)致流水線啟動(dòng)和受限的并發(fā)。類似接下來(lái)要介紹的分塊技術(shù)，本文提出的密鋪算法使得完全并行成為可能，同步中的所有數(shù)據(jù)塊都可以并發(fā)執(zhí)行。

diamond tiling：Bondhugula等[32]首先提出了用于1D Stencil的diamond分塊。之后Bandishti等[3,9]將其擴(kuò)展到了高維Stencil。Grosser等[7]粗化了菱形的尖端，使其在2D成了六邊形，在3D變成了八面體。從本質(zhì)上來(lái)說(shuō)，diamond分塊方法通過(guò)平行六面體切割和旋轉(zhuǎn)空間得到最大并行，并行執(zhí)行的分塊波陣面共同組成了一個(gè)正交向量平行于時(shí)間維度的超平面。這種方法的主要優(yōu)勢(shì)在于它可以實(shí)現(xiàn)同時(shí)啟動(dòng)。

cache oblivious tiling：為了在未知內(nèi)存結(jié)構(gòu)參數(shù)的情況下充分利用數(shù)據(jù)局部性，提出了很多關(guān)于Stencil的cache oblivious（緩存無(wú)關(guān)）算法。Frigo和Strumpen提出了第一個(gè)串行[33]和并行[34]緩存無(wú)關(guān)Stencil算法。緩存無(wú)關(guān)平行四邊形分塊方法同時(shí)切分空間和時(shí)間維度，可以看作是帶波陣面并行的緩存無(wú)關(guān)版本時(shí)間傾斜算法。Pochoir[11]實(shí)現(xiàn)了同時(shí)切割所有空間維度的超平面切割法，與串行算法相比，它在提高并行性的同時(shí)保證了相同緩存復(fù)雜度。

split tiling：為了避免波陣面并行中的流水線啟動(dòng)開(kāi)銷，split tiling方法識(shí)別出每個(gè)分塊中沒(méi)有依賴的子塊并使其并發(fā)執(zhí)行，然后將結(jié)果發(fā)送給繼承者使其他區(qū)域也可以并發(fā)執(zhí)行[8]。Grosser等[10]提出了另外一種類似于超空間切割的緩存無(wú)關(guān)范式，這種循環(huán)切割方法可以遞歸地在所有的空間維度上進(jìn)行切割。

hybrid tiling：Strzodka 等[35]結(jié)合diamond分塊、平行四邊形分塊和流水線處理，提出了CATS（cache accurate time skewing）算法。Grosser等[6]介紹了混合六邊形和平行四邊形分塊算法。這些算法將迭代空間分解為兩部分，在時(shí)間維度和一個(gè)空間維度上使用hexagonal tiling或者diamond分塊，在其他空間維度上使用時(shí)間傾斜分塊。hexagonal tiling相當(dāng)于在空間維度上拉伸分塊的傳統(tǒng)diamond分塊，保證了在高階Stencil中每個(gè)分塊最多依賴于三個(gè)前置，是hybrid tiling的擴(kuò)展。hybrid split-tiling則是nested split-tiling和時(shí)間傾斜分塊的結(jié)合。

2.2 局限性

上述提到的diamond分塊[3]、hyperspace cut（cache oblivious tiling）[11]和 nested split-tiling[10]都是最近發(fā)展出的能夠?qū)崿F(xiàn)最大并發(fā)的方法。下文將說(shuō)明這些方法的一些缺點(diǎn)，值得一提的是，本文的密鋪算法是由一個(gè)能夠展現(xiàn)清晰的圖形分塊算法的數(shù)學(xué)框架得出，能夠克服如下這些問(wèn)題。

diamond分塊是一種基于多面體模型的編譯器轉(zhuǎn)換技術(shù)。為了實(shí)現(xiàn)運(yùn)行參數(shù)化以達(dá)到高效的自動(dòng)調(diào)優(yōu)和可移植性，在編譯運(yùn)行時(shí)分塊大小固定是它的主要限制之一。Bertolacci等[13]發(fā)展出在2D上的參數(shù)化版本，但是這并沒(méi)有解決它在自動(dòng)代碼生成上的問(wèn)題。diamond分塊的另一個(gè)缺點(diǎn)是它需要處理位于diamond尖端的小尺寸數(shù)據(jù)塊。Grosser等[7]提出在2D中diamond可以擴(kuò)展為六邊形分塊，并選擇截?cái)嗟陌嗣骟w作為3D中diamond的擴(kuò)展。但是他們并沒(méi)有給出高維（3D或更高維度）的分塊圖形模型。最后，diamond分塊用skewed hyper-rectangles填充整個(gè)迭代空間，很難找出合適的分塊大小以保證同時(shí)啟動(dòng)。同時(shí)，也很難設(shè)計(jì)出適應(yīng)于多層內(nèi)存的多層次版本。

常見(jiàn)的緩存無(wú)關(guān)實(shí)現(xiàn)總是受限于遞歸的開(kāi)銷并且只能達(dá)到有限的速度提升[36]。此外，很難將諸如SIMDization、數(shù)據(jù)打包、數(shù)據(jù)預(yù)取、數(shù)據(jù)填充、NUMA-aware調(diào)度、指令選擇（緩存旁路存?。┑冗M(jìn)一步的優(yōu)化加入到遞歸控制結(jié)構(gòu)中[37]。同時(shí)，由于緩存無(wú)關(guān)算法中沒(méi)有統(tǒng)一的數(shù)據(jù)塊形狀，粗化也很難實(shí)現(xiàn)。

hyperspace cut算法采用交錯(cuò)的時(shí)間和空間切割，這些算法的遞歸分治實(shí)現(xiàn)引入了子塊之間的人為依賴關(guān)系[37]。Tang等[38-39]提出cache oblivious wavefront可以消減這些錯(cuò)誤的依賴并且提高并行度，縮短算法的關(guān)鍵路徑。然而，對(duì)波陣面數(shù)據(jù)結(jié)構(gòu)的檢查不可避免地會(huì)引發(fā)運(yùn)行時(shí)間開(kāi)銷。

nested split-tiling最大的問(wèn)題就是其過(guò)高的同步開(kāi)銷，以d維的Stencil為例，需要同步2d次。需要指出，雖然Tang等[11]正確地證明了與本文算法相同的同步開(kāi)銷。但是，在仔細(xì)研究了他們的代碼后，發(fā)現(xiàn)這個(gè)基于普通的緩存無(wú)關(guān)遞歸框架的實(shí)現(xiàn)采用了部分和nested split-tiling一樣的策略，同樣將在d維的Stencil中同步2d次。雖然可以采用動(dòng)態(tài)隊(duì)列減少同步開(kāi)銷，但是運(yùn)行時(shí)調(diào)度的成本可能會(huì)將其抵消。

3 算法

本文使用 1 階星型 Stencil（1D 3-point，2D 5-point以及3D 7-point）作為例子，將展示這種方法適用于所有類型的Stencil。由于Gauss-Seidel Stencil強(qiáng)制45°并行啟動(dòng)，且不能通過(guò)分塊方法提升至完全并發(fā)啟動(dòng)[8]，因此本文只討論Jacobi Stencil。

3.1 重構(gòu)1D Stencil

首先介紹一種重構(gòu)1D diamond分塊的方法，這種方法是新的兩層n維Stencil分塊方法的基礎(chǔ)，可以看作是用不同形狀的分塊對(duì)迭代空間進(jìn)行密鋪。

1D 3-point Stencil的diamond分塊方法如圖1（a）所示。每個(gè)diamond被平行于x軸的橫線劃分為兩個(gè)大小相等的子三角形，迭代空間被正三角形和倒三角形填充，整個(gè)執(zhí)行過(guò)程為正三角形和倒三角形的交錯(cuò)計(jì)算。每條橫線代表一次同步，橫線上每個(gè)格點(diǎn)都更新了相同的時(shí)間步長(zhǎng)。兩條橫線之間為一時(shí)間分塊（time tile），其中每個(gè)網(wǎng)格點(diǎn)在相同的時(shí)間維度開(kāi)始計(jì)算并更新相同的時(shí)間步。觀察正三角形和倒三角形開(kāi)始計(jì)算的維度，可以看到，正三角形完全基于起始的7個(gè)格點(diǎn)進(jìn)行計(jì)算，而倒三角形起始于一個(gè)格點(diǎn)，基于這個(gè)格點(diǎn)和正三角形中已經(jīng)更新了一定時(shí)間步的格點(diǎn)進(jìn)行計(jì)算。正三角形和倒三角形起始維度的兩種形狀對(duì)迭代空間完成了密鋪。

每個(gè)正三角形和倒三角形可以通過(guò)各個(gè)格點(diǎn)在時(shí)間維度上更新的次數(shù)來(lái)表示。以圖1（a）為例，對(duì)于包含7個(gè)數(shù)據(jù)格點(diǎn)的數(shù)據(jù)塊（三角形在數(shù)據(jù)空間上的投影），按照各格點(diǎn)在該數(shù)據(jù)塊中更新的時(shí)間步數(shù)可表示為（0,1,2,3,2,1,0），其中心點(diǎn)更新了3次且其鄰居的更新步長(zhǎng)與其到中心的距離成反比。觀察倒三角形的結(jié)束維度，其更新時(shí)間步數(shù)為（3,2,1,0,1,2,3）。這樣，當(dāng)正三角形和倒三角形計(jì)算完畢所有格點(diǎn)到達(dá)橫線位置，所有的格點(diǎn)都更新了3次。整個(gè)執(zhí)行過(guò)程如圖1（b）所示。

Fig.1 Tiling of 1D Stencil圖1 1D Stencil的分塊

當(dāng)擴(kuò)展到更高維時(shí)應(yīng)注意兩個(gè)問(wèn)題。第一，雖然正三角形和倒三角形的投影形式相同，但是其執(zhí)行過(guò)程不同。正三角形中，所有的格點(diǎn)同時(shí)開(kāi)始計(jì)算，而倒三角形中所有的格點(diǎn)同時(shí)結(jié)束計(jì)算。第二，可以通過(guò)重用時(shí)間分塊之間的塊消除分割線上的同步。

3.2 2D Stencil密鋪

首先介紹最大更新的基本思想，最大更新即給定一個(gè)塊，其中每個(gè)點(diǎn)沿著時(shí)間維度向上更新直到無(wú)法滿足Stencil定義的依賴關(guān)系。形式化而言，對(duì)于任意格點(diǎn)a，當(dāng)存在至少一個(gè)鄰居格點(diǎn)a′∈neighbor(a)使得a′所在時(shí)間維度嚴(yán)格小于a，即time[a′]＜time[a]或a′?Bi時(shí)停止更新。當(dāng)給定一個(gè)d維數(shù)據(jù)空間塊Bi，通過(guò)關(guān)聯(lián)每個(gè)格點(diǎn)的更新時(shí)間步數(shù)來(lái)生成其(d+1)維擴(kuò)展塊，記為Bi。

由于Bi之間的轉(zhuǎn)化是高度對(duì)稱的，因此可以在B0的子塊上討論整個(gè)迭代空間的密鋪算法。具體來(lái)說(shuō)，假設(shè)有一B0塊，其中心點(diǎn)位于坐標(biāo)系的原點(diǎn)，則只考慮其中位于第一象限的點(diǎn)。在第i階段最大更新完成后，記這個(gè)子塊為B+i。將每個(gè)Bi塊中第一個(gè)（最后）時(shí)間步更新的格點(diǎn)構(gòu)成的區(qū)域稱為起始（結(jié)束）區(qū)域。

與前述的1D diamond分塊相似，將2D Stencil的迭代空間劃分為若干時(shí)間分塊，所有的格點(diǎn)在一個(gè)時(shí)間分塊開(kāi)始時(shí)位于相同時(shí)間維度，在一個(gè)時(shí)間分塊結(jié)束時(shí)更新相同的時(shí)間步。每個(gè)時(shí)間分塊中包含三個(gè)階段，即三種塊形狀。2D Stencil的三種階段塊形狀以及其轉(zhuǎn)化過(guò)程的俯視圖如圖2所示。其中B0是一個(gè)正方形，每個(gè)B0沿著x維和y維分割（虛線）為4塊，每個(gè)子塊包含B0的一邊，每?jī)蓚€(gè)共享同一邊的子塊合并形成B1。需要注意的是，雖然所有的B1有相同形狀，但分別是從兩個(gè)方向構(gòu)造出來(lái)的，在這里不妨記x軸方向兩B0塊之間的為B11，平行于y軸的兩B0塊之間的為B12。將每個(gè)B1沿著其與合并維度正交的維度劃分，四個(gè)相鄰的子塊結(jié)合在一起形成了B2。

2D 5-point Stencil某一時(shí)間分塊中的更新過(guò)程如表1所示。表中第二、三、四和五行分別為B0塊、B11塊、B12塊和B2塊的更新展示，其中b=3。第二列為各階段塊完成最大更新后的對(duì)應(yīng)更新步數(shù)，之后的第三、四和五列依次為各階段塊在對(duì)應(yīng)時(shí)間步更新的區(qū)域，為了方便理解在最后一列給出各階段塊的形狀。階段塊在各個(gè)時(shí)間步的計(jì)算都是數(shù)據(jù)空間中的一個(gè)矩形，而且當(dāng)前時(shí)間步矩形區(qū)域范圍可由前一時(shí)間步矩形區(qū)域計(jì)算得出，因此在笛卡爾坐標(biāo)系中當(dāng)給定起始區(qū)域坐標(biāo)范圍，即可確定整個(gè)塊各個(gè)時(shí)間步的計(jì)算范圍?？梢钥吹剑珺0塊在給出起始區(qū)域后，其計(jì)算范圍在x和y維度依照Stencil依賴遞減，而B(niǎo)11塊和B12塊在某一維度減小時(shí)另一維度增加，B2塊兩個(gè)維度都在增加，這是由Stencil的依賴關(guān)系決定的，在計(jì)算B0塊時(shí)，只能根據(jù)起始區(qū)域的值進(jìn)行計(jì)算，而計(jì)算B11塊、B12塊和B2塊時(shí)，雖然起始區(qū)域給定，但是B0塊已經(jīng)將其起始區(qū)域外的依賴點(diǎn)計(jì)算完畢，因此可以計(jì)算起始區(qū)域之外的點(diǎn)?？梢钥吹皆谡麄€(gè)計(jì)算過(guò)程中，同類分塊間不存在依賴關(guān)系，任一分塊中不存在無(wú)用格點(diǎn)，即不存在冗余計(jì)算。

Fig.2 Tessellating data space of 2D Stencil圖2 2D Stencil數(shù)據(jù)空間的密鋪

為了形式化描述2維Stencil的這些分塊，如圖3所示，在空間中引入笛卡爾坐標(biāo)系。令一個(gè)B0塊的中心點(diǎn)為原點(diǎn)（0,0），則其他B0塊的中心點(diǎn)為(2k0b，2k1b)，其中ki為整數(shù)。此B0在非負(fù)象限相鄰的B11塊、B12塊和B2塊的中心點(diǎn)坐標(biāo)分別為 (b，0)、(0，b)和(b，b)，將坐標(biāo)為0的維度記為起始維度，把其他維度稱為結(jié)束維度。選取B0中一點(diǎn)a(a0，a1)，其在第i階段的最大更新時(shí)間可以表示為T(mén)i(a0，a1)，點(diǎn)a在B0坐標(biāo)系統(tǒng)里的新坐標(biāo)為(a0，a1)，在B11坐標(biāo)系統(tǒng)里的新坐標(biāo)為 (b-a0，a1)，在B2坐標(biāo)系統(tǒng)里的新坐標(biāo)為 (b-a0，ba1)。在階段i，當(dāng)時(shí)間等于點(diǎn)a在結(jié)束維度到中心點(diǎn)的最大距離值時(shí)，a開(kāi)始更新，當(dāng)時(shí)間等于其在開(kāi)始維度到中心點(diǎn)的最大距離與b的差值時(shí)結(jié)束更新。以點(diǎn)a（2,1）為例，將各塊中心點(diǎn)移動(dòng)至原點(diǎn)，則a點(diǎn)坐標(biāo)取正后為（2,1）、（1,1）和（1,2）。用圓圈圈出點(diǎn)a在各階段所處更新位置，如表1中所示，B0塊中a點(diǎn)在第1步時(shí)更新了1步，在B1塊中y維為起始維度，第2步時(shí)更新了1步，在B2塊中第3步時(shí)更新了1步。

Fig.3 Cartesian coordinate systems圖3 引入笛卡爾坐標(biāo)系

3.3 高維、高階和周期邊界

本節(jié)將密鋪擴(kuò)展至高維和高階Stencil，并討論周期邊界的問(wèn)題。d維Stencil每個(gè)時(shí)間分塊的計(jì)算都包含d+1個(gè)階段。記d+1個(gè)階段中鑲嵌于d維數(shù)據(jù)空間的d+1種數(shù)據(jù)塊為Bi(0≤i≤d)，即在階段i，使用Bi鑲嵌數(shù)據(jù)空間。其中B0為超立方體（在1D時(shí)為線段，2D時(shí)為正方形，3D時(shí)為立方體，4D時(shí)為超立方體，…）。

Table 1 Tessellating data space of 2D Stencil表1 2D Stencil的密鋪迭代

令B0中一點(diǎn)a(a0，a1，…，ad-1)，其在第i階段的最大更新時(shí)間可以表示為T(mén)i(a0，a1，…，ad-1)。不失一般性的，假設(shè)包含a的B0塊的中心點(diǎn)為 (b，…，b，0，…，0)，其中前i個(gè)坐標(biāo)為b。點(diǎn)a在Bi坐標(biāo)系統(tǒng)里的新坐標(biāo)為 (b-a0，b-a1，…，b-ai-1，ai，ai+1，…，ad-1)。則開(kāi)始計(jì)算時(shí)間步和結(jié)束計(jì)算時(shí)間步為：

將其合并可以得出：

特別的：

高階Stencil在實(shí)際應(yīng)用中十分常見(jiàn)。如果Stencil中某一維度上的階數(shù)為m，可以將每m個(gè)點(diǎn)結(jié)合在一起，稱之為超級(jí)節(jié)點(diǎn)。這樣就將超級(jí)節(jié)點(diǎn)中的m-階依賴減小為了1階依賴。在圖4中給出了在x維上m=2的例子。對(duì)高階依賴的所有其他維度進(jìn)行類似操作就可以使用Bi對(duì)這些1階依賴的超級(jí)節(jié)點(diǎn)執(zhí)行密鋪操作，而超級(jí)節(jié)點(diǎn)中所有格點(diǎn)的更新步數(shù)與其在Bi中的更新步數(shù)相同。

Fig.4 Handle high-order Stencil圖4 處理高階Stencil

對(duì)于周期邊界的情況，需要將一側(cè)邊界的一個(gè)塊或高維中的多個(gè)塊進(jìn)行拉伸。輸入尺寸不是塊大小整數(shù)倍的1D Stencil例子如圖5所示，其中灰色點(diǎn)為邊界格點(diǎn)額外需要的數(shù)據(jù)。通過(guò)將一個(gè)diamond塊沿著空間維度拉伸而產(chǎn)生了一個(gè)六邊形塊，保證在兩側(cè)邊界的子塊為相同類型分塊。

4 實(shí)現(xiàn)

基于上述清晰的數(shù)學(xué)描述，可以直接得出對(duì)應(yīng)代碼。在第5章，將提供手工實(shí)現(xiàn)的性能結(jié)果。當(dāng)然，自動(dòng)代碼生成的工具也很容易實(shí)現(xiàn)，這將是下一步的工作內(nèi)容。

Fig.5 Handle periodic boundary condition圖5 處理周期邊界

因?yàn)橥浑A段的所有塊都可以同時(shí)執(zhí)行，可以很容易地實(shí)現(xiàn)密鋪算法的并行化。因此只在共享內(nèi)存的設(shè)備上使用了OpenMP pragma parallel for指令。

4.1 粗化

對(duì)于3D或者高維Stencil，其各個(gè)維度的尺寸通常都小于1 000。為了更好地利用硬件預(yù)取功能，現(xiàn)有的做法通常是將單位跨度的維度切掉。但是這樣就會(huì)限制并發(fā)性。此外，時(shí)間維度的尺寸b（由最大更新產(chǎn)生）通常受限于一個(gè)塊各個(gè)維度中最小尺寸，因此在空間維度不相等的劃分可能會(huì)引發(fā)時(shí)間維度上受限的數(shù)據(jù)重用。

粗化的另一個(gè)驅(qū)動(dòng)因素是本文的密鋪算法可能會(huì)導(dǎo)致無(wú)效的數(shù)據(jù)訪問(wèn)模式。以2D為例，如表1中所示，兩個(gè)四面體B1塊在第一個(gè)或者最后一個(gè)時(shí)間步會(huì)產(chǎn)生行主數(shù)組的列中的元素的更新，這種低效的空間位置使用可能會(huì)降低性能。

密鋪框架是高度可修改的，它允許對(duì)各個(gè)維度塊尺寸的修改，并通過(guò)沿著各個(gè)維度拉伸來(lái)實(shí)現(xiàn)粗粒度化。為了便于解釋，2D Stencil的粗粒度化如圖6所示。為了改善空間局部性，B0的結(jié)束區(qū)域（圖6中的深灰色區(qū)域）和Bd的起始區(qū)域（圖6中的黑色區(qū)域）粗化為相同大小的2D塊。

如上所述，在第一階段和最后一階段的數(shù)據(jù)塊有相同的形狀。在合并時(shí)段之間的Bd和B0可以減少同步并提高數(shù)據(jù)重用。通過(guò)實(shí)驗(yàn)可以得出，如圖6所示，兩個(gè)B0在某一維度上的間距等于B0結(jié)束時(shí)在這個(gè)維度上的大小。這也就保證了B0的結(jié)束塊形和Bd的起始?jí)K形完全一致。

Fig.6 Realization of 2D Stencil圖6 2D Stencil的實(shí)現(xiàn)

4.2 具體實(shí)現(xiàn)

本節(jié)以數(shù)據(jù)規(guī)模為N2的非周期2D Stencil為例介紹算法的具體實(shí)現(xiàn)。由于算法支持任意階數(shù)的Stencil，因此參數(shù)化其在x軸、y軸方向的斜率分別為XSLOPE和YSLOPE，因此數(shù)組計(jì)算范圍為XSLOPE～N+XSLOPE-1。

如圖6所示，在算法中，所有Bi中的每個(gè)時(shí)間步計(jì)算的數(shù)據(jù)空間都是一個(gè)矩形，根據(jù)Stencil計(jì)算本身的依賴關(guān)系和最大更新的思想，當(dāng)前時(shí)間步的矩形可由前一時(shí)間步的數(shù)值以及斜率計(jì)算得出，由前述可知，在同一時(shí)間分塊中，各階段塊規(guī)律分布，因此由某一階段塊起始區(qū)域的(xmin，ymin)即可計(jì)算得出該階段塊各時(shí)間步需要計(jì)算的數(shù)據(jù)空間大小，同時(shí)計(jì)算得出當(dāng)前時(shí)間段所有相同階段塊的所有信息。

Bx和By分別為B0塊中起始區(qū)域在空間維度中的長(zhǎng)度，tb為更新步數(shù)，這三個(gè)參數(shù)也是算法的主要調(diào)優(yōu)選項(xiàng)之一，由這三個(gè)參數(shù)可計(jì)算得出其結(jié)束區(qū)域的區(qū)塊大?。?/p>

如前所述，為了實(shí)現(xiàn)算法粗粒度化，B0以及B2分塊之間在x和y方向的距離分別為Bx+bx和By+by，記為ix、iy，則在x軸、y軸方向各有B0塊個(gè)數(shù)。

其中，┐表示向上取整。對(duì)于B11塊和B12塊，觀察圖6中形狀可得出：

因?yàn)橹灰嬖贐0或B1塊，必定有B2塊使其達(dá)到最大更新，故：

不考慮越界問(wèn)題，由圖可以確定各塊的點(diǎn)坐標(biāo)，選取起始區(qū)域左下角塊的左下角點(diǎn)作為索引點(diǎn)，例如，B0塊的索引點(diǎn)為：

觀察圖6可以看出，B0塊更新后B1塊、B2塊依次更新，完成一次最大更新，之后在B2塊的結(jié)束區(qū)域繼續(xù)向上更新，需要注意的是這時(shí)對(duì)應(yīng)之前B2塊的位置已經(jīng)變?yōu)锽0塊，B11和B12的位置也發(fā)生了互換，在這里引入level，在其奇偶發(fā)生變化時(shí)B0和B2、B11和B12的位置相互轉(zhuǎn)換，為了提高代碼的可讀性，令T從 -tb開(kāi)始，將所有B0塊和B2塊融合在一起記為B02，將其坐標(biāo)進(jìn)行整合：

所有6個(gè)索引點(diǎn)如圖6中黑點(diǎn)所示。

算法的主體代碼如算法1所示，其中包含4層循環(huán)。最外層循環(huán)控制時(shí)間分塊以tb遞增，第一層循環(huán)內(nèi)包含兩個(gè)3層循環(huán)，其中第一個(gè)用于計(jì)算B0塊和B2塊，第二個(gè)用于計(jì)算B1塊。如前文所述，算法將B0塊與B2塊結(jié)合在一起形成B02塊，B02塊從原B2塊的起始區(qū)域開(kāi)始計(jì)算，到其對(duì)應(yīng)的B0塊的結(jié)束區(qū)域?yàn)橹菇Y(jié)束計(jì)算，其計(jì)算的時(shí)間步長(zhǎng)為2×tb，需要計(jì)算的塊個(gè)數(shù)根據(jù)level不同為B0或B2塊個(gè)數(shù)。需要注意的是，雖然令時(shí)間步從-tb開(kāi)始，但是實(shí)際計(jì)算中需要用max(tt，0)限定其計(jì)算范圍。如圖6所示，當(dāng)給定塊編號(hào)n，結(jié)合x(chóng)軸和y軸方向間距ix和iy以及索引點(diǎn)坐標(biāo)(xleft02，ybottom02)，即可確定當(dāng)前塊的索引坐標(biāo) (myxmin，myymin)進(jìn)而確定 (myxmax，myymax)，對(duì)于B2塊，其x軸和y軸方向的計(jì)算范圍隨時(shí)間步都向外擴(kuò)張，B0塊恰好相反，其x軸和y軸方向的計(jì)算范圍隨時(shí)間步都向內(nèi)收縮，由于其變化幅度恰好互負(fù)，可以簡(jiǎn)單地通過(guò)abs(t+1-tt-tb)×SLOPE來(lái)確定當(dāng)前時(shí)間步的計(jì)算區(qū)域。對(duì)于B1塊，其起始時(shí)間步比B02高tb，(myxmin，myymin)和(myxmax，myymax)的計(jì)算方式與B02塊相同，由于B11塊和B12塊形狀不同，其在x軸方向和y軸方向的變化恰好相反，在這里分開(kāi)計(jì)算，可以使用一個(gè)簡(jiǎn)單的判斷語(yǔ)句，令前#B11個(gè)塊去計(jì)算B11塊，剩下的#B12個(gè)塊去算B12。如前文所述，Bi塊的每個(gè)時(shí)間步計(jì)算區(qū)域都是矩形，當(dāng)確定了各時(shí)間步(myxmin，myymin)和(myxmax，myymax)的計(jì)算方法，之后就是簡(jiǎn)單的兩層循環(huán)更新格點(diǎn)了。

算法12D Stencil

5 實(shí)驗(yàn)結(jié)果

5.1 實(shí)驗(yàn)環(huán)境

實(shí)驗(yàn)環(huán)境為配置有兩個(gè)2.70 GHz的Intel Xeon E5-2670處理器的機(jī)器，并在單核到24核上進(jìn)行了測(cè)試。單核L1和L2大小分別為32 KB和256 KB，一個(gè)socket上的12個(gè)核共享30 MB的L3 cache。ICC版本為16.0.1，優(yōu)化選項(xiàng)為“-O3-openmp”。

將密鋪算法和上文提到的三種高度并行方案中的兩種進(jìn)行了比較，即Pluto和Pochoir。其中Pluto采用diamond分塊方法，Pochoir采用hyperspace cut分塊方法。將split tiling代碼和本文方案進(jìn)行比較是不合理的，因?yàn)樗鼈儾捎昧祟~外的有效向量化方法。本文只比較具有相同編譯選項(xiàng)的原始核心代碼，由編譯器本身來(lái)完成向量化。

測(cè)試包括四種星型Stencil（Heat-1D 3-point、Heat-1D 5-point、Heat-2D 5-point和Heat-3D 7-point）以及兩種盒型Stencil（Heat-2D 9-point和Heat-3D 27-point）。各種測(cè)試的用例如表2所示，其中周期邊界與非周期邊界的數(shù)據(jù)規(guī)模相同。Benchmark測(cè)試參數(shù)大小與文獻(xiàn)[3]中的參數(shù)大小一致，保證數(shù)據(jù)規(guī)模大小能夠真實(shí)反映算法性能同時(shí)不會(huì)導(dǎo)致過(guò)長(zhǎng)的測(cè)試時(shí)間。由于密鋪方案的代碼包含更多的可調(diào)參數(shù)，將其他參數(shù)設(shè)置為分塊大小的一半或兩倍，使其與Pluto和Pochoir中的分塊近似。

Table 2 Problem size of Benchmark表2 問(wèn)題規(guī)模

5.2 分析

總體來(lái)說(shuō)，本文算法在星型Stencil上獲得了有競(jìng)爭(zhēng)力的性能數(shù)據(jù)，并且在盒型Stencil上表現(xiàn)出更好的性能。Pluto和Pochoir展示了在文獻(xiàn)[3]中相似的性能趨勢(shì)。Pochoir表現(xiàn)出比Pluto更好的可擴(kuò)展性，而Pluto經(jīng)常獲得更好的性能。密鋪的代碼實(shí)現(xiàn)了更好的加速和可擴(kuò)展性。

1D Stencil的結(jié)果如圖7（a）和圖7（b）所示。所有的這三種方案都顯示出了線性可擴(kuò)展性。新算法的性能和Pluto的性能相似，并優(yōu)于Pochoir的性能。這是因?yàn)槊茕佀惴ê蚉luto產(chǎn)生相同的菱形劃分代碼，而Pochoir利用動(dòng)態(tài)調(diào)度方法生成不同大小的梯形塊。

Fig.7 1D non-periodic boundary condition圖7 1D非周期邊界

Fig.8 2D non-periodic boundary condition圖8 2D非周期邊界

圖8（b）和圖8（a）中顯示的是非周期2D Stencil的性能結(jié)果。對(duì)于Heat-2D 5-point Stencil，本文的代碼和Pochoir展現(xiàn)出相似的性能表現(xiàn)和可擴(kuò)展性。Pluto則如文獻(xiàn)[3]中提到的受限于負(fù)載不均衡。在24核時(shí)，Pluto只比另兩種高不到5%。對(duì)于2D 9-point Stencil，Pluto和Pochoir表現(xiàn)出了和5-point Stencil中相同的趨勢(shì)。但是本文的代碼平均比他們提升了14%和20%。周期邊界2D Stencil的結(jié)果如圖9（b）和圖9（a）所示。最高性能比Pochoir提升了10%和20%。

非周期的3D Stencil的性能結(jié)果如圖10（a）和圖10（b）所示。本文的代碼和Pochoir的代碼相比于Pluto有更好的可擴(kuò)展性。和2D星型例子相似，在超過(guò)21核時(shí)Pluto表現(xiàn)出優(yōu)于Pochoir和本文的性能。但是在3D 27-point Stencil中，本文的代碼性能顯著優(yōu)于Pluto和Pochoir，最高提高74%和100%，平均提升達(dá)30%和99%。周期邊界3D 7-point Stencil和3D 27-point Stencil的結(jié)果如圖11（b）和圖11（a）所示。最高性能比Pochoir提升了25%和40%。

6 結(jié)論

本文提出了一種新的針對(duì)Stencil的高并發(fā)的并行框架。這是一種基于新的兩層時(shí)間分塊的算法。通過(guò)一系列塊來(lái)鑲嵌空間維度，然后這些塊的擴(kuò)展密鋪了整個(gè)迭代空間。密鋪算法提供了簡(jiǎn)潔的并行化框架，并揭示了網(wǎng)格點(diǎn)的坐標(biāo)和更新時(shí)間步長(zhǎng)之間的關(guān)系，該算法可以最大化并發(fā)執(zhí)行，在執(zhí)行過(guò)程中沒(méi)有冗余計(jì)算，有簡(jiǎn)潔的循環(huán)條件并可以適應(yīng)Stencil不同的尺寸、形狀、階數(shù)和邊界條件。實(shí)驗(yàn)結(jié)果表明，該方案在星型Stencil上實(shí)現(xiàn)了相當(dāng)?shù)男阅埽诤行蚐tencil上表現(xiàn)出更好的性能。對(duì)于1D Stencil，由于密鋪算法表現(xiàn)為和Pluto相同的diamond分塊，性能與其持平，對(duì)于2D Stencil和3D Stencil密鋪算法都有更好的性能表現(xiàn)，其中在非周期3D27p Stencil中，相較于Pluto有12%的性能提升，周期邊界3D27p Stencil相較Pochoir的性能最高提升40%。

Fig.9 2D periodic boundary condition圖9 2D周期邊界

Fig.10 3D non-periodic boundary condition圖10 3D非周期邊界

Fig.11 3D periodic boundary condition圖11 3D周期邊界

未來(lái)將基于該框架設(shè)計(jì)一套自動(dòng)代碼生成工具。由于其中包含相較于其他方法更多的參數(shù)，現(xiàn)在的工作正致力于尋找可以自動(dòng)找出最優(yōu)分塊大小的方法。