基于细胞稀疏存贮方案的有限元刚度矩阵组装

1、基于细胞稀疏存贮方案的有限元刚度矩阵组装1陈璞1,陈斌11北京大学力学与工程科学系,北京 (100871E-mail:chenpu摘要:在现代有限元分析的直接解法和迭代解法中都需要组装刚度矩阵。本文在1,2的基础上,叙述了两种基于细胞稀疏存贮方案有限元的符号刚度矩阵的组装方法。在与传统的方法进行了比较的基础上,提出了改进效率的方法。关键词:有限元分析稀疏矩阵高性能计算中图分类号:O241.821.引言近年来,稀疏解法逐渐替代了带宽解法与变带宽解法,成为了有限元分析中首选的直接解法1,2,3。不同于带宽解法或变带宽解存贮方案,稀疏解法要求建立符号矩阵(symbolic matrix,它是总体刚度

2、矩阵对应的Bool值矩阵。刚度矩阵中的非零元在Bool值矩阵为1,其余则为0。在用图讨论矩阵运算时符号矩阵又称为邻接图(adjacent graph。在稀疏解法中符号矩阵涉及几乎所有的步骤,如:符号组装;填充元优化;符号分解;(数值刚度矩阵的组装;数值分解;消元与回带等。因此,减小符号矩阵的大小在稀疏解法中具有十分重要的意义。传统的稀疏符号矩阵的组装方法直接采用了基于方程的数据结构5,符号矩阵占用的内存为数值刚度矩阵占用内存的1/2(4字节整数,8字节实数。对于给定的单元类型,符号矩阵的大小一般正比于节点数。如果分析的问题不太大,符号矩阵可以整个在内存中处理。但模型较大时有限元的符号矩阵的可能

3、超出内存允许的范围。于是就出现了各种符号矩阵的压缩存贮方法,其中最著名的是文献6的方法。在工程有限元分析中,因为求解位移向量场的需要,总体刚度矩阵是由一系列小的子矩阵组成。利用这一性质,文1,2给出了一个称为细胞稀疏存贮方案的符号矩阵压缩方法,其压缩比大于Sherman的方法6。本文从两个方面改进符号矩阵的组装方法:1利用文1,2的超方程和细胞符号矩阵概念直接组装细胞符号矩阵。在工程有限元分析中,它的大小仅是传统符号矩阵的1/36到1/9;2给出了两种不同的组装方案。它们分别适用于不同的内存数据安排。数值试验表明,本文提出的方法可以在仅有28M内存的情形下在内存快速地组装100万阶3维8节点实

4、体有限元方程组的符号矩阵。2.稀疏矩阵的存贮方案为了更清楚地讨论问题,取一个6阶的总体刚度矩阵如下:1 高等学校博士学科点专项科研基金资助项目(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 是非零实数。由于对称性,有限元刚度矩阵采用了上三角行存贮方法(或等价地,下三角列存贮方法。为了有效地对行实施运算,在有限元分析的计算中矩阵A 较常用的稀疏存贮方案是紧凑行向量的顺序组合表示法(又称为带辅助行向量的二元组表示法。这种表示法需要3个数组

5、IA(1:neq, JA(1:nzr 和 PA(1:nzr 来表示矩阵 A 的上三角部分4,5,其中neq 和nzr 分别是线性方程组的阶数和是非零元的个数。以4字节整数和8字节实数计,这种存贮方案占用的内存为 (neq+nzr*4 + nzr*8 字节。表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数组IA(1:neq=6是行索引,它记录了每一紧湊行向量在非零元的列指标数组JA(1:nzr=14 以及它的数值数组 PA(1:nzr=14 中的末

6、地址。第k 行的非零元个数为IA(k-IA(k-1(假设IA(0=0。文1,2提出了一种细胞存贮方案。它的想法是将刚度矩阵中具有相同的非零元位置的行(列所对应的方程定义为超方程,与超方程划分相应的子矩阵为细胞。式(1中A 的第2个表达式就是按细胞方式写的。显然,符号矩阵可以在细胞意义表达,在给定超方程划分的意义下,它与原始的表达是完全等价的。在这一表示法之下,有限元的总体刚度矩阵被视作超行向量的顺序组合。根据超方程的定义,式(1中有5个超方程。记mneq 为超方程数,在超行紧凑方式下,式(1中的第2个矩阵可以用5个数组来描述。表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数组 SUPER (1:mneq=5是超方程编号到原方程编号的映射。数组 ISA 是非零细胞的超列指标(或称细胞指标 JSA 的索引,LSA是数值数组 PSA 以超方程为单位的索引。这种存贮方案占用的内存为 (mneq*3+nzrc*4 + nzr*8 字节。一般的稀疏存贮方案与细胞存贮方案之间的差异在于列指标的数量。细胞存贮方案的指标比传统存贮方案大大减少。对于三维实体分析,细胞指标数nzrc仅仅

8、是传统方案列指标数nzr的约1/9。对于三维壳体分析,则指标数相对更少。使用细胞存贮方案的好处不仅在于指标集合的大幅度减小,而且可以大幅度地减少矩阵运算中的指标操作,从而达到提高符号运算和数值运算速度的目的。结合循环展开技术8,笔者等发展的细胞存贮稀疏解法在速度上达到了很高的水平1,2。在实施稀疏矩阵的数值组装之前,必须知道其对应的符号矩阵 (IA,JA 或等价地,但规模小很多的细胞符号矩阵 (ISA,JSA。传统的有限元稀疏符号矩阵的组装方法直接采用了基于方程的数据结构4,5,符号矩阵占用的内存为数值刚度矩阵占用内存的1/2(4字节整数, 8字节实数。本文建议采用基于超方程的细胞数据结构实施

9、符号矩阵组装,填充元优化和符号分解。在工程有限元分析中可以从单元刚度矩阵的位置向量通过标记单元-方程向量中所有的不连续方程号来做出超方程映射SUPER(1:mneq。在进行细胞符号矩阵组装时,还需要SUPER的逆映射MAEQ(1:neq,它将方程号映射到超方程号。3.组装方案1基于方程的行紧湊存贮方案的符号组装方案1见文献5。它的核心是将单元-方程数组转置成方程-单元数组,由此可方便地找到与某个方程相关的全体方程。此方案要求在内存中存放全体单元的单元-方程数组与方程-单元数组。组装表(1的符号矩阵所需的内存neq+nzr+|LM|+|EQ|个整数,其中|LM|和|EQ|分别为存贮全体单元-方程

10、数组和方程-单元数组所需的整数个数。而组装表(2的细胞符号矩阵所需的内存 mneq+nzrc+|LMC|+|EQC|个整数,其中|LMC|和|EQC|分别为存贮全体单元-超方程数组和超方程-单元数组所需的整数个数。如果我们要在超方程的意义下进行符号组装,只需将方程理解为超方程。但在进行组装之前,须将单元-方程数组变换并压缩成单元-超方程数组。此步骤需要的内存为neq+|LM|个整数。附录中给出的子程序SYMADD1完成单元-方程数组到成方程-单元数组的反转与符号矩阵组装两项功能,可用于通常的行紧凑存贮或超行紧凑存贮。4.组装方案2除了行紧湊方案,还可以用链表存贮方案来表示式(1。下图给出了(1

11、式中矩阵第一行的链表达方式。 - 3 -细胞符号矩阵的链表存贮方案包含一个指针数组 ISP(1:nzrc和一个列指标数组JSB(1:nzrc表示。指针数组ISP的前mneq个数是各超行的指针头,指针为零则表示这一行的结束。表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如果在任何时刻只能有一个“单元-(超方程”数组在内存中,行链表存贮方案可以作为符号组装时的中间存贮方案4。这是因为链表存贮方案允许方便地实施插入操作。为了效率的因素

12、,链表在算法实现时不宜直接用编译语言中的指针变量,而采用数组伪指针。附录中的SYSADD2是链表组装的一种实现。在有限元分析的数据中一般只提供单元-方程数组,如果要实现细胞链接存贮方案,可以在子程序SORTLM中嵌入从方程到超方程的映射MAEQ。在通常意义和细胞意义下实现行链表存贮方案所需的内存为2*nzr+|lm|个整数和neq+2*nzrc+|lm|个整数,其中|lm|是一个单元-方程数组中的整数个数。5.数值例题大量的实际的工程问题被用来检验本文提出的符号组装方案,限于篇幅仅在表4和表5中列出其中的一部分,它们在以前的研究中已多次使用。所有的数值试验都在一台带IDE磁盘,操作系统为中文W

13、indows2000的Pentium IV 850机器上进行,编译器选用的是Compaq Visual Fortran 6.1a。代码采用了标准的Fortran 90语言,它已经被移植到了Sun Sparc 30和Compaq XP计算机平台上。试验的目标是比较各种组装方法所占用的内存。从表5可以看到,细胞稀疏存贮方案显著在减少组装符号方程组的存贮量。当然数据操作量也大为减少。据测试,对于不超过100万阶的三维8节点四面体的有限元方程组,细胞符号矩阵 (ISA,JSA 可以在28 MB内存中处理,处理时间也仅为2分钟左右。而传统的符号矩阵 (IA,JA 的大小约为210MB。处理时间为16分钟

14、左右。表4 试验例题以及它们的说明工程例题说明neq mneq 单元数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节点单元,10x10x10网格16,8605,4001,0001,217,832 127,382 PKUSTK08 21节点单元,11x11x11网格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 组装符号矩阵所需的内存 (KB传统符号矩阵细胞符号矩阵例题组装方案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如果采用基于方程的数据结构,组装方案2明显地比组装方案1占用的存贮量大。在采用了基于超方程的细胞数据结构之后,组装方案2平均比组装方案1占用的存贮量小一些。从计算效率上考虑,紧凑行向量存贮与行链接存贮的查寻时间复杂度均为O(s,其中s为方程一行中非零元素的平均个数4。在采用了细胞概念后,查寻时间复杂度下降为O(s*mneq/neq。与此同时,符号三角分解的时间复杂度之比为1: (mneq/neq3。在数值试验中发现,对于高阶的符号矩阵组装,方案1 所耗的时间明显长于

18、方案2所耗时间。仔细分析算法发现,方案1中的“散射-收集”步骤隐含了一个时间复杂度为O(neq或O(mneq的比较运算5。而方案2组装1行的时间复杂度分为O(s3或O(s*mneq/neq3。改进的办法是将“散射-收集”替换成附录3的“比较收集-排序”,它的时间复杂度也是O(s3或O(s*mneq/neq3。6.结语本文讨论了细胞存贮意义下的符号矩阵以及它的组装。细胞存贮方案在工程有限元分析中显著地减小了符号矩阵的大小,它使得我们有可能在较小的内存中组装相当大规模的有限元方程的符号矩阵。理论和数值试验表明,我们所提出的细胞意义下组装方案可以替代传统的组装方案。如果需要传统意义下的符号矩阵,我们

19、仍然可以在细胞意义下组装,而后再行展开。组装-展开的思想也可以用于有限元中常用的带宽与变带宽存贮方案。我们可以在细胞- 5 -稀疏意义组装有限元刚度矩阵,然后再转换成带宽与变带宽存贮。对大规模的计算,这种方法可以节约可观的刚度矩阵组装时间。最后我们指出,在细胞意义下进行符号组装,数值组装,填充元优化和矩阵分解等运算比在传统意义下的效率高。致谢作者感谢崔俊芝院士的鼓励和郑东博士的讨论。附录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 参数:c mneq: 方程数c ntt: 单元数c IA(: 行紧凑方式的上三角行符号矩阵的索引c JA(: 行紧凑方式的上三角行符号矩阵的列指标c eqp(: "

21、方程-单元"数组c len_eq: 全体"方程-单元"数组长度之和c lmp: "单元-方程"数组c len_lm: 全体"单元-方程"数组长度之和c maxEdge: 调用时为最大的可能列指标数c 返回时为实际的列指标数c -c . 统计各方程相关的单元数目 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指向每个方程-单元数组的起始位置el_ab = mneqdo el = 1, mneqel_end = EQP(elEQP(el = el_abel_ab = el_ab + el_endend doc . 建立方程-单元数组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 ! 对每个方程循环nd_ab = nd_endnd_end = EQP(nn! 以下部分可用比较收集-排序算法替换do nnd = nd_ab+1, nd_end ! 对与nn相关的单元循环 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 参数:c ntape: 单元-方程数组的顺序文件通道号c ntt: 单元总数c ndmax: 单元刚度矩阵的最大自由度数c lm( 单元-方程数组c IP(: 链表示法的中的链c JB(: 链表示法的中的列指标c maxEdge: 调用时为最大的可能列指标数c 返回时为实际的列指标数c mneq: 方程个数c 局部参数:c nn 单元刚度矩阵的自由度数c ndnz 单元刚度矩阵的不重复非零自由度数c accEdge: 已累积的列指标数c -c . 初始化call z

26、ero(lkstru, 2*maxEdgeaccEdge = mneqc . 形成总刚矩阵的链式表达do el=1,ntt ! 对所有单元循环- 8 -read(ntape lm(1:nn ! 每次读入一个单元-方程数组 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相关的单元循环 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陈璞, 孙树立, 袁明武. 有限元快速解法. 力学学报,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谢柏青,佘晓歌. 算法与数据结构. 高等教育出版社, 北京, 2001.5Pissanetzky S. Spar