CN103399841A - 基于gpu的稀疏矩阵lu分解方法 - Google Patents

Abstract

本发明提出一种基于GPU的稀疏矩阵LU分解方法，包括以下步骤：A：在CPU上建立待仿真的电路网单的电路矩阵，并对该电路矩阵进行预处理和稀疏化；B：根据预处理结果选择对稀疏矩阵LU进行处理的计算平台；C：如果使用的计算平台为GPU，则通过GPU对稀疏矩阵LU进行分解；D：通过CPU对稀疏矩阵LU中非零元的值进行回代求解，以完成对稀疏矩阵LU的分解。本发明的实施例可根据分解中浮点运算次数的多少自动选择计算平台，能够合理分配CPU和GPU上的任务，保证GPU上众多线程间的协同运行，从而提高运算速度，减少电路仿真时间。

Description

基于GPU的稀疏矩阵LU分解方法

技术领域

本发明涉及电子设计自动化与并行计算领域，特别涉及一种基于GPU的稀疏矩阵LU分解方法。

背景技术

稀疏矩阵LU分解是线性代数中的一项基本操作，在电路仿真、结构力学、经济建模等许多领域有广泛应用。在电路仿真中，电路被表示为矩阵，电路矩阵极其稀疏，通常每行只有约5个非零元，这是由在电路中（除少量电源、地、时钟节点外）一个节点不能连接太多元件这一原因决定的。求解电路的过程，也就是利用稀疏矩阵LU分解求解稀疏矩阵方程Ax=b的过程。在一次电路仿真中，稀疏矩阵LU分解位于牛顿—拉夫森迭代和瞬态迭代两层循环之中，需要迭代地进行很多次，这一步骤是当前最常用的电路仿真SPICE(SimulationProgram with Integrated Circuit Emphasis)的主要瓶颈。目前对超大规模集成电路的仿真的时间可能长达数周甚至数月。因此稀疏矩阵LU分解的并行化是一个急需解决的问题。

矩阵LU分解的目标是找到上三角矩阵L和下三角矩阵U使得

上述方程解为

$l_{\mathrm{ij}}=a_{\mathrm{ij}}-\underset{k=1}{\overset{i-1}{\mathrm{\Σ}}}{l}_{\mathrm{ik}}\·u_{\mathrm{kj}},j\≤i$

$u_{\mathrm{ij}}=(a_{\mathrm{ij}}-\underset{k=1}{\overset{j-1}{\mathrm{\Σ}}}{l}_{\mathrm{ik}}\·u_{\mathrm{kj}})/l_{\mathrm{jj}},j>i$

上式可以写成

$x_i=a_i-\underset{k=1}{\overset{i-1}{\mathrm{\Σ}}}{l}_{\mathrm{ik}}\·u_k$

l_ij=x_ij,j≤i

$u_{\mathrm{ij}}=\frac{x_{\mathrm{ij}}}{x_{\mathrm{jj}}},j>i$

当矩阵A是稀疏矩阵时，其LU因子矩阵L和U往往也是稀疏矩阵，此时，当且仅当lik≠0时，需要从xi中减去uk。

LU分解完成后，求解Ax=b的问题即变成求解Ly=b和Ux=y两个三角方程的问题，这一步称为回代求解。

理论上，将稀疏矩阵当作普通稠密矩阵进行LU分解，也可得到正确结果的。加州大学伯克利分校（University of California, Berkeley）、田纳西大学诺克斯维尔分校（The Universityof Tennessee, Knoxville）的研究人员已经能在GPU上高效的实现稠密矩阵的LU分解，在GTX 280上达到388十亿次每秒(Gflop/s,Giga Floating point OPerations per second)，相关研究可参见“V.Volkov and J.Demmel,“LU,QR and Cholesky factorizations using vectorcapabilities of gpus,”EECS Department,University of California,Berkeley,Tech.Rep.UCB/EECS-2008-49,May2008”和“S.Tomov,R.Nath,H.Ltaief,and J.Dongarra,“Dense linearalgebra solvers for multicore with gpu accelerators,”IEEE International Symposium on ParallelDistributed Processing Workshops and Phd Forum IPDPSW,pp.1–8,2010”。但是，简单计算可知这样做的效率是极低的。以稀疏矩阵onetone2（来自佛罗里达大学稀疏矩阵库）为例，该矩阵大小为36k*36k，对其进行稠密LU分解，即使性能到达1000Gflop/s，仍需15.5秒，而单核CPU上的稀疏LU分解只需要不到1秒。

