一种基于空间密铺的星型Stencil并行算法

曹 杭1,2 袁 良1 黄 珊1,2 张云泉1 徐勇军1 陆鹏起1,2 张广婷1

1(计算机体系结构国家重点实验室(中国科学院计算技术研究所) 北京 100190)

2(中国科学院大学 北京 100049)

摘 要 Stencil计算(模板计算)是科学工程应用中一类常见的嵌套循环算法.分块方法是提高数据局部性和并行性的高效优化技术之一,目前已有大量针对分块方法的探索,但现有工作往往对不同Stencil形状都采用同一处理方法.首先在空间层面引出“自然块”的概念来区分星型Stencil和盒型Stencil的特征,然后提出一个新的针对星型Stencil的2层密铺方案,此方案中自然块和它的后继块可以密铺数据空间区域,这些分块沿着时间维度扩展,能够密铺整个迭代空间.此外,针对星型Stencil设计了一个新颖的“2次更新”优化技术,改善了核内数据重用模式.理论分析表明:此方案相比现有方法有更低的缓存复杂度,实验结果证实了此方案的有效性.

关键词 Stencil计算;密铺;星型Stencil;盒型Stencil;自然块

Stencil计算(模板计算)是科学工程应用中一类常见的嵌套循环算法,也被称为“结构化网格计算”,是13个伯克利核心计算模式之一.在动态规划、图像处理等领域的很多算法中也包含相似的依赖模式.Stencil定义了一种更新结点所需邻居结点的模式.Stencil计算在时间维度上迭代更新规则的d维网格(也称为数据空间),数据空间沿着时间维度更新,产生的d+1维空间被称为迭代空间.

与其他的数值算法不同,例如稠密或稀疏矩阵计算类型的多样性主要来源于输入数据如稠密或稀疏矩阵的数据量,Stencil的计算类型本质上取决于问题的类型.Stencil可以从不同的角度进行分类,如格点维度(1维,2维,…),邻居数目,即阶数(3点,5点,…),形状(盒型,星型,…),依赖类型(高斯-赛德尔,雅可比,…)和边界条件(常量的,周期的,…)等.

Stencil计算具有计算密集度低的特征,d维Stencil的简单实现由d+1个循环组成,其中最外层循环遍历时间维度,内层循环更新d维空间的所有格点.简单实现的数据重用度较差且典型的带宽受限.

分块是提高多重嵌套循环的数据局部性和并行度最有效的优化技术之一,目前有大量针对Stencil计算分块方案的研究.空间分块算法促进了2维至更高维度Stencil在一个时间步内的数据重用.为了进一步增强数据重用度,往往同时考虑时间维度和空间维度,后文将详细介绍这些技术.另一类工作采用简单的超矩形分块形状[1-4]并引入冗余计算来解决块间数据依赖.这些研究包括自动调优框架[5-12]、编程模型[13-14]、CPU[2,15-18]、GPU[19-23]和Phi[24-25]上的手动调优实现.

不同的Stencil类型也衍生出不同的解决方案.Pluto系统[26]延伸到周期性边界条件Stencil的处理[27].菱形分块最初是针对1维Stencil的分块[28],然后扩展到更高维度的Stencil[29].高速缓存参数无关(cache oblivious)方案从最初的算法[30-31],CORALS(cache oblivious parallelograms)方案[32],发展到能处理任意Stencil并达到最大并发度的Pochoir系统[33].

据我们所知,目前没有针对不同Stencil形状的解决方案,特别是盒型和星型Stencil.星型Stencil只依赖于每个坐标轴上的点,盒型Stencil在此之外还有对角线上的数据依赖.不区分盒型和星型Stencil的根本原因可能是盒型Stencil包含星型Stencil,任何适用于盒型Stencil的算法也能应用到星型Stencil.

本文首先在数据空间上提出“自然块”的概念来确定星型和盒型Stencil之间的区别.然后针对星型Stencil提出了一个新的2层密铺方案,在此方案中,自然块和它的后继块可以密铺整个空间,而他们沿着时间维度扩展后形成对迭代空间的密铺.密铺框架类似于文献[34]中提出的方法.我们将统一阐述基于自然块概念的2种方案并将已有的方案称为盒型密铺而将新提出的方案称为星型密铺.

在密铺方案中区别盒型和星型Stencil有2个主要的优点,分别利用不同层次上的数据局部性.首先,当给定缓存大小时,相比盒型密铺,用星型Stencil的自然块密铺能得到更大的块,这样可以更有效地减少内存系统的压力;其次,我们针对第二后继块中的元素提出了一个新颖的“2次更新”技术,可对1个元素连续进行2次更新,进一步提升了核心数据的重用.

1 相关工作

1.1 分块技术

分块技术[35-39]是探究多层嵌套循环数据局部性和并行度最有效的转换技术之一.Wonnacott和Strout[40]对比了一些现有分块方案的可扩展性,下面我们将总结一些和Stencil计算有关的技术.

1) 重叠分块(overlapped tiling).高性能领域中经常在手动调优的实现中采用超矩形分块[1-4].为了执行多个时间步,往往采用冗余计算[41]来解决分块之间的依赖,这就是重叠分块[28,42].Philips和Fatica[22]提出了GPU上手动调优的手写3.5维分块实现方案,在2.5维分块方案[4]上增加了时间维度上的划分.Demmel等人[43-44]减少了2倍的冗余计算.形状规则的超矩形可以获得更高的并发度和更多细粒度优化.

2) 时间偏移(time skewing).时间偏移分块[45-47](分块形状在2维上是平行四边形,在3维上是平行六面体,在更高维上是超平行体)消除了冗余计算,但是大多数方法只用单一形状的块填充空间,导致流水线启动和有限的并发度[48-49].

3) 菱形分块(diamond tiling).Bondhugula等人[26]对1维Stencil采用菱形分块.Bandishti等人[29,50]将此方法扩展到更高空间维度的Stencil.Grosser等人[51]对2维的菱形尖端粗粒度化形成六边形,在3维上演化成截掉顶端的八面体.

