一种利用时间局部性的格子玻尔兹曼方法并行加速方法与流程

本发明属于计算机技术领域，特别涉及一种利用时间局部性的格子玻尔兹曼方法并行加速方法。

背景技术：

格子玻尔兹曼方法(lbm)是一种应用于计算流体动力学的数值仿真方法。它基于量子统计力学的理论，将流体看做大量有质量无体积的微粒，可以通过统计这些微粒的碰撞规律来描述流体的真实流动。格子玻尔兹曼方法主要由格子(离散速度模型)、平衡态分布函数、分布函数的演化方程组成，采用不同的格子、平衡态分布函数以及演化方程，一般便能解决各式各样的流体问题。

lbgk模型是一种针对普通湍流问题求解的经典单松弛模型，采用d2q9格子的lbgk模型经常用于格子玻尔兹曼方法的加速研究。采用d2q9格子的lbm，每个微粒要和二维空间上相邻的8个微粒进行碰撞和速度交换，演化过程可分为碰撞(collision)和流动(streaming)两步，再加上宏观物理量的更新这一步，从执行角度来看，单时间迭代步内一共具有三个多维嵌套循环，每个多维嵌套循环分别构成一个单独的循环域。

格子玻尔兹曼方法具有显式并行的特点，单时间迭代步计算过程中三个多维嵌套循环均各自满足doall循环的特征。doall循环是指不携带跨迭代依赖的循环(所有的迭代都可以与其他迭代进行完全并行，具有数据无关性)，因此这些循环可以采取并行计算。由于格子玻尔兹曼方法适合在大规模计算集群上实现并行计算，因此广泛应用于多种场景的计算流体动力学相关模拟仿真。虽然在doall循环上直接进行简单并行就能够很大程度地提升计算速度，但是原始的格子玻尔兹曼方法仍然具有诸多优化的途径，能够进一步提高计算速度。

多年来计算机cpu速度的高速增长和内存速度的缓慢增长使得cpu与内存之间的速度差距越来越大，导致了内存数据访问成为性能瓶颈。现代计算机引入多级cache的存储层次来降低内存访问的性能瓶颈，cache的层次越高越接近cpu，cpu对其的访问速度就越快。cache优化利用数据局部性原理为程序的高速运行提供了有效的支持，将访问次数较多的较小工作集放入高层次cache(私有l1/l2cache)，以减少cache失效，保持计算性能。重用距离是度量cache局部性的一个指标，其意义是串行程序运行时连续两次访问同一数据变量之间所访问的不同数据变量的数目，重用距离越大，该数据变量的局部性越差。数据依赖是程序执行数据流所产生的执行次序约束，在动作b发生之前，动作a必须先发生，称为b依赖于a。对于一个二维循环而言，若迭代实例(i0，j0)中语句s1依赖于迭代实例(i1，j1)中语句s2，则定义依赖距离向量为(i0-i1，j0-j1)。合法的依赖距离向量的值必须按词典序大于0。循环变换技术是指在不破坏程序正确执行结果的条件下，改变原始循环的执行次序的一种程序优化技术，是循环偏斜、循环融合、循环分块等一系列技术的合称。循环偏斜是指通过内层循环相对于外层循环的偏斜，来消除某些阻止其他循环变换的依赖条件，循环偏斜并不能改变循环的执行次序，对局部性也没有任何增强。循环融合是将两个及以上的循环合并为一个循环，用以增强循环体大小，增强局部性，消除某些阻止其他循环变换的状态。循环分块是将原始的迭代空间进行划分，使工作集减小，从而增强局部性，减少cache失效。同时，程序的并行化为了增大并行粒度也需要将循环进行分块。循环分块要求循环交换后的新循环所携带的依赖距离向量必须保持词典序为正。

在共享内存架构上，针对lbm的传统优化方法主要依靠对单时间迭代步内的三个doall循环(碰撞、流动、宏观量计算)分别进行循环分块，依次并行计算这三个doall循环，计算完成后进入下一个时间迭代步，重复计算过程，直至计算结束。传统优化方法采用循环分块主要是为了实现数据分布，而且因为三个doall循环(碰撞、流动、宏观量计算)各自的可供开发的局部性很少，循环分块几乎得不到局部性收益。