稀疏矩阵LU分解的方法的研究已经活跃了几十年。目前学术界与工业界已有许多软件。稀疏矩阵LU分解的方法主要可分为两类。一类是基于稠密块的方法，采用这类方法的软件主要有SuperLU和Pardiso（其中Pardiso有GPU版本）。这类方法通过行列调换在矩阵中形成稠密子块，并对这些子块使用稠密矩阵的LU分解方法。但是，在像电路矩阵这样的极稀疏矩阵中，通常无法形成稠密子块。因此，这种方法不适合用于电路仿真。

另一种是G/P left-looking方法，该方法在“J.R.Gilbert and T.Peierls,“Sparse partialpivoting in time proportional to arithmetic operations,”SIAM J.Sci.Statist.Comput.,vol.9,no.5,pp.862–874,1988”一文中提出。该方法较强的数据依赖性使得它的并行化十分困难。据了解，G/P left-looking方法尚未在GPU上实现。同时，数据间的高依赖性还造成并行G/Pleft-looking方法只能高效运行在共享内存的计算设备上，如现场可编程门阵列(FieldProgrammable Gate Array,FPGA)、多核CPU和图形处理器(GPU)。FPGA上实现的稀疏LU分解对大矩阵的可拓展性受限于FPGA很有限的片上资源。由于共享内存的CPU核数通常很有限（大多数商用CPU都不超过6核，如Intel Xeon X5680,AMD Phenom II），稀疏矩阵LU分解在多核CPU上性能无法进一步提升。

发明内容

本发明旨在至少解决上述技术问题之一。

为此，本发明的目的在于提出一种基于GPU的稀疏矩阵LU分解方法，该方法可根据分解中浮点运算次数的多少自动选择计算平台，能够合理分配CPU和GPU上的任务，保证GPU上众多线程间的协同运行，从而提高运算速度，减少电路仿真时间。

为了实现上述目的，本发明的实施例提出了一种基于GPU的稀疏矩阵LU分解方法，包括以下步骤：A：在CPU上建立待仿真的电路网单的电路矩阵，并对所述电路矩阵进行预处理和稀疏化；B：根据预处理结果选择对所述稀疏矩阵LU进行处理的计算平台；C：如果对所述稀疏矩阵LU进行处理所使用的计算平台为GPU，则通过所述GPU对所述稀疏矩阵LU进行分解，包括：C1：将所述电路矩阵的非零元位置、值以及所述稀疏矩阵的LU因子的非零元位置输入所述GPU；C2：确定群并行模式和流水线并行模式的分界点；C3：依次计算属于所述群并行模式的各层直至所有层均计算完成；C4：通过所述流水线并行方式计算剩余各层；C5：将所述稀疏矩阵LU中非零元的值写回到所述CPU；D：通过所述CPU对所述稀疏矩阵LU中非零元的值进行回代求解，以完成对稀疏矩阵LU的分解。

根据本发明实施例的基于GPU的稀疏矩阵LU分解方法，可根据分解中浮点运算次数的多少自动选择计算平台，能够合理分配CPU和GPU上的任务，保证GPU上众多线程间的协同运行，从而提高运算速度，减少电路仿真时间；另外，对于分解中浮点运算次数较多的矩阵，GPU稀疏矩阵LU分解相对单核CPU、4核CPU和8核CPU分别可实现7.9倍、2.6倍和1.5倍的加速比。

另外，根据本发明上述实施例的基于GPU的稀疏矩阵LU分解方法还可以具有如下附加的技术特征：

在本发明的实施例中所述对所述电路矩阵进行预处理进一步包括：对数值进行处理以提高所述数值的稳定性；根据近似最小自由度算法减少分解过程中的非零元写入；根据部分选主元的数值分解计算所述矩阵LU的非零元结构。

在本发明的实施例中，根据HSL_MC64方法对所述数值进行处理以提高所述数值的稳定性。

