CN111079078A - 面向结构网格稀疏矩阵的下三角方程并行求解方法 - Google Patents

Abstract

在高性能并行计算中求解结构化网格稀疏下三角方程中优化处理器间通信和高速缓冲存储器cache空间利用的方法，包括：接收结构化网格的求解向量等，启动多处理器开始并行执行；将求解向量分解为多个子块，将每个子块划分为多个柱，将多处理器视为矩形阵列，每个处理器负责计算一个柱，其中处理器在自身的cache中开辟存储依赖数据的缓冲区，从邻居或者自身内存得到计算所需数据并存入缓冲区，完成自己柱计算，将计算完成的数据发送给依赖该数据的其它处理器；处理器矩形阵列遍历所有子块，完成计算所需数据通信，对每个子块进行计算，直到完成整个求解计算。通过处理器间的快速通信传输每次运算的依赖数据，最小化访存开销，最大化资源利用率。

Description

面向结构网格稀疏矩阵的下三角方程并行求解方法

技术领域

本发明涉及数值计算及高性能计算领域，特别涉及一种高性能并行计算机环境下的求解结构化网格形式的下三角稀疏矩阵的方法。

背景技术

随着计算机技术的发展，大规模科学计算的需求逐渐增大。科学计算的很多领域，如数值天气预报模拟，流体力学计算等，都需要对大规模稀疏线性方程组进行求解。结构化网格是指网格区域内所有的内部点都具有相同的毗邻单元。作为网格的一种特殊形式，其具有数据结构简单，生成速度快，适用于流体计算等优势，因而广泛存在于许多真实应用中。

因此，结构化稀疏下三角矩阵求解(Structured Sparse Triangular Solver)在科学计算领域扮演着十分重要的角色。由于结构化网格的特殊性，专门针对结构化网格稀疏下三角矩阵求解的高效并行算法也就存在可能。

高性能技术的发展使得多处理器协同运行，高效利用并行单元进行计算和访存成为可能。在多处理器提供强大计算能力的基础上，处理器自身的高速缓冲存储器(cache)提供的快速访存模式和系统内存(main memory)提供的大容量，但高延迟的访存模式形成了鲜明对比。如何有效cache利用空间，增大cache命中率，减少内存访问次数，增大内存访问粒度，对程序开发设计者提出了挑战。为了帮助编程者达到这一目标，高性能计算机系统的多个处理器之间往往存在细粒度，低延迟的数据通信模式。以国产申威SW26010处理器为例，单核组上的64个众核存在寄存器通信延迟往往在11个时钟周期以内，最高带宽可达到637GB/s，远低于直接访问主存所带来的延迟。因此，巧妙地利用处理器之间的通信可以快速有效地实现处理器间的同步操作和数据共享，从而为发掘并发性和优化数据访存带来了助力。

目前，一种快速的稀疏下三角矩阵求解方法已经被提出。它通过设计稀疏数据布局和生产者消费者配对的通信方法解决了访存和数据依赖的问题。但通用SpTRSV无法适应结构化网格问题，因为需要极其昂贵的预分析才能发掘并发性并获得缓存友好的数据布局，而对于数据结构较为简单的结构网格SpTRSV这是多余的。结构化网格下三角矩阵求解更加依赖合理的任务划分和数据依赖处理。对于结构网格SpTRSV，目前尚缺少一种以有效利用处理器间通信方法求解结构化稀疏下三角矩阵的算法。

发明内容

本发明的旨在解决上述的技术缺失。

本发明的实施例提出一种在高性能并行计算平台上求解多种形式的结构网格稀疏下三角矩阵的算法，包括以下步骤：

S1：程序接收输入的网格规模，接收输入网格，输入向量和求解向量，并启动多处理器开始并行执行；其中，输入网格是一个规模为M×N×L的三维结构化网格矩阵A，输入向量和求解向量是规模为M×N×L的向量。

S2：将整个求解向量分解为若干P₁×P₂×L的子块。将处理器在逻辑上看作矩形网格阵列。每个处理器分别计算大小为L的连续空间的一柱，处理器在自身的cache中开辟存储依赖数据的缓冲区，从邻居或者内存得到计算所需要的数据并存入缓冲区，完成自己所在柱的计算，并将计算完成的数据发送给计算过程依赖自己所计算数据的处理器。

