專利名稱：電路仿真時(shí)電路稀疏矩陣的基于消去圖的并行分解方法
技術(shù)領(lǐng)域：
本發(fā)明涉及一種針對(duì)電路仿真時(shí)電路稀疏矩陣的基于消去圖的并行LU分解方法，屬于電子設(shè)計(jì)自動(dòng)化(EDA)技術(shù)領(lǐng)域。
背景技術(shù)：
在科學(xué)計(jì)算中，求解線性方程組Ax = b(Ax = b的形式如圖1所示，A為nXn矩陣，b和χ都是η維列向量，A和b已知，χ未知待求)具有舉足輕重的地位。Ax = b的求解方法主要可分為兩大類迭代法和直接法。迭代法如雅可比迭代法、高斯——賽德爾迭代法、超松馳迭代法等，都是先設(shè)定一組解的初值，然后通過(guò)一個(gè)迭代公式進(jìn)行迭代，使得解逐漸向真實(shí)的解靠近，當(dāng)誤差小于給定值時(shí)，迭代收斂。但是迭代法通常對(duì)矩陣A的要求較高，在一般的計(jì)算問(wèn)題中可能難以滿足。不同于迭代法，直接法通過(guò)固定的求解公式直接對(duì) A進(jìn)行有限步驟的計(jì)算，從而求得所需要的解。相比于迭代法，直接法的適用性更強(qiáng)，穩(wěn)定性更好。因此直接法在各類科學(xué)計(jì)算中用得更為普遍。傳統(tǒng)的直接法如高斯消元法、LU分解法(也稱為三角分解法)、喬里斯基分解法等。在一般的問(wèn)題中，矩陣A通常是非常稀疏的。比如對(duì)于電路仿真問(wèn)題，通常電路矩陣中每行的非零元個(gè)數(shù)僅為5個(gè)左右，而且與矩陣規(guī)模無(wú)關(guān)。針對(duì)稀疏矩陣特別設(shè)計(jì)的求解算法，可以大大降低運(yùn)算復(fù)雜度而減少運(yùn)算時(shí)間。在各類直接解法中，運(yùn)用最多的是LU分解法。LU分解法指的是將一個(gè)nXn的方矩陣A分解成一個(gè)下三角矩陣L和一個(gè)上三角矩陣U的乘積(其中U矩陣的對(duì)角線元素都為1)，即A = LU,其中L和U也都是nXn的矩陣。從而求解線性方程組Ax = b的問(wèn)題轉(zhuǎn)化成求解兩個(gè)三角方程Ly = b和te = y (y是η維列向量)。圖2顯示了 LU分解的基本形式，其中L和U未寫數(shù)值的地方都是0。稀疏矩陣的LU分解過(guò)程包括預(yù)處理和數(shù)值分解兩部分。其中預(yù)處理是指使用一定的算法對(duì)矩陣進(jìn)行行列交換以達(dá)到在數(shù)值分解過(guò)程中減少運(yùn)算量的目的，有些預(yù)處理方法中還引入了選主元(選主元指的是將絕對(duì)值大的元素通過(guò)矩陣行列交換操作，交換到對(duì)角線的位置上)的步驟，以保證分解過(guò)程中的數(shù)值穩(wěn)定性。不管采用哪種預(yù)處理方法，預(yù)處理之后的矩陣只是對(duì)原矩陣進(jìn)行了一些行、列的交換，而沒(méi)有其他的變化。本發(fā)明主要針對(duì) LU分解算法的后者，即數(shù)值分解部分，如無(wú)特別說(shuō)明，后文提及的LU分解皆指數(shù)值分解部分。本發(fā)明并不局限于某種特定的預(yù)處理方法，對(duì)于所有預(yù)處理方法都能適用。目前針對(duì)稀疏矩陣的LU分解算法大致可分為向左看算法(Left-looking Algorithm)和向右看算法(Right-looking Algorithm)兩大類，其中向左看算法由于對(duì)稀疏矩陣存儲(chǔ)結(jié)構(gòu)有著很好的適應(yīng)性，在LU分解軟件中得到的廣泛的應(yīng)用。目前加州大學(xué)伯克利分校(University of California, Berkeley) Sherry Li 開發(fā)的 SuperLU 和弗羅里達(dá)大學(xué)(University of Florida) Tim Davis開發(fā)的KLU等軟件的LU分解部分都是以向左看算法為基礎(chǔ)，針對(duì)大規(guī)模稀疏矩陣特點(diǎn)進(jìn)行了優(yōu)化。在向左看LU分解算法中，從第1列到第η列依次分解A的每一列，即每次計(jì)算一個(gè)列向量。對(duì)每一列的分解可以概括為三個(gè)步驟符號(hào)分析、數(shù)值分解和數(shù)值分配，如圖3所示。以正在分解A的第k列為例說(shuō)明。符號(hào)分析指從矩陣A的第k列的非零元結(jié)構(gòu)(非零元結(jié)構(gòu)指的是一個(gè)集合，這個(gè)集合包括第k列上所有非零元素所在的行的行號(hào))計(jì)算出LU分解完成后第k列上的非零元結(jié)構(gòu)。具體的符號(hào)分析所采用的方法，參見文獻(xiàn) J. R. Gilbert and Ε. Ng, Predicting structure in nonsymmetric sparse matrix factorizations, Graph Theory and Sparse Matrix Computation， Springer Verlag， 1993 (該方法可稱為“Gilbert符號(hào)分析方法”)。數(shù)值分解步驟是根據(jù)這一列的符號(hào)分析結(jié)果進(jìn)行數(shù)值計(jì)算，從而獲得第k列的所有非零位置上的數(shù)值結(jié)果。數(shù)值分析所采用的算法， 參見文獻(xiàn) J. R. Gilbert and Τ. Peierls, Sparse partial pivoting in time proportional to arithmetic operations, SIAM J. Sci. Statist. Comput. , vol. 9, pp. 862-874,1988(該算法可稱為Gilbert/Peierls算法，簡(jiǎn)稱G/P算法)。最后一個(gè)步驟是數(shù)值分配，即將第k 列的數(shù)值計(jì)算結(jié)果(是一個(gè)列向量)中行號(hào)小于等于k的部分分配給U矩陣，行號(hào)大于等于k的部分分配給L矩陣。依次對(duì)k=l，2，…，η循環(huán)執(zhí)行上述步驟，即完成了對(duì)整個(gè)矩陣A的LU分解。雖然目前已經(jīng)有很多基于LU分解的軟件包，比如SuperLU，KLU，PARDISO，UMFPack 等共幾十個(gè)軟件，但是在這些軟件中，并行的版本非常少，而針對(duì)通用多核CPU平臺(tái)的并行軟件更是寥寥無(wú)幾。由于LU分解過(guò)程中的高度數(shù)據(jù)依賴性，以及并行時(shí)多線程操作所帶來(lái)的額外代價(jià)，造成了并行的LU分解軟件數(shù)據(jù)如此之少。