4) 缓存无关分块(cache oblivious tiling).缓存无关技术旨在内存层次结构参数未知时充分利用数据局部性.Frigo和Strumpen提出了第1个串行[30]和并行[31]缓存无关Stencil算法.缓存无关平行四边形方法[32]同时分离时间空间维度,但是会导致波阵面并行.Pochior[33]实现了多维空间划分,能够同时划分所有空间维度.

5) 分裂分块(split tiling).分裂分块方法发掘每个块中的独立子块,然后将这些子块发送给他们的后继块,使得剩下的区域同样可以并发执行[28,52-54].Grosser等人[55]提出了另一个类似于具有缓存无关范式的多维空间划分的分裂分块方法.嵌套分裂分块(the nested split-tiling)[56]能在所有空间维度递归分裂.

6) 混合分块(hybrid tiling).Strzodka等人[57]提出的CATS(cache accurate time skewing)算法结合了菱形分块、平行四边形分块和流水线处理.Grosser等人[58]采用混合六边形和平行四边形的分块算法.这些算法取某一维空间和时间维度组成平面,用六边形分块和菱形分块进行划分,在剩下的空间维度中采用时间偏移分块,从而分解迭代空间.混合分裂分块方法[56]结合了嵌套分裂分块和时间偏移分块.

7) 密铺分块(tessellating).密铺方案[34]适用于包括星型和盒型等不同Stencil类型.空间被一系列分块密铺,分块之间可以无冗余并行执行.相应地,这些分块在时间维度上扩展后得到的扩展块可以在迭代空间中形成密铺.简单明了的数学特性使得这种方法可以与细粒度优化方法相结合.

1.2 优化方法

为了使核心程序达到更高的性能,核心内的优化技术如计算重用[59-60]、计算重组[61]、向量化[15,62-64]和数据布局转换[65-66]等十分关键.但是这些方法大多没有对不同的Stencil形状进行区分.本文设计实现了一个新颖的技术来提升星型Stencil核心程序级别上的数据重用.这种技术也可以应用到其他的分块方法上.

2 算法描述

本节首先介绍一个新的概念,即自然块,来区分不同的Stencil形状特征.然后围绕自然块和后继块这2个方面,详细描述新提出的星型Stencil密铺方案以及已有的盒型密铺方案[34].最后将密铺方案应用于高阶Stencil和不同的边界条件.

Gauss-Seidel Stencil需要执行逐点45°流水线启动,使用分块技术也不能使其完全并发启动[28].因此我们仅考虑Jacobi Stencil.此外,星型Stencil只能应用到2维数据空间中,我们只讨论2维Stencil.在第5节的实验中,对于3维星型Stencil将保留1维数据空间不进行切分.

密铺3维迭代空间对应于计算2维Stencil,我们首先将其划分成多个相同的时间片.在一个时间片开始时所有格点处于同一时间维度,结束时所有格点都更新了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的自然块.

对于1维Stencil,盒型和星型形状是相同的,时间步长为4的自然块是0,1,2,3,4,3,2,1,0.自然块中的数字代表相应位置上格点的时间更新步数.因为开始时所有格点处于同一时间维度,所以它们在第1个时间步同时开始更新.代表数据空间中物理格点的集合时,自然块用B0表示(之后将只会用到自然块的形状,即盒型Stencil的自然块是正方形的,星型Stencil的自然块是菱形的).自然块B0为Stencil定义的扩展,将一个2D9P Stencil的网格点(x,y),更新t个时间步使之从(x,y,0)迭代到(x,y,t),需要(2t+1)2个网格点,具体范围为:(x+a,y+b,0),其中a∈[-t,t],b∈[-t,t];而对于2D5P Stencil,范围为2t2+2t+1个网格点,其中a∈[-t,t],b∈[-t,t],且a+bt.而当自然块定义中t=1时,自然块与Stencil定义等价.在B0基础上赋予每个格点1个新的时间维度来表示其时间更新步数,此时自然块用0表示,用于解释迭代空间的密铺.

最大更新方法[34]与自然块概念相近.具体定义为:在数据空间中给定B,沿着时间维度最大更新每个格点,直到不满足Stencil定义的依赖关系为止.形式上而言,对所有格点bB,满足以下2个条件中任一个时,停止更新:

1) time[b]=tt是给定的最大时间更新步数;

2) 存在至少一个b的邻居b′,b′的时间维度严格小于b(time[b′]<time[b])或b′∉B.

自然块和最大更新在某种程度上可以看作是对偶的概念.给定同样的格点集B0,且所有格点更新前都处于同一时间维度,则这2种机制实际上会产生同样的0.此外,因为分块之间没有依赖关系,用这2种方法产生的所有分块都能并发执行.

但是这2种概念并不完全等价,实际上它们互为补充.一方面,最大更新没有给出密铺数据空间块的具体形状.盒型密铺[34]采用盒型Stencil的自然块(超立方体)来密铺盒型和星型Stencil,盒型Stencil包含星型Stencil,因此可以保证结果的正确性.在本工作中,我们将采用不同类型Stencil各自的自然块.

另一方面,单一的自然块0不能密铺迭代空间.我们需要其他2个后继块B1B2.但是它们包含的格点在更新前初始时间维度并不完全相同,这使我们不能采用自然块的概念.最大更新方法可以无视格点的时间维度,因此迭代空间中后继块的扩展块12可以由最大更新方法产生.

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

盒型Stencil自然块B0的形状是正方形,星型Stencil自然块的形状是菱形,这2种形状都能密铺2维数据空间.图2(a)和图2(b)(忽略灰色区域)分别表示盒型Stencil的4个自然块和星型Stencil的12个自然块.注意这些自然块可以并行执行.

Fig. 2 Tessellating data space with B0
图2 用B0密铺数据空间