S3：处理器阵列依次遍历所有的子块，并完成必要的数据通信，对求解向量中的每个子块进行计算，直到完成整个结构化网格的求解计算。

上述方法中，步骤S1中，启动多处理器开始并行执行包括：启动处理器组中所有可用处理器资源，每个线程绑定一个处理器资源。这里通常启动的运算资源为全部P₁×P₂个处理器，对应P₁×P₂个线程。

上述方法中，步骤S1中，输入网格，输入向量和求解向量进一步包括：输入网格是一个规模为M×N×L的三维结构化网格矩阵A，用一个M×N×L×rnz的一维数组表示。其中，M×N×L表示网格空间的大小，代表求解的矩阵共有M×N×L行和列，每一行有rnz个非零元，其中有一个非零元分布在矩阵的对角线上。由于输入网格是结构化网格，每一行的非零元分布有固定的位置关系。输入向量是一个规模为M×N×L的向量B，求解向量是一个规模为M×N×L的向量X，。令A(i，j，k，r)(0≤r＜rnz)依赖的数据位置为(DI_r，DJ_r，DK_r)(假定0号非零元代表矩阵对角线上的元素)则求解计算关系如下：

上述方法中，每一行的非零元分布有一定的位置关系包括：对于网格中任意两个不同的位置(i，j，k)和(i′，j′，k′)，有

DI_r-i＝DI′_r-i′，

DJ_r-j＝DJ′_r-j′，DK_r-k＝DK′_r-k′.

同时，非零元的相对位置分布满足一定的条件限制。

上述方法中，非零元的相对位置分布满足一定的条件限制进一步包括：求解向量用一维数组表示，三个维度的大小分别为M，N和L，其在网格中的下标位置x与三个维度分别对应的坐标I，J，K表示如下：

同时输入满足条件A：

条件A：

i≥j&I(i)-I(j)≤δ_I&J(i)-J(j)≤δ_J&K(i)-K(j)≤δ_K；

其中l_ij≠0且i≥j表明对网格下标i位置数值的求解依赖于网格下标j位置的求解结果。为针对不同情况进行求解，对条件A进行不同程度的简化和分离。

上述方法中，对条件A进行不同程度的简化和分离包括：将条件A分离和简化分别得到条件B，条件C和条件D，描述如下：

条件B：

i≥j&I(i)-I(j)+J(i)-J(j)+K(i)-K(j)≤1；

条件C：

i≥j&I(i)-I(j)+J(i)-J(j)≤1&K(i)-K(j)≤δ_K；

条件D：

i≥j&I(i)-I(j)≤δ_I&J(i)-J(j)≤δ_J&K(i)＝＝K(j)；

条件B，条件C和条件D可以通过组合得到条件A。

上述方法中，条件B，条件C和条件D可以通过组合得到条件A包括：条件B是三维结构网格中最简单的问题。条件C是在条件B的基础上，保持i，j方向的依赖关系不变，放松k方向的限制。条件D则是只考虑了i，j方向的依赖处理。最后只需将处理条件C和条件D的算法合并就可以得到处理条件A的算法。

上述方法中，步骤S2中，将处理器看作矩形阵列包括：假设多个处理器以P₁行P₂列的矩形阵列方式排列，只有每一行和每一列之间的处理器可以相互通信。在网格上位于第i行j列的处理器用CPE(i，j)表示，并把在阵列位置上与其相邻的处理器称为邻居。例如CPE(i，j)和CPE(i，j+1)，CPE(i，j)和CPE(i+1，j)之间为邻居关系。该发明的实施过程中，对处理器阵列的构造需要考虑不同处理器之间进行通信的代价，根据高性能计算系统的处理器个数和通信关系对矩阵阵列进行构造。

上述方法中，根据高性能计算系统的处理器个数和通信关系对矩阵阵列进行构造进一步包括：对于部分处理器架构，如国产处理器SW26010，其处理器本身按照8×8的矩阵阵列形式排列并提供行/列的快速通信，因此可以直接以处理器间拓扑关系作为矩阵阵列。对于其他提供处理器间快速通信的处理器架构，根据处理器的数目构造，使得使用的处理器规模P₁×P₂尽可能接近总处理器数目，从而最大化资源利用率。如果P₁×P₂不等于总处理器数目，多出的处理器会闲置。同时使矩阵阵列尽可能接近正方形。构造原则应是：尽量将通信代价小的处理器放在同一行/列。

