CN110704023B

CN110704023B - 一种基于拓扑排序的矩阵分块划分方法及装置

Info

Publication number: CN110704023B
Application number: CN201910916488.5A
Authority: CN
Inventors: 郭根华; 周振亚; 王晓光; 徐启迪
Original assignee: Beijing Empyrean Technology Co Ltd
Current assignee: Beijing Empyrean Technology Co Ltd
Priority date: 2019-09-26
Filing date: 2019-09-26
Publication date: 2021-10-22
Anticipated expiration: 2039-09-26
Also published as: CN110704023A

Abstract

一种基于拓扑排序的矩阵分块划分方法，包括：1）根据LU分解规则，求出矩阵各列之间的依赖关系；2）对所述依赖关系进行拓扑排序，将所有列分成至少一层；3）根据层次关系，对所述矩阵列的位置进行重排；4）对位置重排后的新矩阵进行分块划分；5）按层级顺序对同一层级的所述分块进行并行计算。本发明的一种基于拓扑排序的矩阵分块划分方法，可以更加有效地利用计算机的CPU/GPU资源，提高矩阵并行计算的能力。

Description

一种基于拓扑排序的矩阵分块划分方法及装置

技术领域

本发明涉及集成电路设计技术领域，特别是涉及一种基于拓扑排序的矩阵分块划分方法。

背景技术

随着超大规模集成电路的发展，在电路方程的求解过程中，最耗时的是线性方程组的求解。因此，对矩阵进行分块和并行计算，成为了一种广泛应用的仿真加速方法。

常用的线性方程组的求解方法是直接法，即对系数矩阵进行LU分解。而电路仿真产生的矩阵是稀疏矩阵。传统的处理稀疏矩阵LU分解的方法是对矩阵进行划分，或者说对矩阵进行分块，然后按照分块结果，对矩阵进行压缩。压缩之后的矩阵的一个元素对应原始矩阵中的一个分块小矩阵。这种划分只考虑了矩阵的稀疏性，划分出来的各个子矩阵大小不一，在并行计算的时候都要分别占用一个单独的任务，即耗费一个CPU/GPU计算单元，因此并行效率不高。

发明内容

为了解决现有技术存在的不足，本发明的目的在于提供一种基于拓扑排序的矩阵分块划分方法，可以更加有效地利用计算机的CPU/GPU资源，提高矩阵并行计算的能力。

为实现上述目的，本发明提供的基于拓扑排序的矩阵分块划分方法，包括以下步骤：

1）根据LU分解规则，求出矩阵各列之间的依赖关系；

2）对所述依赖关系进行拓扑排序，将所有列分成至少一层；

3）根据层次关系，对所述矩阵列的位置进行重排；

4）对位置重排后的新矩阵进行分块划分；

5）按层级顺序对同一层级的所述分块进行并行计算。

进一步地，所述步骤1）进一步包括，

对尺寸为n的矩阵进行LU分解，分解为下三角矩阵L和上三角矩阵U；

数组a表示存储的元素，a[i][j]表示第i+1行第j+1实列上的元素，其中，0<=i<n，0<=j<n；

所述矩阵列号为m，其中0<=m<n，“->”或“<-”表示依赖的方向；

当m=0时，第1列不依赖任何列；

当m>0时，求第m+1实列所依赖实列的列号，其步骤如下：

（1001）将U矩阵中位于第m+1实列中非零元素的行号组成行范围，所述行范围为[r1，r2]，其中，0<=r1<=r2<m；

（1002）设i为自变量，从r1开始循环，到r2结束，称为外循环；

（1003）若a[i][m]=0，将i的值增加1，进行下一次外循环；若a[i][m]=1，跳到步骤（1004）；

（1004）将L矩阵中位于第i+1列非零元素的行号组成行范围，所述行范围为[r3，r4]，其中，0<=r3<=r4<m；

（1005）若r4<m，将i的值增加1，进行下一次外循环；否则，跳到步骤（1006）；

