基于細(xì)胞稀疏存貯方案的有限元?jiǎng)偠染仃嚱M裝

1、基于細(xì)胞稀疏存貯方案的有限元?jiǎng)偠染仃嚱M裝1陳璞1,陳斌11北京大學(xué)力學(xué)與工程科學(xué)系,北京 (100871E-mail:chenpu摘要:在現(xiàn)代有限元分析的直接解法和迭代解法中都需要組裝剛度矩陣。本文在1,2的基礎(chǔ)上,敘述了兩種基于細(xì)胞稀疏存貯方案有限元的符號(hào)剛度矩陣的組裝方法。在與傳統(tǒng)的方法進(jìn)行了比較的基礎(chǔ)上,提出了改進(jìn)效率的方法。關(guān)鍵詞:有限元分析稀疏矩陣高性能計(jì)算中圖分類號(hào):O241.821.引言近年來(lái),稀疏解法逐漸替代了帶寬解法與變帶寬解法,成為了有限元分析中首選的直接解法1,2,3。不同于帶寬解法或變帶寬解存貯方案,稀疏解法要求建立符號(hào)矩陣(symbolic matrix,它是總體剛度

2、矩陣對(duì)應(yīng)的Bool值矩陣。剛度矩陣中的非零元在Bool值矩陣為1,其余則為0。在用圖討論矩陣運(yùn)算時(shí)符號(hào)矩陣又稱為鄰接圖(adjacent graph。在稀疏解法中符號(hào)矩陣涉及幾乎所有的步驟,如:符號(hào)組裝;填充元優(yōu)化;符號(hào)分解;(數(shù)值剛度矩陣的組裝;數(shù)值分解;消元與回帶等。因此,減小符號(hào)矩陣的大小在稀疏解法中具有十分重要的意義。傳統(tǒng)的稀疏符號(hào)矩陣的組裝方法直接采用了基于方程的數(shù)據(jù)結(jié)構(gòu)5,符號(hào)矩陣占用的內(nèi)存為數(shù)值剛度矩陣占用內(nèi)存的1/2(4字節(jié)整數(shù),8字節(jié)實(shí)數(shù)。對(duì)于給定的單元類型,符號(hào)矩陣的大小一般正比于節(jié)點(diǎn)數(shù)。如果分析的問(wèn)題不太大,符號(hào)矩陣可以整個(gè)在內(nèi)存中處理。但模型較大時(shí)有限元的符號(hào)矩陣的可能

3、超出內(nèi)存允許的范圍。于是就出現(xiàn)了各種符號(hào)矩陣的壓縮存貯方法,其中最著名的是文獻(xiàn)6的方法。在工程有限元分析中,因?yàn)榍蠼馕灰葡蛄繄?chǎng)的需要,總體剛度矩陣是由一系列小的子矩陣組成。利用這一性質(zhì),文1,2給出了一個(gè)稱為細(xì)胞稀疏存貯方案的符號(hào)矩陣壓縮方法,其壓縮比大于Sherman的方法6。本文從兩個(gè)方面改進(jìn)符號(hào)矩陣的組裝方法:1利用文1,2的超方程和細(xì)胞符號(hào)矩陣概念直接組裝細(xì)胞符號(hào)矩陣。在工程有限元分析中,它的大小僅是傳統(tǒng)符號(hào)矩陣的1/36到1/9;2給出了兩種不同的組裝方案。它們分別適用于不同的內(nèi)存數(shù)據(jù)安排。數(shù)值試驗(yàn)表明,本文提出的方法可以在僅有28M內(nèi)存的情形下在內(nèi)存快速地組裝100萬(wàn)階3維8節(jié)點(diǎn)實(shí)

4、體有限元方程組的符號(hào)矩陣。2.稀疏矩陣的存貯方案為了更清楚地討論問(wèn)題,取一個(gè)6階的總體剛度矩陣如下:1 高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金資助項(xiàng)目(20030001112- 1 - 2 -=665544332211h SYM g f ed c ba A (11(22(3344(55(66a b c d e f g SYMh =(1 其中a ,b ,h 是非零實(shí)數(shù)。由于對(duì)稱性,有限元?jiǎng)偠染仃嚥捎昧松先切写尜A方法(或等價(jià)地,下三角列存貯方法。為了有效地對(duì)行實(shí)施運(yùn)算,在有限元分析的計(jì)算中矩陣A 較常用的稀疏存貯方案是緊湊行向量的順序組合表示法(又稱為帶輔助行向量的二元組表示法。這種表示法需要3個(gè)數(shù)組

5、IA(1:neq, JA(1:nzr 和 PA(1:nzr 來(lái)表示矩陣 A 的上三角部分4,5,其中neq 和nzr 分別是線性方程組的階數(shù)和是非零元的個(gè)數(shù)。以4字節(jié)整數(shù)和8字節(jié)實(shí)數(shù)計(jì),這種存貯方案占用的內(nèi)存為 (neq+nzr*4 + nzr*8 字節(jié)。表1 緊湊行向量的順序組合表示法方程號(hào)12 3 4 5 6IA 4 6 9 11 13 14JA1 3 4 62 5345 4 5 56 6PA 11a b c 22d33e f 44g55h66數(shù)組IA(1:neq=6是行索引,它記錄了每一緊湊行向量在非零元的列指標(biāo)數(shù)組JA(1:nzr=14 以及它的數(shù)值數(shù)組 PA(1:nzr=14 中的末

6、地址。第k 行的非零元個(gè)數(shù)為IA(k-IA(k-1(假設(shè)IA(0=0。文1,2提出了一種細(xì)胞存貯方案。它的想法是將剛度矩陣中具有相同的非零元位置的行(列所對(duì)應(yīng)的方程定義為超方程,與超方程劃分相應(yīng)的子矩陣為細(xì)胞。式(1中A 的第2個(gè)表達(dá)式就是按細(xì)胞方式寫的。顯然,符號(hào)矩陣可以在細(xì)胞意義表達(dá),在給定超方程劃分的意義下,它與原始的表達(dá)是完全等價(jià)的。在這一表示法之下,有限元的總體剛度矩陣被視作超行向量的順序組合。根據(jù)超方程的定義,式(1中有5個(gè)超方程。記mneq 為超方程數(shù),在超行緊湊方式下,式(1中的第2個(gè)矩陣可以用5個(gè)數(shù)組來(lái)描述。表2 緊湊行向量的細(xì)胞稀疏存貯方案超方程1 2 3 4 5SUPER

7、 1 2 4 5 6 ISA 3 5 7 9 10 JSA 1 3 5 2 4 3 4 5 6 6 LSA 4 6 11 13 14PSA 11(b a c 22 d4433eg f 55 h 66數(shù)組 SUPER (1:mneq=5是超方程編號(hào)到原方程編號(hào)的映射。數(shù)組 ISA 是非零細(xì)胞的超列指標(biāo)(或稱細(xì)胞指標(biāo) JSA 的索引,LSA是數(shù)值數(shù)組 PSA 以超方程為單位的索引。這種存貯方案占用的內(nèi)存為 (mneq*3+nzrc*4 + nzr*8 字節(jié)。一般的稀疏存貯方案與細(xì)胞存貯方案之間的差異在于列指標(biāo)的數(shù)量。細(xì)胞存貯方案的指標(biāo)比傳統(tǒng)存貯方案大大減少。對(duì)于三維實(shí)體分析,細(xì)胞指標(biāo)數(shù)nzrc僅僅

8、是傳統(tǒng)方案列指標(biāo)數(shù)nzr的約1/9。對(duì)于三維殼體分析,則指標(biāo)數(shù)相對(duì)更少。使用細(xì)胞存貯方案的好處不僅在于指標(biāo)集合的大幅度減小,而且可以大幅度地減少矩陣運(yùn)算中的指標(biāo)操作,從而達(dá)到提高符號(hào)運(yùn)算和數(shù)值運(yùn)算速度的目的。結(jié)合循環(huán)展開(kāi)技術(shù)8,筆者等發(fā)展的細(xì)胞存貯稀疏解法在速度上達(dá)到了很高的水平1,2。在實(shí)施稀疏矩陣的數(shù)值組裝之前,必須知道其對(duì)應(yīng)的符號(hào)矩陣 (IA,JA 或等價(jià)地,但規(guī)模小很多的細(xì)胞符號(hào)矩陣 (ISA,JSA。傳統(tǒng)的有限元稀疏符號(hào)矩陣的組裝方法直接采用了基于方程的數(shù)據(jù)結(jié)構(gòu)4,5,符號(hào)矩陣占用的內(nèi)存為數(shù)值剛度矩陣占用內(nèi)存的1/2(4字節(jié)整數(shù), 8字節(jié)實(shí)數(shù)。本文建議采用基于超方程的細(xì)胞數(shù)據(jù)結(jié)構(gòu)實(shí)施

9、符號(hào)矩陣組裝,填充元優(yōu)化和符號(hào)分解。在工程有限元分析中可以從單元?jiǎng)偠染仃嚨奈恢孟蛄客ㄟ^(guò)標(biāo)記單元-方程向量中所有的不連續(xù)方程號(hào)來(lái)做出超方程映射SUPER(1:mneq。在進(jìn)行細(xì)胞符號(hào)矩陣組裝時(shí),還需要SUPER的逆映射MAEQ(1:neq,它將方程號(hào)映射到超方程號(hào)。3.組裝方案1基于方程的行緊湊存貯方案的符號(hào)組裝方案1見(jiàn)文獻(xiàn)5。它的核心是將單元-方程數(shù)組轉(zhuǎn)置成方程-單元數(shù)組,由此可方便地找到與某個(gè)方程相關(guān)的全體方程。此方案要求在內(nèi)存中存放全體單元的單元-方程數(shù)組與方程-單元數(shù)組。組裝表(1的符號(hào)矩陣所需的內(nèi)存neq+nzr+|LM|+|EQ|個(gè)整數(shù),其中|LM|和|EQ|分別為存貯全體單元-方程

10、數(shù)組和方程-單元數(shù)組所需的整數(shù)個(gè)數(shù)。而組裝表(2的細(xì)胞符號(hào)矩陣所需的內(nèi)存 mneq+nzrc+|LMC|+|EQC|個(gè)整數(shù),其中|LMC|和|EQC|分別為存貯全體單元-超方程數(shù)組和超方程-單元數(shù)組所需的整數(shù)個(gè)數(shù)。如果我們要在超方程的意義下進(jìn)行符號(hào)組裝,只需將方程理解為超方程。但在進(jìn)行組裝之前,須將單元-方程數(shù)組變換并壓縮成單元-超方程數(shù)組。此步驟需要的內(nèi)存為neq+|LM|個(gè)整數(shù)。附錄中給出的子程序SYMADD1完成單元-方程數(shù)組到成方程-單元數(shù)組的反轉(zhuǎn)與符號(hào)矩陣組裝兩項(xiàng)功能,可用于通常的行緊湊存貯或超行緊湊存貯。4.組裝方案2除了行緊湊方案,還可以用鏈表存貯方案來(lái)表示式(1。下圖給出了(1

11、式中矩陣第一行的鏈表達(dá)方式。 - 3 -細(xì)胞符號(hào)矩陣的鏈表存貯方案包含一個(gè)指針數(shù)組 ISP(1:nzrc和一個(gè)列指標(biāo)數(shù)組JSB(1:nzrc表示。指針數(shù)組ISP的前mneq個(gè)數(shù)是各超行的指針頭,指針為零則表示這一行的結(jié)束。表3是式(1的細(xì)胞符號(hào)矩陣鏈表存貯方案的一個(gè)例子。表3 細(xì)胞稀疏超行鏈表存貯方案1 2 3 4 5 6 7 8 9 10ISP 10 7 9 8 0 0 0 0 0 6JSB 1 2 3 4 5 5 4 5 4 3如果在任何時(shí)刻只能有一個(gè)“單元-(超方程”數(shù)組在內(nèi)存中,行鏈表存貯方案可以作為符號(hào)組裝時(shí)的中間存貯方案4。這是因?yàn)殒湵泶尜A方案允許方便地實(shí)施插入操作。為了效率的因素

12、,鏈表在算法實(shí)現(xiàn)時(shí)不宜直接用編譯語(yǔ)言中的指針變量,而采用數(shù)組偽指針。附錄中的SYSADD2是鏈表組裝的一種實(shí)現(xiàn)。在有限元分析的數(shù)據(jù)中一般只提供單元-方程數(shù)組,如果要實(shí)現(xiàn)細(xì)胞鏈接存貯方案,可以在子程序SORTLM中嵌入從方程到超方程的映射MAEQ。在通常意義和細(xì)胞意義下實(shí)現(xiàn)行鏈表存貯方案所需的內(nèi)存為2*nzr+|lm|個(gè)整數(shù)和neq+2*nzrc+|lm|個(gè)整數(shù),其中|lm|是一個(gè)單元-方程數(shù)組中的整數(shù)個(gè)數(shù)。5.數(shù)值例題大量的實(shí)際的工程問(wèn)題被用來(lái)檢驗(yàn)本文提出的符號(hào)組裝方案,限于篇幅僅在表4和表5中列出其中的一部分,它們?cè)谝郧暗难芯恐幸讯啻问褂谩Ｋ械臄?shù)值試驗(yàn)都在一臺(tái)帶IDE磁盤,操作系統(tǒng)為中文W

13、indows2000的Pentium IV 850機(jī)器上進(jìn)行,編譯器選用的是Compaq Visual Fortran 6.1a。代碼采用了標(biāo)準(zhǔn)的Fortran 90語(yǔ)言,它已經(jīng)被移植到了Sun Sparc 30和Compaq XP計(jì)算機(jī)平臺(tái)上。試驗(yàn)的目標(biāo)是比較各種組裝方法所占用的內(nèi)存。從表5可以看到,細(xì)胞稀疏存貯方案顯著在減少組裝符號(hào)方程組的存貯量。當(dāng)然數(shù)據(jù)操作量也大為減少。據(jù)測(cè)試,對(duì)于不超過(guò)100萬(wàn)階的三維8節(jié)點(diǎn)四面體的有限元方程組,細(xì)胞符號(hào)矩陣 (ISA,JSA 可以在28 MB內(nèi)存中處理,處理時(shí)間也僅為2分鐘左右。而傳統(tǒng)的符號(hào)矩陣 (IA,JA 的大小約為210MB。處理時(shí)間為16分鐘

14、左右。表4 試驗(yàn)例題以及它們的說(shuō)明工程例題說(shuō)明neq mneq 單元數(shù)nzr nzrc PKUSTK01 北京植物園溫室22,0443,74110,149500,712 15,740 PKUSTK02 飛躍雙塔建筑10,8001,8002,149410,400 12,150 PKUSTK03 大連群倉(cāng)63,33610,55622,1671,596,876 48,756 PKUSTK07 21節(jié)點(diǎn)單元,10x10x10網(wǎng)格16,8605,4001,0001,217,832 127,382 PKUSTK08 21節(jié)點(diǎn)單元,11x11x11網(wǎng)格22,2097,1391,3311,624,440 17

15、1,051- 4 -PKUSTK09 群倉(cāng)33,9605,66014,811808,800 24825 PKUSTK10 4 筒群倉(cāng)80,67613,44614,4002,194,830 66570 PKUSTK11 圍堰87,80414,63435,7502,652,858 79,788 PKUSTK12 吉建大廈94,65322,82211,5923,803,485 310,247 PKUSTK13 機(jī)械零件94,89331,58016,4603,355,860 389,607表5 組裝符號(hào)矩陣所需的內(nèi)存 (KB傳統(tǒng)符號(hào)矩陣細(xì)胞符號(hào)矩陣?yán)}組裝方案1組裝方案2組裝方案1組裝方案2PKUST

16、K01 3194.173912.00388.79209.27PKUSTK02 2032.743206.44168.25137.30PKUSTK03 9720.5712475.781090.48628.50PKUSTK07 5366.229514.50761.411061.22PKUSTK08 7156.5912691.121019.481423.28PKUSTK09 5367.736318.94652.44326.79PKUSTK10 11834.0917147.301165.48835.41PKUSTK11 15530.5820725.641632.27966.52PKUSTK12 1763

17、6.3629714.912414.992793.73PKUSTK13 17734.3726217.843489.233414.67如果采用基于方程的數(shù)據(jù)結(jié)構(gòu),組裝方案2明顯地比組裝方案1占用的存貯量大。在采用了基于超方程的細(xì)胞數(shù)據(jù)結(jié)構(gòu)之后,組裝方案2平均比組裝方案1占用的存貯量小一些。從計(jì)算效率上考慮,緊湊行向量存貯與行鏈接存貯的查尋時(shí)間復(fù)雜度均為O(s,其中s為方程一行中非零元素的平均個(gè)數(shù)4。在采用了細(xì)胞概念后,查尋時(shí)間復(fù)雜度下降為O(s*mneq/neq。與此同時(shí),符號(hào)三角分解的時(shí)間復(fù)雜度之比為1: (mneq/neq3。在數(shù)值試驗(yàn)中發(fā)現(xiàn),對(duì)于高階的符號(hào)矩陣組裝,方案1 所耗的時(shí)間明顯長(zhǎng)于

18、方案2所耗時(shí)間。仔細(xì)分析算法發(fā)現(xiàn),方案1中的“散射-收集”步驟隱含了一個(gè)時(shí)間復(fù)雜度為O(neq或O(mneq的比較運(yùn)算5。而方案2組裝1行的時(shí)間復(fù)雜度分為O(s3或O(s*mneq/neq3。改進(jìn)的辦法是將“散射-收集”替換成附錄3的“比較收集-排序”,它的時(shí)間復(fù)雜度也是O(s3或O(s*mneq/neq3。6.結(jié)語(yǔ)本文討論了細(xì)胞存貯意義下的符號(hào)矩陣以及它的組裝。細(xì)胞存貯方案在工程有限元分析中顯著地減小了符號(hào)矩陣的大小,它使得我們有可能在較小的內(nèi)存中組裝相當(dāng)大規(guī)模的有限元方程的符號(hào)矩陣。理論和數(shù)值試驗(yàn)表明,我們所提出的細(xì)胞意義下組裝方案可以替代傳統(tǒng)的組裝方案。如果需要傳統(tǒng)意義下的符號(hào)矩陣,我們

19、仍然可以在細(xì)胞意義下組裝,而后再行展開(kāi)。組裝-展開(kāi)的思想也可以用于有限元中常用的帶寬與變帶寬存貯方案。我們可以在細(xì)胞- 5 -稀疏意義組裝有限元?jiǎng)偠染仃?然后再轉(zhuǎn)換成帶寬與變帶寬存貯。對(duì)大規(guī)模的計(jì)算,這種方法可以節(jié)約可觀的剛度矩陣組裝時(shí)間。最后我們指出,在細(xì)胞意義下進(jìn)行符號(hào)組裝,數(shù)值組裝,填充元優(yōu)化和矩陣分解等運(yùn)算比在傳統(tǒng)意義下的效率高。致謝作者感謝崔俊芝院士的鼓勵(lì)和鄭東博士的討論。附錄1 稀疏矩陣的行緊湊組裝算法c - SYMADD1subroutine SYMADD1(EQP,LMP,LEN_EQ,LEN_LM,ntt,mneq,IA,JA,maxEdgeinteger(4 mneq,nt

20、t,maxEdgeinteger(4 LEN_EQ,EQP(0:LEN_EQ,LEN_LM,LMP(0:LEN_LMinteger(4 IA(0:mneq, JA(maxEdgeinteger(4 nd_ab, nd_end, el_ab, el_endinteger(4 nnd, n, nd, el, accEdge, jc -c 程序功能: 建立行緊湊方式的上三角行符號(hào)矩陣c 作者:陳璞c 日期: 04/06/98c 參數(shù):c mneq: 方程數(shù)c ntt: 單元數(shù)c IA(: 行緊湊方式的上三角行符號(hào)矩陣的索引c JA(: 行緊湊方式的上三角行符號(hào)矩陣的列指標(biāo)c eqp(: "

21、方程-單元"數(shù)組c len_eq: 全體"方程-單元"數(shù)組長(zhǎng)度之和c lmp: "單元-方程"數(shù)組c len_lm: 全體"單元-方程"數(shù)組長(zhǎng)度之和c maxEdge: 調(diào)用時(shí)為最大的可能列指標(biāo)數(shù)c 返回時(shí)為實(shí)際的列指標(biāo)數(shù)c -c . 統(tǒng)計(jì)各方程相關(guān)的單元數(shù)目 eqp(1:mneqcall zeroi (EQP, mneq+1el_end = nttdo el = 1, nttel_ab = el_endel_end = LMP(eldo j = el_ab+1, el_endEQP(LMP(j = EQP(LMP(j +

22、1- 6 -end doend doc . 將EQP(1:mneq指向每個(gè)方程-單元數(shù)組的起始位置el_ab = mneqdo el = 1, mneqel_end = EQP(elEQP(el = el_abel_ab = el_ab + el_endend doc . 建立方程-單元數(shù)組EQP(0 = mneqel_end = nttdo el = 1, nttel_ab = el_endel_end = LMP(eldo j = el_ab+1, el_endEQP(LMP(j = EQP(LMP(j + 1n = EQP(LMP(jEQP(n = elend doend doc . 形

23、成符號(hào)矩陣call zeroi (IA,mneq+1accEdge = 0nd_end = mneqdo nn = 1, mneq ! 對(duì)每個(gè)方程循環(huán)nd_ab = nd_endnd_end = EQP(nn! 以下部分可用比較收集-排序算法替換do nnd = nd_ab+1, nd_end ! 對(duì)與nn相關(guān)的單元循環(huán) el = EQP(nndel_ab = lmp(el-1el_end = lmp(eldo n = el_ab+1, el_end ! 散射if (lmp(n.ge.nn ia(lmp(n = nnend doend dodo n = nn, mneq ! 收集if (ia(

24、n.eq.nn thenaccEdge = accEdge + 1ja(accEdge = n- 7 -end ifend doia(nn = accEdgeend domaxEdge = accEdgereturnend附錄2 稀疏矩陣的鏈表組裝算法C - SYMADD2 subroutine SYMADD2 (ntape,ntt,ndmax,LM,IP,JB,maxEdge,mneq integer(4 ntape,ntt,ndmax,maxEdge,mneqinteger(4 IP(maxEdge, JB(maxEdge, lm(ndmaxinteger(4 el,nn,ndnz,acc

25、Edge,i,j,ij,ijpc -c 程序功能:建立鏈方式的上三角行符號(hào)矩陣c 作者:陳璞,鄭東c 日期: 12/28/97c 參數(shù):c ntape: 單元-方程數(shù)組的順序文件通道號(hào)c ntt: 單元總數(shù)c ndmax: 單元?jiǎng)偠染仃嚨淖畲笞杂啥葦?shù)c lm( 單元-方程數(shù)組c IP(: 鏈表示法的中的鏈c JB(: 鏈表示法的中的列指標(biāo)c maxEdge: 調(diào)用時(shí)為最大的可能列指標(biāo)數(shù)c 返回時(shí)為實(shí)際的列指標(biāo)數(shù)c mneq: 方程個(gè)數(shù)c 局部參數(shù):c nn 單元?jiǎng)偠染仃嚨淖杂啥葦?shù)c ndnz 單元?jiǎng)偠染仃嚨牟恢貜?fù)非零自由度數(shù)c accEdge: 已累積的列指標(biāo)數(shù)c -c . 初始化call z

26、ero(lkstru, 2*maxEdgeaccEdge = mneqc . 形成總剛矩陣的鏈?zhǔn)奖磉_(dá)do el=1,ntt ! 對(duì)所有單元循環(huán)- 8 -read(ntape lm(1:nn ! 每次讀入一個(gè)單元-方程數(shù)組 call sortlm (lm,nn,ndnz ! 將lm的非零元素排序并壓縮 do i = 1, ndnz-1do j = i+1, nnc . 將(lm(i,lm(j插入鏈表,找到插入的位置ijp = lm(iij = IP(lm(ido while (ij.gt.0 .and. JB(ij.lt.lm(jijp = ijij = IP(ijend doif ( ij .

27、eq. 0 thenc . 在鏈表末插入(lm(i,lm(jaccEdge = accEdge + 1JB(accEdge = lm(jIP(ijp = accEdgeelse if (JB(ij .gt. lm(j thenc . 在鏈表中間插入(lm(i,lm(jaccEdge = accEdge + 1JB(accEdge = lm(jIP(ijp = accEdgeIP(accEdge = ijendifend doend do ! do i = n, ndnz-1end do ! do el = 1, nttmaxEdge = accEdgereturnend附錄3 比較收集-排序算

28、法accEdgeLast = accEdgedo nnd = nd_ab+1, nd_end ! 對(duì)與nn相關(guān)的單元循環(huán) el = EQP(nndel_ab = lmp(el-1el_end = lmp(eldo n = el_ab+1, el_end ! 比較收集if (lmp(n.lt.nn goto 100do j = accEdgeLast+1, accEdgeif (lmp(n.eq.ja(j goto 100end do- 9 -accEdge = accEdge + 1ja(accEdge = lmp(n100 continueend doend docall sortkm (j

29、a(accEdgelast+1,accEdge-accEdgeLast !排序ia(nn = accEdge參考文獻(xiàn)1陳璞, 孫樹(shù)立, 袁明武. 有限元快速解法. 力學(xué)學(xué)報(bào),2002, 34(2: 216-222.2Chen P, Zheng D, Sun SL, Yuan MW. High Performance Sparse Static Solver in Finite Element Analysiswith Loop-unrolling. Advances in Engineering Software, 2003, 34: 203-215.3Damhaug AC, Reid J, Bergseth A. The Impact of an Efficient Linear Solver on Finite Element Analysis.Computers & Structures, 1999; 72: 594-604.4謝柏青,佘曉歌. 算法與數(shù)據(jù)結(jié)構(gòu). 高等教育出版社, 北京, 2001.5Pissanetzky S. Spar