上述方法中，使矩阵阵列尽可能接近正方形进一步包括：在不违背上述构造矩阵阵列方法原则的情况下，使|P₁-P₂|尽可能小，使得阵列形状更接近于正方形。这是由于，在一个P₁×P₂的矩阵阵列中，对于大部分下三角稀疏矩阵求解，处理器计算过程中相互依赖的层级关系近似于用从左下到右上，与水平方向成45度角的斜线从下往上进行多次切分，每一条线代表一个层级。在面积相同的所有矩形中，正方形的这种相互依赖的层级数目最少。

上述方法中，步骤S2中，处理器分别计算大小为L的连续空间的一柱进一步包括：对于规模为M×N×L的求解向量，每个处理器一次计算X(i，j，0∶L)，其中0≤i＜M，0≤j＜N，表示要计算的柱在水平面上的位置。

上述方法中，步骤S2中，处理器从邻居或者内存得到计算所需要的数据并存入缓冲区包括：当处理器位于处理器阵列边缘时，其依赖某一柱数据可能不在目前处理器阵列遍历到的子块范围内。此时，这样计算依赖的数据要么为0，对应整个网格矩阵的边界部分；要么已经在之前遍历到的子块中计算完成并存放在内存中。此时，需要将对应位置的缓冲区置为全0，或通过DMA将读入到缓冲区中。

当处理器不位于处理器阵列前部边缘时，或位于处理器阵列边缘但其依赖的某一柱数据在子块范围内。此时，该柱数据由处理器阵列的另一处理器负责计算，则需要接收来自该处理器的通信数据并存入缓冲区中。

上述方法中，步骤S2中，将计算完成的数据发送给计算过程依赖自己所计算数据的处理器包括：处理器完成一柱中某一部分数据X(i，j，z₁∶z₂)的计算，若根据网格结构，若该处理器某一邻居计算X(i′，j′，z′)依赖该数据，则算法规定该处理器会在计算完成后将数据X(i，j，z)发送给该邻居；若某一非邻居处理器计算X(i′，j′，z′)依赖该数据，则算法规定该处理器会经过一次或几次中间结点的将数据传送给该非邻居处理器。

上述方法中，处理器会经过一次或几次中间结点的将数据传送给该非邻居处理器包括：此时网格输入满足条件D，每一个处理器的计算将依赖其相邻的其他多个处理器，并且有的处理器与它并不在同一行或者同一列。这里规定处理器阵列不能做对角方向的通信，这时需要通过另外一个或几个处理器做跳转，因为非直线的通信会显著增大通信和同步的开销。此时，中间节点相应的计算一定也依赖于该处理器X(i，j，z₁∶z₂)的数据。若不谨慎地处理以上通信行为，会使通信过程变得复杂，并且有可能会使两个或者更多的处理器之间形成通信环，影响性能并造成死锁。因此需要进行通信路径规约和对角划分。

上述方法中，通信路径规约包括：条件D中，此时δ_I和δ_J为不大于P的任意自然数(大于P的情况超出了处理器阵列的边长，不作考虑)。如果处理器B的计算依赖处理器A，那么对于处理器A，有一条出边，而对于处理器B，则有一条入边。对于满足条件D的问题，每个处理器都有2δ_Iδ_J+δ_I+δ_J条入边和出边。可以根据不同边的重要程度将边分为关键边和非关键边。实际通信时，只要做关键边上的通信即可完成整个通信过程。

上述方法中，将边分为关键边和非关键边包括：如果处理器A到处理器B有一条边，并且可以发现其他一条以处理器A为起点，且以处理器B为终点的路径存在，那么该边为非关键边，反之为关键边。可以发现，若CPE(i，j)表示行编号为i，列编号为j的处理器，对其所有的2δ_Iδ_J+δ_I+δ_J条出边，只有两条关键边：CPE(i，j)→CPE(i，j+1)和CPE(i，j)→CPE(i+1，j-δ_J)。其他的通信过程，都可由这样的关键边得到。