發(fā)明內(nèi)容
本發(fā)明的目的是提供一種針對(duì)電路仿真中電路矩陣的并行LU分解方法。本發(fā)明提出使用消去圖來(lái)表示LU分解過(guò)程中的數(shù)據(jù)依賴性并用來(lái)進(jìn)行并行的任務(wù)調(diào)度。該算法從對(duì)電路矩陣的符號(hào)分析結(jié)果中提取出消去圖，然后基于該消去圖的不同結(jié)構(gòu)，提出兩種不同的并行方式群并行與流水線并行。在LU數(shù)值分解過(guò)程中這兩種并行方式動(dòng)態(tài)的調(diào)用，以適應(yīng)不同的數(shù)據(jù)依賴性。電路仿真時(shí)電路稀疏矩陣的基于消去圖的并行分解方法，其特征在于，是在計(jì)算機(jī)中按照以下步驟實(shí)現(xiàn)的步驟(1)，輸入要仿真解析的電路的網(wǎng)單；步驟(2)，建立nXn 的電路稀疏矩陣，包括 add20、circuit_l、circuit_2、 add32、rajat03、coupled、circuit_3、onetonel、onetone2、cktll752_dc_l、circuit_4、 ASIC_100ks、ASIC_100k、del、trans4、G2_circuit、transient、ASIC_320ks、ASIC_320k、 rajat30>ASIC_680ks>ASIC_680k>G3_circuit>FreescaleUrajat31 \)JsR circuit5M^|^ 的任何一個(gè)電路稀疏矩陣；步驟(3)，選擇對(duì)角塊模式和非對(duì)角塊模式中的任何一種模式，對(duì)步驟O)中所建立的電路稀疏矩陣進(jìn)行預(yù)處理，得到預(yù)處理之后的電路稀疏矩陣A ；步驟，根據(jù)John R. GiIbert和Timothy Peierls提出的非零元符號(hào)分析方法， 對(duì)步驟C3)得到的所述nXn的電路稀疏矩陣A進(jìn)行符號(hào)分析，從第1列到第η列依次完成對(duì)所述電路稀疏矩陣A的下三角矩陣L和上三角矩陣U的η列內(nèi)的非零結(jié)構(gòu)的計(jì)算，A = LU ；
步驟(5)，利用步驟(4)得到的U矩陣的非零結(jié)構(gòu)，計(jì)算消去圖消去圖共有η個(gè)頂點(diǎn)，每個(gè)分立的頂點(diǎn)i唯一地對(duì)應(yīng)于所述U矩陣中的第i列，頂點(diǎn)之間不存在任何連線， 但每個(gè)頂點(diǎn)都有一個(gè)“所處的層”屬性值，用于表示每個(gè)頂點(diǎn)以及該頂點(diǎn)在所述U矩陣中對(duì)應(yīng)的列在消去圖中處于哪一層，該屬性值記為level (i)，建立所述消去圖的步驟如下步驟(5. 1)，從第1列到第η列循環(huán)所有列，對(duì)其中的每一列k，取出所述U矩陣中第k列中不包括對(duì)角線的上三角部分的所有非零元，共ζ個(gè)，所述ζ個(gè)非零元所在行的行號(hào)分另U標(biāo)記為j1，j2，…，jz；步驟(5.2)，若ζ = 0，則 level (k) = 1 ；步驟(5.3)，若ζ > 0，則取出所述各行號(hào)對(duì)應(yīng)的相同序號(hào)的列所處的層的屬性值 level α)，level (J2)，…，level (jz)，求出其中的最大值 m = max (level (J1)， Ievel(J2),…，level (jz))，于是，所述第k列在所述消去圖中所處的層的屬性值為 level (k) = m+1 ；步驟(5. 4)，得到所述η列中每一列所處的層的屬性值之后，就得到了整個(gè)消去圖，從該消去圖中所有頂點(diǎn)所處的層的屬性值中求出最大值，即為整個(gè)消去圖的總層數(shù)，記
為 LtOtal ；步驟(6)，根據(jù)步驟(4)得到的L矩陣和U矩陣的非零結(jié)構(gòu)以及步驟(5)得到的消去圖，按以下步驟進(jìn)行并行的LU數(shù)值分解步驟(6. 1)，輸入線程數(shù)目Thread，并創(chuàng)建每個(gè)線程，所述線程數(shù)目Thread的值應(yīng)大于或等于1，小于或等于所使用的計(jì)算機(jī)上的CPU的核心數(shù)；步驟(6. 2)，輸入一個(gè)用于考察所述消去圖中各層頂點(diǎn)總數(shù)V(Iv)的閾值Vth，該閾值Vth在所述線程數(shù)目Thread的1倍-10倍之間取值；步驟(6. 3)，從第1層到第Ltotal層依次考察所述消去圖中每一層Iv的頂點(diǎn)總數(shù) V(Iv)若V(lv) < Vth，則對(duì)第Iv層采用流水線并行方法，若V(lv)≥Vth，則對(duì)第Iv層采用群并行方法，把當(dāng)前任何一種所述并行方法的開始層記為L(zhǎng)start,結(jié)束層記為L(zhǎng)st。p，Lstart = lv,按以下方法從開始層Lstot開始向后考察各層，第lv+1層，…，直到第Lt。al層為止，尋找所述的結(jié)束層Lst。p 在采用流水線并行方法時(shí)，若有一層c滿足V(C)≥Vth，則結(jié)束層Lst。p = c-1，在采用群并行的方法時(shí)，若有一層C’滿足V(c’ ) < Vth，則結(jié)束層Lst。p = c’ -1，步驟(6. 4)，把從步驟(6.3)中得到的從所述消去圖中的第Lstart到第Lst。p層這些層中的所有頂點(diǎn)號(hào)，根據(jù)步驟(6.3)中得到的相應(yīng)的一種并行方法，按以下步驟進(jìn)行多線程方式的并行LU數(shù)值分解步驟(6. 4. 1)，當(dāng)從開始層Lstart到結(jié)束層Lst。p采用群并行方法時(shí)，在循環(huán)進(jìn)行到其中任何一層時(shí)，把該層上的所有頂點(diǎn)平均分配給所述各個(gè)線程，所述各個(gè)線程都采用 John R. Gilbert和Timothy Peierls提出的LU數(shù)值分解方法對(duì)分配到的各列進(jìn)行并行的 LU數(shù)值分解，等待所有線程都計(jì)算完成后再進(jìn)入下一層，循環(huán)往復(fù)，遍歷從開始層Lstart到結(jié)束層Lst。p的所有層，步驟(6. 4. 2)，當(dāng)從開始層Lstot到結(jié)束層Lst。p采用流水線并行方法時(shí)，把從開始層Lstart到結(jié)束層Lst。p的所有層中的列號(hào)按層的順序排成一個(gè)序列，并給這個(gè)序列中每一個(gè)列都置一個(gè)標(biāo)志，初始時(shí)所有標(biāo)志都是“未完成”，然后每個(gè)線程從該序列的頭部開始依次取出一個(gè)列號(hào)，并從序列中刪去所取到的列號(hào)，用John R. Gilbert和Timothy Peierls 提出的LU數(shù)值分解方法對(duì)所取到的列號(hào)進(jìn)行LU數(shù)值分解，在線程完成一列的LU數(shù)值分解后，把該列對(duì)應(yīng)的標(biāo)志寫成“已完成”，每個(gè)線程都重復(fù)這一過(guò)程直到所述序列為空為止；步驟(7)，在滿足收斂條件并完成瞬態(tài)仿真后，打印仿真結(jié)果。在所述步驟(6.4.2)中，在進(jìn)行某一列χ的LU數(shù)值分解時(shí)，要判斷所述U矩陣的第y行、第X列上的值u(y，x)是否為0:若，為0，則對(duì)第χ列繼續(xù)進(jìn)行LU數(shù)值分解，若非0，則第χ列進(jìn)行數(shù)值分解時(shí)依賴于第y列的數(shù)值分解結(jié)果，若第y列的標(biāo)志為“未完成”，則分解第X列的線程需要等待，直到所述第y列的標(biāo)志變成“已完成”后再繼續(xù)進(jìn)行下一步計(jì)算。利用本發(fā)明提出的并行LU分解算法，可以在通用多核CPU平臺(tái)上進(jìn)行并行的LU 分解，從而加速電路仿真的速度。在一系列測(cè)試電路矩陣上的測(cè)試結(jié)果表明本發(fā)明的LU分解時(shí)間在并行線程數(shù)為1 8個(gè)時(shí)比LU分解軟件KLU快1. 66 7. 72倍。本發(fā)明所提出的并行LU分解算法與矩陣中的數(shù)據(jù)依賴性高度協(xié)調(diào)，從而可以保證在計(jì)算過(guò)程中的高度并行化。


