专利名称：一种基于迭代空间条块的并行有限差分模版方法
技术领域：
本发明属于工程计算领域，涉及一种基于迭代空间条块的并行有限差分模版方法。
背景技术：
随着科学计算需要解决的问题渐趋复杂和并行计算机的迅速发展，并行计算已成为解决大规模并行科学计算问题必不可少的手段。在计算数学和计算物理等科学与工程计算领域，很多问题最终都归结为求解稀疏线性代数方程组。因此，在并行计算机上高效求解大规模稀疏线性代数方程组已成为当前科学计算领域的一项重要任务。由于受舍入误差、计算机内存和计算复杂度的限制，对大规模问题，直接求解该类方程组几乎是不可能的，通常采用有限差分模版计算方法，迭代法的主要思想是通过构造有效的迭代格式，在有限步数内收敛于方程的精确解。有限差分模版并行方法的实现一直是数值计算的重要研究对象。传统的做法是通过使用基于区域分解的多色排序方法来实现面向集群的并行GS方法，但是当数据量增大时，数据局部性成下降趋势，并且在每次迭代计算过程，都需要通信和同步。另一方面，循环分块技术可以提高循环并行度和数据局部性优化。为此，研究人员对迭代空间分块做了很多的研究。这些研究主要集中在两个方面(1)数据局部性优化研究；(2)提高并行效率的分块尺寸/形状研究。通过数据依赖向量实现数据分块使其能够最大化复用，但是分块技术会引入数据块的空间依赖性，导致其基于空间条块的流水线执行方式不利于并行执行。 但模版迭代方法并行化还有4个问题需要进一步的优化
数据局部性问题。大多数方法能够有效的提高迭代内数据局部性，但对迭代间的数据局部性优化效果有限。可扩展性问题。传统的并行化迭代方法在迭代内和迭代间都需要同步操作以维护数据依赖关系。处理机进行全局的同步，会增加开销时间，当处理机的台数增多时，全局同步的代价变得更加重要，并且影响方法的可扩展性。通信和同步开销问题。由于传统的并行化方法需要在每次迭代过程中通过通信操作得到边界数据，通信开销制约了并行方法的效率。并且当问题给定时，随着处理机台数的增大，并行纯计算时间在减少，而通信时间在不断增加，这必将影响并行方法的可扩展性。依赖关系问题。迭代空间划分所引入的条块依赖关系严重阻碍了循环的并行化。
发明内容
本发明的目的在于鉴于现有方法存在的上述问题，提出基于迭代空间条块的并行有限模版方法。在不增加数据通信量的情况下，通过负超平面对迭代空间进行划分，减少通信启动开销和同步时间，提高数据局部性。并且通过对空间迭代块的重新排序，改变块间数据依赖性，有效地实现迭代方法的并行化。本发明解决技术问题所采用的技术方案是
本发明首先通过对迭代空间进行时间轴方向划分成网格条块，实现对同一网格块进行递归式多次迭代步更新，从而在不改变串行模版迭代方法性质的同时，提高条块内数据局部性。然后通过对网格条块重排序，实现交错条块有限差分模版方法的并行化。步骤1根据处理器拓扑进行空间计算区域划分
依据处理器数目P及处理器的拓扑结构，将空间计算区域划分为Pl XP2个子空间 domain，其中 Pl XP2=P 且 P1=P2 ； 步骤2指定迭代方向
交错条块迭代方法分为奇数k次迭代和偶数k次迭代，设定奇数k次迭代与偶数k次迭代执行方向相反；
步骤3对子空间进行网格条块划分
通过时间轴将各子空间划分为多个网格条块，同串行迭代方法，并且奇数k次和偶数k 次划分方向相反；
步骤4对网格条块重新排序
根据处理器数目P按照CydiC(P)排列方法对所有的网格条块进行重排序； 步骤5以网格条块为单位执行奇数k次迭代数据更新
执行奇数k次迭代更新，更新顺序参照重排序之后的网格条块序；当更新奇数k次中的发送条块和混合条块后，将边界数据发送给相应的接收条块和混合条块，而接收数据的条块在接收完数据后再进行迭代更新；
步骤6以网格条块为单位执行偶数k次迭代数据更新
执行偶数k次迭代更新，更新顺序参照重排序之后的网格条块序；当更新偶数k次中的发送条块和混合条块后，将边界数据发送给相应的接收条块和混合条块，而接收数据的条块在接收完数据后再进行迭代更新。本发明具有的有益效果是
1、本发明方法通过网格条块重排序减少了 cache缺失率、通信启动和同步次数。2、本发明方法比传统的区域分解方法和红黑排序并行方法具有更好的数据局部性，并行效率和可扩展性。
图1为迭代空间交错条块串行模版方法，⑴为从1到k次的正向执行，(2)为从 k+Ι到业次的反向执行。其中不同线条表示不同时间的迭代范围。图2为本发明方法的系统框图。图3为三维迭代空间的两类条块划分方法，(1)为对称遍历条块串行执行的前半部分对称遍历，其中迭代步是从第1到k次迭代，(2)为对称遍历条块串行执行的后半部分对称遍历，其中迭代步是从第k+l到业次迭代，(3)为对称遍历条块重排后的前半部分对称遍历，其中迭代步是从第1到k次迭代，(4)为对称遍历条块重排后后半部分对称遍历其中迭代步是从第k+l到业次迭代。图4为条块超平面划分俯视图。图5为图3所示两类条块划分方法的条块依赖图，(1)为图3(1)和图3(2)中划分方法的条块依赖图，( 为图3C3)和图3(4)中划分方法的条块依赖图。其中数字前面的标号表示不同的迭代方向f表示奇数次迭代，b表示为偶数次迭代。
图6为四个子区域的并行化，(1)为对称遍历过程的正向部分，(2)为对称遍历过程反向部分。
具体实施例方式下面结合附图和实施例，对本发明的具体实施方案作进一步详细描述 传统的迭代方法执行序过程如下
在一次迭代内部依据网格点的顺序，依次对所有网格点进行迭代更新操作 for ( t=l; t<=T; t++ ) for ( i=l; i<N-l; i++ ) for ( j=l; j<N-l; j++ )
A[t+1] [i] [j] = (A[t+l] [i-1] [j]+A[t] [i+1] [j]+A[t] [i] [j] + A[t+l][i][j-l] + A[t][i][j+l])/5
其中T表示迭代次数，NXN表示矩阵的大小。在上述方法中，数据复用的方式共有以下四种
行优先存储格式下，在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_l] [j]、A[i] [j]和 A[i+1] [j]过程中复用 A[i] [j]。在时间步中，复用A[i][j]。由于数据存取位置相近，因此出现第一类和第二类数据复用。当cache的容量C 大于2N (N为一维边界长度)时，A[i+l][j]与A[i-l][j]在行优先存储格式下相邻2N个数据，因此会出现第三类数据复用。当N2小于C时，所有数据都会一次性读取到cache中， 因此会出现第四类数据复用即提高了数据空间局部性。但这种“串行”执行序不仅导致其数据依赖关系阻碍了迭代方法并行化执行，而且当数据量大时，其数据局部性较差。其原因在于当数组大于cache容量时，本次更新的数据在下一次更新之前已经被写回内存。而且当数据量大时，传统迭代方法中多次迭代会导致数据cache的容量缺失。此外，多次迭代使数据的地址转换信息会周期性地在TLB中进行存取，也是影响其性能的一个重要因素。当矩阵规模增加时，大量的TLB容量缺失会严重影响程序的性能。因此，本发明采取迭代空间条块串行迭代执行序，其核心思想是改变以往传统的以迭代次序为执行序的特点，将迭代空间分块引入执行序中，如图1所示。本发明还通过网格条块重排序，来实现交错条块有限差分模版方法的并行化。其具体步骤描述如图2所示
步骤1根据处理器拓扑进行空间计算区域划分
依据处理器数目P及处理器的拓扑结构，将空间计算区域划分为Pl XP2个子空间 i/offlai/ ，其中 Ρ1ΧΡ2=Ρ 且 Ρ1=Ρ2，如图 3(3)禾Π图 3 )所示； 步骤2指定迭代方向
交错条块迭代方法分为奇数k次迭代和偶数k次迭代，设定奇数k次迭代与偶数k次迭代执行方向相反；例如，使用LB-RT顺序(左下到右上)执行奇数k次迭代，使用RT-LB 顺序(右上到左下)执行偶数k次迭代。步骤3对子空间进行网格条块划分将各子空间划分为多个网格条块，同串行迭代方法，并且奇数k次和偶数k次划分方向相反；
以二维Possion方程为例，迭代空间交错条块串行迭代方法执行过程如下 1)网格条块划分
m维网格点x(il，…，im)组成m维数据空间data_space (m)。 n_l维数据空间data_Space (n-1)和迭代维T的组合可以被看作η维的迭代空间iter_ space (I1，…，Ilri，Τ)，其中的每个点都可由一个ηΧ 1维的列向量来表示，即
F = Ci1,...,^f，其中.,U从左至右分别代表网格点在数据空间中的I1,…，Ilri维的维度
坐标。数据空间中网格点Wi1，…，im)在k次迭代的值，在迭代空间中可表示为!^(‘… ，‘)。例如，三维迭代空间iter_SpaCe(I，J，T)由二维网格点x(i，j)和时间维T组成。 uk(i, j)表示网格点x(i，j)在k次迭代的值。用区域分解方法将迭代空间iter_SpaCe(I，J, Τ)在T=O处进行数据划分。如图 1所示，定义Λ为子空间Sub_iter_Space(p，q，0)的行数，定义及为子空间sub_iter_ space (ρ, q, 0)的列数，划分后使得每个子空间sub_iter_space (p，q, 0)中的网格点数为R= I1 χ Λ，A与及满足式(1)
I1 >k Π I2 >k (1) 其中，k为单向迭代次数。2)在网格条块的基础上，沿时间轴对迭代空间进行划分
划分方法采用时滞技术，对每层迭代的子空间修正边界，图1中虚边界表示修正后的边界。定义有向图G(V，E)存放相邻网格条块的关系.若网格条块Vi与网格条块Vj边界相
连且Vi<\，则<Vi，、> e Ε。定义Pmfe (Vi, Vj, k)为在第k次迭代中属于条块Vi但与条块Vj
相邻的边界数据。算法描述如下 /*条块边界修正算法*/ for (n=0;n<N;n+=2)
for( k= n*K+l ； k<=(n+l)*K; k++) { //奇数 K 次迭代修正边界 Foreach Ky1, v」> E
sub_iter_space (Vi, k+1) = sub_iter_space (Vi, k) _ Pnode (Vi, Vj, k)
sub_iter_space (v」，k+1) = sub_iter_space (v」，k) + Pnode (Vi, Vj, k)
}
for (k=(n+l)^K+l; k〈=(n+l)*K; k++) { //偶数 K 次迭代修正边界 Foreach〈ν” ν』> E
sub_iter_space (Vjj k+1) = sub_iter_space (Vij k) - Pnode (Vjj Vij k) sub_iter_space (Vij k+1) = sub_iter_space (Vij k) + Pnode (Vjj Vij k)
ι
3)按空间网格条块顺序执行模版迭代方法
网格条块生成后，按条块顺序执行迭代方法。以条块内部网格点层作为内部执行序，条块内部迭代次数作为中间执行序，条块序作为外部执行序，更新每个网格点的值。执行过程如图1所示，a. 1-a. 4为条块正向执行顺序，b. l_b. 4为条块反向执行顺序，其中数据块中3种不同的线条框分别代表连续的3次迭代计算。步骤4对网格条块重新排序
根据处理器个数P按照CyCliC(P)排列对所有的网格条块进行重排序。传统的多面体模型通过迭代空间多面体矩阵B，依赖矩阵D，超平面矩阵H，条块依赖关系矩阵S描述迭代空间的数据划分及依赖关系。但传统多面体模型中引入超平面执行顺序会降低迭代空间并行化的效率，因此需要改变超平面方向消除引入的数据依赖关系， 提高迭代计算的并行性。在标准条块迭代方法中超平面由一组法向量组成。法向量定义了条块在同一超平面方向的两面，其中Iii (乃定义为I方向第/条块区域tile(A I)在I方向超平面的后平面，^(7+1)定义为第/条块区域的I方向超平面的前平面。在条块中前平面不会引入I方向其他条块的依赖，而后平面中的数据在计算时需要I方向相邻条块的数据，因此引入了 I方向数据依赖。例如V <Ρ'·0 , V fik(q,rj€ suhdomain(m)
，且ftKiM)与iAfeJ)相邻，在计算时，会使用fife(ftZ)的边界数据。在边界条块中通过建立I方向的负超平面可以减少条块之间在I方向的依赖关系。将后平面、(乃替换为负超平面、、满足
权利要求
1.一种基于迭代空间条块的并行有限差分模版方法，其特征在于该方法包括如下步骤步骤1根据处理器拓扑进行空间计算区域划分，具体是依据处理器数目P及处理器的拓扑结构，将空间计算区域划分为P1XP2个子空间SizliZofflai/ ，其中P1XP2=P且P1=P2 ； 步骤2指定迭代方向，具体是交错条块迭代方法分为奇数k次迭代和偶数k次迭代， 设定奇数k次迭代与偶数k次迭代执行方向相反；步骤3对子空间进行网格条块划分，具体是将各子空间划分为多个网格条块，并且奇数k次和偶数k次划分方向相反；步骤4对网格条块重新排序，具体是根据处理器数目P按照cyclic (P)排列方法对所有的网格条块进行重排序；步骤5以网格条块为单位执行奇数k次迭代数据更新，具体是执行奇数k次迭代更新，更新顺序参照重排序之后的网格条块序；当更新奇数k次中的发送条块和混合条块后， 将边界数据发送给相应的接收条块和混合条块，而接收数据的条块在接收完数据后再进行迭代更新；步骤6以网格条块为单位执行偶数k次迭代数据更新，具体是执行偶数k次迭代更新，更新顺序参照重排序之后的网格条块序；当更新偶数k次中的发送条块和混合条块后， 将边界数据发送给相应的接收条块和混合条块，而接收数据的条块在接收完数据后再进行迭代更新。
全文摘要
本发明涉及一种基于迭代空间条块的并行有限差分模版方法。传统的并行化迭代方法在迭代内和迭代间都需要同步操作以维护数据依赖关系。本发明首先通过对迭代空间进行时间轴方向划分成网格条块，实现对同一网格块进行递归式多次迭代步更新，从而在不改变串行模版迭代方法性质的同时，提高条块内数据局部性。然后通过对网格条块重排序，实现交错条块有限差分模版方法的并行化。本发明方法比传统的区域分解方法和红黑排序并行方法具有更好的数据局部性，并行效率和可扩展性。
文档编号G06F17/00GK102200962SQ20111005128
公开日2011年9月28日 申请日期2011年7月25日 优先权日2011年7月25日
发明者万健, 任永坚, 张伟, 张纪林, 徐向华, 蒋从锋 申请人:杭州电子科技大学