（1006）若a[m][i]=1，则依赖关系成立，即m->i，然后将i的值增加1，进行下一次外循环；否则，跳到步骤（1007）；

（1007）将L矩阵中位于第m+1列非零元素的行号组成行范围，所述行范围为[r5，r6]，其中，m<=r5<=r6<n；

（1008）若r4<r5，将i的值增加1，进行下一次外循环；否则，跳到步骤（1009）；

（1009）设自变量为j，从r5开始循环，到r4结束，称为内循环1；设另一自变量为k，从r5开始循环，到r6结束，称为内循环2；内循环1和内循环2同时开始，跳到步骤（1010）；

（1010）若j==k，a[j][i]=1且a[k][m]=1，则内循环1和内循环2都终止，依赖关系成立，即m->i，然后将i的值增加1，进行下一次外循环；否则，跳到步骤（1011）；

（1011）若j>k，则k=k+1；若j<k，则j=j+1；若j==k，且a[j][i]!=1或 a[k][m] !=1，则k=k+1且j=j+1，跳到步骤（1012）

（1012）若j<=r4且k<=r6，返回步骤（1010）继续比较，否则内循环1和内循环2都终止，将i的值增加1，进行下一次外循环。

进一步地，当第m+1实列不依赖于第i+1实列，而第m+1行依赖于第i+1行时，所述第m+1列依赖于所述第i+1列，其中，0<=i<m。

进一步地，所述步骤2）进一步包括，将具有依赖关系的列分为N层，层号从下至上分别为1，2，3，...，N，其中，

第1层各个列不依赖任何列；

第2层各个列，至少依赖第1层中的一个列；

第n层各个列依赖第1层至第（n-1）层中至少一列。

进一步地，所述步骤3）进一步包括以下步骤：

（501）将第1层第1列置于新矩阵的第1列；

（502）将第1层的第n列置于新矩阵的第n列；

（503）通过行列变换调整第n列中非零元素的位置使第n列与第1列至第（n-1）列不存在依赖关系；

（504）对第n层按（501）至（503）中的步骤进行重排。

进一步地，所述新矩阵从左到右存在层次关系，每个层次内各列之间不存在依赖关系，第1层不依赖任何列，第n层各列依赖第1层至第（n-1）层中至少一列。

进一步地，所述步骤4）进一步包括，

从第1层开始按层级顺序和设定的粒度大小对所述新矩阵进行分块划分。

进一步地，所述步骤5）进一步包括，除第1层外，当前层需要前一层所有分块计算完成之后开始当前层的计算。

为实现上述目的，本发明还提供一种计算机可读存储介质，其上存储有计算机指令，所述计算机指令运行时执行如上文所述的基于拓扑排序的矩阵分块划分方法步骤。

为实现上述目的，本发明还提供一种基于拓扑排序的矩阵分块划分装置，包括存储器和处理器，所述存储器上储存有在所述处理器上运行的计算机指令，所述处理器运行所述计算机指令时执行如上文所述的基于拓扑排序的矩阵分块划分方法步骤。

本发明的一种基于拓扑排序的矩阵分块划分方法，具有以下有益效果：

1）求出了各列之间的依赖关系，将它们进行拓扑排序和位置重排后，按层进行分块划分，最后得到的各个分块，绝大部分都比较均匀。

2）位置重排后，新矩阵中列数极少的层都分布于尾部，划分后对应的碎块也就都排在了后面，根据流水线的工作原理，将时间短的任务排在最后可以缩短总的运行时间。

3）不仅大大减少了碎块产生的数量，还将大部分碎块的计算任务移到了最后，从而使得前期运行的大部分并行的任务都处于饱满或接近饱满的状态，最大限度地利用了CPU/GPU的计算资源，提高了并行计算的效率。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，并与本发明的实施例一起，用于解释本发明，并不构成对本发明的限制。在附图中：

图1为根据本发明的基于拓扑排序的矩阵分块划分方法流程图；

图2为根据本发明的矩阵的LU表示法以及对应的L矩阵和U矩阵；