上述方法中，其他的通信过程，都可由这样的关键边得到包括：

对于边：CPE(i，j)→CPE(i+a，j+b)，1＜a≤δ_I，-δ_J≤b≤δ_J，

可以发现一条路径：CPE(i，j)→CPE(i+1，j)→CPE(i+a，j+b)；

对于边：CPE(i，j)→CPE(i+1，j+b)，0＜b≤δ_J，

可以发现一条路径：CPE(i，j)→CPE(i，j+b)→CPE(i+1，j+b)；

对于边：CPE(i，j)→CPE(i+1，j+b)，-δ_J≤b＜0，

可以发现一条路径：CPE(i，j)→CPE(i，j+1)→CPE(i+1，j+b)；

对于边：CPE(i，j)→CPE(i，j+b)，1＜b≤δ_J，

可以发现一条路径：CPE(i，j)→CPE(i，j+1)→CPE(i，j+b).

因而以上2δ_Iδ_J+δ_I+δ_J条出边都是非关键边，入边的推导同理。

上述方法中，对角划分包括：经过上面的通信路径规约后，每个处理器只剩下两条出边需要进行通信。对于边CPE(i，j)→CPE(i，j+1)，可以直接通过行或者列通信完成。而边CPE(i，j)→CPE(i+1，j-δ_J)仍然需要对角通信。为了避免寻找跳转节点从而形成通信环进而导致死锁的情况发生，需要采用对角划分的方法。例如：在一个8×8的子块中，柱(0，0)由CPE(0，0)处理，柱(1，0)由CPE(1，δ_J)处理，柱(0，2)由CPE(2，δ_J+1)处理，依此类推，列方向上负责的处理器顺序排列。在这样的划分中，会有一些柱超出了长方体的问题范围，这时仍然要分配处理器去处理它们，只负责通信而没有计算任务。

上述方法中，步骤S3中，完成必要的数据通信包括：当某一处理器位于处理器边缘时，其依赖的数据很可能是前一个子块中某个处理器的计算结果。如果简单地让该计算结果存入内存再进行DMA读入，效率较低。因此，位于处理器阵列最后一行的处理器，如考虑δ_I＝δ_J＝1的情况，CPE(P₁-1，0)-CPE(P₁-1，P₂-1)，会在计算完成后将自己和邻居的部分计算数据发给位于处理器阵列第一行的处理器，即CPE(0，0)-CPE(0，P₂-1)，以便其开展下一个子块的计算。

根据本发明另一方面，提供了一种计算机可读存储介质，其上存储有计算机指令，所述指令当被高性能并行计算机执行时，执行上述方法。

根据本发明实施例的方法，通过处理器间的快速通信传输每次运算的依赖数据，最小化访存开销，最大化计算资源利用率，从而对各种形式的结构网格稀疏下三角矩阵方程进行高效求解。

附图说明

图1为根据本发明实施例的利用高性能并行计算机对结构化网格稀疏下三角矩阵求解中优化处理器通信和cashe利用的方法的总体流程图。

图2为根据本发明一个实施例的对满足条件C的结构化网格进行求解的算法流程图。

图3为根据本发明一个实施例的对满足条件D的结构化网格进行求解的算法流程图。

图4为根据本发明的一个实施例在求解满足条件B和条件C的结构化网格过程示意图。

图5(a)-(d)为根据本发明的一个实施例在处理满足条件D的结构化网格过程中提供的通信路径规约和对角划分方法示意图。

图6为根据本发明的一个实施例在神威太湖之光超级计算机上运行的性能测试结果示意图。

具体实施方式

以下将通过参考附图对示例性实施例进行详细描述，附图意在描绘示例性实施例而不应被解释为对权利要求的预期范围加以限制。

图1为根据本发明实施例的利用高性能并行计算机对结构化网格稀疏下三角矩阵求解中优化处理器通信和cashe利用的方法的总体流程图。

S1：接收结构化网格的网格规模，接收输入网格，输入向量和求解向量，并启动多处理器开始并行执行；

