CN102200962A

CN102200962A - 一种基于迭代空间条块的并行有限差分模版方法

Info

Publication number: CN102200962A
Application number: CN2011100512858A
Authority: CN
Inventors: 张纪林; 徐向华; 万健; 蒋从锋; 张伟; 任永坚
Original assignee: Hangzhou Dianzi University
Current assignee: Hangzhou Dianzi University
Priority date: 2011-07-25
Filing date: 2011-07-25
Publication date: 2011-09-28

Abstract

本发明涉及一种基于迭代空间条块的并行有限差分模版方法。传统的并行化迭代方法在迭代内和迭代间都需要同步操作以维护数据依赖关系。本发明首先通过对迭代空间进行时间轴方向划分成网格条块，实现对同一网格块进行递归式多次迭代步更新，从而在不改变串行模版迭代方法性质的同时，提高条块内数据局部性。然后通过对网格条块重排序，实现交错条块有限差分模版方法的并行化。本发明方法比传统的区域分解方法和红黑排序并行方法具有更好的数据局部性，并行效率和可扩展性。

Description

一种基于迭代空间条块的并行有限差分模版方法

技术领域

本发明属于工程计算领域，涉及一种基于迭代空间条块的并行有限差分模版方法。

背景技术

随着科学计算需要解决的问题渐趋复杂和并行计算机的迅速发展，并行计算已成为解决大规模并行科学计算问题必不可少的手段。在计算数学和计算物理等科学与工程计算领域，很多问题最终都归结为求解稀疏线性代数方程组。因此，在并行计算机上高效求解大规模稀疏线性代数方程组已成为当前科学计算领域的一项重要任务。

由于受舍入误差、计算机内存和计算复杂度的限制，对大规模问题，直接求解该类方程组几乎是不可能的，通常采用有限差分模版计算方法，迭代法的主要思想是通过构造有效的迭代格式，在有限步数内收敛于方程的精确解。

有限差分模版并行方法的实现一直是数值计算的重要研究对象。传统的做法是通过使用基于区域分解的多色排序方法来实现面向集群的并行GS方法，但是当数据量增大时，数据局部性成下降趋势，并且在每次迭代计算过程，都需要通信和同步。另一方面，循环分块技术可以提高循环并行度和数据局部性优化。为此，研究人员对迭代空间分块做了很多的研究。这些研究主要集中在两个方面：（1）数据局部性优化研究；（2）提高并行效率的分块尺寸/形状研究。通过数据依赖向量实现数据分块使其能够最大化复用，但是分块技术会引入数据块的空间依赖性，导致其基于空间条块的流水线执行方式不利于并行执行。但模版迭代方法并行化还有4个问题需要进一步的优化：

数据局部性问题。大多数方法能够有效的提高迭代内数据局部性，但对迭代间的数据局部性优化效果有限。

可扩展性问题。传统的并行化迭代方法在迭代内和迭代间都需要同步操作以维护数据依赖关系。处理机进行全局的同步，会增加开销时间，当处理机的台数增多时，全局同步的代价变得更加重要，并且影响方法的可扩展性。

通信和同步开销问题。由于传统的并行化方法需要在每次迭代过程中通过通信操作得到边界数据，通信开销制约了并行方法的效率。并且当问题给定时，随着处理机台数的增大，并行纯计算时间在减少，而通信时间在不断增加，这必将影响并行方法的可扩展性。

依赖关系问题。迭代空间划分所引入的条块依赖关系严重阻碍了循环的并行化。

发明内容

本发明的目的在于：鉴于现有方法存在的上述问题，提出基于迭代空间条块的并行有限模版方法。在不增加数据通信量的情况下，通过负超平面对迭代空间进行划分，减少通信启动开销和同步时间，提高数据局部性。并且通过对空间迭代块的重新排序，改变块间数据依赖性，有效地实现迭代方法的并行化。

本发明解决技术问题所采用的技术方案是：

本发明首先通过对迭代空间进行时间轴方向划分成网格条块，实现对同一网格块进行递归式多次迭代步更新，从而在不改变串行模版迭代方法性质的同时，提高条块内数据局部性。然后通过对网格条块重排序，实现交错条块有限差分模版方法的并行化。

步骤1根据处理器拓扑进行空间计算区域划分：

