电路仿真时电路稀疏矩阵的基于消去图的并行分解方法

摘要

电路仿真时电路稀疏矩阵的基于消去图的并行分解方法，属于电子设计自动化(EDA)领域，其特征在于，从电路矩阵的符号分析结果中提取出消去图，用来表示LU分解过程中的数据依赖性，根据消去图的不同结构使用两种不同的并行方法：基于群的并行和流水线的并行，从而降低LU分解的运算时间、加快电路仿真的速度，在一系列测试电路矩阵上的测试结果表明本发明的LU分解时间在并行线程数为1～8个时比LU分解软件KLU快1.66～7.72倍。

说明书

技术领域

本发明涉及一种针对电路仿真时电路稀疏矩阵的基于消去图的并行LU分解方法，属于电子设计自动化(EDA)技术领域。

背景技术

在科学计算中，求解线性方程组Ax＝b(Ax＝b的形式如图1所示，A为n×n矩阵，b和x都是n维列向量，A和b已知，x未知待求)具有举足轻重的地位。Ax＝b的求解方法主要可分为两大类：迭代法和直接法。迭代法如雅可比迭代法、高斯——赛德尔迭代法、超松驰迭代法等，都是先设定一组解的初值，然后通过一个迭代公式进行迭代，使得解逐渐向真实的解靠近，当误差小于给定值时，迭代收敛。但是迭代法通常对矩阵A的要求较高，在一般的计算问题中可能难以满足。不同于迭代法，直接法通过固定的求解公式直接对A进行有限步骤的计算，从而求得所需要的解。相比于迭代法，直接法的适用性更强，稳定性更好。因此直接法在各类科学计算中用得更为普遍。传统的直接法如高斯消元法、LU分解法(也称为三角分解法)、乔里斯基分解法等。

在一般的问题中，矩阵A通常是非常稀疏的。比如对于电路仿真问题，通常电路矩阵中每行的非零元个数仅为5个左右，而且与矩阵规模无关。针对稀疏矩阵特别设计的求解算法，可以大大降低运算复杂度而减少运算时间。

在各类直接解法中，运用最多的是LU分解法。LU分解法指的是将一个n×n的方矩阵A分解成一个下三角矩阵L和一个上三角矩阵U的乘积(其中U矩阵的对角线元素都为1)，即A＝LU，其中L和U也都是n×n的矩阵。从而求解线性方程组Ax＝b的问题转化成求解两个三角方程Ly＝b和Ux＝y(y是n维列向量)。图2显示了LU分解的基本形式，其中L和U未写数值的地方都是0。

稀疏矩阵的LU分解过程包括预处理和数值分解两部分。其中预处理是指使用一定的算法对矩阵进行行列交换以达到在数值分解过程中减少运算量的目的，有些预处理方法中还引入了选主元(选主元指的是将绝对值大的元素通过矩阵行列交换操作，交换到对角线的位置上)的步骤，以保证分解过程中的数值稳定性。不管采用哪种预处理方法，预处理之后的矩阵只是对原矩阵进行了一些行、列的交换，而没有其他的变化。本发明主要针对LU分解算法的后者，即数值分解部分，如无特别说明，后文提及的LU分解皆指数值分解部分。本发明并不局限于某种特定的预处理方法，对于所有预处理方法都能适用。

目前针对稀疏矩阵的LU分解算法大致可分为向左看算法(Left-looking Algorithm)和向右看算法(Right-looking Algorithm)两大类，其中向左看算法由于对稀疏矩阵存储结构有着很好的适应性，在LU分解软件中得到的广泛的应用。目前加州大学伯克利分校(University of California，Berkeley)Sherry Li开发的SuperLU和弗罗里达大学(University of Florida)Tim Davis开发的KLU等软件的LU分解部分都是以向左看算法为基础，针对大规模稀疏矩阵特点进行了优化。