S2：将整个求解向量分解为多个子块，将每个子块划分为多个柱，每个柱代表求解向量在内存地址空间中连续的一段，将多处理器逻辑上视为矩形阵列，每个处理器负责计算一个柱，其中每个处理器在自身的cache中开辟用于存储依赖数据的缓冲区，从邻居处理器或者自身内存得到计算所需要的数据并存入所述缓冲区，完成自己所在柱的计算，并将计算完成的数据发送给计算过程依赖所述计算完成的数据的其它处理器。

S3：所述处理器矩形阵列依次遍历所有的子块，并完成计算所需的数据通信，对求解向量中的每个子块进行计算，直到完成整个结构化网格的求解计算。图2为根据本发明对满足条件C的结构化网格进行求解的算法流程图。下文参考图2，以条件C中δ_K＝1为例，对本发明实施例提供的求解结构化下三角稀疏矩阵方法进行详细阐释，包括以下步骤，

步骤S1：程序接受网格矩阵A，输入向量B，求解向量X。将规模为M×N×L的三维结构化网格矩阵A划分为多个P₁×P₂×L大小的矩阵，从而可以调度全部P₁×P₂个处理器进行计算。计算每个子块时，单个处理器负责计算一柱规模为L的数据。应当注意，当输入网格矩阵的规模不满足P₁或P₂的整数倍时，处理器阵列将在最后一行/列的子块中启用部分计算资源进行计算，资源数刚好等于最后不完整子块包含的柱数。

步骤S2：处理器计算行编号rid和列编号cid，并开启一个大小为3×L的缓冲区，令其为buf[3×L]。规定下文中a[index：size]表示a数组中从下标index开始，长度为size的部分数据，处理器的行列标号为i，j。buf缓冲区包含三个长度为L的数据，buf[0∶L]表示行方向i上的依赖数据，buf[L∶L]表示列方向j上的依赖数据，他们分别用来存放邻居X[i-1，j，0∶L]和X[i，j-1，0∶L]已经计算完成的数据。buf[2×L∶L]表示自身负责计算的柱的数据，其初始值和计算结果都存在其中。

同时在开辟缓冲区时，要将右端向量B[i，j，0∶L]读入并存储在buf[2×L∶L]中，并将其余部分初始化为0。

步骤S3_1：如果计算位置k＝0时需要进行判断，如果当前处理器位于处理器阵列的边界，但是其处理的柱并不在边界，这种情况在第一个子块后面的其他子块的处理中都会出现。这时，目标柱所依赖的数据在前一个子块中已经被计算完成并且写入到主存中，需要重新读入cache中。而DMA操作只需进行一次，在k＝0时将L长度的内存的依赖数据读入即可。

在开始描述通信部分前，规定putc(a，i)为每个处理器往第i行列通信发送数据a，putr(a，j)为往第j列行通信发送数据a。getc(i)为从第i行列通信接收得到的数据，getr(j)为从第j列行通信接收得到的数据。

步骤S3_2：当处理器不位于阵列边缘时，在计算每一个未知数时，需要从相邻的处理器获得所需的数据。例如(考虑δ_K＝1的情况)，当CPE(i，j)计算X(i，j，k)时，可能会依赖CPE(i，j-1)计算的X[i，j-1，k-1]，X[i，j-1，k]和X[i，j-1，k+1]以及CPE(i-1，j)计算的X[i-1，j，k-1]，X[i-1，j，k]和X[i-1，j，k+1]。而对于一个CPE的计算而言，相邻两个未知数依赖的大部分数据是可以复用的。例如：CPE(i，j)在计算X[i，j，k]会依赖X[i-1，j，k-1]，X[i-1，j，k]和X[i-1，j，k+1]，计算X[i，j，k+1]会依赖X[i-1，j，k]，X[i-1，j，k+1]和X[i-1，j，k+2]，这个过程中X[i-1，j，k]

和X[i-1，j，k+1]都是可以复用的，复用数量为2×δ_K。因此，每个非阵列边缘的处理器在计算一个柱的第一个未知数(k＝0)时，需要等待相邻的每个处理器发送(1+δ_K)个数据，而计算第二个及之后的未知数时，只需要等待相邻处理器发送一个数据即可。

在每一个柱开始计算之前，每个处理器先等待每个相邻的处理器发送δ_K个数据。然后在k逐渐增大的过程中，每一个处理器从它每一个邻居处理器各多获得一个数据即可开始计算当前的X[i，j，k]。计算最后δ_K个未知数时，由于所有的依赖数据都已经获得，这时不再需要进行通信。