依据处理器数目P及处理器的拓扑结构，将空间计算区域划分为P1×P2个子空间sub_domain，其中P1×P2=P且P1=P2；

步骤2指定迭代方向：

交错条块迭代方法分为奇数k次迭代和偶数k次迭代，设定奇数k次迭代与偶数k次迭代执行方向相反；

步骤3对子空间进行网格条块划分：

通过时间轴将各子空间划分为多个网格条块，同串行迭代方法，并且奇数k次和偶数k次划分方向相反；

步骤4对网格条块重新排序：

根据处理器数目P按照cyclic(P)排列方法对所有的网格条块进行重排序；

步骤5以网格条块为单位执行奇数k次迭代数据更新：

执行奇数k次迭代更新，更新顺序参照重排序之后的网格条块序；当更新奇数k次中的发送条块和混合条块后，将边界数据发送给相应的接收条块和混合条块，而接收数据的条块在接收完数据后再进行迭代更新；

步骤6以网格条块为单位执行偶数k次迭代数据更新：

执行偶数k次迭代更新，更新顺序参照重排序之后的网格条块序；当更新偶数k次中的发送条块和混合条块后，将边界数据发送给相应的接收条块和混合条块，而接收数据的条块在接收完数据后再进行迭代更新。

本发明具有的有益效果是：

1、本发明方法通过网格条块重排序减少了cache缺失率、通信启动和同步次数。

2、本发明方法比传统的区域分解方法和红黑排序并行方法具有更好的数据局部性，并行效率和可扩展性。

附图说明

图1为迭代空间交错条块串行模版方法，(1)为从1到k次的正向执行，(2)为从k+1到2k次的反向执行。其中不同线条表示不同时间的迭代范围。

图2为本发明方法的系统框图。

图3为三维迭代空间的两类条块划分方法，(1)为对称遍历条块串行执行的前半部分对称遍历，其中迭代步是从第1到k次迭代，(2)为对称遍历条块串行执行的后半部分对称遍历，其中迭代步是从第:k+1到2k次迭代，(3)为对称遍历条块重排后的前半部分对称遍历，其中迭代步是从第1到k次迭代， (4)为对称遍历条块重排后后半部分对称遍历其中迭代步是从第:k+1到2k次迭代。

图4为条块超平面划分俯视图。

图5为图3所示两类条块划分方法的条块依赖图，(1)为图3(1)和图3(2)中划分方法的条块依赖图，(2)为图3(3)和图3(4)中划分方法的条块依赖图。其中数字前面的标号表示不同的迭代方向：f表示奇数次迭代，b表示为偶数次迭代。

图6为四个子区域的并行化，(1)为对称遍历过程的正向部分，(2)为对称遍历过程反向部分。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方案作进一步详细描述：

传统的迭代方法执行序过程如下：

在一次迭代内部依据网格点的顺序，依次对所有网格点进行迭代更新操作：

for ( t=1; t<=T; t++ )

for ( i=1; i<N-1; i++ )

for ( j=1; j<N-1; j++ )

A[t+1][i][j]=(A[t+1][i-1][j]+A[t][i+1][j]+A[t][i][j] + A[t+1][i][j-1] + A[t][i][j+1])/5

其中T表示迭代次数，N×N表示矩阵的大小。

在上述方法中，数据复用的方式共有以下四种：

行优先存储格式下，在j循环中复用A[i][j]，A[i][j+1]……

在j循环中，更新A[i][j-1]、A[i][j]和A[i][j+1]过程中复用A[i][j]。

在i循环中，更新A[i-1][j]、A[i][j]和A[i+1][j]过程中复用A[i][j]。

在时间步中，复用A[i][j]。

由于数据存取位置相近，因此出现第一类和第二类数据复用。当cache的容量C大于2N（N为一维边界长度）时，A[i+1][j]与A[i-1][j]在行优先存储格式下相邻2N个数据，因此会出现第三类数据复用。当N²小于C时，所有数据都会一次性读取到cache中，因此会出现第四类数据复用即提高了数据空间局部性。但这种“串行”执行序不仅导致其数据依赖关系阻碍了迭代方法并行化执行，而且当数据量大时，其数据局部性较差。其原因在于当数组大于cache容量时，本次更新的数据在下一次更新之前已经被写回内存。而且当数据量大时，传统迭代方法中多次迭代会导致数据cache的容量缺失。此外，多次迭代使数据的地址转换信息会周期性地在TLB中进行存取，也是影响其性能的一个重要因素。当矩阵规模增加时，大量的TLB容量缺失会严重影响程序的性能。因此，本发明采取迭代空间条块串行迭代执行序，其核心思想是：改变以往传统的以迭代次序为执行序的特点，将迭代空间分块引入执行序中，如图1所示。本发明还通过网格条块重排序，来实现交错条块有限差分模版方法的并行化。