下面给出数学公式将Bi扩展为i(i=1,2,3),即赋予每个点(x,y)在i中的时间更新步数,时间更新步数在盒型Stencil中用TBi(x,y)表示,在星型Stencil中用TSi(x,y)表示.密铺方案的特性使得我们只需推导一个自然块B0中的公式.我们将这个自然块的中心点(更新t=4个时间步的格点)在笛卡儿坐标系中的坐标设为(0,0),2个坐标轴的正方向分别是向右和向上.此外,由Bi密铺到Bi+1密铺的转变允许我们只需关注满足条件xy≥0的点(x,y),因此我们只对每个Bi密铺图形中的灰色三角区域进行说明,并给出i图形中灰色三角形区域的公式.需要注意的是,图2~8中的灰色三角形区域表示同一格点集合,这个集合是总能属于一个Bi块的最大格点集合.如图3所示,根据自然块的定义,1维盒型和星型Stencil中0的计算为

(1)

Fig. 3 Determining the center points of B1
图3 由0确定B1块的中心点

由图3可知,这些自然块中只有中心点被更新了t次,单独由0并不能密铺更新步数为t的时间片,更不能密铺整个3维迭代空间.我们需要引入另外2个后继块B1B2来密铺数据空间,并且利用其对应的12完整地密铺迭代空间,即更新所有格点t个时间步.

后继块B1,B2等是对自然块B0进行分裂而后重组格点所产生的.分裂和重组的关键在于确定这些B1的中心点,之后每个格点可以简单地划分给离它最近的中心点.这种方法保证除边界上的点外,每个格点只属于某一个B1,并且B1可以密铺数据空间.

自然块的一个性质是中心点在时间维度更新步数最多,因此最显然的处理方式是将更新次数为0的格点凸集的中心点作为后继块的中心点.需要注意的是,后继块的中心点总在前驱块的边界上,只用考虑边界上标为0的格点凸集.

由图3可知,盒型和星型Stencil中所有标为0的格点分别在正方形和菱形的边界上.根据凸集的要求,4条边被视为4个凸集,因此方框中的4个格点0即为对应4个B1的中心点.

图4(a)和图4(b)中展示了由B0B1密铺数据空间的转变过程,其中虚线代表对B0的划分,可以看到每个B0被划分成4个子块.图4(c)和图4(d)展示了用B1密铺数据空间的方法,其中共享同一个中心点的B0子块合并组成1个B1.经过此转变过程,盒型Stencil由4个B0转换为12个B1,星型Stencil由12个B0转换为16个B1.

Fig. 4 Tessellating data space with B1
图4 用B1块密铺数据空间

通过最大更新策略,相应的1如图5所示.1阶盒型和星型Stencil的1可以计算得出:

(2)

Fig. 5 1
图5 1

结合01密铺迭代空间相当于将对应格点的值相加,即TB0(x,y)+TB1(x,y)=t-yTS0(x,y)+TS1(x,y)=t-x+y.图6展示了累加后的值.我们只用检查灰色三角形区域的结果,将01中相应灰色三角形区域中的值相加可以验证其结果正确性.

Fig. 6 Determining the center points of B2
图6 由01确定B2的中心点

通过相似的步骤确定图6(a)和图6(b)中所示的2个0点是B2的中心点.从B1B2的转换过程如图7(a)和图7(b)中虚线所示,由B2密铺数据空间的方式如图7(c)和图7(d)所示.盒型Stencil的密铺方式由12个B1转换为9个B2,星型Stencil由16个B1转换为13个B2.

Fig. 7 Tessellating data space with B2
图7 用B2块密铺数据空间

相似地,采用最大更新策略,相应的2如图8所示.1阶盒型和星型Stencil中2可以计算得出:

(3)

Fig. 8 Natural block B2(2)
图8 自然块B2(2)

结合0,1,2块密铺迭代空间相当于将对应格点值相加,由式(1)~(3)相加可得:

每个格点都更新t步,因此i(i=1,2,3)可以密铺盒型和星型Stencil迭代空间.

需要注意的是存在2种不同类型的1.以图4(a)(c)和图7(a)中的盒型Stencil为例,尽管所有的B1都是菱形的,但它们从2个不同方向产生,即△+▽→◇-和◁+▷→◇|,并且将被沿着与产生方向相垂直的方向划分,即◇-→◁+▷和◇|→△+▽.这2种类型分别对应于满足xy≥0或yx≥0的三角形块.第2种1的生成公式可以由交换xy得到.

2.3 高阶Stencil和边界条件

2.1和2.2节中针对1阶Stencil进行了讨论,2阶甚至更高阶Stencil计算在实际应用中也很常见.假设Stencil在1个维度的阶数为m,可以将沿着这一维度的m个格点组合成1个超格点.对于所有高阶依赖关系,我们都可以采用这种组合方式将m阶依赖减少为超格点间的1阶依赖.采用Bi密铺的方法可以应用到由超格点组成的1阶Stencil中,在每个i中,同一超格点中的所有格点更新相同的时间步数.

解决周期性边界条件问题时,我们可能需要在1维Stencil中对1个块进行拉伸,或者在高维Stencil中拉伸1组块.采用拉伸过的六边形块来保证在目标维度2边的2个子块可以组成1个自然块.

3 实现方法

基于第2节描述的数学框架可以手工实现代码.本节将介绍手工实现代码的具体细节.自动代码生成工具的实现将作为未来的工作.

3.1 粗粒度化

粗粒度化的实现有2个动机.首先,密铺方案需要在i的尖端处理Bi中较小的数据区域,特别是B0的开始和B2的结束.较差的硬件预取利用率或其他因素可能会限制硬件能力的发挥.其次,相邻的块可能会造成共享缓存上的伪共享.