在向左看LU分解算法中，从第1列到第n列依次分解A的每一列，即每次计算一个列向量。对每一列的分解可以概括为三个步骤：符号分析、数值分解和数值分配，如图3所示。以正在分解A的第k列为例说明。符号分析指从矩阵A的第k列的非零元结构(非零元结构指的是一个集合，这个集合包括第k列上所有非零元素所在的行的行号)计算出LU分解完成后第k列上的非零元结构。具体的符号分析所采用的方法，参见文献J.R.Gilbert and E.Ng，Predicting structure in nonsymmetric sparse matrix factorizations，Graph Theory and Sparse Matrix Computation，Springer Verlag，1993(该方法可称为“Gilbert符号分析方法”)。数值分解步骤是根据这一列的符号分析结果进行数值计算，从而获得第k列的所有非零位置上的数值结果。数值分析所采用的算法，参见文献J.R.Gilbert and T.Peierls，Sparse partial pivoting in time proportional to arithmetic operations，SIAM J.Sci.Statist.Comput.，vol.9，pp.862-874，1988(该算法可称为Gilbert/Peierls算法，简称G/P算法)。最后一个步骤是数值分配，即将第k列的数值计算结果(是一个列向量)中行号小于等于k的部分分配给U矩阵，行号大于等于k的部分分配给L矩阵。依次对k＝1，2，…，n循环执行上述步骤，即完成了对整个矩阵A的LU分解。

虽然目前已经有很多基于LU分解的软件包，比如SuperLU，KLU，PARDISO，UMFPack等共几十个软件，但是在这些软件中，并行的版本非常少，而针对通用多核CPU平台的并行软件更是寥寥无几。由于LU分解过程中的高度数据依赖性，以及并行时多线程操作所带来的额外代价，造成了并行的LU分解软件数据如此之少。

发明内容

本发明的目的是提供一种针对电路仿真中电路矩阵的并行LU分解方法。本发明提出使用消去图来表示LU分解过程中的数据依赖性并用来进行并行的任务调度。该算法从对电路矩阵的符号分析结果中提取出消去图，然后基于该消去图的不同结构，提出两种不同的并行方式：群并行与流水线并行。在LU数值分解过程中这两种并行方式动态的调用，以适应不同的数据依赖性。

电路仿真时电路稀疏矩阵的基于消去图的并行分解方法，其特征在于，是在计算机中按照以下步骤实现的：

步骤(1)，输入要仿真解析的电路的网单；

步骤(2)，建立n×n的电路稀疏矩阵，包括add20、circuit_1、circuit_2、add32、rajat03、coupled、circuit_3、onetone1、onetone2、ckt11752_dc_1、circuit_4、ASIC_100ks、ASIC_100k、dc1、trans4、G2_circuit、transient、ASIC_320ks、ASIC_320k、rajat30、ASIC_680ks、ASIC_680k、G3_circuit、Freescale1、rajat31以及circuit5M在内的任何一个电路稀疏矩阵；

步骤(3)，选择对角块模式和非对角块模式中的任何一种模式，对步骤(2)中所建立的电路稀疏矩阵进行预处理，得到预处理之后的电路稀疏矩阵A；

步骤(4)，根据John R.Gilbert和Timothy Peierls提出的非零元符号分析方法，对步骤(3)得到的所述n×n的电路稀疏矩阵A进行符号分析，从第1列到第n列依次完成对所述电路稀疏矩阵A的下三角矩阵L和上三角矩阵U的n列内的非零结构的计算，A＝LU；

步骤(5)，利用步骤(4)得到的U矩阵的非零结构，计算消去图：消去图共有n个顶点，每个分立的顶点i唯一地对应于所述U矩阵中的第i列，顶点之间不存在任何连线，但每个顶点都有一个“所处的层”属性值，用于表示每个顶点以及该顶点在所述U矩阵中对应的列在消去图中处于哪一层，该属性值记为level(i)，建立所述消去图的步骤如下：

步骤(5.1)，从第1列到第n列循环所有列，对其中的每一列k，取出所述U矩阵中第k列中不包括对角线的上三角部分的所有非零元，共z个，所述z个非零元所在行的行号分别标记为j₁，j₂，…，j_z；

步骤(5.2)，若：z＝0，则level(k)＝1；

步骤(5.3)，若：z＞O，则取出所述各行号对应的相同序号的列所处的层的属性值level(j₁)，level(j₂)，…，level(j_z)，求出其中的最大值m＝max(level(j₁)，level(j₂)，…，level(j_z))，于是，所述第k列在所述消去图中所处的层的属性值为level(k)＝m+1；