其具体步骤描述如图2所示：

步骤1根据处理器拓扑进行空间计算区域划分：

依据处理器数目P及处理器的拓扑结构，将空间计算区域划分为P1×P2个子空间sub_domain，其中P1×P2=P且P1=P2，如图3(3)和图3(4)所示；

步骤2指定迭代方向：

交错条块迭代方法分为奇数k次迭代和偶数k次迭代，设定奇数k次迭代与偶数k次迭代执行方向相反；例如，使用LB-RT顺序(左下到右上)执行奇数k次迭代，使用RT-LB顺序(右上到左下)执行偶数k次迭代。

步骤3对子空间进行网格条块划分：

将各子空间划分为多个网格条块，同串行迭代方法，并且奇数k次和偶数k次划分方向相反；

以二维Possion方程为例，迭代空间交错条块串行迭代方法执行过程如下：

1) 网格条块划分：

m维网格点x(i1,…,im)组成m维数据空间data_space(m)。 n-1维数据空间data_space(n-1)和迭代维T的组合可以被看作n维的迭代空间iter_space(I₁,…,I_n-1,T)，其中的每个点都可由一个n×1维的列向量来表示，即

Figure 2011100512858100002DEST_PATH_IMAGE002

，其中

Figure 2011100512858100002DEST_PATH_IMAGE004

从左至右分别代表网格点在数据空间中的I₁,…,I_n-1维的维度坐标。数据空间中网格点x(i₁,…,i_m)在k次迭代的值，在迭代空间中可表示为u^k(i₁,…,i_n-1)。例如，三维迭代空间iter_space(I,J,T)由二维网格点x(i,j)和时间维T组成。u^k(i,j)表示网格点x(i,j)在k次迭代的值。

用区域分解方法将迭代空间iter_space(I,J,T)在T=0处进行数据划分。如图1所示，定义l ₁为子空间sub_iter_space(p,q,0)的行数，定义 l ₂为子空间sub_iter_space(p,q,0)的列数，划分后使得每个子空间sub_iter_space(p,q,0)中的网格点数为R= l ₁ x l ₂，l ₁与l ₂满足式(1)：

l ₁ >k ∩ l ₂ >k (1)

其中，k为单向迭代次数。

2) 在网格条块的基础上，沿时间轴对迭代空间进行划分：

划分方法采用时滞技术，对每层迭代的子空间修正边界，图1中虚边界表示修正后的边界。定义有向图G(V,E)存放相邻网格条块的关系.若网格条块v_i与网格条块v_j边界相连且v_i<v_j，则<v_i,v_j>

E。定义(v_i,v_j,k)为在第k次迭代中属于条块v_i但与条块v_j相邻的边界数据。算法描述如下：

/*条块边界修正算法*/

for(n=0;n<N;n+=2)

for( k= n*K+1; k<=(n+1)*K; k++) { //奇数K次迭代修正边界

Foreach <v_i,v_j> E

sub_iter_space(v_i,k+1)= sub_iter_space(v_i,k)- Pnode(v_i,v_j,k)

sub_iter_space(v_j,k+1)= sub_iter_space(v_j,k)+ Pnode(v_i,v_j,k)

}

for(k=(n+1)*K+1; k<=(n+1)*K; k++){ //偶数K次迭代修正边界

Foreach <v_i,v_j> E

sub_iter_space(v_j,k+1)= sub_iter_space(v_i,k)- Pnode(v_j,v_i,k)

sub_iter_space(v_i,k+1)= sub_iter_space(v_i,k)+ Pnode(v_j,v_i,k)

}

3) 按空间网格条块顺序执行模版迭代方法：