传统优化方法未曾考虑对最外层的时间维进行利用。因为lbm的最外层循环是时间迭代步，也就意味着每次时间迭代内访问到内层循环的内存地址序列是一样的。如果一个变量在时间迭代步tn内被访问到，它将被存储在cache中，当到达时间迭代步tn+1时，如果该变量还没有被替换出cache，则处理器可以直接从cache中访问到该变量，形成一次重用。但是，因为一次时间迭代步内将访问到内层循环的整个内存地址序列，访问的不同变量的数量很多，导致重用距离大于cache容量，变量在被时间迭代步tn+1重用之前就会被替换出cache，导致cache失效。循环分块可以用于改变循环的内存地址访问次序，从而缩小重用距离。将循环分块应用在lbm上，一次时间迭代内将不再访问到整个内层的内存地址序列，而是一个分块内的内存地址序列，选择合适分块大小，可以使得一个分块内的所有的变量都可以一直存储在cache中，这样的话，下一次时间迭代再次访问这个分块的内存地址序列都不会产生cache失效。利用这种方法可以开发程序的时间局部性优势，获得加速。

技术实现要素：

本发明的目的在于提供一种利用时间局部性的格子玻尔兹曼方法并行加速方法，提高格子玻尔兹曼方法的计算速度。

本发明是通过以下技术方案来实现：

一种利用时间局部性的格子玻尔兹曼方法并行加速方法，包括步骤：

1)将单时间迭代步内的三个空间维的doall循环融合成一个doacross循环；

2)对融合后的doacross循环执行循环偏斜，消除与时间维相关的负依赖，形成融合时间维的doacross循环；

3)使用循环分块技术，对融合时间维的doacross循环执行循环分块，形成分块大小为a×a×t的多个分块；

4)对分块实现波阵面并行。

优选地，在步骤1)中，取偏斜因子为1对存在负依赖的三个空间维的doall循环执行循环偏斜以消除负依赖，并实现循环融合。

优选地，在步骤2)中，取偏斜因子为2对融合后的doacross循环执行循环偏斜，消除与时间维相关的负依赖，形成融合时间维的doacross循环。

优选地，在步骤3)中，执行循环分块时，分块大小为a×a×t，其中，a为空间维分块因子，t为时间维分块因子；

空间维分块因子a可以根据以下公式进行计算并取整：

cache_capacity：私有l1/l2cache的容量，

num_of_data：最内层循环单次迭代所访问的变量的数目，

datatype：变量所占字节数。

优选地，在步骤3)中，融合时间维的doacross循环采用矩形分块的方式执行循环分块，得到多个分块。

优选地，在步骤4)中，步骤3)所得分块采用post/wait操作实现波阵面并行。

优选地，在步骤4)中，post/wait操作的方式为：线程t2被分配给子块n2时，先执行wait操作，等待接收执行上一个波阵面对应子块n1的线程t1发出的信号；待接收到线程t1发出的信号之后，线程t2立即执行子块n2，计算结束后，线程t2执行post操作，发送信号给下一个波阵面的对应子块n3，以启动处于等待的线程t3执行子块n3。

优选地，所述格子玻尔兹曼方法为采用d2q9模型的格子玻尔兹曼方法。

优选地，所述方法应用在共享内存架构上。

与现有技术相比，本发明具有以下有益的技术效果：

本发明提供的利用时间局部性的格子玻尔兹曼方法并行加速方法，其先通过循环变换技术将原有的单时间迭代步内三个空间维的doall循环融合成单时间步迭代的单个doacross循环，然后通过执行循环偏斜消除与时间维相关的负依赖，形成融合时间维的doacross循环，然后通过循环分块技术，形成多个分块，并对这些分块实现波阵面并行。本发明通过循环变换技术消除负依赖、合并迭代空间、实现循环分块，将原有的单时间迭代步内三个doall循环合并成多时间步迭代的单个doacross循环，减少了循环迭代次数，增加了数据的时间局部性，并实现了分块间的可并行执行。