步骤(5.4)，得到所述n列中每一列所处的层的属性值之后，就得到了整个消去图，从该消去图中所有顶点所处的层的属性值中求出最大值，即为整个消去图的总层数，记为L_total；

步骤(6)，根据步骤(4)得到的L矩阵和U矩阵的非零结构以及步骤(5)得到的消去图，按以下步骤进行并行的LU数值分解：

步骤(6.1)，输入线程数目Thread，并创建每个线程，所述线程数目Thread的值应大于或等于1，小于或等于所使用的计算机上的CPU的核心数；

步骤(6.2)，输入一个用于考察所述消去图中各层顶点总数V(lv)的阈值V_th，该阈值V_th在所述线程数目Thread的1倍-10倍之间取值；

步骤(6.3)，从第1层到第L_total层依次考察所述消去图中每一层lv的顶点总数V(lv)：

若：V(lv)＜V_th，则对第lv层采用流水线并行方法，

若：V(lv)≥V_th，则对第lv层采用群并行方法，

把当前任何一种所述并行方法的开始层记为L_start，结束层记为L_stop，L_start＝lv，按以下方法从开始层L_start开始向后考察各层，第lv+1层，…，直到第L_toal层为止，寻找所述的结束层L_stop：

在采用流水线并行方法时，若有一层c满足V(c)≥V_th，则结束层L_stop＝c-1，

在采用群并行的方法时，若有一层c’满足V(c’)＜V_th，则结束层L_stop＝c’-1，

步骤(6.4)，把从步骤(6.3)中得到的从所述消去图中的第L_start到第L_stop层这些层中的所有顶点号，根据步骤(6.3)中得到的相应的一种并行方法，按以下步骤进行多线程方式的并行LU数值分解：

步骤(6.4.1)，当从开始层L_start到结束层L_stop采用群并行方法时，在循环进行到其中任何一层时，把该层上的所有顶点平均分配给所述各个线程，所述各个线程都采用John R.Gilbert和Timothy Peierls提出的LU数值分解方法对分配到的各列进行并行的LU数值分解，等待所有线程都计算完成后再进入下一层，循环往复，遍历从开始层L_start到结束层L_stop的所有层，

步骤(6.4.2)，当从开始层L_start到结束层L_stop采用流水线并行方法时，把从开始层L_start到结束层L_stop的所有层中的列号按层的顺序排成一个序列，并给这个序列中每一个列都置一个标志，初始时所有标志都是“未完成”，然后每个线程从该序列的头部开始依次取出一个列号，并从序列中删去所取到的列号，用John R.Gilbert和Timothy Peierls提出的LU数值分解方法对所取到的列号进行LU数值分解，在线程完成一列的LU数值分解后，把该列对应的标志写成“已完成”，每个线程都重复这一过程直到所述序列为空为止；

步骤(7)，在满足收敛条件并完成瞬态仿真后，打印仿真结果。

在所述步骤(6.4.2)中，在进行某一列x的LU数值分解时，要判断所述U矩阵的第y行、第x列上的值U(y，x)是否为0：

若，为0，则对第x列继续进行LU数值分解，

若：非0，则第x列进行数值分解时依赖于第y列的数值分解结果，若：第y列的标志为“未完成”，则分解第x列的线程需要等待，直到所述第y列的标志变成“已完成”后再继续进行下一步计算。

利用本发明提出的并行LU分解算法，可以在通用多核CPU平台上进行并行的LU分解，从而加速电路仿真的速度。在一系列测试电路矩阵上的测试结果表明本发明的LU分解时间在并行线程数为1～8个时比LU分解软件KLU快1.66～7.72倍。本发明所提出的并行LU分解算法与矩阵中的数据依赖性高度协调，从而可以保证在计算过程中的高度并行化。

附图说明

图1为Ax＝b的基本形式。

图2为LU分解的基本形式。

图3为G/P算法的总体流程图。

图4为非零元符号分析中的有向无环图举例：图4.a为下三角矩阵L的一个例子，其中第9列已计算完成，图4.b为L的前9列对应的有向无环图。

图5为上三角矩阵U与其对应的消去图举例：图5.a为上三角矩阵U的一个例子，图5.b为与图5.a对应的消去图。

图6为针对电路仿真中电路矩阵的并行LU分解方法的总体流程。

图7为并行LU分解方法的一个例子，对应于图5：表示线程1，表示线程2。