星型密铺的粗粒度化和盒型密铺类似,但是菱形的大小只取决于1个参数即菱形边长,因此星型密铺中我们只需另外1个参数来将02的顶点扩展成1个小的菱形.图9展示了2维Stencil的粗粒度化方法.i中第1个(最后1个)时间步更新的格点区域被称为开始(结束)区域.为了提高空间局部性,0的结束区域(图9(a)中深灰色区域)和开始区域(图9(c)中黑色区域)被粗化为1个同样大小的2维块.粗粒度化使得每个阶段的每1个时间步都能处理1个3维块.例如,1块的结束位置是一个区域而不是1条对角线.

Fig. 9 Coarsening star tessellation
图9 粗粒度化星型密铺

3.2 02块间重用

如图9所示,B0B2块在数据空间中具有同样的形状.因此可以将相邻时间片的B0B2块合并,减少一次同步并增加数据重用.

此外,从图9中可以看出,2个B0在某一维度的距离等于B0结束区域的大小,这使得B0的开始区域和B2的结束区域完全相同.

3.3 1中2次更新技术

星型密铺一个重要的优势是可以在1中连续2次更新元素.从图5(b)中可以观察到,每个元素都将更新偶数次(TS1(x,y)=2y).基于这一特征,我们对每个元素连续更新2次,这种实现技术称为2次更新技术,下面将从1的更新过程中引出这个技术.1常规更新方法为

1==++
+.

等号右边的块中依次包含了1在连续4个时间步中每个时间步内更新的格点,即右边每个块中的格点在同一个时间步内更新.由于同一块中格点的更新次序是任意的,我们可以进一步将中间2个1的块分解成为

+,
+.

最后,将表示中间2个时间步更新格点的1替换为分解后的形式,并将最左块(表示第1个时间步更新格点)和分解后的第1个块结合,将最右块(表示第4个时间步更新格点)和分解后的最后1个块结合,得到的执行过程如下所示,其中每个格点连续更新2次,有效提升核内执行效率.

1==++.

3.4 并行化和向量化

在同一阶段的所有块可以并发执行,因此很容易实现并行化密铺方案.因此,我们只在共享内存机器上简单使用OpenMP编译指示parallel for.对于分布式内存计算机,清晰的密铺方案使得我们可以实现一个简单的数据/计算分布以及高效的数据通信,留做下一步工作.

星型密铺方案只适用于2维Stencil,而且B1的计算很难向量化,因此我们只对此方案在3维星型Stencil上的计算性能进行评估.因为3维Stencil中每个维度上的大小通常小于1 000,我们选择不划分单位跨度维度来利用硬件预取的能力,采用pragma simd对最内层循环进行向量化.

3.5 代码实现

代码实现主体包含5层循环,如算法1所示.最外层循环控制时间维度上时间片的切换,与分块相关的时间片大小为tBt0即为当前时间片的开始时间.内部包含2个4层循环,第1个循环用来进行B0块和B2块的更新,第2个循环更新B1块.

在第1个循环中,第1层循环根据时间片层数的奇偶遍历所有B0块或B2块.2结束时间维度即为下一个时间片0开始时间维度,且共用相同的数据空间,可以将02合并为02一起更新,除了迭代空间边界上的块,合并后的更新时间步长一般为2×tB.第2层循环遍历02所在的时间维度,t代表时间维度值,跨度为2×tB,但t的值会被时间维度的下界0和上界T所限定.求出02块在特定时间上的数据空间区域位置即菱形边界xmin,xmax,ymin,ymax后,在第3层和第4层循环中遍历更新02对应数据空间区域中的格点.

第2个循环的第1层循环遍历所有B1块,在星型密铺中,B1块是沿2个不同的对角线方向合并产生的,可以按照对角线方向的不同将B1块进一步细分为B11块和B12块,对应的dy值分别为1和-1,同时根据B1块的位置可以确定对应数据空间区域边界xmin,ymax.第2层循环遍历B1所在的时间维度,跨度为tB,并根据t计算特定时间维度上对角线区域边界xmax,ymin.第3层循环遍历对角线区域中每条对角线,确定每条对角线的起始坐标x,y和结束横坐标xr.第4层循环遍历某一条对角线上的所有格点,每个格点被连续更新2次.

算法1. 2D星型Stencil算法.

输入:时间步长tB;迭代空间上下界0,T;数据空间菱形边界xmin,xmax,ymin,ymax;数据空间(t,x,y);

输出:更新后的数据空间(t,x,y).

① for t0= -tBto T step tBdo

② #pragma omp parallel for private(xmin,

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

③ for n=0 to B0or B2do

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

T) do

set(xmin,xmax,ymin,ymax);

update(t,x,y);

/*在边界内更新B0B2*/

⑦ end for

⑧ end for

⑨ #pragma omp parallel for private(xmin,

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

⑩ for n=0 to B11 or B12 do

set(xmin,xmax,ymin,ymax,x,y,xr,dy);

while x<xrdo

update(t,x,y) twice;

/*2次更新B1块*/

end while

end for

end for

4 实验结果

4.1 实验环境

Fig. 10 Performance of size 2563
图10 问题规模为2563的性能

我们的实验运行在Intel Xeon E5-2670处理器上,处理器时钟频率为2.70 GHz,实验规模由单核扩展到所有12核.每个核独享32 KB L1 cache和256 KB L2 cache,12个核共享1个统一的30 MB L3 cache.使用版本号为16.0.1的ICC编译器编译程序并使用优化标志‘-O3 -openmp’.

将我们的方案与1.1节提及的2种高并发方案Pluto和Girih进行比较.盒型密铺、Pochoir和SDSL(Stencil domain specific language)[56]因其较差的性能结果(特别是3维Stencil)不作考虑.测试用例包括3种3维星型Stencil:1阶(3D7P)、2阶(3D13P)和4阶(3D25P)的具有非周期性边界条件的情况.测试的问题规模为2563×500,3843×800,5123×1 000.

Girih提供一种自动调优组件来选择不同线程数时的最优分块大小.对于Pluto和星型密铺方案,我们将测试所有可能的分块大小,然后选择12核时性能最优的分块大小并应用于所有1~12核测试中.