进一步地，采用post/wait操作对这些分块实现波阵面并行，避免了隐式栅栏同步时线程的等待，从而加速了并行计算的过程。计算速度快于传统方法，最终达到了并行计算加速的目的。

附图说明

图1为跨时间迭代步引起的负依赖距离向量示意图。

图2为波阵面并行执行次序示意图

图3为post/wait操作的信号传递示意图

具体实施方式

下面结合具体的实施例对本发明做进一步的详细说明，所述是对本发明的解释而不是限定。

一种利用时间局部性的格子玻尔兹曼方法并行加速方法，包括以下步骤：

1)将单时间迭代步内的三个空间维doall循环融合成一个doacross循环。原始代码结构如表1所示，单时间迭代步中存在三个空间维的嵌套doall循环，由于这三个循环的迭代空间都是一样的，因此可以将三个循环的循环体合并执行，这样可以减少两倍的循环迭代次数，同时为后续循环优化提供唯一的循环域。

表1原始代码结构

其中，第二个循环体要消费第一个循环体的计算结果，第三个循环体要消费第二个循环体的计算结果，因此三个循环之间存在数据依赖。对存在数据依赖的多个循环进行融合，需要满足后一个循环体在某个迭代实例(i0，j0)时所消费的数据是前一个循环体已经计算得到的数据这一条件。由于d2q9格子玻尔兹曼方法一个迭代实例的变量要消费相邻共9个迭代实例的变量，考虑最坏情况，在相同时间迭代步下，第二个循环体的迭代实例为(i0，j0)时需要消费第一个循环体的迭代实例为(i0+1，j0+1)所产生的数据。因此要实现循环融合，第二个循环体的迭代实例(i0，j0)必须依赖于第一个循环体本应该随后才执行的迭代实例(i0+1，j0+1)，从而产生最大为(-1，-1)的依赖距离向量。要消除负依赖(即保证依赖距离向量满足词典序为正，消除分量的负值)才能正确进行循环融合，因此取偏斜因子(skewfactor)为1分别对两个空间维进行循环偏斜。消除负依赖之后，三个doall循环被融合成一个doacross循环，其代码结构如表2所示：

表2循环融合后的代码结构

2)对融合了三个doall循环的doacross循环执行循环偏斜，消除与时间维相关的负依赖。该doacross循环仍然限制在单时间迭代步内，要融合多个时间迭代步来分块执行doacross循环，首先需要解决各时间迭代步之间的负依赖。由于融合时间迭代步的循环进行循环分块后，整个循环的执行顺序并非按照最外层的时间迭代步的顺序来执行，而是对包含多个时间迭代步的分块按照分块间的数据依赖关系先后执行。如图1所示最坏情况，如果某个迭代实例(t0，i0-2，j0-2)与(t0-1，i0，j0)分别从属于先执行的分块n1和后执行的分块n2，那么分块n1在t0时对f[i0-2][j0-2]进行的写操作理论上要早于分块n2在t0-1时对f[i0-2][j0-2]进行的读操作，从而覆盖掉完成分块n2的计算所依赖的数据。从程序的计算逻辑来说，这个读写顺序将子块n2还没来得及读的数据覆盖掉了，出现“写后读(read-after-write)”的错误，计算结果会受影响而全部出错。为了避免这种现象，这两个不同时间迭代步上的迭代需要放在同一个分块中。经过依赖分析可知，此时循环的两个空间维相对于时间维存在最坏为(-2，-2)的依赖距离向量。为了消除与时间维相关的负依赖，取偏斜因子为2对此时的doacross循环执行循环偏斜，得到代码结构如表3所示的融合了时间维的doacross循环。

表3循环偏斜后的代码结构