具体实施方式

为了实现上述目的，本发明采用如下技术方案，其实施步骤为：

步骤1：对输入的n×n电路矩阵A进行符号分析，从第1列到第n列依次完成对L和U的n列的非零元结构的计算；

步骤2：利用步骤1得到的U矩阵的非零结构，通过最早开始算法，获得消去图；

步骤3：基于步骤1得到的L和U的非零结构以及步骤2得到的消去图，进行并行的LU数值分解。在数值分解过程中，采用两种不同的并行方法：群并行和流水线并行。这两种并行方法在数值分解过程中动态地根据消去图的结构进行选择。从而完成并行的LU分解过程。

本发明提出的针对电路仿真中电路矩阵的并行LU分解算法，结合附图详细说明如下。本文中所有的n均是指矩阵的维度，即矩阵A是n×n阶方阵，L和U分别是n维下三角矩阵和n维上三角矩阵，见图2。

步骤1：对输入的n×n电路矩阵A进行符号分析，从第1列到第n列依次完成对L和U的n列的非零元结构的计算；

非零元符号分析采用了J.R.Gilbert and E.Ng，Predicting structure in nonsymmetric sparse matrix factorizations，Graph Theory and Sparse Matrix Computation，Springer-Verlag，1993中符号分析的方法(该方法可称为“Gilbert符号分析方法”)。对矩阵的第k列进行非零元的符号分析通过以下1.1～1.4完成：

1.1：取A的第k列中的所有非零元所在行的行号为集合B；

1.2：将对第1～k 1列的符号分析结果中的L矩阵看做一张有向无环图，矩阵L的每一列对应图中一个节点，即第j列对应图中节点j。如果第j列上的第i行为非零元，即L_ij≠0，则图中存在一条从j指向i的边。如图4所示的例子，该例子是已完成第1～9列符号分析的L以及L所对应的有向无环图；

1.3：从集合B中的每一个元素为起点在L矩阵对应的有向无环图中进行深度优先搜索。深度优先搜索(Depth-First Search)的过程可表述为：假设顶点v为深度优先遍历的起点(源点)，那么首先访问出发点v，并将其标记为“已访问过”；然后依次从v出发搜索v的每个邻接点w。若w没有被访问过，则以w为新的出发点继续进行深度优先搜索，否则顶点w已被访问过，则跳过它继续寻找v的所有邻接点中没有被访问过的顶点。因此深度优先遍历是一个递归的过程，直到图中所有的点都被访问过为止。在遍历的同时将B中每一个元素为起点进行深度优先遍历所访问到的顶点号按顺序存入一个对应的序列，这样集合B有几个元素，就获得了几个遍历序列；

1.4：将所有深度优先搜索出的节点按图的拓扑结构排列，所获得的序列就是第k列的符号分析结果，这一步实际上实质上是将1.3中获得的深度优先遍历的序列按反序排列所得。假设B中有3个非零元{a1，a2，a3}，通过1.3的方法深度优先搜索得到的序列分别为R_a1，R_a2和R_a3，其中R_a3中不包含R_a1和R_a2中已搜索的节点，那么将这三个序列反向排序，即{R_a3，R_a2，R_a1}便是按图的拓扑结构的排列。

通过对第1～n列均执行1.1～1.4这4步，就可以实现对整个L和U矩阵的非零元的符号分析。

以图4中的有向无环图为例，第1～9列的非零元结构已计算完成，现需计算第10列分解后的非零元结构。假设矩阵A的第10列中第4，6，10行的值为非零，即集合B＝{4，6，10}。则从B中的第一个元素4开始获得深度优先遍历序列R₄＝{4，9，12，13}，将这些节点标记为已访问过。从B的第二个元素6开始深度优先遍历可得到R6＝{6，10}，同样，从B的第三个元素10开始深度优先遍历可获得(空集)。因此非零元的集合{R10，R6，R4}＝{6，10，4，9，12，13}就是预测出的L和U在第10列中的非零元所在的行号。

步骤2：利用步骤1得到的U矩阵的非零结构，通过最早开始算法，获得消去图；