图3为根据本发明一实施例中重排后的新矩阵A’与原始矩阵A的对比示意图；

图4为根据本发明一实施例中新矩阵A’和原始矩阵A进行分块划分后的结果对比示意图。

具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

图1为根据本发明的基于拓扑排序的矩阵分块划分方法流程图，下面将参考图1，对本发明的基于拓扑排序的矩阵分块划分方法进行详细描述。

为了便于叙述，本发明中的“矩阵”都指代方阵，即行和列数目是相等的。另外，为了避免重复描述，除非特别说明，本发明中的“列”都指代“行或者列”，而“实列”都特指狭隘意义上与行对应的列。

首先，在步骤110，根据LU分解规则，求出矩阵各列之间的依赖关系。该步骤中，对矩阵A进行LU分解，得到单位下三角矩阵L（即L的对角线元素为1）和上三角矩阵U。为了便于描述，我们把L矩阵和U矩阵可以合并成一个方阵，下三角是L，上三角是U。上述表示方法我们称为LU表示法。本发明所指的依赖关系，就是在这种表示法的基础上来建立的。

优选地，各列之间的依赖关系，是根据LU分解法的规则，得出的各列之间的运算先后顺序。该步骤中，假设一个尺寸大小为n的方阵A，其对应的列号分别为1，2，3，…，n。用数组a表示方阵中元素的存储，其中a[i][j]（0<=i<n，0<=j<n）表示第i+1行第j+1实列上的元素，它的值为1表明是非零元素，为0表示是零元素。对于A中的任意一列，假设列号为m(0<=m<n)，其依赖的所有列的列号必定会小于m。在本发明的叙述中，用箭头“->”或“<-”表示依赖的方向，例如j<-m或m->j都表示第m+1列依赖于第j+1列。

优选地，当m=0时，第1列不依赖其他任何列；当m>0时，求第m+1实列所依赖的所有实列的列号，其步骤如下：

步骤1101，找到位于第m+1实列对角线元素以上位置的所有非零元素（都位于U矩阵中），得到一个由这些非零元素的行号组成的行范围，假设这个范围为[r1，r2]，其中0<=r1<=r2<m。

步骤1102，从r1开始，到r2时结束，假设i为循环的自变量，开始进行外循环。

步骤1103，如果a[i][m]=0，则本次循环结束，将i的值增加1，进行下一次循环；如果a[i][m]=1，跳到步骤1104。

步骤1104，找到位于第i+1列对角线元素以下位置的所有非零元素(都位于在L矩阵中)，这些非零元素的行号将组成一个行范围，假设为[r3，r4]，其中i<r3<=r4<n。

步骤1105，如果r4<m，则本次循环结束，将i的值增加1，进行下一次外循环；否则，跳到步骤1106。

步骤1106，如果a[m][i]=1，则本次循环结束，且依赖关系成立，即m->i，将i的值增加1并进行下一次外循环；否则，跳到步骤1107。

步骤1107，找到位于第m+1列对角线元素以下位置的所有非零元素（都位于L矩阵），得到一个由这些非零元素的行号组成的范围，假设这个范围为[r5，r6]，其中m<r5<=r6<n。

步骤1108，如果r4<r5，则本次循环结束，将i的值增加1，进行下一次循环；否则，跳到步骤1109。

步骤1109，同时开始两个循环，为了与步骤1102中的循环区分，分别称为内循环1和内循环2。假设内循环1的自变量为j，从r5开始，到r4结束；内循环2的自变量为k，从r5开始，到r6结束。如果循环过程中发现满足j=k且a[j][i]=1且a[k][m]=1，则内循环1和内循环2都终止，本次外循环也结束，依赖关系成立，即m->i，进行下一次外循环；否则，跳到步骤1110。

步骤1110，在不满足步骤1108的前提下，如果j>k，则k=k+1;如果j<k，则j=j+1，如此下去，直到内循环1和内循环2都结束。

优选地，求第m行所依赖的所有行的方法与上述求第m列所依赖的所有列的步骤类似，因此，不重复阐述。