Fig. 11 Performance of size 3843
图11 问题规模为3843的性能

4.2 结果分析

图10~12分别显示了星型密铺(深色柱)、Girih(白色柱)和Pluto(浅色柱)在3D7P,3D13P,3D25P星型Stencil上大小为2563,3843,5123时的性能.其中性能GStencil的定义为:以最内层循环计算量看为一个Stencil,单位时间完成给定的Stencil迭代更新计算的次数为GStencil.由图10~12中可以明显看出,Pluto对应的性能条柱性能都比我们的代码和Girih低,且其结果不是一条直线,特别是在6核或更多核时.这和文献[29]中的结果一致.因此,下面主要分析星型密铺方案和Girih.

Fig. 12 Performance of size 5123
图12 问题规模为5123的性能

在问题规模大小为2563的3D7P星型Stencil中,星型密铺方案性能比Girih高,相比于Girih,星型密铺方案在2核时最高加速1.24倍,平均加速1.19倍.对于图10中更高阶的Stencil类型,我们的方案仍保持较高的性能,但随着阶数增加,这3种方法对应的性能柱变得更加接近,因为对于给定的自然块大小,更高阶的Stencil会产生更小的时间维度分块,造成更大的内存数据传输,后面将对此进行定量分析.

类似地,问题规模变大时3条性能曲线也会相互接近,也是因为Pluto和我们的密铺方案不对单位跨度的维度进行分块,导致时间维度分块大小受限.我们的方案在核数较少时仍取得了性能提升,在8~12核时虽然性能提升速度比Girih慢,但仍达到了可相比的性能.例如,在大小为3843的3D13P星型Stencil的实验结果中,星型密铺和Girih在12核上的性能分别为4.2,4.4 GStencil.

Fig. 13 Memory transfers of size 2563
图13 问题规模为2563的数据传输量

图13~15中显示了对应的最后一级cache和内存之间数据传输量.Girih利用整个L3 cache来最大限度开发数据局部性,因此具有最小的数据传输量.Girih在9个测试中的数据传输量都不规则变化,与密铺或Pluto的变化曲线不同,因为自动调优方案在不同核数时确定的分块大小不一样,例如,大小为2563的3D7P Stencil的Girih测量值在4核上是10.2 GBps,6核上是19.8 GBps,相应的菱形分块宽度最优值分别为32和16.Girih事实上是菱形分块和波前并行的混合方法,并采用细粒度并行,所有线程同时处理一个块,因此我们不对它进行理论分析.

Fig. 14 Memory transfers of size 3843
图14 问题规模为3843的数据传输量

Fig. 15 Memory transfers of size 5123
图15 问题规模为5123的数据传输量

对于3D25P Stencil,星型密铺方案在3个问题规模上的最优分块大小分别是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)中星型密铺方案的数据传输量开始增长的位置.尽管核数更大时超出了L3 cache大小,但由于高阶依赖,数据块将快速收缩.此外,用时间维度块大小为2的分块密铺方案等价于原始算法,因此它们产生同样的缓存-内存传输量即T×N2.Pluto在3D25P Stencil上也有类似地变化趋势和原因.

对于另外2个低阶Stencil,L3 cache能保留所有12个核的自然块.例如大小为2563的3D7P Stencil中自然块占用内存为2.5 MB,12核将占用所有30 MB L3 cache.因此,星型密铺的带宽变化接近于水平.但是性能在核数变多时增长变缓,这可能是由于有限的L3带宽.

最后比较盒型密铺的性能.总体上看盒型密铺和Girih性能相近.例如在大小为2563的3D7P Stencil中,盒型密铺的最高性能是6.9 GStencils而Girih是7.0 GStencil.盒型密铺对应的块大小是23×11,约为星型密铺最优块大小33×16的倍.分块大小使得这2种方案都在12核时完全占用L3 cache.星型和盒型密铺数据传输量的测量值分别为19.3 GBps和26.4 GBps,同样有的关系.

4 总 结

本文提出了一个新的并行框架,能够高并发执行星型Stencil,并以“自然块”的新概念来标识不同形状Stencil的特征.自然块及其后继块能够密铺数据空间,并且它们沿着时间维度扩展后能形成迭代空间的密铺.我们的方案能够导出简洁的并行框架,揭示一个格点的坐标和它更新时间步之间的关系.此外,专为星型Stencil研发了一个新颖的实现方式即“2次更新”,能够在第1个后继块的执行过程中对每个元素更新2次,从而改进核内数据重用模式.实验结果证明了此方法的有效性.

未来我们将着重于自动调优方法来高效寻找最优块大小,基于提出的框架设计一个工具来自动生成Stencil代码.

参考文献

[1]Ding C, He Yun. A ghost cell expansion method for reducing communications in solving pde problems[C] //Proc of the ACM/IEEE SC’01. New York: ACM, 2001: 55-55

[2]Nguyen A, Satish N, Chhugani J, et al. 3.5-d blocking optimization for Stencil computations on modern CPUs and GPUs[C] //Proc of the ACM/IEEE SC’10. New York: ACM, 2010: 1-13

[3]Rastello F, Dauxois T. Efficient tiling for an ode discrete integration program: Redundant tasks instead of trapezoidal shaped-tiles[C] //Proc of the IEEE IPDPS’02. Piscataway, NJ: IEEE, 2002: 1-8

[4]Rivera G, Tseng C W. Tiling optimizations for 3d scientific computations[C] //Proc of the 2000 ACM/IEEE Conf on Supercomputing. Piscataway, NJ: IEEE, 2000: 32-32

[5]Christen M, Schenk O, Burkhart H. Patus: A code generation and autotuning framework for parallel iterative Stencil computations on modern microarchitectures[C] //Proc of the 25th IEEE Int Symp on Parallel and Distributed Processing. Piscataway, NJ: IEEE, 2011: 676-687