在消去图中，共有n个顶点，每个顶点表示矩阵中对应的列，即消去图中顶点i唯一对应矩阵中的第i列，下文中凡是涉及到消去图中的“顶点”的说法，同时也是指矩阵中对应的列，两者为一一对应的关系。消去图的物理意义是矩阵中每一列可以进行LU数值分解的最早时间。在LU数值分解过程中，列与列之间存在着数据相关性，比如第x列依赖于第y列，即第x列的数值分解会用到第y列的数值分解结果，那么第x列必须等待第y列完成数值分解后才能进行数值分解。消去图中不存在边(即顶点与顶点之间不存在线连)，但每个顶点都有一个“所处的层”属性值，用于记录每个顶点在消去图中处于哪一层，也就是该顶点在矩阵中对应的列可以进行数值分解的最早时间。对于第i列在消去图中所处的层，记为level(i)。

步骤2的计算过程如下：

从第1列到第n列循环所有列，对其中每一列(记为第k列)，取出U矩阵第k列的上三角部分(不包括对角线)的所有非零元所在行的行号，共z个，它们分别为j₁，j₂，…，j_z。如果z大于0，则求出所有非零元行号对应的相同序列的列所处的层的最大值为m，即m＝max(level(j₁)，level(j₂)，…，level(j_z))，于是第k列所处的层为level(k)＝m+1。如果z等于0，那么level(k)＝1。

如图5.a的例子所示，这是一个已完成符号分解的U矩阵的例子。首先取出第1列的上三角部分(不包括对角线)的所有非零元，为0个，那么level(1)＝1。对第2列，同样也是level(2)＝1。对第3列，有一个非零元，行号为1，所以level(3)＝max(level(1))+1＝2。对第4列，有一个非零元，行号为2，所以level(4)＝max(level(2))+1＝2。对第5列，有一个非零元，行号为3，所以level(5)＝max(level(3))+1＝3。同样的方法获得level(6)＝1，level(7)＝2。又如第8列，有两个非零元，行号分别为3、5，所以level(8)＝max(level(3)，level(5))+1＝max(2，3)+1＝4，即第8列处于第4层。以此类推，不再赘述。

获得每一列所处的层之后，即得到了整个消去图，从所有顶点所处的层中求出最大值，即得到整个消去图的总层数，记为L_total。图5.a对应的消去图如图5.b所示，共有8层。

步骤3：基于步骤1得到的L和U的非零结构以及步骤2得到的消去图，进行并行的LU数值分解。在数值分解过程中，采用两种不同的并行方法：群并行和流水线并行。这两种并行方法在数值分解过程中动态地根据消去图的结构进行选择。从而完成并行的LU分解过程。

步骤3的方法见图6所示，分以下四小步进行。

3.1输入线程数目Thread，并创建Thread个线程。Thread不能超过所使用的计算机上的CPU的核心数目，也不能小于1。

3.2输入一个阈值V_th用于区分两种并行方法。V_th是一个经验常数，没有特定的范围限制。一般可取Thread的1～10倍。

3.3从消去图的第1层到第L_total层依次考察每一层，当考察到第lv层时，检查消去图中这一层上的顶点总数(记为V(lv))，如果V(lv)＜V_th，则这一层将采用流水线并行方法，否则采用群并行方法。记L_start和L_stop为当前并行方法的开始层和结束层，置L_start等于lv(记为“当前层”)。从第lv层开始向后面的层(即第lv+1层、第lv+2层……直到第L_total层)寻找L_stop，方法为，如果当前层采用流水线并行，那么在第lv层之后的层中找到第一个层c，满足V(c)≥V_th；如果当前层采用群并行，那么就以同样的方法在当前层之后的层中找满足V(c’)＜V_th的第一个层c’。那么L_stop＝c’-1。从L_start层到L_stop层将采用相同的并行方法。

3.4将从L_start到L_stop这些层中的所有顶点号，根据在3.3中确定的并行方法，进行多线程方式的并行的LU数值分解。群并行和流水线并行两种方法分别解释如下。

3.4.1如果从L_start层到L_stop层采用群并行，那么循环从L_start到L_stop的所有层，当循环进行到第h层时，将这一层上的所有列号平均分给各个线程，然后各个线程对所分配到的列进行并行的LU数值分解，LU数值分解的方法采用文献J.R.Gilbert and T.Peierls，Sparse partial pivoting in time proportional to arithmetic operations，SIAM J.Sci.Statist.Comput.，vol.9，pp.862-874，1988中的方法，即Gilbert/Peierls算法，简称G/P算法。每个线程会分配到1个或多个列号，每个线程对它们逐列进行LU数值分解。等待所有线程都计算完成后，循环进入第h+1层进行同样的任务分配和计算，这一过程直到L_stop层计算完成为止。