优选地，根据矩阵运算的规则，要得到第m+1列所依赖的其他所有列的列号，既要看第m+1实列的依赖关系，也要看第m+1行的依赖关系。假设第i+1列是第m+1列之前的一列，即0<=i<m，如果第m+1实列不依赖于第i+1实列，而第m+1行依赖于第i+1行，那么依赖关系m->i仍然成立，即第m+1列仍然依赖于第i+1列。也就是说，在求某一列的依赖关系时，最后的结果要将对应实列的依赖关系与行的依赖关系的进行与操作。

在步骤120，对依赖关系进行拓扑排序，将所有列分为多个层。该步骤中，矩阵各列之间的依赖关系建立之后，下一步就是进行拓扑排序。拓扑排序将具有依赖关系的各列分为多个层。假设共N层，层号从下往上分别为1，2，3，...，N。第1层中各个列不依赖其他任何列，是可以直接进行计算的。第2层中的各列，肯定依赖第1层中某列或某些列，因此，必须等其依赖的第1层中的各列运算全部完成之后才能开始本层的运算。第3层的各列可能直接依赖第1层中的某列或某些列，或者依赖第2层中的某列或某些列（因此也间接地依赖第1层中的某列或某些列），因此，必须等其直接和间接依赖的各列运算全部完成之后，才能开始本层的运算。后续各层的依赖关系以此类推。

在步骤130，根据层次关系，对矩阵的某些列的位置进行重排。该步骤中，为了与原始矩阵A区分，将位置重排后的矩阵称为新矩阵A’。重排的基本步骤如下：

步骤1301，将第1层的第1列排在新矩阵A’的第1列。

步骤1302，如果第1层不止一列，则将第1层的第2列置于新矩阵A’的第2列。由于第1层的第2列并不一定是原始矩阵A的第2列，因此，为了保证与新矩阵A’的第1列没有依赖关系，可能需要通过多次行列变换来调整该列中某些非零元素的位置。由于通过行列变换进行元素调整是一种基本的矩阵操作，因此，本发明不作具体阐述。

步骤1303，如果第1层还有更多的列，继续按步骤1302进行处理。

步骤1304，将后面的各层按第1层的方法进行类似处理。

优选地，重排后的新矩阵，位置调整过的列仍然对其所属层中前面的列没有任何依赖关系。该步骤中，位置重排后的新矩阵A’，从左到右依然存在层次关系。每个层次内的各列之间没有依赖关系，但是，除了第1层之外，任一层都一定依赖于其前面一层或几层的某些列。

在步骤140，对位置重排后的新矩阵按层进行分块划分。

优选地，分块划分必须按层级顺序进行，即第一层开始，依次往后；每层经过划分之后，至多有一个碎块，也就是该层的最后一个分块，分块划分完成后。该步骤中，经过位置重排后得到的新矩阵A’，从左到右依次是拓扑排序后第1层的列号，第2层的列号，第3层的列号，…，等等。当然，这些列的元素位置，经过多次行列变换之后，可能与原始列不完全相同。但是，各层之间的依赖关系没有改变。由于除了第1层之外，每一层的列都依赖于前面一层或几层的列，因此，分块划分必须按层进行，并且从新矩阵A’的左侧开始，也就是第1层开始。

优选地，根据粒度大小对层进行分块划分，即从第一列开始，每次拆分出一个粒度大小的块，剩余的继续类似处理，直到该层所有的列都划分完毕，大部分块的大小都等于或接近粒度大小。该步骤中，假设对第i层进行分块划分，这个层有m列，分块的粒度大小为g。如果m%g==0，则可以分为m/g个块，否则，块的数量为m/g+1，其中最后一个块的大小为m%g，其列数肯定小于粒度大小。如果最后一个块的列数非常少，则可能成为碎块。可以看出，新矩阵A’经过划分之后，每一层至多有一个碎块。由于经过拓扑排序之后，一般列数很少的层主要分布于尾部，因此，产生的大部分碎块也都位于最后。