[6]Gysi T, Grosser T, Hoefler T. Modesto: Data-centric analytic optimization of complex Stencil programs on heterogeneous architectures[C] //Proc of the 29th ACM Int Conf on Supercomputing. New York: ACM, 2015: 177-186

[7]Kamil S, Chan C, Oliker L, et al. An auto-tuning framework for parallel multicore Stencil computations[C] //Proc of 2010 IEEE Int Symp on Parallel & Distributed Processing. Piscataway, NJ: IEEE, 2010: 1-12

[8]Luo Yulong, Tan Guangming, Mo Zeyao, et al. Fast: A fast Stencil autotuning framework based on an optimal-solution space model[C] //Proc of the 29th ACM on Int Conf on Supercomputing. New York: ACM, 2015: 187-196

[9]Lutz T, Fensch C, Cole M. Partans: An autotuning framework for Stencil computation on multi-GPU systems[J]. ACM Transactions on Architecture and Code Optimization, 2013, 9(4): No.59

[10]Mametjanov A, Lowell D, Ma C C, et al. Autotuning Stencil-based computations on GPUs[C] //Proc of 2012 IEEE Int Conf on Cluster Computing. Piscataway, NJ: IEEE, 2012: 266-274

[11]Rocha R C O, Pereira A D, Ramos L, et al. Toast: Automatic tiling for iterative Stencil computations on GPUs[J]. Concurrency and Computation: Practice and Experience, 2017 , 29(8): No.e4053

[12]Zhang Yongpeng, Mueller F. Auto-generation and auto-tuning of 3d Stencil codes on GPU clusters[C] //Proc of the 10th Int Symp on Code Generation and Optimization. Piscataway, NJ: IEEE, 2012: 155-164

[13]Maruyama N, Nomura T, Sato K, et al. Physis: An implicitly parallel programming model for Stencil computations on large-scale GPU-accelerated supercomputers[C] //Proc of the 2011 Conf for High Performance Computing, Networking, Storage and Analysis. New York: ACM, 2011: No.11

[14]Unat D, Cai Xing, Baden S B. Mint: Realizing cuda performance in 3d Stencil methods with annotated c[C] //Proc of the 25th Int Conf on Supercomputing. New York: ACM, 2011: 214-224

[15]Datta K, Murphy M, Volkov V, et al. Stencil computation optimization and auto-tuning on state-of-the-art multicore architectures[C] //Proc of the 2008 ACM/IEEE Conf on Supercomputing. Piscataway, NJ: IEEE, 2008: No.4

[16]Prieto M, Llorente I M, Tirado F. Data locality exploitation in the decomposition of regular domain problems[J]. IEEE Transactions on Parallel and Distributed Systems, 2000, 11(11): 1141-1150

[17]Venkatasubramanian S, Vuduc R W. Tuned and wildly asynchronous Stencil kernels for hybrid CPU/GPU systems[C] //Proc of the 23rd Int Conf on Supercomputing. New York: ACM, 2009: 244-255

[18]Williams S, Oliker L, Carter J, et al. Extracting ultra-scale lattice boltzmann performance via hierarchical and distributed auto-tuning[C] //Proc of the ACM/IEEE SC’11. New York: ACM, 2011: No.55

[19]Christen M, Schenk O, Neufeld E, et al. Parallel data-locality aware Stencil computations on modern micro-architectures[C] //Proc of the IEEE IPDPS’09. Piscataway, NJ: IEEE, 2009

[20]Krotkiewski M, Dabrowski M. Efficient 3d Stencil computations using CUDA[J]. Parallel Computing, 2013, 39(10): 533-548

[21]Olschanowsky C, Strout M M, Guzik S, et al. A study on balancing parallelism, data locality, and recomputation in existing pde solvers[C] //Proc of the 2014 Int Conf for High Performance Computing, Networking, Storage and Analysis. Piscataway, NJ: IEEE, 2014: 793-804

[22]Phillips E H, Fatica M. Implementing the himeno benchmark with CUDA on GPU clusters[C] //Proc of the IEEE IPDPS’10. Piscataway, NJ: IEEE, 2010

[23]Rawat P S, Hong Changwan, Ravishankar M, et al. Effective resource management for enhancing performance of 2d and 3d Stencils on GPUs[C] //Proc of the 9th Annual Workshop on General Purpose Processing Using Graphics Processing Unit. New York: ACM, 2016: 92-102

[24]Heybrock S, Joo B, Kalamkar D D, et al. Lattice qcd with domain decomposition on Intel xeon phi co-processors[C] //Proc of the ACM/IEEE SC’14. New York: ACM, 2014: 69-80

[25]Yount C, Duran A. Effective use of large high-bandwidth memory caches in HPC Stencil computation via temporal wave-front tiling[C] //Proc of the 7th Int Workshop on Performance Modeling, Benchmarking and Simulation of High Performance Computer Systems. Piscataway, NJ: IEEE, 2016: 65-75

[26]Bondhugula U, Hartono A, Ramanujam J, et al. A practical automatic polyhedral parallelizer and locality optimizer[C] //Proc of the 29th ACM SIGPLAN Conf on Programming Language Design and Implementation. New York: ACM, 2008: 101-113

[27]Bondhugula U, Bandishti V, Cohen A, et al. Tiling and optimizing time-iterated computations on periodic domains[C] //Proc of the 23rd Int Conf on Parallel Architecture and Compilation Techniques. New York: ACM, 2014: 39-50

[28]Krishnamoorthy S, Baskaran M, Bondhugula U, et al. Effective automatic parallelization of Stencil computations[C] //Proc of the ACM SIGPLAN 2007 Conf on Programming Language Design and Implementation. New York: ACM, 2007: 235-244

[29]Bandishti V, Pananilath I, Bondhugula U. Tiling Stencil computations to maximize parallelism[C] //Proc of the Int Conf on High Performance Computing Networking, Storage and Analysis. Piscataway, NJ: IEEE, 2012: No.40