网格条块生成后，按条块顺序执行迭代方法。以条块内部网格点层作为内部执行序，条块内部迭代次数作为中间执行序，条块序作为外部执行序，更新每个网格点的值。执行过程如图1所示，a.1-a.4为条块正向执行顺序，b.1-b.4为条块反向执行顺序，其中数据块中3种不同的线条框分别代表连续的3次迭代计算。

步骤4对网格条块重新排序：

根据处理器个数P按照cyclic(P)排列对所有的网格条块进行重排序。

传统的多面体模型通过迭代空间多面体矩阵B，依赖矩阵D，超平面矩阵H，条块依赖关系矩阵S描述迭代空间的数据划分及依赖关系。但传统多面体模型中引入超平面执行顺序会降低迭代空间并行化的效率，因此需要改变超平面方向消除引入的数据依赖关系，提高迭代计算的并行性。在标准条块迭代方法中超平面由一组法向量组成。法向量定义了条块在同一超平面方向的两面，其中h_i(l)定义为I方向第l条块区域tile(l，I)在I方向超平面的后平面，h_i(l+1)定义为第l条块区域的I方向超平面的前平面。在条块中前平面不会引入I方向其他条块的依赖，而后平面中的数据在计算时需要I方向相邻条块的数据，因此引入了I方向数据依赖。例如：

，

，且

与

相邻，在计算

时，会使用