优选地，根据一个CPU/GPU核的计算处理单元数来设置分块的粒度大小。该步骤中，分块的粒度大小，它指的是一个分块所包含的列数。这个值理论上可以自由设定，但一般会根据一个CPU/GPU核的计算处理单元数来设置，例如16或32或64等，这样可以使得一个任务正好占满一个CPU/GPU核的所有计算处理单元，最大限度地利用到计算机的硬件资源，提高并行的效率。

在步骤150，按层级顺序对同一层级的所述分块进行并行计算。该步骤中，每层按设定的粒度大小进行划分后，得到的所有分块都可以同时进行计算；除了第一层之外，当前层必须等前一层的所有分块都计算完成之后，本层的分块才可以开始计算。

下面结合一具体实施例对本发明的一种基于拓扑排序的矩阵分块划分方法做进一步的说明。

首先，为了阐述的方便，我们以{x₁，x₂，x₃，…，x_n}的形式表示一个集合，元素之间用逗号隔开，且其位置代表先后顺序。{}则表示集合为空。下面以n=6的一个方阵A为例，对A矩阵基于拓扑排序进行分块划分的过程进行举例说明。

图2为根据本发明的矩阵的LU表示法以及对应的L矩阵和U矩阵，如图2所示，划分的具体步骤如下：

步骤一，求LU分解后矩阵的依赖关系，过程如下：

a) 第1列肯定不依赖任何列，因此，既没有实列依赖关系，也没有行依赖关系，综合而得：0->{}。

b) 第2列，先看实列依赖关系：第2实列对角线上面的元素只有1个，为{0}，因此，没有实列依赖关系。再看行依赖关系：第2行对角线上左边的元素只有一个，为{0}，因此，也没有行依赖关系。综合而得：1->{}。

c) 第3列，先看实列依赖关系：第3实列对角线上面的元素有2个，为{0，0}，因此，没有实列依赖关系。再看行依赖关系：第3行对角线左边的元素有两个，为{1，0}，继续看是否依赖第1行。因为a[0][2]=0，继续看U矩阵中第1行和第3行的元素中，是否有列号相等且都为非零的元素。U矩阵中，第1行的元素中列号大于2的所有元素为{1，0，0}，第3行的所有元素为{0，0，1}，不符合依赖要求，因此，没有行依赖关系。综合而得：2->{}。

d) 第4列，先看实列依赖关系：第4实列对角线上面的元素有3个，为{1，1，0}，继续看是否依赖第1实列和第2实列。对于第1实列，因为a[3][0]=1，所以，第4实列依赖于第1实列，即3->0。再看第2实列，因为a[3][1]=0，继续看L矩阵中第2实列和第4实列的元素中，是否有行号相等且都为非零的元素。L矩阵中，第2实列的元素中列号大于3的所有元素为{0，0}，第4实列的所有元素为{0，1}，因此，不符合依赖要求。所以，最终的实列依赖关系只有3->{0}。再看行依赖关系：第4行对角线左边的元素有3个，为{1，0，1}，需要看是否依赖第1行和第3行。因为已知第4实列依赖于第1实列，所以，这里只需要看第4行是否依赖第3行就行了。因为a[2][3]=0，继续看U矩阵中第3行和第4行的元素中，是否有列号相等且都为非零的元素。U矩阵中，第3行的元素中列号大于3的所有元素为{0，1}，第4行的所有元素为{1，1}，因为a[3][5]=1且a[4][5]=1，符合依赖要求，所以，行依赖关系存在，即3->2。综合而得：3->{0，2}。

e) 第5列，先看实列依赖关系：第5实列对角线上面的元素有4个，为{0，0，0，1}，因此，只需要看是否依赖于第4实列，因为a[4][3]=0，继续看L矩阵中第4实列和第5实列的元素中，是否有行号相等且都为非零的元素。L矩阵中，第4实列的元素中列号大于4的所有元素为{0}，第5实列的所有元素为{1}，不符合依赖要求。所以，没有实列依赖关系；再看行依赖关系：第5行对角线左边的元素有4个，为{0，0，0，0}，因此，没有行依赖关系。综合而得：4->{}。

