曹 杭 袁 良 黃 珊 張?jiān)迫?徐勇軍 陸鵬起 張廣婷

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

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

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

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

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

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

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

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

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

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

1 相關(guān)工作

1.1 分塊技術(shù)

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

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

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

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

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

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

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

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

1.2 優(yōu)化方法

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

2 算法描述

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

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

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

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

2.1 自然塊B0(0)

當(dāng)所有格點(diǎn)在時(shí)間維度上的值相同時(shí)，更新某一格點(diǎn)t個(gè)時(shí)間步所需的最小鄰居點(diǎn)集被稱為時(shí)間步為t的自然塊，簡(jiǎn)稱自然塊.

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

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

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

2) 存在至少一個(gè)b的鄰居b′，b′的時(shí)間維度嚴(yán)格小于b(time[b′]

但是這2種概念并不完全等價(jià)，實(shí)際上它們互為補(bǔ)充.一方面，最大更新沒有給出密鋪數(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個(gè)自然塊和星型Stencil的12個(gè)自然塊.注意這些自然塊可以并行執(zhí)行.

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

(1)

Fig. 3 Determining the center points of B1圖3 由0確定B1塊的中心點(diǎn)

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

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

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

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

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

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

3 實(shí)現(xiàn)方法

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

3.1 粗粒度化

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

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

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

3.4 并行化和向量化

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

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

3.5 代碼實(shí)現(xiàn)

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

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

算法1.2D星型Stencil算法.

輸入：時(shí)間步長(zhǎng)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 實(shí)驗(yàn)結(jié)果

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

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

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

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

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

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

4.2 結(jié)果分析

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

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

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

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

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

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

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

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

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

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

4 總 結(jié)

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

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