CN102142052B - 一种针对电路仿真中电路稀疏矩阵的快速lu分解方法 - Google Patents

Abstract

一种针对电路仿真中电路稀疏矩阵的快速LU分解方法，属于EDA技术领域。本发明所述方法包括：符号分析和LU分解计算；其中，所述符号分析是对预处理后的矩阵进行分析预测出该矩阵在完成LU分解后的矩阵L和U的各列非零元结构；所述LU分解计算是基于所述符号分析得出的矩阵L和U的非零元结构，依次对矩阵的各列进行数值求解和数值分配处理，得出LU分解后的矩阵L和矩阵U。本发明的特征在于，将符号分析从LU分解循环中分离出来，在电路仿真过程中可将符号分析放在牛顿-拉夫逊迭代的循环外执行，只执行一次符号分析，能有效降低电路仿真的复杂度，提高LU分解的速度，进而加速整个电路仿真的速度。

Description

一种针对电路仿真中电路稀疏矩阵的快速LU分解方法

技术领域

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

背景技术

在传统的电路设计过程中，为了验证集成电路的连接和功能的完整性，以及预测电路的行为，通常采用在面包板和印刷电路板上做实验来验证设计结果；但是对于集成电路的验证与传统电路不同，在集成电路成品生产之前，无法使用实验的方式来验证其性能。为了提高集成电路正式生产时的良品率并降低集成电路设计成本，有必要在实际生产阶段开始之前，使用计算机辅助集成电路设计(CAD)工具对其电路特性进行检查，以确保该电路的性能。

集成电路从诞生至今，经历了小规模(SSI)、中规模(MSI)、大规模(LSI)、超大规模(VLSI)等多个发展阶段。随着微电子工业的发展，IC设计规模越来越庞大，系统结构越来越复杂，常用的CAD仿真工具(如SPICE、Spectre等)在面对规模和复杂度不断增大的VLSI芯片设计仿真任务时遇到了许多问题，其中之一就是仿真周期随着问题电路的规模呈现超线性的增长，对于常见的AD，PLL等数模混合电路，其前仿周期已经长达几天，后仿周期甚至数月以至于无法进行。IC设计制造业界提出并被大家广泛接受和认可的对仿真软件的要求是：“理想的电路仿真应该在12小时内完成”，这样可以最大程度的实现人机交互运作，提高芯片设计工作的效率，缩短设计周期，大大缩减设计成本。如何提高常用仿真软件的仿真效率，缩短其仿真VLSI电路的周期，以适应现代集成电路设计的要求，是EDA领域一个备受关注的课题。图1中给出电路仿真中瞬态分析的程序流程图，图1中虚线框内的操作为牛顿-拉夫逊迭代的过程，在电路仿真中需要多次执行(对于非线性电路，可能需要执行上万次甚至更多)。这些操作中电路矩阵LU的分解尤为耗时，提高LU分解的速度对整个电路仿真的加速起着举足轻重的作用。

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

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

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

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

发明内容

本发明的目的是提供一种针对电路仿真的快速LU分解算法。该算法调整了向左看算法中不同步骤的执行顺序，如图4所示，该算法将符号分析从循环中分离出来，单独执行。这样利用电路仿真过程中矩阵非零元结构不变的特点，可以将从LU分解中分离出来的符号分析这一步放在牛顿-拉夫逊迭代的循环外执行，因此只计算一次即可，而在牛顿-拉夫逊迭代的循环中多次执行的只有数值分解和数值分配这两步。通过这种调整降低了电路仿真的整体复杂度。修改后的电路仿真瞬态分析的程序流程图如图5所示。

本发明的特征在于，是在计算机中按以下的步骤实现的：

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

步骤(2)，建立n×n的电路稀疏矩阵，包括rajat05，rajat19，add20，circuit_1，circuit_2，add32，meg4，rajat03，coupled，circuit_3，onetone2，ckt11752_dc_1，circuit_4，ASIC_100k，ASIC_100ks，dc1，trans4，G2_circuit，transient，ASIC_320ks，ASIC_680k，ASIC_680ks在内的任何电路稀疏矩阵，

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

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

步骤(5)，按照设定的瞬态仿真迭代次数，依次按以下步骤进行牛顿-拉夫逊迭代：

步骤(5.1)，按John R.Gilbert和Timothy Peierls提出的数值计算方法，根据步骤(4)得到的所述非零结构，对步骤(3)得到的电路稀疏矩阵A进行LU数值分解和数值分配，A＝LU，获得三角方程，Ly＝b和Ux＝y，x和b是电路线性方程组中的n维列向量，y是求解过程中的n维列向量，

步骤(5.2)，解步骤(5.1)得到的三角方程Ly＝b和Ux＝y，并判断解是否收敛：

若：不收敛，则按照杨华中在《电子电路的计算机辅助分析与设计方法》中介绍的方法，更新电路中的器件模型参数并根据更新后的器件模型参数更新所述电路稀疏矩阵A的元素数值，返回步骤(5.1)，