在本发明的实施例中，所述根据预处理结果选择对所述稀疏矩阵LU进行处理的计算平台，进一步包括：对所述矩阵LU中每行的非零元按列号进行排序；建立消去图，并对所述消去图进行分层；根据浮点运算次数确定要使用的平台；如果所述浮点运算次数大于或等于预设次数，则使用GPU平台；如果所述浮点运算次数小于所述预设次数，则使用CPU平台。

在本发明的实施例中，所述建立消去图，并对所述消去图进行分层，进一步包括：建立消去图，其中，所述消去图中包括节点，且各节点表示各行。对所述消去图进行分层，其中，所述节点的层号为所述节点所依赖的行所处的最大层号加1。

在本发明的实施例中，所述对所述消去图进行分层之后，还包括：依次处理所述消去图的每一层，所述每一层内部的各行可按照任意顺序进行计算。其中，处理所述消去图的方式包括：群并行模式和流水线并行模式，所述群并行模式指所述每一层有足够多的行被并行分解；所述流水线并行模式指所述每一层只有很少的行，且来自不同层的行重叠执行分解。

在本发明的实施例中，所述通过所述流水线并行方式计算剩余各层，进一步包括：如果检测到未计算完成的行，则设为第k行，并计算所述第k行；依次检测所述第k行所依赖的行，并设为第j行，并判断所述第j行是否计算完成；如果所述第j行计算完成，则从第k行中减去第j行乘以其相应系数；如果第j行未完成，则等待预设时间后重新检测第j行是否计算完成；依次确认所述第k行所依赖的每一行均计算完成后，从所述第k行中减去该行乘以其相应系数；将第k行的计算结果写回到所述GPU中，并判定第k行已计算完成，并进一步处理下一个未完成的行，直至所有行均已完成计算。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1为根据本发明一个实施例的基于GPU的稀疏矩阵LU分解方法的流程图；

图2为根据本发明一个实施例的基于GPU的稀疏矩阵LU分解方法的一个消去图示意图；

图3为根据本发明另一个实施例的基于GPU的稀疏矩阵LU分解方法的流程图；

图4为根据本发明一个实施例的基于GPU的稀疏矩阵LU分解方法的对非零元访问情况示意图；

图5为根据本发明一个实施例的基于GPU的稀疏矩阵LU分解方法的不活跃warp引起死锁的示意图；

图6为根据本发明一个实施例的基于GPU的稀疏矩阵LU分解方法的基于GTX580GPU的稀疏矩阵LU分解实现的带宽与活跃warp数之间的关系示意图；和

图7为根据本发明一个实施例的基于GPU的稀疏矩阵LU分解方法的各矩阵的浮点运算次数及到达的带宽示意图。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。

在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性。

在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。

本发明描述中的“warp”指的是GPU上的线程组，一个线程组中包含的线程数量可以通过GPU的产品手册查到。在本发明的一个实施例中，使用的NVIDIA GTX580GPU上一个warp包含32个线程。在本发明的描述中，“warp”与“线程组”是等价的术语。

以下结合附图详细描述根据本发明实施例的基于GPU的稀疏矩阵LU分解方法。

图1为根据本发明一个实施例的基于GPU的稀疏矩阵LU分解方法的流程图。如图1所示，根据本发明一个实施例的基于GPU的稀疏矩阵LU分解方法，包括以下步骤：

步骤S101，在CPU上建立待仿真的电路网单的电路矩阵，并对电路矩阵进行预处理和稀疏化。

具体地，对电路矩阵进行预处理时，首先对数值进行处理以提高数值的稳定性，在本发明一个实施例中，通过HSL_MC64方法对数值进行处理以提高数值的稳定性，然后根据近似最小自由度算法减少分解过程中的非零元写入，最后根据部分选主元的数值分解计算矩阵LU的非零元结构。

步骤S102，根据预处理结果选择对稀疏矩阵LU进行处理的计算平台。

具体地，对稀疏矩阵LU中每行的非零元按列号进行排序，建立消去图，并对消去图进行分层，然后根据浮点运算次数确定要使用的平台，具体而言，如果浮点运算次数大于或等于预设次数，则使用GPU平台，如果浮点运算次数小于预设次数，则使用CPU平台。其中，预设次数根据具体情形预先设定。