的边界数据。在边界条块中通过建立I方向的负超平面可以减少条块之间在I方向的依赖关系。将后平面h_i(l)替换为负超平面h_i`满足：

(2)

举例说明如图4所示，该条块区域为区域分解后在各处理器上的第一个条块(如图3(3)和图3 (4)中标注4的数据块)。当引入J向负超平面和I向负超平面，由超平面组成的迭代条块可以独立的执行迭代计算，而不需要引入其他依赖关系。从而在各处理器执行过程中，此类条块可以被同时执行更新迭代。

在维度j中可以独立执行的条块区域用

表示，其余非独立执行部分用表示。在每一维度j中，条块区域均可被分为两类：独立执行的条块区域和非独立执行的条块区域。因此在n维空间，迭代空间中共有2n类条块区域。由于我们仅考虑在迭代空间的并行执行过程，因此不考虑时间维。例如在图1中2维迭代空间(I,J)中条块区域可分为以下四类：，

，

，

。基于条块顺序的迭代执行的基本准则是：在每一维中先执行独立条块，然后将相应维度的边界依赖数据转递给非独立条块；当非独立条块得到所需要的依赖数据后，再执行非独立条块数据的迭代更新。例如在图3(3)中，属于

类别的条块为1，2，3，4；属于

类别的条块是5，6，7，8；属于

类别的条块是9，10，11，12；属于

类别的条块是13，14，15，16。而在图3(4)中，属于

类别的条块为16，15，14，13；属于

类别的条块是12，11，10，9；属于

类别的条块是8，7，6，5；属于

类别的条块是4，3，2，1。其执行顺序和通信顺序如图5(2)所示。

图5(1)显示网格条块按条块顺序串行执行迭代更新，从正向第一块数据块f_tile(1)（含标号f的块1）到反向最后一块数据块b_tile(1)（含标号b的块1）长度为32，因此没有并行度。图5(2)显示通过增加边界条块的负超平面进而改变网格条块的执行顺序，可以实现四个进程并行执行迭代更新，因此并行度为4。为了描述方便，本文在边界条块增加负超平面之后将图5(1)的条块顺序标号改为如图5(2)的顺序标号，进程1依次执行f_tile(1)、f_tile(5)、f_tile(9)、f_tile(13)、b_tile(13)、b_tile(9)、b_tile(5)、b_tile(1)条块，进程2依次执行f_tile(2)、f_tile(6)、f_tile(10)、f_tile(14)、b_tile(14)、b_tile(10)、b_tile(6)、b_tile(2)条块。进程3依次执行f_tile(3)、f_tile(7)、f_tile(11)、f_tile(15)、b_tile(15)、b_tile(11)、b_tile(7)、b_tile(3)条块，进程4依次执行f_tile(4)、f_tile(8)、f_tile(12)、f_tile(16)、b_tile(16)、b_tile(12)、b_tile(8)、b_tile(4)条块。其中，四个进程在正向执行条块迭代更新过程中，分别在执行所属第一个条块及第三个条块后执行同步操作以维护条块间的数据依赖关系。同样四个进行在反向执行条块更新过程中，也分别在执行所属第一个条块及第三个条块后执行同步操作以维护条块间的数据依赖关系。显然，第2种划分方式的执行速度是第1种划分方式的四倍。这两种方式最大的不同在于初始化时的排序方式不同。因为排序方式的不同代表了相邻条块间的执行顺序的差异，因此条块排序影响了条块间的数据依赖关系。通过改变网格条块间数据依赖关系，提高条块执行的并行度。如图5(2)所示。区域分解方法是实现分布式内存并行化的主要方法，但是传统的区域分解方法仅在空间维度上实现求解空间的分解，并没有考虑时间维度。为了实现并行类似于GS迭代本身有串行性质的迭代方法，交错条块并行方法改进了区域分解方法，使其沿时间维度对求解空间进行划分，形成不同的网格条块。为了方便说明，将迭代空间iter_spact(I₁,I₂,T)划分为四个子空间sub_domain1，sub_domain2，sub_domain3，sub_domain4。如图6所示。

空间条块的边界网格点必须发送给其它相邻子空间以维护数据依赖关系。子空间中的条块按通信类别分为四类：发送条块、接收条块、混合条块、非通信条块。发送条块在执行更新后将自身边界网格数据发送给接收条块。接收条块必须在接受其它子空间条块发送的边界网格数据后才执行更新。混合条块在执行前后需要接收数据和发送数据。其余的条块为非通信条块，其自身计算不需要其它处理器中网格条块的边界值。在交错条块方法中，按条块的执行序可以减少处理器之间的通信开销。当迭代执行2K次时，通信执行两次，通信数据量为2V，通信时间满足式(3)：

(3)

(4)

交错条块并行GS方法中需要正反方向各一次通信，每一次通信需要K个边界值，Ts是通信的启动时间。如式(3)、式(4)所示，T_comm是交错条块方法的通信时间，

是传统的区域分解并行GS方法中的通信时间。显然，在交错条块方法的通信开销中，启动时间比区域分解方法减少了2×(k-1)T_s。

步骤5以网格条块为单位执行奇数k次迭代数据更新：

执行奇数k次迭代更新，更新顺序参照重排序之后的网格条块序；过程如步骤4中所描述，当更新奇数k次中的发送条块和混合条块后，将边界数据发送给相应的接收条块和混合条块，而接收数据的条块在接收完数据后再进行迭代更新；

步骤6以网格条块为单位执行偶数k次迭代数据更新：

执行偶数k次迭代更新，更新顺序参照重排序之后的网格条块序；过程如步骤4中所描述，当更新偶数k次中的发送条块和混合条块后，将边界数据发送给相应的接收条块和混合条块，而接收数据的条块在接收完数据后再进行迭代更新。

Claims

1.一种基于迭代空间条块的并行有限差分模版方法，其特征在于该方法包括如下步骤：

步骤1根据处理器拓扑进行空间计算区域划分，具体是：依据处理器数目P及处理器的拓扑结构，将空间计算区域划分为P1×P2个子空间sub_domain，其中P1×P2=P且P1=P2；

步骤2指定迭代方向，具体是：交错条块迭代方法分为奇数k次迭代和偶数k次迭代，设定奇数k次迭代与偶数k次迭代执行方向相反；

步骤3对子空间进行网格条块划分，具体是：将各子空间划分为多个网格条块，并且奇数k次和偶数k次划分方向相反；

步骤4对网格条块重新排序，具体是：根据处理器数目P按照cyclic(P)排列方法对所有的网格条块进行重排序；

步骤5以网格条块为单位执行奇数k次迭代数据更新，具体是：执行奇数k次迭代更新，更新顺序参照重排序之后的网格条块序；当更新奇数k次中的发送条块和混合条块后，将边界数据发送给相应的接收条块和混合条块，而接收数据的条块在接收完数据后再进行迭代更新；

步骤6以网格条块为单位执行偶数k次迭代数据更新，具体是：执行偶数k次迭代更新，更新顺序参照重排序之后的网格条块序；当更新偶数k次中的发送条块和混合条块后，将边界数据发送给相应的接收条块和混合条块，而接收数据的条块在接收完数据后再进行迭代更新。