步骤S3通信的具体形式为：

buf[k]←getc(rid-1)，

buf[L+k]←getr(cid-1).

步骤S4：使用已经获取的依赖数据对X[i，j，k]进行计算，这里不同网格的计算方式和依赖关系不同，对计算过程不作描述。

步骤S5：X[i，j，k]的计算完成后，将计算得到的数据通过行/列通信发给邻居处理器。

步骤S5通信的具体形式为：

putc(X[2×L+k]，rid+1)，

putr(X[2×L+k]，cid+1).

步骤S6：处理器完成整个柱的计算后，将所计算的柱传入内存的X数组中，具体传递关系为为X[i，j，0∶L]←buf[2×L∶L]。

图3为根据本发明对满足条件D的结构化网格进行求解的算法流程图。下文参考图3，对本发明实施例提供的求解结构化下三角稀疏矩阵方法进行详细阐释，包括以下步骤，与图2说明中重复的部分不再说明。

步骤S1：程序接受网格矩阵A，输入向量B，求解向量X。将规模为M×N×L的三维结构化网格矩阵A划分为多个P₁×P₂×L大小的矩阵，从而可以调度全部P₁×P₂个处理器进行计算。与条件C计算不同的是，这里需要进行对角划分。令CPE(rid，cid)在各个子块负责的柱的行列编号在迭代中用变量i，j表示，用a∶b∶c表示首项为a，公差为b，末项大小不大于c的等差数列，则

i＝rid∶P∶M-1，

j＝cid-rid×δ_J∶P∶N-1+δ_J.

注意若j的值超过了矩阵规模N的范围，此时该处理器在这一柱上不作计算和写操作，仅进行数据的接收和发送。

步骤S2：开辟大小为buf[(2δ_Iδ_J+δ_I+δ_J+1)×L]大小的缓冲区，分别用来存放2δ_Iδ_J+δ_I+δ_J个依赖柱数据以及自身的数据。

同时在开辟缓冲区时，要将右端向量B[i，j，0∶L]读入并存储在buf[(2δ_Iδ_J+δ_I+δ_J)×L∶L]中，并将其余部分初始化为0。

步骤S3：当处理器位于处理器阵列的第一行时，需要用DMA方式将前一行的某个依赖柱X[i′，j′，k]读入，其中x′，j′满足条件C的限制。

当处理器当前计算的柱位于网格矩阵边缘时，其依赖的前一行或几行在网格矩阵之外，此时将缓冲区对应部分置0即可。

步骤S4：通过处理器间通信的方式获取部分或全部依赖数据。值得注意的是，即使当前处理器位于处理器阵列第一列，也可能会接收前一个子块计算完成后阵列中最后一列处理器发送的依赖数据。这样，阵列第一列的处理器获取其前一列或前几列的依赖数据时无需进行DMA操作。具体形式如下，这里无需详细说明发送数据的先后关系。

buf[0∶(2δ_Iδ_J+δ_I+δ_J)L]←getr/getc.

步骤S5：使用已经获取的依赖数据对X[i，j，k]进行计算，这里不同网格的计算方式和依赖关系不同，对计算过程不作描述。

步骤S6：通过处理器间通信的方式将计算得到的数据发送给邻居。具体形式如下，这里无需详细说明发送数据的先后关系。

putr(buf[L∶(2δ_Iδ_J+δ_I+δ_J)L]，(cid+1)mod P)，

putc(buf[L∶(2δ_Iδ_J+δ_I+δ_J)L]，(rid+1)mod P).

当处理器位于阵列的第一列时，不作行通信；当j＝＝N-1+δ_J，不作列通信。因为以上没有对应的接收情况。

步骤S7：处理器完成整个柱的计算后，将所计算的柱传入内存的X数组中，具体传递关系为为X[i，j，0∶L]←buf[(2δ_Iδ_J+δ_I+δ_J)L∶L]。这里默认i，j在网格规定的范围内。