另外，在上述步骤S102中，建立消去图，并对消去图进行分层进一步包括：首先建立消去图，其中，该消去图中包括节点，且各节点表示各行，其次对消去图进行分层，且各节点的层号为该节点所依赖的行所处的最大层号加1。

进一步地，在对消去图进行分层之后，依次处理消去图中的每一层，且每一层内部的各行可按照任意顺序进行计算，其中，处理消去图的方式包括：群并行模式和流水线并行模式。群并行模式指每一层有足够多的行被并行分解；流水线并行模式指每一层只有很少的行，且来自不同层的行重叠执行分解。

作为具体的示例，在上述步骤S102中，换言之，即首先对矩阵L和U中每行的非零元按照列号排序，以提高数据局域性。在GPU程序中，如果连续的线程访问连续的内存空间，可以实现较高的带宽。但预处理后，LU因子中的非零元是乱序的，这影响了访存的连续性。因此对L和U中每行的非零元按列号排序，以提高主存中的数据的局域性。如图4所示，排序后，一行中相邻的非零元更有可能被相邻的线程处理。从而，无需对GPU程序做任何修改，就获得了1.7倍的加速。排序过程只需进行一次，并将其融合到预处理中。实验证明，对非零元进行快速排序的代价很小。

其次，建立消去图，并对消去图分层。消去图中各节点表示行，当且仅lik≠0时，即第i行依赖于第k行，也即需要从xi中减去uk，生成从节点k指向节点i的边。如图2所示，为一个消去图的示意图。其中，任何符合消去图的更新顺序都是可行。具体来说，只要保证对任意的两行xi和uk（xi依赖于uk），从xi中减去uk时，uk已经计算完成，这个顺序就是可行的。如图2所示，是一个消去图的例子，消去图可表示left-looking方法中的并行性和要求的时序关系。假设图中虚线圆圈表示的行已经分解完成，而行8,9,10正在被处理。实线箭头代表的操作是目前可以执行的，这些操作中是存在并行性的。同时要注意，对某一未完成行的更新操作需要按照严格的顺序执行。例如，在图2所示的时刻，行9可以依次减去行4,6,7，对应图2中的实线箭头。而且，行9减去行4,6,7后，还必须等待行8完成，才能减去行8。对行10也存在类似的情况。进一步地，对消去图分层，某个节点的层号，是它所依赖的行所处的最大层号加1。依次处理每一层，每层内部的各行可以按照任意顺序计算。这种顺序是可行的，因为某一行所依赖的行必定都处于它之前的层，都会先于该行被计算。消去图中，属于同一层的行可以被独立地并行地处理。然而，除了最前面的几层以外，后面每层往往只含有一两个节点，不具备足够的并行性。此时，需要将不同层的计算重叠起来。为此，将每层有足够多的行被并行的分解称为群并行模式，将每层只有很少的行，来自不同层的行重叠执行称作流水线并行模式。

最后，根据分解过程中涉及的浮点运算次数确定使用的平台。具体地，如果浮点运算次数大于或等于预设次数F0，则使用GPU平台；反之，使用CPU平台。其中，预设次数F0的具体值与所使用的CPU和GPU型号有关。在本发明一个优选实施例中，使用的是IntelE5680CPU和NVIDIA GTX580GPU，相应的F0=200M。

步骤S103，如果对稀疏矩阵LU进行处理所使用的计算平台为CPU，则采用中国发明专利“201110337789.6，陈晓明，汪玉，杨华中，针对电路仿真的自适应并行LU分解方法”中的方法进行稀疏LU分解；如果对稀疏矩阵LU进行处理所使用的计算平台为GPU，则通过GPU对稀疏矩阵LU进行分解，具体包括以下步骤:

步骤S1031，将电路矩阵的非零元位置、值以及稀疏矩阵的LU因子的非零元位置输入GPU。

步骤S1032，确定群并行模式和流水线并行模式的分界点。具体地，从消去树的第一层开始逐层检查每层上的点的个数，直到某一层上的点数小于预设的warp（线程组）数量时，这一层之前的所有层采用群并行模式，这一层以及这一层之后的所有层采用流水线模式，所述分界点确定。