3.4.2如果从L_start层到L_stop层采用流水线并行，将从L_start到L_stop的所有层中的列号按层的顺序排成一个序列，并给这个序列中每一个列都置一个标志，初始时所有标志都是“未完成”。之后每个线程从序列的头部开始依次取一个列号并进行LU数值分解计算，并将所取的列号从序列中删去。这种并行方式下LU数值分解仍然采用G/P算法。每个线程完成一列的LU数值分解计算后，将这一列对应的标志写成“已完成”状态。每个线程在计算时，如果需要用到其他线程的LU数值分解计算结果(即其他列的LU数值分解结果)，并且如果所需要的列的标志是“未完成”，则这个线程需要等待，直到所需要的列的标志变成“已完成”后方可以继续运算。判断一列(记为第x列)的数值分解时是否要用到另一列(记为第y列)的数值分解结果的方法是检查U(y，x)(U矩阵的第y行、第x列上的值)是否为0，如果该值不为0，那么x列进行数值分解时依赖于y列的数值分解结果，否则不依赖。每个线程完成一列的LU数值分解后，将这一列的标志改成“已完成”，再从序列中取一个列号，并将其从序列中删去。这个过程一直重复到序列为空为止。

重复执行3.3和3.4两步骤，对从消去图的第1层到第L_toal层循环完所有层之后，就完成了对矩阵A的并行LU数值分解。

以图7的例子进行举例说明，这个例子是与图5对应的。在这个例子中假设线程数Thread＝2，阈值V_th＝2。从消去图的第一层开始考察，第一层有3个列号，那么采用群并行方法，并找到L_start＝1和L_stop＝3。于是从第1层到第3层都采用群并行方法，即L_start＝1和L_stop＝3。首先数值分解第1层，将第1层的列号平均分给2个线程，线程1分配到第1列和第2列，线程2分配到第6列，于是这两个线程各自独立地对分配到的列进行LU数值分解计算。等两个线程都完成后，同样的方法进入第2层的计算，依此类推第3层。之后考察第4层，只有一个列号，于是采用流水线并行方法，并得到L_start＝4和L_stop＝8。将从第4层到第8层中所有列号排成一个序列得{8，11，12，13，14}，线程1从中取8，并将8从序列中删去，线程2从中取出11，并将11从序列中删去，然后这两个线程各自独立地进行LU数值分解。由于线程2所分解的第11列依赖于线程1所分解的第8列的LU数值分解结果(因为图5.a中U(8，11)不为0)，所以当线程2的LU数值分解计算过程中需要用到第8列的LU数值分解结果时，线程2需要等待直到第8列的状态变成“已完成”后才可以继续分解第11列。当线程1完成第8列的LU数值分解后，将第8列标志成“已完成”，从序列中取出12，并将12从序列中删去，然后开始对第12列进行LU数值分解，同样的道理，第12列的LU数值分解可能需要等待第11列的LU数值分解结果，对后面的列以此类推。

本发明提出的针对电路仿真中电路矩阵的并行LU分解方法拥有良好的加速性能。测试结果举例如下表，所用测试矩阵均是从美国佛罗里达大学矩阵收集网站所下载获得。该表显示了本发明提出的并行LU分解方法与KLU进行分解时间比较的结果。表中所统计的时间单位均为秒，“KLU时间”这一列表示KLU所用的LU分解时间，P＝1，2，4，8这四列分别表示本发明在线程数Thread分别为1，2，4，8时的LU分解时间和相对于KLU的加速比。从该表可以看出本文的并行方法在1～8个线程并行时可以取得比KLU快1.66～7.72倍的加速比。

注：表中，所有时间列的单位为“秒”，“KLU时间”表示KLU进行LU分解所用的时间，P＝1、2、4、8四列表示本发明提出的方法在1线程、2线程、4线程以及8线程时进行LU分解所用的时间，以及与KLU所用时间相比所获得的加速比。