若：已经收敛，则执行步骤(6)，

步骤(6)，判断瞬态仿真是否结束：

若：尚未结束，则计算下一个时间节点，重新依次更新器件模型和矩阵A的元素数值，返回步骤(5.1)，

若：已结束，则打印解的信息，结束。

利用本发明中将符号分析和数值分解分离的方法，可以实现对非零元结构不变，但矩阵内数值变化的矩阵分解问题的优化，尤其是对于电路仿真这种问题，可将原本需执行多次(电路仿真迭代的次数)的符号分析的时间缩短到只执行一次。在整体上，很大程度的降低了电路仿真中LU分解的时间，进而降低电路仿真的时间。实验结果表明本发明中将符号分析和数值分解分离之后，LU分解的速度提高到了原来的1.1-20.1倍(平均加速3.29倍)。

附图说明

图1为电路仿真瞬态分析的程序流程图。

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

图3为向左看LU分解算法的流程。

图4为本发明提出的LU分解方法流程。

图5为采用本发明中LU算法后的电路仿真瞬态分析的程序流程图。

图6为本发明中所采用的符号分析方法中使用到的有向无环图示例，其中，图6.1为前9列

已完成符号分析的下三角矩阵L，图6.2为图6.1中L矩阵前9列对应的有向无环图。

具体实施方式

虽然很多学者在编写LU分解软件时考虑到了电路矩阵的特点，在向左看算法的基础上进行了相应的改进(例如KLU中抛弃了SuperLU中的超节点(Supernode)而使用块三角(BlockTriangular Form，BTF)的方式)，但还没有学者结合电路仿真流程的特点，将LU分解中的符号分析与数值分解分离开执行，以减少相同的符号分析的执行次数，进而降低整个电路仿真过程中LU分解的复杂度，减少处理时间。

本发明提出的降低电路仿真整体复杂度的LU分解方法，结合附图说明如下，其实施步骤依次为：

步骤1：对预处理(交换矩阵的行、列)之后的矩阵进行符号分析，从第1列到第n列依次完成对L和U的n列内非零元结构的计算；

步骤2：利用步骤1中获得的非零结构，对A进行A＝LU的数值分解，得到L和U矩阵。这一步是通过从第1列到第n列依次完成对L和U的n列中所有数值的计算来完成。

详细的步骤说明如下：

步骤1：对预处理(交换矩阵的行、列)之后的矩阵进行符号分析，从第1列到第n列依次完成对L和U的n列内非零元结构的计算；

非零元符号分析采用了John R.Gilbert，Timothy Peierls的文献Sparse PartialPivoting in Time Proportional to Arithmetic Operations中符号分析方法的方法。对矩阵的第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的边。如图6所示的例子，图6为已完成第1-9列符号分析的L；

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

1.4：将所有深度优先搜索出的节点按图的拓扑结构排列，所获得的序列就是第k列的符号分析结果，这一步实际上实质上是将1.3中获得的深度优先遍历的序列按反序排列所得。假设B中有3个非零元{a1，a2，a 3}，通过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矩阵的非零元的符号分析。

以图6中的有向无环图为例，第1-9列的非零元结构已计算完成，现需计算第10列分解后的非零元结构。假设矩阵A的第10列中第4，6，10行的值为非零，即集合B＝{4，6，10}。则从B中的第一个元素4开始获得深度优先遍历序列R₄＝{4，9，12，13}，将这些节点标记为已访问过。从B的第二个元素6开始深度优先遍历可得到R₆＝{6，10}，同样，从B的第三个元素10开始深度优先遍历可获得

