基于細胞稀疏存貯方案的有限元剛度矩陣組裝

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

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

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

4、體有限元方程組的符號矩陣。2.稀疏矩陣的存貯方案為了更清楚地討論問題,取一個6階的總體剛度矩陣如下:1 高等學(xué)校博士學(xué)科點專項科研基金資助項目(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ù)。由于對稱性,有限元剛度矩陣采用了上三角行存貯方法(或等價地,下三角列存貯方法。為了有效地對行實施運算,在有限元分析的計算中矩陣A 較常用的稀疏存貯方案是緊湊行向量的順序組合表示法(又稱為帶輔助行向量的二元組表示法。這種表示法需要3個數(shù)組

5、IA(1:neq, JA(1:nzr 和 PA(1:nzr 來表示矩陣 A 的上三角部分4,5,其中neq 和nzr 分別是線性方程組的階數(shù)和是非零元的個數(shù)。以4字節(jié)整數(shù)和8字節(jié)實數(shù)計,這種存貯方案占用的內(nèi)存為 (neq+nzr*4 + nzr*8 字節(jié)。表1 緊湊行向量的順序組合表示法方程號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 行的非零元個數(shù)為IA(k-IA(k-1(假設(shè)IA(0=0。文1,2提出了一種細胞存貯方案。它的想法是將剛度矩陣中具有相同的非零元位置的行(列所對應(yīng)的方程定義為超方程,與超方程劃分相應(yīng)的子矩陣為細胞。式(1中A 的第2個表達式就是按細胞方式寫的。顯然,符號矩陣可以在細胞意義表達,在給定超方程劃分的意義下,它與原始的表達是完全等價的。在這一表示法之下,有限元的總體剛度矩陣被視作超行向量的順序組合。根據(jù)超方程的定義,式(1中有5個超方程。記mneq 為超方程數(shù),在超行緊湊方式下,式(1中的第2個矩陣可以用5個數(shù)組來描述。表2 緊湊行向量的細胞稀疏存貯方案超方程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是超方程編號到原方程編號的映射。數(shù)組 ISA 是非零細胞的超列指標(biāo)(或稱細胞指標(biāo) JSA 的索引,LSA是數(shù)值數(shù)組 PSA 以超方程為單位的索引。這種存貯方案占用的內(nèi)存為 (mneq*3+nzrc*4 + nzr*8 字節(jié)。一般的稀疏存貯方案與細胞存貯方案之間的差異在于列指標(biāo)的數(shù)量。細胞存貯方案的指標(biāo)比傳統(tǒng)存貯方案大大減少。對于三維實體分析,細胞指標(biāo)數(shù)nzrc僅僅

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

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

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

11、式中矩陣第一行的鏈表達方式。 - 3 -細胞符號矩陣的鏈表存貯方案包含一個指針數(shù)組 ISP(1:nzrc和一個列指標(biāo)數(shù)組JSB(1:nzrc表示。指針數(shù)組ISP的前mneq個數(shù)是各超行的指針頭,指針為零則表示這一行的結(jié)束。表3是式(1的細胞符號矩陣鏈表存貯方案的一個例子。表3 細胞稀疏超行鏈表存貯方案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ù)組在內(nèi)存中,行鏈表存貯方案可以作為符號組裝時的中間存貯方案4。這是因為鏈表存貯方案允許方便地實施插入操作。為了效率的因素

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

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

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

15、1,051- 4 -PKUSTK09 群倉33,9605,66014,811808,800 24825 PKUSTK10 4 筒群倉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 機械零件94,89331,58016,4603,355,860 389,607表5 組裝符號矩陣所需的內(nèi)存 (KB傳統(tǒng)符號矩陣細胞符號矩陣例題組裝方案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占用的存貯量大。在采用了基于超方程的細胞數(shù)據(jù)結(jié)構(gòu)之后,組裝方案2平均比組裝方案1占用的存貯量小一些。從計算效率上考慮,緊湊行向量存貯與行鏈接存貯的查尋時間復(fù)雜度均為O(s,其中s為方程一行中非零元素的平均個數(shù)4。在采用了細胞概念后,查尋時間復(fù)雜度下降為O(s*mneq/neq。與此同時,符號三角分解的時間復(fù)雜度之比為1: (mneq/neq3。在數(shù)值試驗中發(fā)現(xiàn),對于高階的符號矩陣組裝,方案1 所耗的時間明顯長于

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

19、仍然可以在細胞意義下組裝,而后再行展開。組裝-展開的思想也可以用于有限元中常用的帶寬與變帶寬存貯方案。我們可以在細胞- 5 -稀疏意義組裝有限元剛度矩陣,然后再轉(zhuǎn)換成帶寬與變帶寬存貯。對大規(guī)模的計算,這種方法可以節(jié)約可觀的剛度矩陣組裝時間。最后我們指出,在細胞意義下進行符號組裝,數(shù)值組裝,填充元優(yōu)化和矩陣分解等運算比在傳統(tǒng)意義下的效率高。致謝作者感謝崔俊芝院士的鼓勵和鄭東博士的討論。附錄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 程序功能: 建立行緊湊方式的上三角行符號矩陣c 作者:陳璞c 日期: 04/06/98c 參數(shù):c mneq: 方程數(shù)c ntt: 單元數(shù)c IA(: 行緊湊方式的上三角行符號矩陣的索引c JA(: 行緊湊方式的上三角行符號矩陣的列指標(biāo)c eqp(: "

21、方程-單元"數(shù)組c len_eq: 全體"方程-單元"數(shù)組長度之和c lmp: "單元-方程"數(shù)組c len_lm: 全體"單元-方程"數(shù)組長度之和c maxEdge: 調(diào)用時為最大的可能列指標(biāo)數(shù)c 返回時為實際的列指標(biāo)數(shù)c -c . 統(tǒng)計各方程相關(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指向每個方程-單元數(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、成符號矩陣call zeroi (IA,mneq+1accEdge = 0nd_end = mneqdo nn = 1, mneq ! 對每個方程循環(huán)nd_ab = nd_endnd_end = EQP(nn! 以下部分可用比較收集-排序算法替換do nnd = nd_ab+1, nd_end ! 對與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 程序功能:建立鏈方式的上三角行符號矩陣c 作者:陳璞,鄭東c 日期: 12/28/97c 參數(shù):c ntape: 單元-方程數(shù)組的順序文件通道號c ntt: 單元總數(shù)c ndmax: 單元剛度矩陣的最大自由度數(shù)c lm( 單元-方程數(shù)組c IP(: 鏈表示法的中的鏈c JB(: 鏈表示法的中的列指標(biāo)c maxEdge: 調(diào)用時為最大的可能列指標(biāo)數(shù)c 返回時為實際的列指標(biāo)數(shù)c mneq: 方程個數(shù)c 局部參數(shù):c nn 單元剛度矩陣的自由度數(shù)c ndnz 單元剛度矩陣的不重復(fù)非零自由度數(shù)c accEdge: 已累積的列指標(biāo)數(shù)c -c . 初始化call z

26、ero(lkstru, 2*maxEdgeaccEdge = mneqc . 形成總剛矩陣的鏈?zhǔn)奖磉_do el=1,ntt ! 對所有單元循環(huán)- 8 -read(ntape lm(1:nn ! 每次讀入一個單元-方程數(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 ! 對與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參考文獻1陳璞, 孫樹立, 袁明武. 有限元快速解法. 力學(xué)學(xué)報,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