步骤S1033，依次计算属于群并行模式的各层直至所有层均计算完成。具体地，依次计算属于群并行模式的各层，GPU上的每个线程组负责计算一行，一个工作组内部的线程操作同一行的不同非零元。具体过程即，从正在被计算的行（设为第k行）中依次减去它所依赖的行（设为第j行）乘以其相应系数（即lkj）。通过维护一个标记数组，记录每行是否已完成计算，初始化为未完成，某行计算结束后，标记该行为已完成，最后重复上述过程直到所有属于群并行模式的层计算完成。

步骤S1034，通过流水线并行方式计算剩余各层。具体地，如果检测到未计算完成的行，则将该行设为第k行，并计算第k行，依次检测第k行所依赖的行，并设为第j行，并判断第j行是否计算完成，如果第j行计算完成，则从第k行中减去第j行乘以其相应系数（ljk），如果第j行未计算完成，则等待预设时间后重新检测第j行是否计算完成（此时一定有某个线程组正在并行地计算第j行，第j行最终一点会被完成，因而此处不会形成死锁），重复上述过程，即依次确认第k行所依赖的每一行均计算完成后，则从第k行中减去该行乘以其相应系数，最后将第k行的计算结果写回到GPU中，并判定第k行已计算完成，并继续处理下一个未完成的行，直至所有行均已完成计算，即矩阵的LU分解已完成。

步骤S1035，将稀疏矩阵LU中非零元的值写回到CPU。即将GPU中的分解结果写回到CPU。

步骤S104，通过CPU对稀疏矩阵LU中非零元的值进行回代求解，已完成对稀疏矩阵LU的分解。

在上述的示例中，涉及到了GPU众多线程之间的协同，以及保证GPU不同线程组之间的时序关系，这是基于GPU的稀疏矩阵LU分解方法中的一个难点。在GPU上保证时序关系必须精确控制发布的线程组的个数。普通GPU程序中允许一些warp在开始时不活跃，而等其他warp完成后再活跃。然而，在并行left-looking算法的流水线并行模式中，我们必须确保发布的所有warp从一开始就活跃，否则将会造成死锁。如图5所示，说明了这一情况。假设在一个只支持两个活跃warp的GPU上发布了三个warp，群并行模式可以顺利完成，因为warp1和2最终总会结束运行，于是warp3开始运行。但在流水线并行模式中，行9和行10依赖于行8，而行8被分配给了不活跃的warp3，于是活跃的warp（1和2）一直等待，陷入死循环。反过来又使warp3没有机会活跃，从而形成死锁。因此流水线并行模式中能被并行分解的最大行数就是GPU程序执行时活跃warp的数目。

对许多GPU程序而言，更多的活跃warp数目通常意味着更高的性能，但如果发布了太多的活跃warp，同时被分解的行数过多，可能会出现如下情况：矩阵下方的行依赖于一些之前的行，由于这些行没有完成，下方的这些行只能等待。在这里等待使用死循环实现的，这些行的分解无法进行，却仍然会占用运算单元。这反而会造成性能的下降。这个转折点不会在CPU上出现，但在GTX580GPU上可以观察到这个现象。如图6所示，在图6中，给出了在GTX580总分解其中四个矩阵，在发布不同warp数时达到的带宽。最高性能并不是在活跃warp数最多时达到的，而是在大约每个流处理器有24个warp时达到的。这意味着在GTX580中已经充分利用了稀疏矩阵LU分解中流水线并行模式中的并行性。采用一个自动调整策略来确定GPU上最佳的warp数，即针对某个输入矩阵，使用不同的warp数目，多次执行上述GPU稀疏矩阵LU分解方法，找到使性能最佳的warp数，用于后续迭代之中。

在GTX580上，最高可实现74%的峰值带宽。考虑到访存并不完全连续，而且warp有时候需要等待，所以这个带宽已经接近峰值的192GB/s。如果今后GPU峰值带宽能够进一步提升，则稀疏矩阵LU分解的性能也将进一步提升。

图3为根据本发明另一个实施例的基于GPU的稀疏矩阵LU分解方法的流程图。

如图3所示，根据本发明另一个实施例的基于GPU的稀疏矩阵LU分解方法，包括以下步骤：