(空集)。因此非零元的集合{R₁₀，R₆，R₄＝{6，10，4，9，12，13}就是预测出的L和U在第10列中的非零元的位置。

步骤2：利用步骤1中获得的非零结构，对A进行A＝LU的数值分解，得到L和U矩阵。这一步是通过从第1列到第n列依次完成对L和U的n列中所有数值的计算来完成。数值分析采用文献J.R.Gilbert and T.Peierls，Sparse partial pivoting in time proportionalto arithmetic operations，SIAM J.Sci.Statist.Comput.，vol.9，pp.862-874，1988中的方法。分解过程从第1列到第n列按列依次进行，对每一列的分解可通过如下3步实现：

2.1：展开：将矩阵A的第k列数值存入数组，记这个数组为D_k(是一个列向量)；

2.2：更新：根据步骤1中得到的第k列的非零元结构，从第1个非零元开始，如果这个非零元在对角线上方(即它的行号i＜k)，那么就通过这个非零元来对当前数组中的元素进行更新，更新方法为D_k＝D_k-A_ik*L_i，其中A_ik为矩阵A中第i行第k列位置上的值，L_i为矩阵L中第i列的值(列向量)。这一过程直至第k列中所有非零元都已搜索过。

2.3数值分配：将2.2完成过后的D_k中所有的非零元素，如果非零元在对角线上方(包括对角线，即行号i≤k)，则将这个非零元的行号以及数组中对应的值写入U矩阵，否则将这个非零元的行号以及数组中对应的值与对角元的比值写入L矩阵。

依次对矩阵A的第1-n列执行2.1-2.3这3步操作，便可实现对该矩阵的LU数值分解。由于在电路仿真问题中，步骤1只需执行一次，需要多次循环执行的只有步骤2，因此电路仿真中LU分解的时间主要取决于步骤2。下表显示了对弗罗里达大学稀疏矩阵集(University of Florida Sparse Matrix Collection)中所收集的部分电路矩阵，使用本发明中的方法和LU分解软件KLU测试的分解时间以及本方法对KLU的加速比，测试本方法和KLU的测试环境均为Intel i7 930 CPU，6G DDR3 1066MHz。由于KLU可以配置成使用对角块(BTF)和不使用对角块的模式，为了保证比较的公平性，比较了本方法相比使用对角块和不使用对角块模式下的KLU的性能。

从下表结果可以看出，在使用对角块的模式下，本方法相比KLU能取得1.1-7.3倍的加速，平均加速2.78倍；当不使用对角块时，本方法相比KLU能去的1.1-20.1倍的加速，平均加速3.80倍。综合考虑两种模式下的性能，加速比为3.29倍。

1中预处理使用对角块模式，本方法和KLU使用相同预处理，KLU使用了对角块模式；

2中预处理不使用对角块模式，本方法和KLU使用相同预处理，KLU不使用对角块模式。

Claims (1)

1.一种针对电路仿真中电路稀疏矩阵的快速LU分解方法，其特征在于，是在计算机中按以下的步骤实现的：

步骤（1），输入要解析的电路的网单，

步骤（2），建立n×n的电路稀疏矩阵，包括rajat05，rajat19，add20，circuit_1，circuit_2，add32，meg4，rajat03，coupled，circuit_3，onetone2，ckt11752_dc_1，circuit_4，ASIC_100k，ASIC_100ks，dc1，trans4，G2_circuit，transient，ASIC_320ks，ASIC_680k，ASIC_680ks在内的任何电路稀疏矩阵，

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

步骤（4），对步骤（3）得到的所述n×n的电路稀疏矩阵A从第1列到第n列依次完成对所述电路稀疏矩阵A的下三角矩阵L和上三角矩阵U的n列内的非零结构的计算，其步骤是依次对k=1,2,…,n执行以下步骤：

步骤（4.1），所述电路稀疏矩阵A的第k列的所有非零元的行号构成集合B，B中的元素个数为z，

步骤（4.2），根据所述下三角矩阵L的第1列到第k-1列的非零元，构造有向无环图G，使所述下三角矩阵L的每一列对应所述图G中的每一个顶点，所述下三角矩阵L中第i行第j列的非零元对应所述图G中从顶点j指向顶点i的边，

步骤（4.3），从所述集合B中的每一个元素开始，在所述有向无环图G中做深度优先遍历，获得遍历序列S₁,S₂…,S_z，

步骤（4.4），将所述遍历序列S₁,S₂…,S_z中的所有元素按所述有向无环图G的拓扑结构排序，获得排好序的序列S，所述S中的所有元素即为所述下三角矩阵L和上三角矩阵U的第k列的符号结构，

步骤（5），按照设定的瞬态仿真迭代次数，依次按以下步骤进行牛顿-拉夫逊迭代：

步骤（5.1），根据步骤（4）得到的所述非零结构，对步骤（3）得到的电路稀疏矩阵A进行LU数值分解和数值分配，A=LU，获得三角方程，Ly=b和Ux=y，x和b是电路线性方程组中的n维列向量，y是求解过程中的n维列向量，其步骤是依次对k=1,2,…,n执行以下步骤：

步骤（5.1.1），把所述矩阵A的第k列的数值存入一个列向量D_k，

步骤（5.1.2），对所述上三角矩阵U的第k列中的每一个非零元行号i，计算D_k=D_k-A_ik*L_i，其中A_ik是所述矩阵A中第i行第k列位置上的值，L_i是所述下三角矩阵L第i列的值，

步骤（5.1.3），把所述列向量D_k中行号小于等于k的元素数值存入所述上三角矩阵U的第k列，把所述列向量D_k中行号大于k的元素数值存入所述下三角矩阵L的第k列，

步骤（5.2），解步骤（5.1）得到的三角方程Ly=b和Ux=y，并判断解是否收敛：

若：不收敛，则更新电路中的器件模型参数并根据更新后的器件模型参数更新所述电路稀疏矩阵A的元素数值，返回步骤（5.1），

若：已经收敛，则执行步骤（6），

步骤（6），判断瞬态仿真是否结束：

若：尚未结束，则计算下一个时间节点，重新依次更新器件模型和矩阵A的元素数值，返回步骤（5.1），

若：已结束，则打印解的信息，结束。

Images