f) 第6列，先看实列依赖关系：第6实列对角线上面的元素有5个，为{0，0，1，1，0}，需要看是否依赖第3实列和第4实列。先看第3实列，因为a[5][2]=1，所以第6实列依赖于第3实列，即5->2。再看第4实列，因为a[5][3]=0，继续看L矩阵中第4实列和第6实列的元素中，是否有行号相等且都为非零的元素。L矩阵中，第4实列的元素中列号大于5的所有元素为{}，第6实列的所有元素为{}，不符合依赖要求。所以，实列依赖关系为：5->{2}。再看行依赖关系：第6行对角线左边的元素有5个，为{0，0，1，1，0}，需要看是否依赖第3行和第4行。因为已知第6实列依赖于第3实列，因此，只需要看是否依赖于第4行即可。因为a[3][5]=1，所以第6行依赖于第4行，即5->3。综合而得，5->{2，3}。

g) 最后，综合所有列的依赖关系，得到一个依赖关系表，如下：

0->{}

1->{}

2->{}

3->{0，2}

4->{}

5->{2，3}

步骤二，进行拓扑排序，得到多个层。拓扑排序的过程如下：

a) 遍历第一步中的矩阵依赖关系表，把没有依赖关系的列号都提出来，放在第一层，即{0，1，2，4}，同时更新依赖关系列表，把第1层中的列号都去掉，此时依赖关系表就变成了：3->{}，5->{3}。

b) 遍历更新后的矩阵关系表，将没有依赖关系的元素提出来，放在下一层，因此，第2层的元素就产生了：{3}，同时进一步更新依赖关系表，把第2层中的列号都去掉，此时依赖关系表就变成了：5->{}。

c) 遍历更新后的矩阵关系表，将没有依赖关系的元素提出来，放在下一层，因此，第3层的元素就产生了：{5}。同时进一步更新依赖关系表，把第3层中的元素都去掉，此时依赖关系表变为了：{}。

d) 因为依赖关系表已经为空，拓扑排序结束。

排序后，层数共为3层，具体结果如下：

第1层：{0，1，2，4}

第2层：{3}

第3层：{5}

步骤三，对A矩阵进行重排，得到一个新矩阵A’。

新矩阵A’初始状态没有任何列。重排的步骤如下：

a) 将第1层的第1列加入到A’中。

b) 将第1层的第2列加入到A’中，因为第1层的第2列也是原始矩阵A中的第2列，因此，不需要做行列变换或位置调整。

c) 将第1层的第3列加入到A’中，因为第1层的第3列也是原始矩阵A中的第3列，因此，不需要做行列变换或位置调整。

d) 将第1层的第4列加入到A’中，因为第1层的第4列是原始矩阵A中的第5列，因此，可能需要做行列变换或位置调整。

e) 将第2层的第1列加入到A’中。

f) 将第3层的第1列加入到A’中。

图3为根据本发明一具体实施例中重排后的新矩阵A’与原始矩阵A的对比示意图。如图3所示，新矩阵A’中的列号旁边括起来的数字表示在原始矩阵A中对应的列号，即重排之前的列号。

步骤四，对重排后的新矩阵A’按层进行分块划分。

图4为根据本发明一具体实施例中新矩阵A’和原始矩阵A进行分块划分后的结果对比示意图，如图4所示，因为矩阵A的规模比较小（n=6），为了便于说明，将分块的粒度大小设为2。新矩阵A’的第1层可以拆分为两个块，即第1列和第2列为一个分块，第3列和第4列为一个分块。第2层和第3层都只有一列，分别作为一个单独的分块，因此，总共4个分块，并行运行时的任务数分别为2、1、1。而如果按传统方法，A矩阵的第一层为{0，1，2}，第二层为{3}，第三层为{4}，第四层为{5}。第一层划分为两块，即第1列和第2列为一个分块，第3列为一个分块，接下来第2、3、4层各为一个分块。虽然第一层的两个分块和第三层的那个分块可以并行计算，但是需要占用3个任务，所以，总共需要5个任务，并行运行时的任务数分别为3、1、1。