[30]Frigo M, Strumpen V. Cache oblivious Stencil computations[C] //Proc of the 19th Annual Int Conf on Supercomputing. New York: ACM, 2005: 361-366

[31]Frigo M, Strumpen V. The cache complexity of multithreaded cache oblivious algorithms[C] //Proc of the 18th Annual ACM Symp on Parallelism in Algorithms and Architectures. New York: ACM, 2006: 271-280

[32]Strzodka R, Shaheen M, Pajak D, et al. Cache oblivious parallelograms in iterative Stencil computations[C] //Proc of the 24th ACM Int Conf on Supercomputing. New York: ACM, 2010: 49-59

[33]Tang Yuan, Chowdhury R A, Kuszmaul B C, et al. The pochoir Stencil compiler[C] //Proc of the 23rd Annual ACM Symp on Parallelism in Algorithms and Architectures. New York: ACM, 2011: 117-128

[34]Yuan Liang, Zhang Yunquan, Guo Peng, et al. Tessellating Stencils[C] //Proc of the Int Conf for High Performance Computing, Networking, Storage and Analysis. New York: ACM, 2017: No.49

[35]Irigoin F, Triolet R. Supernode partitioning[C] //Proc of the 15th Annual ACM Symp on Principles of Programming Languages. New York: ACM, 1988: 319-329

[36]Lam M D, Rothberg E E, Wolf M E. The cache performance and optimizations of blocked algorithms[C] //Proc of the 4th Int Conf on Architectural Support for Programming Languages and Operating Systems. New York: ACM, 1991: 63-74

[37]McKellar A C, Coffman E G Jr. Organizing matrices and matrix operations for paged memory systems[J]. Communications of the ACM, 1969, 12(3): 153-165

[38]Wolf M E, Lam M S. A data locality optimizing algorithm[C] //Proc of the ACM SIGPLAN Conf on Programming Language Design and Implementation. New York: ACM, 1991: 30-44

[39]Wolfe M. More iteration space tiling[C] //Proc of the 1989 ACM/IEEE Conf on Supercomputing. New York: ACM, 1989: 655-664

[40]Wonnacott D G, Strout M M. On the scalability of loop tiling techniques[C] //Proc of the 3rd Int Workshop on Polyhedral Compilation Techniques. New York: ACM, 2013: 3-11

[41]Meng J, Skadron K. Performance modeling and automatic ghost zone optimization for iterative Stencil loops on GPUs[C] //Proc of the 23rd Int Conf on Supercomputing. New York: ACM, 2009: 256-265

[42]Holewinski J, Pouchet L N, Sadayappan P. High-performance code generation for Stencil computations on GPU architectures[C] //Proc of the 26th ACM Int Conf on Supercomputing. New York: ACM, 2012: 311-320

[43]Demmel J, Hoemmen M F, Mohiyuddin M, et al. Avoiding communication in computing krylov subspaces[OL].[2020-02-10]. http://citeseerx.ist.psu.edu/

[44]Mohiyuddin M, Hoemmen M, Demmel J, et al. Minimizing communication in sparse matrix solvers[C] //Proc of the Conf on High Performance Computing Networking, Storage and Analysis. New York: ACM, 2009: No.36

[45]Jin Guohua, Mellor-Crummey J, Fowler R. Increasing temporal locality with skewing and recursive blocking[C] //Proc of the 2001 ACM/IEEE Conf on Supercomputing. New York: ACM, 2001: No.43

[46]Song Yonghong, Li Zhiyuan. New tiling techniques to improve cache temporal locality[C] //Proc of the 1999 ACM SIGPLAN Conf on Programming Language Design and Implementation. New York: ACM, 1999: 215-228

[47]Wonnacott D. Achieving scalable locality with time skewing[J]. International Journal of Parallel Programming, 2002, 30(3): 181-221

[48]Malas T, Hager G, Ltaief H, et al. Multicore-optimized wavefront diamond blocking for optimizing Stencil updates[J]. SIAM Journal on Scientific Computing, 2015, 37(4): C439-C464

[49]Wellein G, Hager G, Zeiser T, et al. Efficient temporal blocking for Stencil computations by multicore-aware wavefront parallelization[C] //Proc of the 33rd Annual IEEE Int Computer Software and Applications Conf. Piscataway, NJ: IEEE, 2009: 579-586

[50]Pananilath I, Acharya A, Vasista V, et al. An optimizing code generator for a class of lattice-boltzmann computations[J]. ACM Transactions on Architecture & Code Optimization, 2015, 12(2): No.14

[51]Grosser T, Verdoolaege S, Cohen A, et al. The relation between diamond tiling and hexagonal tiling[J]. Parallel Processing Letters, 2014, 24(3): No.1441002

[52]Cui Huimin, Wang Lei, Fan Dongrui, et al. Landing Stencil code on godson-t[J]. Journal of Computer Science and Technology, 2010, 25(4): 886-894

[53]Wonnacott D. Using time skewing to eliminate idle time due to memory bandwidth and network limitations[C] //Proc of the 14th Int Parallel and Distributed Processing Symposium. Piscataway, NJ: IEEE, 2000: 171-180

[54]Zhou Xing, Giacalone J P, Garzarán M J, et al. Hierarchical overlapped tiling[C] //Proc of the 10th Int Symp on Code Generation and Optimization. New York: ACM, 2012: 207-218

[55]Grosser T, Cohen A, Kelly P H J, et al. Split tiling for GPUs: Automatic parallelization using trapezoidal tiles[C] //Proc of the 6th Workshop on General Purpose Processor Using Graphics Processing Units. New York: ACM, 2013: 24-31

[56]Henretty T, Veras R, Franchetti F, et al. A Stencil compiler for short-vector simd architectures[C] //Proc of the 27th Int ACM Conf on Int Conf on Supercomputing. New York: ACM, 2013: 13-24

[57]Strzodka R, Shaheen M, Pajak D, et al. Cache accurate time skewing in iterative Stencil computations[C] //Proc of 2011 Int Conf on Parallel Processing. Piscataway, NJ: IEEE, 2011: 571-581