步骤S301，在CPU上进行预处理。其中，预处理包括：1）使用HSL_MC64方法提高数值稳定性；2）使用近似最小自由度算法减少分解过程中的非零元填入；3）使用部分选主元的数值分解计算L和U的非零元结构；4）对矩阵L和U中每行的非零元按列号排序以提高数据局域性。预处理在CPU上进行，后续的迭代的数值分解过程在GPU上进行。在电路仿真的迭代过程中，L和U的非零元结构保持不变。将预处理后的矩阵记为矩阵A。

步骤S302，对矩阵L和U中的非零元进行排序。具体地，对矩阵L和U中每行的非零元按列号排序，可提高数据局域性。

进一步地，通过对矩阵进行图的深度优先遍历，得到消去图，并将消去图分层。并根据分解过程中涉及的浮点次数确定是否使用GPU平台。具体地，如果浮点运算次数少于F0，则使用CPU平台；如果浮点运算次数多于F0，则使用GPU进行稀疏矩阵LU分解，并进一步执行步骤S303。

步骤S303，将矩阵A的L和U中的非零元位置写入GPU内存。

步骤S304，将矩阵A的非零元的值写入GPU内存。

步骤S305，判断是否处于群并行模式，如果是则执行步骤S306，否则执行步骤S307。

步骤S306，在GPU上并行地分解该层中的各行。完成后返回继续执行步骤S305。具体而言，在群并行模式中，每层有足够多的行，可以独立地并行地被分解。依次并行地处理每层中的各行，各行被分配给GPU的各个工作组，一个工作组内部的线程操作同一行的不同非零元。维护一个标记数组，记录每行是否已完成计算。初始化为未完成，处理完某行后，标记该行已完成。

步骤S307，在GPU上按流水线并行模式分解剩余所有行。具体地，在流水线并行模式中，每层的行很少，其并行模式与群并行模式有重要区别。在流水线并行模式中，需要保证时序关系。每次更新前，必须检查被依赖的行的计算是否已经完成。如果尚未完成，则必须等待一段时间之后重新检查，用这种方法可保证所要求的时序关系。

步骤S308，从GPU中读取矩阵LU中非零元的值，进一步地，将该值写回到CPU，并回代求解。

需要说明的是，在上述图3所示的流程图的步骤中，步骤S301至步骤S303，只需执行一次，而步骤S304至步骤S308需要多次循环执行。

作为具体的示例，下表1显示了根据本发明一个实施例的基于GPU的稀疏矩阵LU分解的性能，以及与CPU版本实现的加速效果。表1中所列测试矩阵来自佛罗里达大学矩阵收集（http://www.cise.ufl.edu/research/sparse/matrices/）。

表1.GPULU分解性能及加速比¹

1除KLU之外，其他CPU结果均指我们自己的CPU实现

2矩阵大小

3矩阵A中的非零元

4分解中的浮点运算次数（兆为单位）

5所有的平均值均为几何平均

表1

图7为根据本发明一个实施例的基于GPU的稀疏矩阵LU分解方法的各矩阵的浮点运算次数及到达的带宽示意图。

如图7所示，给出了各个矩阵分解中的浮点运算次数以及分解该矩阵时在GPU上实现的带宽。基于GPU的稀疏矩阵LU分解，对于分解过程中浮点运算次数较多（本发明的实施例中大于200M，即B组）的矩阵非常有效。在这些矩阵上，基于GPU的LU分解相对单核CPU、4核CPU、8核CPU分解具有7.9倍、2.6倍和1.5倍的加速。由于大多数商用CPU都不超过6核，因此对于分解中运算量较大的矩阵使用GPU进行LU分解，能带来可观的性能提升。且应用在电路仿真领域时，还能提高SPICE电路仿真器的性能。

需要说明的是，在分解C组矩阵时，性能大大超过CPU。这是因为，分解C组矩阵的过程中会出现许多反常浮点数（denormal floating point numbers）。其中，反常浮点数用来表示极小的实数。而CPU处理这种数据很慢，因此分解这些矩阵时性能很差。

另外，GPU上的稀疏矩阵LU分解本身也具有重要意义。目前，GPU上的基本线性代数程序库(Basic Linear Algebra Subprogram，BLAS)中已有三角矩阵求解的模块，但尚没有稀疏矩阵LU分解模块。而LU分解的作用可将一个一般矩阵方程转换为两个三角矩阵方程。