3)使用循环分块技术，对融合时间维的doacross循环执行循环分块。我们采用矩形分块的方式，将三维循环的迭代实例划分为多个长方体分块(迭代实例的标准数量为a×a×t，其中a为空间维分块因子，即空间维上迭代次数，t为时间维分块因子，即时间维上迭代次数)。需要注意的是，由于代码转换过程执行了多次循环偏斜，循环边界上的分块的迭代实例数量并没有达到标准值。执行循环分块后，循环代码结构如表4所示。

表4循环分块后的代码结构

由于在不同时间迭代步的相同空间上迭代实例的内存地址相同，因此分块的工作集大小只和分块在两个空间维的分块因子a相关，时间维分块因子t不会影响工作集大小，没有强制选择要求。选择私有l1/l2cache的容量作为子块的数据容量极限，因此分块的工作集大小为a×a，分块因子a可以根据以下公式(1)进行计算并取整：

其中，

cache_capacity：私有l1/l2cache的容量，

num_of_data：最内层循环单次迭代所访问的变量的数目，

datatype：变量所占字节数。

4)采用post/wait操作对循环分块后的分块实现波阵面并行。由于各分块之间存在一定的数据依赖，需要分块按照一定的先后执行顺序才能实现有限的并行，分块内部的迭代实例则串行执行。由依赖关系可知传统波阵面并行的执行次序，如图2所示，箭头表示依赖关系和执行次序，正方形表示分块，数字表示执行此分块的线程号，虚线表示隐式同步栅栏。由于步骤3)完成循环分块时已经相对于时间维进行了偏斜，消除了时间维的影响，且工作集只与空间维迭代相关，因此示意图只需要表示分块在空间维方向的执行次序情况便能体现整个融合时间维的doacross循环分块的执行情况。假设只能提供4个线程，最开始的时候，只有一个分块可以执行，该分块执行完之后，能为空间维的迭代推进方向上的两个分块提供依赖数据，这两个分块可以并行执行。当这两个分块执行结束，两个分块可以一起为空间维的推进方向上的三个分块提供依赖数据，这三个分块也可以并行执行。以此类推，直至达到最大线程数的并行，每一次分块的并行执行形成一个波阵面，这种并行叫做波阵面并行。但是传统的波阵面并行存在一定的隐式栅栏同步开销，处于当前波阵面的线程执行完一个分块之后不能立即执行下一个波阵面的分块，而是需要等待当前波阵面所有分块都被执行完成后才能开始对下一个波阵面分块的执行，因此采用post/wait操作来对传统波阵面并行进行优化。post/wait操作的信号传送如图3所示，箭头代表信号量的传送，正方形表示分块，数字表示执行此分块的线程号。假设只能提供四个线程，post/wait操作要求一个线程t2被分配给分块n2时，先执行wait操作，等待接收执行上一个波阵面对应分块n1的线程t1发出的信号。待接收到这个信号之后，线程t2立即执行分块n2，计算结束后，线程t2执行post操作，发送信号给下一个波阵面的对应分块n3，以启动处于等待的线程t3执行分块n3。由于post/wait操作允许线程执行完当前波阵面一个分块之后立即被分配给下一个波阵面一个已经接收到信号的未执行分块，避免了隐式栅栏同步时线程的等待，从而加速了并行计算的过程。表5给出了对经过步骤1)、2)、3)获得的分块应用基于post/wait操作的波阵面并行同步以实现d2q9格子玻尔兹曼方法加速并行的具体操作。

表5基于post/wait操作的波阵面并行同步策略

经过实验发现，在intel服务器上(2×6核xeone5645，2.4ghz，12gb内存，32kb的l1cache，256kb的l2cache，使用openmp多线程)，采用a＝40，t＝200的分块因子，本发发明提供的并行加速方法与原并行方法相比，在执行不同线程数(4,6,8,10,12)和不同问题规模(二维网格大小500×500,1000×1000,2000×2000,3000×3000,4000×4000,5000×5000)时，计算速度性能均能获得提升，平均提升11％的计算速度性能。尤其是当二维网格大小4000×4000、10线程)时，可以获得最大性能提升，达到25％。