图4为本发明的一个实施例在求解满足条件B和条件C的结构化网格过程示意图。该实施例所在的并行计算机系统中存在8×8＝64个处理器。本发明提供的方法将整个求解向量分解为若干个子块，每一个子块的大小为8×8×L，在每一个子块中有8+8+L-2＝14+L个层次集合。每一个处理器计算大小为L的一柱，则一个子块正好可以分配给64个处理器，然后处理器阵列依次遍历所有的子块。在每一个子块的计算中，由CPE(0，0)开始计算第0层(仅有一个未知数)，然后CPE(0，0)、CPE(0，1)和CPE(1，0)分别处理第1层的三个未知数，再之后CPE(0，0)、CPE(0，1)、CPE(1，0)、CPE(0，2)、CPE(1，1)和CPE(2，0)开始处理第2层。整个计算过程就像是在立方体的顶角推倒了一个多米诺骨牌。

图5(a)-(d)为本发明的一个实施例在处理满足条件D的结构化网格过程中提供的通信路径规约和对角划分方法示意图。该实施例所在的并行计算机系统中存在4×4＝16个处理器，本发明提供的方法将整个求解向量分解为若干个子块，每一个子块的大小为4×4×L。求解的结构化网格矩阵满足条件D，且δ_I＝δ_J＝1。参考图5(a)，这里展示M-N平面的示意图，通信是发生在L维的未知数计算时。每一个处理器至多会从相邻的4个处理器获得数据，并且向其它4个相邻处理器发送数据。通过仔细分析可以发现，这里存在了一些冗余的通信路径。例如，CPE(1，1)依赖CPE(0，0)、CPE(0，1)、CPE(0，2)和CPE(1，0)，而CPE(0，1)依赖CPE(0，0)，CPE(0，2)依赖CPE(0，1)。如果去除CPE(0，0)→CPE(1，1)和CPE(0，1)→CPE(1，1)的通信，整个计算的依赖关系仍然可以满足。同理可以去除所有冗余的通信路径，而只保留剩下的关键路径。如图5(c)所示，红色的箭头表示需要保留的关键路径，而黑色的箭头表示可以去除的冗余路径。此时，每个处理器只需从其他至多两个处理器获得数据，并向其它至多两个处理器发送数据即可。图5(b)示出了不进行对角划分时每个从核负责计算的柱，图5(d)示出了对角划分后每个从核负责计算的柱，当处理器分配到超出矩阵规模范围的柱时，只进行通信，不进行计算

根据本发明实施例的方法，通过处理器间的快速通信传输每次运算的依赖数据，最小化访存开销，最大化计算资源利用率，从而对各种形式的结构网格稀疏下三角矩阵方程进行高效求解。

图6为本发明的一个实施例在神威太湖之光超级计算机上运行的性能测试结果示意图，其中MPE代表单处理器运行，CPE代表多处理器并行运行。根据本发明提供的算法，对3D7和3D19结构网格稀疏下三角矩阵进行求解。测试网格数从40×40×40增大到200×200×200，按理论数据访存量计算带宽：

对于3D7矩阵，平均带宽可达19.1GB/s。在最大网格数的测试下，带宽达21.3GB/s。与StreamTriad测试基准相比，平均带宽利用率达84.4％。与单机的串行算法相比，平均加速了15.5倍。对于3D19矩阵，平均带宽可达21.2GB/s。在最大网格数的测试下，带宽达22.1GB/s。与Stream Triad测试基准相比，平均带宽利用率达93.7％。与单机的串行算法相比，平均加速了17.7倍。本文提出的结构网格稀疏下三角矩阵求解算法没有额外的数据访问，性能的主要损失在于算法启动与结束阶段时，只有一部分处理器在工作。

根据本发明实施例的方法，通过处理器间的快速通信传输每次运算的依赖数据，最小化访存开销，最大化计算资源利用率，从而对各种形式的结构网格稀疏下三角矩阵方程进行高效求解。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (10)

1.一种在高性能并行计算环境中求解结构化网格稀疏下三角方程的计算过程中优化处理器间通信和高速缓冲存储器cache空间利用的方法，其特征在于，包括以下步骤：

接收结构化网格的网格规模，接收输入网格，输入向量和求解向量，并启动多处理器开始并行执行；