本发明提出了一种基于拓扑排序的矩阵分块划分方法，其基本原理为：根据LU分解规则，求出矩阵各列之间的依赖关系，接下来对依赖关系进行拓扑排序，将列划分为多个层，然后根据层次关系，对矩阵中某些列的位置进行重排，形成一个新的矩阵，最后对这个新矩阵按层进行分块划分，得到大量尺寸比较均匀的分块。可以更加有效地利用计算机的CPU/GPU资源，提高矩阵并行计算的能力。

本领域普通技术人员可以理解：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims

1.一种基于拓扑排序的矩阵分块划分方法，应用于电路仿真领域，包括以下步骤：

1）根据LU分解规则，求出矩阵各列之间的依赖关系；

所述步骤1）进一步包括，

当m=0时，第1列不依赖任何列；

当m>0时，求第m+1实列所依赖实列的列号，其步骤如下：

（1012）若j<=r4且k<=r6，返回步骤（1010）继续比较，否则内循环1和内循环2都终止，将i的值增加1，进行下一次外循环；

2）对所述依赖关系进行拓扑排序，将所有列分成至少一层；

3）根据层次关系，对所述矩阵列的位置进行重排；

4）对位置重排后的新矩阵进行分块划分；

所述步骤4）进一步包括，根据一个CPU/GPU核的计算处理单元数来设置分块的粒度大小，分块的粒度大小指的是一个分块所包含的列数，使得一个任务正好占满一个CPU/GPU核的所有计算处理单元，最大限度地利用到计算机的硬件资源；

5）按层级顺序对同一层级的所述分块进行并行计算。

2.根据权利要求1所述的基于拓扑排序的矩阵分块划分方法，其特征在于，当第m+1实列不依赖于第i+1实列，而第m+1行依赖于第i+1行时，所述第m+1列依赖于所述第i+1列，其中，0<=i<m。

3.根据权利要求1所述的基于拓扑排序的矩阵分块划分方法，其特征在于，所述步骤2）进一步包括，将具有依赖关系的列分为N层，层号从下至上分别为1，2，3，...，N，其中，

第1层各个列不依赖任何列；

第2层各个列，至少依赖第1层中的一个列；

第n层各个列依赖第1层至第（n-1）层中至少一列。

4.根据权利要求1所述的基于拓扑排序的矩阵分块划分方法，其特征在于，所述步骤3）进一步包括以下步骤：

（501）将第1层第1列置于新矩阵的第1列；

（502）将第1层的第n列置于新矩阵的第n列；

（504）对第n层按（501）至（503）中的步骤进行重排。

5.根据权利要求4所述的基于拓扑排序的矩阵分块划分方法，其特征在于，所述新矩阵从左到右存在层次关系，每个层次内各列之间不存在依赖关系，第1层不依赖任何列，第n层各列依赖第1层至第（n-1）层中至少一列。

6.根据权利要求1所述的基于拓扑排序的矩阵分块划分方法，其特征在于，所述步骤4）进一步包括，

7.根据权利要求1所述的基于拓扑排序的矩阵分块划分方法，其特征在于，所述步骤5）进一步包括，除第1层外，当前层需要前一层所有分块计算完成之后开始当前层的计算。

8.一种基于拓扑排序的矩阵分块划分装置，其特征在于，包括存储器和处理器，所述存储器上储存有在所述处理器上运行的计算机指令，所述处理器运行所述计算机指令时执行权利要求1至7任一项所述的基于拓扑排序的矩阵分块划分方法步骤。

9.一种计算机可读存储介质，其上存储有计算机指令，其特征在于，所述计算机指令运行时执行权利要求1至7任一项所述的基于拓扑排序的矩阵分块划分方法步骤。