[58]Grosser T, Cohen A, Holewinski J, et al. Hybrid hexagonal/classical tiling for GPUs[C] //Proc of Annual IEEE/ACM Int Symp on Code Generation and Optimization. New York: ACM, 2014: No.66

[59]Basu P, Hall M, Williams S, et al. Compiler-directed transformation for higher-order Stencils[C] //Proc of 2015 IEEE Int Parallel and Distributed Processing Symp. Piscataway, NJ: IEEE, 2015: 313-323

[60]Deitz S J, Chamberlain B L, Snyder L. Eliminating redundancies in sum-of-product array computations[C] //Proc of the 15th Int Conf on Supercomputing. New York: ACM, 2001: 65-77

[61]de la Cruz R, Araya-Polo M. Algorithm 942: Semi-Stencil[J]. ACM Transactions on Mathematical Software, 2014, 40(3): 23:1-23:39

[62]Jin Mengyao, Fu Haohuan, Lü Zihong, et al. Libra: An automated code generation and tuning framework for register-limited Stencils on GPUs[C] //Proc of the ACM Int Conf on Computing Frontiers. New York: ACM, 2016: 92-99

[63]Sedaghati N, Thomas R, Pouchet L N, et al. Stvec: A vector instruction extension for high performance Stencil computation[C] //Proc of 2011 Int Conf on Parallel Architectures and Compilation Techniques. Piscataway, NJ: IEEE. 2011: 276-287

[64]Zumbusch G. Vectorized higher order finite difference kernels[C] //Proc of 2012 Int Workshop on Applied Parallel Computing. Berlin: Springer, 2012: 343-357

[65]Hadade I, Mare L. Exploiting simd and thread-level parallelism in multiblock CFD[C] //Proc of the 29th Int Supercomputing Conf. Berlin: Springer, 2014: 410-419

[66]Henretty T, Stock K, Pouchet L N, et al. Data layout transformation for Stencil computations on short-vector simd architectures[C] //Proc of 2011 Int Conf on Compiler Construction. Berlin: Springer, 2011: 225-245

A Parallel Star Stencil Algorithm Based on Tessellating

Cao Hang1,2, Yuan Liang1, Huang Shan1,2, Zhang Yunquan1, Xu Yongjun1, Lu Pengqi1,2, and Zhang Guangting1

1(State Key Laboratory of Computer Architecture (Institute of Computing Technology, Chinese Academy of Sciences), Beijing 100190)

2(University of Chinese Academy of Sciences, Beijing 100049)

Abstract The Stencil computation, i.e. the structured grid computing, is a very common kind of loop nesting algorithm in scientific and engineering applications. The exhaustively studied tiling method is of great effectiveness as one of the transformation techniques to exploit the data locality and parallelism of Stencil computations. However, the state-of-the-art work of tiling often uniformly handles different Stencil shapes. We first present a concept called natural block to identify the difference between the star and box Stencils. Then we propose a new two-level tessellation scheme for star Stencils, where the natural block, as well as its successive blocks can tessellate the spatial space and their extensions along the time dimension are able to form a tessellation of the iteration space. Furthermore, a novel implementation technique called double updating is developed for star Stencils specifically, which updates each element twice continuously and improves the in-core data reuse pattern. In addition, we adopt coarsening and block reuse to enhance the parallelization performance. Theoretical analysis shows that our scheme achieves a better cache complexity than existing methods such as Girih and Pluto. The experiments on performance and bandwidth are conducted on a multicore system. The results demonstrate the effectiveness of our approach.

Key words Stencil computation; tessellation; star Stencil; box Stencil; natural block

(caohang@ict.ac.cn)

中图法分类号 TP301

DOI:10.7544/issn1000-1239.2020.20190734

收稿日期2019-10-12;修回日期:2020-04-07

基金项目国家重点研发计划项目(2016YFB0200803);中国科学院战略性先导科技专项(C类)(XDC01040100); 国家自然科学基金项目(61972376,62072431,61432018);北京市自然科学基金项目(L182053)

This work was supported by the National Key Research and Development Program of China (2016YFB0200803), the Strategic Priority Research Program of Chinese Academy of Sciences (C) (XDC01040100), the National Natural Science Foundation of China (61972376, 62072431, 61432018), and the Beijing Natural Science Foundation (L182053).

通信作者袁良(yuanliang@ict.ac.cn)

Cao Hang, born in 1992. PhD from the Institute of Computing Technology, Chinese Academy of Sciences, Beijing. His main research interests include large-scale parallel computing and heterogeneous parallel computing performance optimization.

Yuan Liang, born in 1984. PhD, assistant professor with the State Key Laboratory of Computer Architecture, Institute of Computing Technology, Chinese Academy of Sciences, Beijing, China. Member of CCF. His main research interests include large-scale parallel computing and heterogeneous computing.

Huang Shan, born in 1994. MSc degree in the State Key Laboratory of Computer Architecture, Institute of Computing Technology, Chinese Academy of Sciences, Beijing, China. Her main research interests include parallel computing model and parallel computing algorithm performance optimization.

Zhang Yunquan, born in 1973. PhD, professor of computer science with the Institute of Computing Technology, Chinese Academy of Sciences, Beijing, China.Senior member of CCF. His main research interests include high performance parallel computing, big data processing.

Xu Yongjun, born in 1979. PhD, PhD supervisor, professor level senior engineer with the Institute of Computing Technology, Chinese Academy of Sciences, Beijing, China. His main research interests include big data processing, artificial intelligence, IoT.

Lu Pengqi, born in 1995. Master candidate in the Institute of Computing Technology, Chinese Academy of Sciences, Beijing. His main research interests is high performance computing.

Zhang Guangting, born in 1987. Master, assistant professor with the State Key Laboratory of Computer Architecture, Institute of Computing Technology, Chinese Academy of Sciences, Beijing, China. Her main research interest is high performance computing.