将整个求解向量分解为多个子块，将每个子块划分为多个柱，每个柱代表求解向量在内存地址空间中连续的一段，将多处理器逻辑上视为矩形阵列，每个处理器负责计算一个柱，其中每个处理器在自身的cache中开辟用于存储依赖数据的缓冲区，从邻居处理器或者自身内存得到计算所需要的数据并存入所述缓冲区，完成自己所在柱的计算，并将计算完成的数据发送给计算过程依赖所述计算完成的数据的其它处理器。

所述处理器矩形阵列依次遍历所有的子块，并完成计算所需的数据通信，对求解向量中的每个子块进行计算，直到完成整个结构化网格的求解计算。

2.根据权利要求1所述的方法，所述结构化网格矩阵是规模为M×N×L的三维结构化网格矩阵，将整个求解向量分解为若干P₁×P₂×L的子块，其中M、N、L为大于等于1的整数，0<P₁≤M,0<P₂≤N,且P₁×P₂不超过并行系统中处理器的个数，每个处理器分别计算大小为L的连续空间的一柱，且只有每一行和每一列之间的处理器可以相互通信。

3.根据权利要求2所述的方法，处理器个数为P₁×P₂个，在不违背构造矩阵阵列方法原则的情况下，使|P₁-P₂|尽可能小，使得阵列形状更接近于正方形。

4.根据权利要求3所述的方法，矩阵阵列构造中，将通信代价小的处理器放在同一行/列。

5.根据权利要求3的方法，启动多处理器开始并行执行包括：启动处理器组中所有可用处理器资源，每个线程绑定一个处理器资源。

6.根据权利要求1到5任一项所述的方法，其中

当处理器位于处理器阵列前部边缘时，将对应计算依赖数据的缓冲区置为全0，或通过DMA将位于当前子块外的依赖数据读入到缓冲区中；当处理器不位于处理器阵列前部边缘时，接受先行完成计算的邻居处理器的通信数据，作为其计算的依赖数据存入缓冲区中；以及

一个处理器完成计算后，通过处理器间的通信将计算得到的数据传送给依赖此数据的邻居或者非邻居处理器，并将结果写入内存。

7.如权利要求6所述的方法，其特征在于，一个处理器将数据传送给非邻居处理器进一步包括：

所述处理器会经过一次或几次中间结点地将数据传送给该非邻居处理器，以及对通信节点之间的关系进行通信路径规约和对角划分过程，以避免死锁。

8.如权利要求7所述的方法，其特征在于，所述通信路径规约和对角划分过程进一步包括：

根据通信网络的拓扑结构将通信边划分为关键边和非关键边，仅保留关键边规定的通信过程；

规约后每个处理器只与两个邻居进行通信。根据对角划分的方法确定各处理器负责计算的柱在求解向量中的位置。划分完成后，处理器阵列中同一列处理器负责的柱会在水平面上沿一条斜向直线排列，从而避免出现通信环路。这种位置计算方式会使部分处理器分配到不在输入矩阵规模范围内的柱，此时处理器负责通信而没有计算任务。

9.根据权利要求1到8任一项所述的方法，所述输入网格，输入向量和求解向量进一步包括：输入网格是一个规模为M×N×L的三维结构化网格矩阵A，用一个M×N×L×rnz的一维数组表示，其中，M×N×L表示网格空间的大小，代表求解的矩阵共有M×N×L行和列，每一行有rnz个非零元，其中有一个非零元分布在矩阵的对角线上。由于输入网格是结构化网格，每一行的非零元分布有固定的位置关系。输入向量是一个规模为M×N×L的向量B，求解向量是一个规模为M×N×L的向量X，设A(i,j,k,r)(0≤r<rnz)依赖的数据位置为(DI_r,DJ_r,DK_r)，其中假定0号非零元代表矩阵对角线上的元素，则求解计算关系如下：

其中，每一行的非零元分布有一定的位置关系包括：对于网格中任意两个不同的位置(i,j,k)和(i′,j′,k′)，有

DJ_r-j＝DJ′_r-j′,DK_r-k＝DK′_r-k′.

同时，非零元的相对位置分布满足预定的条件限制。

10.一种计算机可读存储介质，其上存储有计算机指令，所述指令当被高性能并行计算机执行时，执行权利要求1到9任一项的方法。

Images