根据本发明实施例的基于GPU的稀疏矩阵LU分解方法，可根据分解中浮点运算次数的多少自动选择计算平台，能够合理分配CPU和GPU上的任务，保证GPU上众多线程间的协同运行，从而提高运算速度，减少电路仿真时间；另外，对于分解中浮点运算次数较多的矩阵，GPU稀疏矩阵LU分解相对单核CPU、4核CPU和8核CPU分别可实现7.9倍、2.6倍和1.5倍的加速比。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同限定。

Claims (7)

1.一种基于GPU的稀疏矩阵LU分解方法，其特征在于，包括以下步骤：

A：在CPU上建立待仿真的电路网单的电路矩阵，并对所述电路矩阵进行预处理和稀疏化；

B：根据预处理结果选择对所述稀疏矩阵LU进行处理的计算平台；

C：如果对所述稀疏矩阵LU进行处理所使用的计算平台为GPU，则通过所述GPU对所述稀疏矩阵LU进行分解，包括：

C1：将所述电路矩阵的非零元位置、值以及所述稀疏矩阵的LU因子的非零元位置输入所述GPU；

C2：确定群并行模式和流水线并行模式的分界点；

C3：依次计算属于所述群并行模式的各层直至所有层均计算完成；

C4：通过所述流水线并行方式计算剩余各层；

C5：将所述稀疏矩阵LU中非零元的值写回到所述CPU；

D：通过所述CPU对所述稀疏矩阵LU中非零元的值进行回代求解，以完成对稀疏矩阵LU的分解。

2.如权利要求1所述的基于GPU的稀疏矩阵LU分解方法，其特征在于，所述对所述电路矩阵进行预处理，进一步包括：

对数值进行处理以提高所述数值的稳定性；

根据近似最小自由度算法减少分解过程中的非零元写入；

根据部分选主元的数值分解计算所述矩阵LU的非零元结构。

3.如权利要求2所述的基于GPU的稀疏矩阵LU分解方法，其特征在于，根据HSL_MC64方法对所述数值进行处理以提高所述数值的稳定性。

4.如权利要求1所述的基于GPU的稀疏矩阵LU分解方法，其特征在于，所述根据预处理结果选择对所述稀疏矩阵LU进行处理的计算平台，进一步包括：

对所述矩阵LU中每行的非零元按列号进行排序；

建立消去图，并对所述消去图进行分层；

根据浮点运算次数确定要使用的平台；

如果所述浮点运算次数大于或等于预设次数，则使用GPU平台；

如果所述浮点运算次数小于所述预设次数，则使用CPU平台。

5.如权利要求4所述的基于GPU的稀疏矩阵LU分解方法，其特征在于，所述建立消去图，并对所述消去图进行分层，进一步包括：

建立消去图，其中，所述消去图中包括节点，且各节点表示各行。

对所述消去图进行分层，其中，所述节点的层号为所述节点所依赖的行所处的最大层号加1。

6.如权利要求5所述的基于GPU的稀疏矩阵LU分解方法，其特征在于，所述对所述消去图进行分层之后，还包括：

依次处理所述消去图的每一层，所述每一层内部的各行可按照任意顺序进行计算。其中，

处理所述消去图的方式包括：群并行模式和流水线并行模式，所述群并行模式指所述每一层有足够多的行被并行分解；所述流水线并行模式指所述每一层只有很少的行，且来自不同层的行重叠执行分解。

7.如权利要求1所述的基于GPU的稀疏矩阵LU分解方法，其特征在于，所述通过所述流水线并行方式计算剩余各层，进一步包括：

如果检测到未计算完成的行，则设为第k行，并计算所述第k行；

依次检测所述第k行所依赖的行，并设为第j行，并判断所述第j行是否计算完成；

如果所述第j行计算完成，则从第K行中减去第j行乘以其相应系数；

如果第j行未完成，则等待预设时间后重新检测第j行是否计算完成；

依次确认所述第k行所依赖的每一行均计算完成后，从所述第k行中减去该行乘以其相应系数；

将第k行的计算结果写回到所述GPU中，并判定第k行已计算完成，并进一步处理下一个未完成的行，直至所有行均已完成计算。

Images