圖1為Ax = b的基本形式。圖2為L(zhǎng)U分解的基本形式。圖3為G/P算法的總體流程圖。圖4為非零元符號(hào)分析中的有向無(wú)環(huán)圖舉例圖4. a為下三角矩陣L的一個(gè)例子， 其中第9列已計(jì)算完成，圖4. b為L(zhǎng)的前9列對(duì)應(yīng)的有向無(wú)環(huán)圖。 圖5為上三角矩陣U與其對(duì)應(yīng)的消去圖舉例圖5. a為上三角矩陣U的一個(gè)例子， 圖5. b為與圖5. a對(duì)應(yīng)的消去圖。圖6為針對(duì)電路仿真中電路矩陣的并行LU分解方法的總體流程。圖7為并行LU分解方法的一個(gè)例子，對(duì)應(yīng)于圖5 □表示線程1， Ξ]表示線程2。
具體實(shí)施例方式為了實(shí)現(xiàn)上述目的，本發(fā)明采用如下技術(shù)方案，其實(shí)施步驟為步驟1 對(duì)輸入的ηΧη電路矩陣A進(jìn)行符號(hào)分析，從第1列到第η列依次完成對(duì)L 和U的η列的非零元結(jié)構(gòu)的計(jì)算；步驟2 利用步驟1得到的U矩陣的非零結(jié)構(gòu)，通過(guò)最早開始算法，獲得消去圖；步驟3 基于步驟1得到的L和U的非零結(jié)構(gòu)以及步驟2得到的消去圖，進(jìn)行并行的LU數(shù)值分解。在數(shù)值分解過(guò)程中，采用兩種不同的并行方法群并行和流水線并行。這兩種并行方法在數(shù)值分解過(guò)程中動(dòng)態(tài)地根據(jù)消去圖的結(jié)構(gòu)進(jìn)行選擇。從而完成并行的LU 分解過(guò)程。本發(fā)明提出的針對(duì)電路仿真中電路矩陣的并行LU分解算法，結(jié)合附圖詳細(xì)說(shuō)明如下。本文中所有的η均是指矩陣的維度，即矩陣A是nXn階方陣，L和U分別是η維下三角矩陣和η維上三角矩陣，見圖2。步驟1 對(duì)輸入的nXn電路矩陣A進(jìn)行符號(hào)分析，從第1列到第η列依次完成對(duì)L 和U的η列的非零元結(jié)構(gòu)的計(jì)算；-iJ7 J. R. Gilbert and Ε. Ng, Predicting structure in nonsymmetric sparse matrix factorizations, Graph Theory and Sparse Matrix Computation, Springer-Verlag，1993中符號(hào)分析的方法(該方法可稱為“Gilbert符號(hào)分析方法”)。對(duì)矩陣的第k列進(jìn)行非零元的符號(hào)分析通過(guò)以下1. 1 1. 4完成1. 1 取A的第k列中的所有非零元所在行的行號(hào)為集合B ；1. 2 將對(duì)第1 k 1列的符號(hào)分析結(jié)果中的L矩陣看做一張有向無(wú)環(huán)圖，矩陣L 的每一列對(duì)應(yīng)圖中一個(gè)節(jié)點(diǎn)，即第j列對(duì)應(yīng)圖中節(jié)點(diǎn)j。如果第j列上的第i行為非零元， 即Lu Φ 0，則圖中存在一條從j指向i的邊。如圖4所示的例子，該例子是已完成第1 9 列符號(hào)分析的L以及L所對(duì)應(yīng)的有向無(wú)環(huán)圖；1. 3 從集合B中的每一個(gè)元素為起點(diǎn)在L矩陣對(duì)應(yīng)的有向無(wú)環(huán)圖中進(jìn)行深度優(yōu)先搜索。深度優(yōu)先搜索(D^th-First Search)的過(guò)程可表述為假設(shè)頂點(diǎn)ν為深度優(yōu)先遍歷的起點(diǎn)(源點(diǎn))，那么首先訪問(wèn)出發(fā)點(diǎn)V，并將其標(biāo)記為“已訪問(wèn)過(guò)”；然后依次從ν出發(fā)搜索ν的每個(gè)鄰接點(diǎn)W。若w沒(méi)有被訪問(wèn)過(guò)，則以w為新的出發(fā)點(diǎn)繼續(xù)進(jìn)行深度優(yōu)先搜索，否則頂點(diǎn)w已被訪問(wèn)過(guò)，則跳過(guò)它繼續(xù)尋找ν的所有鄰接點(diǎn)中沒(méi)有被訪問(wèn)過(guò)的頂點(diǎn)。因此深度優(yōu)先遍歷是一個(gè)遞歸的過(guò)程，直到圖中所有的點(diǎn)都被訪問(wèn)過(guò)為止。在遍歷的同時(shí)將B中每一個(gè)元素為起點(diǎn)進(jìn)行深度優(yōu)先遍歷所訪問(wèn)到的頂點(diǎn)號(hào)按順序存入一個(gè)對(duì)應(yīng)的序列，這樣集合B有幾個(gè)元素，就獲得了幾個(gè)遍歷序列；1. 4 將所有深度優(yōu)先搜索出的節(jié)點(diǎn)按圖的拓?fù)浣Y(jié)構(gòu)排列，所獲得的序列就是第k 列的符號(hào)分析結(jié)果，這一步實(shí)際上實(shí)質(zhì)上是將1. 3中獲得的深度優(yōu)先遍歷的序列按反序排列所得。假設(shè)B中有3個(gè)非零元{al，a2，a3}，通過(guò)1. 3的方法深度優(yōu)先搜索得到的序列分別為Ral，Ra2和Ra3，其中Ra3中不包含Ral和Ra2中已搜索的節(jié)點(diǎn)，那么將這三個(gè)序列反向排序，即IRa3UaJ便是按圖的拓?fù)浣Y(jié)構(gòu)的排列。通過(guò)對(duì)第1 η列均執(zhí)行1. 1 1. 4這4步，就可以實(shí)現(xiàn)對(duì)整個(gè)L和U矩陣的非零元的符號(hào)分析。以圖4中的有向無(wú)環(huán)圖為例，第1 9列的非零元結(jié)構(gòu)已計(jì)算完成，現(xiàn)需計(jì)算第10 列分解后的非零元結(jié)構(gòu)。假設(shè)矩陣A的第10列中第4，6，10行的值為非零，即集合B= {4， 6，10}。則從B中的第一個(gè)元素4開始獲得深度優(yōu)先遍歷序列R4= {4，9，12，13}，將這些節(jié)點(diǎn)標(biāo)記為已訪問(wèn)過(guò)。從B的第二個(gè)元素6開始深度優(yōu)先遍歷可得到R6= {6，10}，同樣，從 B的第三個(gè)元素10開始深度優(yōu)先遍歷可獲得RlO= {0}(空集)。因此非零元的集合{R10， R6，R4} = {6,10,4,9,12,13}就是預(yù)測(cè)出的L和U在第10列中的非零元所在的行號(hào)。步驟2 利用步驟1得到的U矩陣的非零結(jié)構(gòu)，通過(guò)最早開始算法，獲得消去圖；在消去圖中，共有η個(gè)頂點(diǎn)，每個(gè)頂點(diǎn)表示矩陣中對(duì)應(yīng)的列，即消去圖中頂點(diǎn)i唯一對(duì)應(yīng)矩陣中的第i列，下文中凡是涉及到消去圖中的“頂點(diǎn)”的說(shuō)法，同時(shí)也是指矩陣中對(duì)應(yīng)的列，兩者為一一對(duì)應(yīng)的關(guān)系。消去圖的物理意義是矩陣中每一列可以進(jìn)行LU數(shù)值分解的最早時(shí)間。在LU數(shù)值分解過(guò)程中，列與列之間存在著數(shù)據(jù)相關(guān)性，比如第χ列依賴于第y列，即第X列的數(shù)值分解會(huì)用到第ι列的數(shù)值分解結(jié)果，那么第X列必須等待第ι列完成數(shù)值分解后才能進(jìn)行數(shù)值分解。消去圖中不存在邊(即頂點(diǎn)與頂點(diǎn)之間不存在線連)， 但每個(gè)頂點(diǎn)都有一個(gè)“所處的層”屬性值，用于記錄每個(gè)頂點(diǎn)在消去圖中處于哪一層，也就是該頂點(diǎn)在矩陣中對(duì)應(yīng)的列可以進(jìn)行數(shù)值分解的最早時(shí)間。對(duì)于第i列在消去圖中所處的層，記為level⑴O步驟2的計(jì)算過(guò)程如下從第1列到第η列循環(huán)所有列，對(duì)其中每一列(記為第k列)，取出U矩陣第k列的上三角部分(不包括對(duì)角線)的所有非零元所在行的行號(hào)，共ζ個(gè)，它們分別為j1;j2，…， jz。如果ζ大于0，則求出所有非零元行號(hào)對(duì)應(yīng)的相同序列的列所處的層的最大值為m，即 m = max (level (J1)，level (J2)，...，level (jz))，于是第 k 列所處的層為 level (k) = m+1。 如果ζ等于0，那么level (k) = 1。如圖5. a的例子所示，這是一個(gè)已完成符號(hào)分解的U矩陣的例子。首先取出第1列的上三角部分(不包括對(duì)角線)的所有非零元，為0個(gè)，那么Ievel(I) =1。對(duì)第2列，同樣也是level (2) = 1。對(duì)第3列，有一個(gè)非零元，行號(hào)為1，所以level (3) = max (level (1))+1 =2。對(duì)第4列，有一個(gè)非零元，行號(hào)為2，所以level (4) = max (level (2))+1 = 2 對(duì)第 5列，有一個(gè)非零元，行號(hào)為3，所以level (5) = max (level (3))+1 = 30同樣的方法獲得 level (6) = Llevel (7) =2。又如第8列，有兩個(gè)非零元，行號(hào)分別為3、5，所以level (8) =max (level (3), level (5))+1 = max (2, 3)+1 = 4，即第 8 列處于第 4 層。以此類推，不再贅述。獲得每一列所處的層之后，即得到了整個(gè)消去圖，從所有頂點(diǎn)所處的層中求出最大值，即得到整個(gè)消去圖的總層數(shù)，記為L(zhǎng)t。tal。圖5. a對(duì)應(yīng)的消去圖如圖5. b所示，共有8層。步驟3 基于步驟1得到的L和U的非零結(jié)構(gòu)以及步驟2得到的消去圖，進(jìn)行并行的LU數(shù)值分解。在數(shù)值分解過(guò)程中，采用兩種不同的并行方法群并行和流水線并行。這兩種并行方法在數(shù)值分解過(guò)程中動(dòng)態(tài)地根據(jù)消去圖的結(jié)構(gòu)進(jìn)行選擇。從而完成并行的LU 分解過(guò)程。步驟3的方法見圖6所示，分以下四小步進(jìn)行。3. 1輸入線程數(shù)目Thread，并創(chuàng)建Thread個(gè)線程。Thread不能超過(guò)所使用的計(jì)算機(jī)上的CPU的核心數(shù)目，也不能小于1。3. 2輸入一個(gè)閾值Vth用于區(qū)分兩種并行方法。Vth是一個(gè)經(jīng)驗(yàn)常數(shù)，沒(méi)有特定的范圍限制。一般可取^Thread的1 10倍。3. 3從消去圖的第1層到第Lt。tal層依次考察每一層，當(dāng)考察到第Iv層時(shí)，檢查消去圖中這一層上的頂點(diǎn)總數(shù)(記為V(lv))，如果V(Iv) <Vth，則這一層將采用流水線并行方法，否則采用群并行方法。記Lstart和Lst。p為當(dāng)前并行方法的開始層和結(jié)束層，置Lstart等于Iv(記為“當(dāng)前層”)。從第Iv層開始向后面的層(即第lv+1層、第lv+2層……直到第 Ltotal層)尋找Lst。p，方法為，如果當(dāng)前層采用流水線并行，那么在第Iv層之后的層中找到第一個(gè)層c，滿足V(C) >Vth;如果當(dāng)前層采用群并行，那么就以同樣的方法在當(dāng)前層之后的層中找滿足V(c’ ) < Vth的第一個(gè)層C’。那么Lst。p = c’ -1。從Lstart層到Lst。p層將采用相同的并行方法。
3. 4將從Lstart到Lst。p這些層中的所有頂點(diǎn)號(hào)，根據(jù)在3. 3中確定的并行方法，進(jìn)行多線程方式的并行的LU數(shù)值分解。群并行和流水線并行兩種方法分別解釋如下。3. 4. 1如果從Lstot層到Lst。p層采用群并行，那么循環(huán)從Lstart到Lst。p的所有層，當(dāng)循環(huán)進(jìn)行到第h層時(shí)，將這一層上的所有列號(hào)平均分給各個(gè)線程，然后各個(gè)線程對(duì)所分配到的列進(jìn)行并行的LU數(shù)值分解，LU數(shù)值分解的方法采用文獻(xiàn)J. R. Gilbert and Τ.Peierls, Sparse partial pivoting in time proportional to arithmetic operations, SIAM J. Sci. Statist. Comput. , vol. 9, pp. 862-874,1988 ψ 的方法，即 Gilbert/Peierls算法，簡(jiǎn)稱G/P算法。每個(gè)線程會(huì)分配到1個(gè)或多個(gè)列號(hào)，每個(gè)線程對(duì)它們逐列進(jìn)行LU數(shù)值分解。等待所有線程都計(jì)算完成后，循環(huán)進(jìn)入第h+Ι層進(jìn)行同樣的任務(wù)分配和計(jì)算，這一過(guò)程直到Lst。p層計(jì)算完成為止。3. 4. 2如果從Lstart層到Lst。p層采用流水線并行，將從Lstmt到Lst。p的所有層中的列號(hào)按層的順序排成一個(gè)序列，并給這個(gè)序列中每一個(gè)列都置一個(gè)標(biāo)志，初始時(shí)所有標(biāo)志都是“未完成”。之后每個(gè)線程從序列的頭部開始依次取一個(gè)列號(hào)并進(jìn)行LU數(shù)值分解計(jì)算， 并將所取的列號(hào)從序列中刪去。這種并行方式下LU數(shù)值分解仍然采用G/P算法。每個(gè)線程完成一列的LU數(shù)值分解計(jì)算后，將這一列對(duì)應(yīng)的標(biāo)志寫成“已完成”狀態(tài)。每個(gè)線程在計(jì)算時(shí)，如果需要用到其他線程的LU數(shù)值分解計(jì)算結(jié)果(即其他列的LU數(shù)值分解結(jié)果)， 并且如果所需要的列的標(biāo)志是“未完成”，則這個(gè)線程需要等待，直到所需要的列的標(biāo)志變成“已完成”后方可以繼續(xù)運(yùn)算。判斷一列(記為第χ列)的數(shù)值分解時(shí)是否要用到另一列(記為第y列)的數(shù)值分解結(jié)果的方法是檢查U (y，x) (U矩陣的第y行、第χ列上的值) 是否為0，如果該值不為0，那么χ列進(jìn)行數(shù)值分解時(shí)依賴于y列的數(shù)值分解結(jié)果，否則不依賴。每個(gè)線程完成一列的LU數(shù)值分解后，將這一列的標(biāo)志改成“已完成”，再?gòu)男蛄兄腥∫粋€(gè)列號(hào)，并將其從序列中刪去。這個(gè)過(guò)程一直重復(fù)到序列為空為止。重復(fù)執(zhí)行3. 3和3. 4兩步驟，對(duì)從消去圖的第1層到第Lt。al層循環(huán)完所有層之后， 就完成了對(duì)矩陣A的并行LU數(shù)值分解。以圖7的例子進(jìn)行舉例說(shuō)明，這個(gè)例子是與圖5對(duì)應(yīng)的。在這個(gè)例子中假設(shè)線程數(shù)Thread = 2，閾值Vth = 2。從消去圖的第一層開始考察，第一層有3個(gè)列號(hào)，那么采用群并行方法，并找到Lstart = 1和Lst。p = 3。于是從第1層到第3層都采用群并行方法，即 Lstart = 1和Lst。p = 3。首先數(shù)值分解第1層，將第1層的列號(hào)平均分給2個(gè)線程，線程1分配到第1列和第2列，線程2分配到第6列，于是這兩個(gè)線程各自獨(dú)立地對(duì)分配到的列進(jìn)行 LU數(shù)值分解計(jì)算。等兩個(gè)線程都完成后，同樣的方法進(jìn)入第2層的計(jì)算，依此類推第3層。 之后考察第4層，只有一個(gè)列號(hào)，于是采用流水線并行方法，并得到Lstart = 4和Lst。p = 8。 將從第4層到第8層中所有列號(hào)排成一個(gè)序列得{8，11，12，13，14}，線程1從中取8，并將 8從序列中刪去，線程2從中取出11，并將11從序列中刪去，然后這兩個(gè)線程各自獨(dú)立地進(jìn)行LU數(shù)值分解。由于線程2所分解的第11列依賴于線程1所分解的第8列的LU數(shù)值分解結(jié)果(因?yàn)閳D5. a中U(8，11)不為0)，所以當(dāng)線程2的LU數(shù)值分解計(jì)算過(guò)程中需要用到第8列的LU數(shù)值分解結(jié)果時(shí)，線程2需要等待直到第8列的狀態(tài)變成“已完成”后才可以繼續(xù)分解第11列。當(dāng)線程1完成第8列的LU數(shù)值分解后，將第8列標(biāo)志成“已完成”，從序列中取出12，并將12從序列中刪去，然后開始對(duì)第12列進(jìn)行LU數(shù)值分解，同樣的道理，第 12列的LU數(shù)值分解可能需要等待第11列的LU數(shù)值分解結(jié)果，對(duì)后面的列以此類推。
本發(fā)明提出的針對(duì)電路仿真中電路矩陣的并行LU分解方法擁有良好的加速性能。測(cè)試結(jié)果舉例如下表，所用測(cè)試矩陣均是從美國(guó)佛羅里達(dá)大學(xué)矩陣收集網(wǎng)站所下載獲得。該表顯示了本發(fā)明提出的并行LU分解方法與KLU進(jìn)行分解時(shí)間比較的結(jié)果。表中所統(tǒng)計(jì)的時(shí)間單位均為秒，“KLU時(shí)間”這一列表示KLU所用的LU分解時(shí)間，P = 1，2，4，8這四列分別表示本發(fā)明在線程數(shù)Thread分別為1，2，4，8時(shí)的LU分解時(shí)間和相對(duì)于KLU的加速比。從該表可以看出本文的并行方法在1 8個(gè)線程并行時(shí)可以取得比KLU快1.66 7. 72倍的加速比。
權(quán)利要求
1.電路仿真時(shí)電路稀疏矩陣的基于消去圖的并行分解方法，其特征在于，是在計(jì)算機(jī)中按照以下步驟實(shí)現(xiàn)的步驟(1)，輸入要仿真解析的電路的網(wǎng)單；步驟(2)，建立nXn的電路稀疏矩陣，包括add20、circuit」、circuit_2、add32、 rajat03、 coupled、 circuit_3、 onetonel、 onetone2、 ckt11752_dc_l、 circuit_4、 ASIC_100ks、ASIC_100k、del、trans4、G2_circuit、transient、ASIC_320ks、ASIC_320k、 rajat30>ASIC_680ks>ASIC_680k>G3_circuit>FreescaleUrajat31 \)JsR circuit5M^|^ 的任何一個(gè)電路稀疏矩陣；步驟(3)，選擇對(duì)角塊模式和非對(duì)角塊模式中的任何一種模式，對(duì)步驟O)中所建立的電路稀疏矩陣進(jìn)行預(yù)處理，得到預(yù)處理之后的電路稀疏矩陣A ；步驟G)，根據(jù)John R. GiIbert和Timothy Peierls提出的非零元符號(hào)分析方法，對(duì)步驟(3)得到的所述nXn的電路稀疏矩陣A進(jìn)行符號(hào)分析，從第1列到第η列依次完成對(duì)所述電路稀疏矩陣A的下三角矩陣L和上三角矩陣U的η列內(nèi)的非零結(jié)構(gòu)的計(jì)算，A = LU5 步驟(5)，利用步驟(4)得到的U矩陣的非零結(jié)構(gòu)，計(jì)算消去圖消去圖共有η個(gè)頂點(diǎn)， 每個(gè)分立的頂點(diǎn)i唯一地對(duì)應(yīng)于所述U矩陣中的第i列，頂點(diǎn)之間不存在任何連線，但每個(gè)頂點(diǎn)都有一個(gè)“所處的層”屬性值，用于表示每個(gè)頂點(diǎn)以及該頂點(diǎn)在所述U矩陣中對(duì)應(yīng)的列在消去圖中處于哪一層，該屬性值記為level (i)，建立所述消去圖的步驟如下步驟(5. 1)，從第1列到第η列循環(huán)所有列，對(duì)其中的每一列k，取出所述U矩陣中第k 列中不包括對(duì)角線的上三角部分的所有非零元，共ζ個(gè)，所述ζ個(gè)非零元所在行的行號(hào)分別標(biāo)記為丄，j2，...，jz;步驟(5. 2)，若ζ = 0，則 level (k) = 1 ；步驟(5. 3)，若z > 0，則取出所述各行號(hào)對(duì)應(yīng)的相同序號(hào)的列所處的層的屬性值 level (丄)，level (j2)，…，Ievel(Jz),求出其中的最大值 m = max (level (J1), Ievel(J2),…，level (jz))，于是，所述第k列在所述消去圖中所處的層的屬性值為 level (k) = m+1 ；步驟(5. 4)，得到所述η列中每一列所處的層的屬性值之后，就得到了整個(gè)消去圖， 從該消去圖中所有頂點(diǎn)所處的層的屬性值中求出最大值，即為整個(gè)消去圖的總層數(shù)，記為L(zhǎng)total ‘步驟(6)，根據(jù)步驟⑷得到的L矩陣和U矩陣的非零結(jié)構(gòu)以及步驟(5)得到的消去圖，按以下步驟進(jìn)行并行的LU數(shù)值分解步驟(6. 1)，輸入線程數(shù)目Thread，并創(chuàng)建每個(gè)線程，所述線程數(shù)目Thread的值應(yīng)大于或等于1，小于或等于所使用的計(jì)算機(jī)上的CPU的核心數(shù)；步驟(6. 2)，輸入一個(gè)用于考察所述消去圖中各層頂點(diǎn)總數(shù)V(Iv)的閾值Vth，該閾值 Vth在所述線程數(shù)目Thread的1倍-10倍之間取值；步驟(6.3)，從第1層到第Lt。tal層依次考察所述消去圖中每一層Iv的頂點(diǎn)總數(shù)V(Iv) 若V(lv) < Vth，則對(duì)第Iv層采用流水線并行方法， 若V(lv) ^ Vth，則對(duì)第Iv層采用群并行方法，把當(dāng)前任何一種所述并行方法的開始層記為L(zhǎng)start,結(jié)束層記為L(zhǎng)st。p，Lstart = lv,按以下方法從開始層Lstart開始向后考察各層，第lv+1層，…，直到第Ltotal層為止，尋找所述的結(jié)束層Lstop 在采用流水線并行方法時(shí)，若有一層C滿足V(C)彡Vth，則結(jié)束層Lst。p = C-1, 在采用群并行的方法時(shí)，若有一層c’滿足V(c’ ) < Vth，則結(jié)束層Lst。p = c’ -1， 步驟(6. 4)，把從步驟(6.3)中得到的從所述消去圖中的第Lstart到第Lst。p層這些層中的所有頂點(diǎn)號(hào)，根據(jù)步驟(6.3)中得到的相應(yīng)的一種并行方法，按以下步驟進(jìn)行多線程方式的并行LU數(shù)值分解步驟(6. 4. 1)，當(dāng)從開始層Lstart到結(jié)束層Lst。p采用群并行方法時(shí)，在循環(huán)進(jìn)行到其中任何一層時(shí)，把該層上的所有頂點(diǎn)平均分配給所述各個(gè)線程，所述各個(gè)線程都采用John R. Gilbert和Timothy Peierls提出的LU數(shù)值分解方法對(duì)分配到的各列進(jìn)行并行的LU數(shù)值分解，等待所有線程都計(jì)算完成后再進(jìn)入下一層，循環(huán)往復(fù)，遍歷從開始層Lstart到結(jié)束層 Lstop的所有層，步驟(6. 4. 2)，當(dāng)從開始層Lstart到結(jié)束層Lst。p采用流水線并行方法時(shí)，把從開始層 Lstart到結(jié)束層Lst。p的所有層中的列號(hào)按層的順序排成一個(gè)序列，并給這個(gè)序列中每一個(gè)列都置一個(gè)標(biāo)志，初始時(shí)所有標(biāo)志都是“未完成”，然后每個(gè)線程從該序列的頭部開始依次取出一個(gè)列號(hào)，并從序列中刪去所取到的列號(hào)，用John R. Gilbert和Timothy Peierls提出的LU數(shù)值分解方法對(duì)所取到的列號(hào)進(jìn)行LU數(shù)值分解，在線程完成一列的LU數(shù)值分解后， 把該列對(duì)應(yīng)的標(biāo)志寫成“已完成”，每個(gè)線程都重復(fù)這一過(guò)程直到所述序列為空為止； 步驟(7)，在滿足收斂條件并完成瞬態(tài)仿真后，打印仿真結(jié)果。
2.根據(jù)權(quán)利要求書1所述的電路仿真時(shí)電路稀疏矩陣的基于消去圖的并行分解方法， 其特征在于，在所述步驟(6.4. 2)中，在進(jìn)行某一列χ的LU數(shù)值分解時(shí)，要判斷所述U矩陣的第y行、第χ列上的值U (y，χ)是否為0 若，為0，則對(duì)第χ列繼續(xù)進(jìn)行LU數(shù)值分解，若非0，則第χ列進(jìn)行數(shù)值分解時(shí)依賴于第y列的數(shù)值分解結(jié)果，若第y列的標(biāo)志為 “未完成”，則分解第χ列的線程需要等待，直到所述第y列的標(biāo)志變成“已完成”后再繼續(xù)進(jìn)行下一步計(jì)算。
全文摘要
電路仿真時(shí)電路稀疏矩陣的基于消去圖的并行分解方法，屬于電子設(shè)計(jì)自動(dòng)化(EDA)領(lǐng)域，其特征在于，從電路矩陣的符號(hào)分析結(jié)果中提取出消去圖，用來(lái)表示LU分解過(guò)程中的數(shù)據(jù)依賴性，根據(jù)消去圖的不同結(jié)構(gòu)使用兩種不同的并行方法基于群的并行和流水線的并行，從而降低LU分解的運(yùn)算時(shí)間、加快電路仿真的速度，在一系列測(cè)試電路矩陣上的測(cè)試結(jié)果表明本發(fā)明的LU分解時(shí)間在并行線程數(shù)為1～8個(gè)時(shí)比LU分解軟件KLU快1.66～7.72倍。
文檔編號(hào)G06F17/50GK102156777SQ20111008808
公開日2011年8月17日 申請(qǐng)日期2011年4月8日 優(yōu)先權(quán)日2011年4月8日
發(fā)明者楊華中, 武偉, 汪玉, 陳曉明 申請(qǐng)人:清華大學(xué)