专利名称：一种针对电路仿真中电路稀疏矩阵的快速lu分解方法
技术领域：
本发明涉及一种针对电路仿真中电路稀疏矩阵的快速LU分解方法，属于电子设 计自动化(EDA)技术领域。
背景技术：
在传统的电路设计过程中，为了验证集成电路的连接和功能的完整性，以及预测 电路的行为，通常采用在面包板和印刷电路板上做实验来验证设计结果；但是对于集成电 路的验证与传统电路不同，在集成电路成品生产之前，无法使用实验的方式来验证其性能。 为了提高集成电路正式生产时的良品率并降低集成电路设计成本，有必要在实际生产阶段 开始之前，使用计算机辅助集成电路设计(CAD)工具对其电路特性进行检查，以确保该电 路的性能。集成电路从诞生至今，经历了小规模(SSI)、中规模(MSI)、大规模(LSI)、超大规 模(VLSI)等多个发展阶段。随着微电子工业的发展，IC设计规模越来越庞大，系统结构越 来越复杂，常用的CAD仿真工具(如SPICE、Spectre等)在面对规模和复杂度不断增大的 VLSI芯片设计仿真任务时遇到了许多问题，其中之一就是仿真周期随着问题电路的规模呈 现超线性的增长，对于常见的AD，PLL等数模混合电路，其前仿周期已经长达几天，后仿周 期甚至数月以至于无法进行。IC设计制造业界提出并被大家广泛接受和认可的对仿真软件 的要求是“理想的电路仿真应该在12小时内完成”，这样可以最大程度的实现人机交互运 作，提高芯片设计工作的效率，缩短设计周期，大大缩减设计成本。如何提高常用仿真软件 的仿真效率，缩短其仿真VLSI电路的周期，以适应现代集成电路设计的要求，是EDA领域一 个备受关注的课题。图1中给出电路仿真中瞬态分析的程序流程图，图1中虚线框内的操 作为牛顿-拉夫逊迭代的过程，在电路仿真中需要多次执行(对于非线性电路，可能需要执 行上万次甚至更多)。这些操作中电路矩阵LU的分解尤为耗时，提高LU分解的速度对整个 电路仿真的加速起着举足轻重的作用。LU分解法指的是将一个nXn的方矩阵A分解成一个下三角矩阵L和一个上三角 矩阵U的乘积，即A = LU，其中L和U也都是nXn的矩阵。从而求解线性方程组Ax = b (χ 和b均为η维列向量)的问题转化成求解两个三角方程Ly = b和te = y (y是η维列向 量)。图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列到第η列依次分解每一列，即每次计算一个列 向量。对每一列的分解可以概括为三个步骤符号分析、数值分解和数值分配，如图3所 示。以正在计算第k列为例说明。符号分析指从矩阵A的第k列的非零元结构(非零元 结构指的是一个集合，这个集合包括第k列上所有非零元素的行号)计算出LU分解完成 后第k列上的非零元结构。具体的符号分析所采用的算法，参见文献J.R.Gi lbert and Ε. Ng, Predictingstructure in nonsymmetric sparse matrix factorizations, Graph Theory and SparseMatrix Computation，Springer-Verlag，1993。数值分析步骤是根据 这一列的符号分析结果进行数值计算，从而获得第k列的所有非零位置上的数值结果。数 值分析所采用的算法，参见文献 J.R. Gilbert and Τ. 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，. . .，η循环执行上述步骤，即完成了对整个矩阵A的LU分解。
发明内容
本发明的目的是提供一种针对电路仿真的快速LU分解算法。该算法调整了向左 看算法中不同步骤的执行顺序，如图4所示，该算法将符号分析从循环中分离出来，单独执 行。这样利用电路仿真过程中矩阵非零元结构不变的特点，可以将从LU分解中分离出来的 符号分析这一步放在牛顿-拉夫逊迭代的循环外执行，因此只计算一次即可，而在牛顿-拉 夫逊迭代的循环中多次执行的只有数值分解和数值分配这两步。通过这种调整降低了电路 仿真的整体复杂度。修改后的电路仿真瞬态分析的程序流程图如图5所示。本发明的特征在于，是在计算机中按以下的步骤实现的步骤(1)，输入要解析的电路的网单，步骤(2)，建立ηΧη 的电路稀疏矩阵，包括 rajat05，rajatl9，add20，circuit_l， circuit_2, add32, meg4, rajat03, coupled, circuit_3, onetone2, cktl1752_dc_l, circuit_4, ASIC一100k， ASIC一100ks， del, trans4, G2_circuit, transient, ASIC_320ks, ASIC_680k, ASIC_680ks在内的任何电路稀疏矩阵，步骤(3)，选择对角块模式和非对角块模式中的任何一种模式，对步骤O)中所建 立的电路稀疏矩阵进行预处理，得到预处理之后的电路稀疏矩阵A，步骤，根据John R. Gilbert和Timothy Peierls提出的非零元符号分析方法， 对步骤C3)得到的所述ηΧη的电路稀疏矩阵A从第1列到第η列依次完成对所述电路稀 疏矩阵A的下三角矩阵L和上三角矩阵U的η列内的非零结构的计算，步骤( ，按照设定的瞬态仿真迭代次数，依次按以下步骤进行牛顿-拉夫逊迭 代
步骤(5. 1)，按John R. Gilbert和Timothy Peierls提出的数值计算方法，根据步 骤(4)得到的所述非零结构，对步骤C3)得到的电路稀疏矩阵A进行LU数值分解和数值分 配，A = LU，获得三角方程，Ly = b和te = y，χ和b是电路线性方程组中的η维列向量，y 是求解过程中的η维列向量，步骤(5. 2)，解步骤(5. 1)得到的三角方程Ly = b和te = 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列到第η 列依次完成对L和U的η列内非零元结构的计算；步骤2 利用步骤1中获得的非零结构，对A进行A = LU的数值分解，得到L和U 矩阵。这一步是通过从第1列到第η列依次完成对L和U的η列中所有数值的计算来完成。详细的步骤说明如下步骤1 对预处理(交换矩阵的行、列)之后的矩阵进行符号分析，从第1列到第η 列依次完成对L和U的η列内非零元结构的计算；非零元符号分析采用了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 将对第l-k-Ι列的符号分析结果中的L矩阵看做一张有向无环图，矩阵L的 每一列对应图中一个节点，即第j列对应图中节点j。如果第j列上的第i行为非零元，即 Lij Φ 0，则图中存在一条从j指向i的边。如图6所示的例子，图6为已完成第1-9列符号 分析的L ；1. 3 从集合B中的每一个元素为起点在L矩阵对应的有向无环图中进行深度优先 搜索。深度优先搜索(D^th-First Search)的过程可表述为假设顶点ν为深度优先遍历 的起点(源点)，那么首先访问出发点V，并将其标记为“已访问过”；然后依次从ν出发搜 索ν的每个邻接点W。若w没有被访问过，则以w为新的出发点继续进行深度优先搜索，否 则顶点w已被访问过，则跳过它继续寻找下一个没有被访问过的顶点。因此深度优先遍历 是一个递归的过程，直到图中所有的点都被访问过为止。在遍历的同时将B中每一个元素 为起点进行深度优先遍历所访问到的顶点号按顺序存入一个对应的序列，这样集合B有几 个元素，就获得了几个遍历序列；1. 4 将所有深度优先搜索出的节点按图的拓扑结构排列，所获得的序列就是第k 列的符号分析结果，这一步实际上实质上是将1. 3中获得的深度优先遍历的序列按反序排 列所得。假设B中有3个非零元{al，a2, a 3}，通过1. 3的方法深度优先搜索得到的序列 分别为Ral，Ra2和Ra3，其中Ra3中不包含Ral和Ra2中已搜索的节点，那么将这三个序列反向 排序，即IRa3UaJ便是按图的拓扑结构的排列。通过对第1-Π列均执行1. 1-1. 4这4步，就可以实现对整个L和U矩阵的非零元 的符号分析。以图6中的有向无环图为例，第1-9列的非零元结构已计算完成，现需计算第10 列分解后的非零元结构。假设矩阵A的第10列中第4，6，10行的值为非零，即集合B= {4， 6，10}。则从B中的第一个元素4开始获得深度优先遍历序列R4= {4，9，12，13}，将这些节 点标记为已访问过。从B的第二个元素6开始深度优先遍历可得到{6，10}，同样，从 B的第三个元素10开始深度优先遍历可获得兄(空集)。因此非零元的集合{R1(l，R6, R4= {6,10,4,9,12,13}就是预测出的L和U在第10列中的非零元的位置。步骤2 利用步骤1中获得的非零结构，对A进行A = LU的数值分解，得到L和 U矩阵。这一步是通过从第1列到第η列依次完成对L和U的η列中所有数值的计算来 完成。数值分析米用文献 J. R. Gilbert and Τ. Peierls，Sparse partial pivoting intime proportionalto arithmetic operations, SIAM J. Sci. Statist. Comput. , vol.9, pp. 862-874，1988中的方法。分解过程从第1列到第η列按列依次进行，对每一列的分解可 通过如下3步实现2. 1 展开将矩阵A的第k列数值存入数组，记这个数组为Dk (是一个列向量)；2. 2 更新根据步骤1中得到的第k列的非零元结构，从第1个非零元开始，如果 这个非零元在对角线上方(即它的行号i < k)，那么就通过这个非零元来对当前数组中的 元素进行更新，更新方法为Dk = Dk-AikIi,其中Aik为矩阵A中第i行第k列位置上的值，Li 为矩阵L中第i列的值(列向量)。这一过程直至第k列中所有非零元都已搜索过。2. 3数值分配将2. 2完成过后的Dk中所有的非零元素，如果非零元在对角线上方 (包括对角线，即行号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的测试环境均为htel 7 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. 一种针对电路仿真中电路稀疏矩阵的快速LU分解方法，其特征在于，是在计算机中 按以下的步骤实现的步骤(1)，输入要解析的电路的网单，步骤(2)，建立 nXn 的电路稀疏矩阵，包括 rajat05，rajatl9，add20, circuit_l, circuit_2, add32, meg4, rajat03, coupled, circuit_3, onetone2, cktll752_dc_l, circuit_4, ASIC一100k， ASIC一lOOks， del, trans4, G2_circuit, transient, ASIC_320ks, ASIC_680k, ASIC_680ks在内的任何电路稀疏矩阵，步骤(3)，选择对角块模式和非对角块模式中的任何一种模式，对步骤O)中所建立的 电路稀疏矩阵进行预处理，得到预处理之后的电路稀疏矩阵A，步骤G)，根据John R. Gilbert和Timothy Peierls提出的非零元符号分析方法，对步 骤C3)得到的所述nXn的电路稀疏矩阵A从第1列到第η列依次完成对所述电路稀疏矩 阵A的下三角矩阵L和上三角矩阵U的η列内的非零结构的计算，步骤(5)，按照设定的瞬态仿真迭代次数，依次按以下步骤进行牛顿-拉夫逊迭代 步骤(5. 1)，按John R. Gilbert和Timothy Peierls提出的数值计算方法，根据步骤 (4)得到的所述非零结构，对步骤C3)得到的电路稀疏矩阵A进行LU数值分解和数值分配， A = LU，获得三角方程，Ly = b和te = y，χ和b是电路线性方程组中的η维列向量，y是 求解过程中的η维列向量，步骤(5. 2)，解步骤(5. 1)得到的三角方程Ly = b和te = y，并判断解是否收敛 若不收敛，则按照杨华中在《电子电路的计算机辅助分析与设计方法》中介绍的方法， 更新电路中的器件模型参数并根据更新后的器件模型参数更新所述电路稀疏矩阵A的元 素数值，返回步骤(5. 1)，若已经收敛，则执行步骤(6)， 步骤(6)，判断瞬态仿真是否结束若尚未结束，则计算下一个时间节点，重新依次更新器件模型和矩阵A的元素数值， 返回步骤(5. 1)，若已结束，则打印解的信息，结束。
全文摘要
一种针对电路仿真中电路稀疏矩阵的快速LU分解方法，属于EDA技术领域。本发明所述方法包括符号分析和LU分解计算；其中，所述符号分析是对预处理后的矩阵进行分析预测出该矩阵在完成LU分解后的矩阵L和U的各列非零元结构；所述LU分解计算是基于所述符号分析得出的矩阵L和U的非零元结构，依次对矩阵的各列进行数值求解和数值分配处理，得出LU分解后的矩阵L和矩阵U。本发明的特征在于，将符号分析从LU分解循环中分离出来，在电路仿真过程中可将符号分析放在牛顿-拉夫逊迭代的循环外执行，只执行一次符号分析，能有效降低电路仿真的复杂度，提高LU分解的速度，进而加速整个电路仿真的速度。
文档编号G06F17/50GK102142052SQ20111007602
公开日2011年8月3日 申请日期2011年3月28日 优先权日2011年3月28日
发明者杨华中, 武伟, 汪玉, 陈晓明 申请人:清华大学