《带有双重循环展开的非对称稀疏线性方程组快速直接解法》由会员分享，可在线阅读，更多相关《带有双重循环展开的非对称稀疏线性方程组快速直接解法（9页珍藏版）》请在金锄头文库上搜索。
1、带有双重循环展开的非对称稀疏线性方程组快速直接解法苑维然1 ，陈璞1 ，刘凯欣k2( 1 北京大学工学院北京；2 北京大学工程研究院北京)警要：本文介绍了一种新的直接解法来求解科学与工程计算中生成的大型非对称稀疏线性方程组。该解法从现有的对称解法中演变出来，其分解过程在矩阵的上、下三角阵中对称行进。该解法中的L D U分解算法利用了双重循环展开技术，并且由于其对称行进的求解方式，可以通过修改已有的对称矩阵分解算法的代码来实现，这提供了从对称解法到非对称解法的快捷转换。在数值测试中求解了若干无网格局部P c 订o v G a l e r k i n 法( M L P G ) 法生成的矩阵，结果表
2、明本文的方法可以大幅度提高了大型非对称稀疏线性方程组的求解速度。关键词：稀疏矩阵线性方程组无网格法高性能计算1 引言在上个世纪，随着有限元法的发展中涌现了大量有效的对称稀疏矩阵的解法，它们适用于在不同计算机环境下进行有限元方程组的求解【l 圳。在计算力学领域，近年来发展了无网格法等新方法，与有限元法的方程组相比，以无网格局部P e 仃o v G a l e r k i n 法( M L P G ) 法1 1 5 , 1 6 1 为代表的无网格法方程组有以下特点：( 1 ) 没有单元刚度矩阵的概念；( 2 ) 系数矩阵是非对称的，但它的非零元分布几乎是对称的：( 3 ) 各阶顺序主子式都是正的。
3、本文利用有限元解法的研究成果来发展这一类非对称稀疏矩阵的高效求解器。在科学与工程计算中，M L P G 法产生的线性方程组的未知量数目随着研究的深入不断变大。文献中解中小例题的方法，例如：满矩阵或带宽矩阵的三角分解已不能适应求解大型应用问题的需求。这一因难主要反映在求解时间和系数矩阵的大小上。快速解法是指非并行计算机上比传统的带宽解法快数倍至数十倍的线性代数方程组的解法。尽管计算机的速度、内存、外存容量等在提高，但随着工程实际问题复杂程度的增加和分析要求的提高，计算机性能的提高并不能完全满足大规模计算的需要。另一方面，更快更节省存贮空间的算法也一直是计算科学所追求的目标之一。快速解法是计算力学
4、的核心技术之一，也是计算力学和计算数学领域中的一个非常重要的课题。本文研究并实现了一个直接非对称矩阵求解器，它能够快速求解大型非对称稀疏线性方程组。在求解的过程中，系数矩阵的三角分解是求解稀疏线性方程组最重要的一步，本文将着重讨论三角分解的策略。在二十世纪7 0 年代，不少有限元软件采用了变带宽存储的解法。最近十多年，总体刚度矩阵的存贮逐渐变为稀疏存贮法。国际上一些著名的有限元分析软件已经将稀疏存储的三角分解作为默认求解方案。与传统的方案比较，稀疏存储解法仅仅需要为刚度矩阵及其因子的非零元分配存储空间，矩阵的运算都在索引形式下进行，只有非零元才被处理。稀疏存储的解法比传统的变带宽存储更快速，需
5、要更少的磁盘空间。当然，稀疏解法的编程远比传统的解法困难。考虑非对称的线性方程 A x = b( 1 )其中，矩阵A = ( 珥，) ，a i ，R 是一个大型稀疏矩阵，向量x 和b 分别是待求的未知向量和已知的国家自然科学基金资助项目( 1 0 2 3 2 0 4 0 ，1 0 5 7 2 0 0 2 )1 2 7输入向量。在我们的问题中，系数矩阵A 在数值上是非对称的，但是A 的上三角和下三角部分在非零元的分布上具有几乎对称的特征。对这样的方程组，可以在矩阵A 中少量结构非对称的位置填入零，使矩阵A 在结构上对称，然后利用对称矩阵已有的成熟方法进行计算。 M L P G 法生成的系数矩阵A
6、 是各阶顺序主子式全为正的，这种特性使矩阵A 存在唯一的三角分解如下： A = L D U( 2 )其中，L 和U 分别是单位下三角和单位上三角矩阵，D 是对角矩阵。在完成三角分解之后，式( 1 ) 的求解可以分为三个步骤：L y = b ，D z = y ，U x = z 嘲。当A 为对称矩阵，U = o 。没有任何疑问，L 和U 的计算是求解线性方程组( 1 ) 中计算量最高的步骤。如果矩阵A 是稀疏的，则一般L 和U 也是稀疏的，但是呸，0 中的( f ，歹) 集是非零元结构L + D + U 的子集。2 存储方案在数值计算中，稀疏矩阵的存储一般可以选用行主元存储或列主元存储。在行主元存
7、储中，矩阵被视作紧凑行向量的顺序组合。在本文中，系数矩阵使用一种混合的存储方案。在式( 3 ) 中，a ，b ，f是上三角部分的非零元，口，矿，f + 是和上三角部分具有相同稀疏结构的下三角部分的非零元。表1中的数组I A ，J C A ，P A 和Q A 为矩阵A 的行列主元存储方案，其中上三角部分在行主元形式下需要 3 个数组I A ，J C A 和P A 来表示，而它的下三角部分可以在列主元形式下用I A ，J C A 和Q A 表示。数组从是非零元的列，行指标数组J C A 以及它的数值数组P A Q A 的行列索引。第k 行的非零元个数为认国一I A ( k - 1 ) ( 设认(
8、0 ) = 0 ) 。一般称( I A ，J C A ) 为符号矩阵，( P A ，Q A ) 为数值矩阵。A =2 2CC + 3 3d f 口egb abdefg4 4hh 5 5il 6 6( 3 )在系数矩阵A 中，一个节点相关的方程经常具有连续的编号，这样在A 中相应的行和列具有相 同的稀疏结构。这样一组连续的方程定义为A 的超方程。超方程的第l 行称为主方程( m a s t e r - e q u a t i o n ) ，其他行称为从方程( s l a v e - e q u a t i o n ) 。由于超方程技术在有限元求解中成功的应用【5 1 ，很自然的考虑到应用超方程来
9、提高非对称情形下解法的性能。超方程的原理可以结合在传统的稀疏存储策略中，只需添加一个整型数组M A S T E R A ( n e q ) 。值为正的M A S T E R A ( O 表示从第f 个方程开始的超方程的大小，而值为0 的M A S T E R A ( 0 表示第f 个方程是一个从方程。例如，根据超方程的定义，式( 3 ) 中矩阵A 有4 个超方程，分别为 1 ) ， 2 ，3 ，4 ) ， 5 ) 和 6 ) 。这种表示为超方程的稀疏矩阵A 可以用S h e r m a n 压缩存储方案t g l 来减少传统存储方案中列，行指标数组J C A 的存储量。修正后的列，行指标数组为
10、J A 。为了使用J A ，还需要一个辅助数组K A ，来指出每一列行的指标修正在J A 中的开始位置。稀疏矩阵A 的存储方案如表1 ，其中n e q 是系数矩阵A 的阶数，咖是压缩的J A 的个数，r l g r 是刚度矩阵上三角部分非零元个数，即矩阵A 中总的非零元个数为，2 z 芦2 。通常n j a n z r 。表l 矩阵存储方案E q u a t i o nl23456I A ( n e q )371 01 21 41 5J C A ( n z r )156234534545566P A ( n z r )1 1ab2 2Cde3 3厂g4 4办5 5f6 6Q A ( n z r
11、 )a b c + d e f g h ZK A ( n e q )l45681 0J A ( n j a )1562345566M A S T E R A ( n e q )l300ll对于A 的因子矩阵L 和U ，可以采用相同的数据结构I u 、J U 、M A S T E R U 、P u 和Q u 来表示，但要注意它们的非零元个数与A 的上下三角部分不同，这是因为在矩阵的三角分解过程中将产生许多新的非零元，一般称之为填充元。与A 的超方程定义不同的是，填充元的产生将使得U 与L 中的不少相邻的行列的对角块是满上厂F 三角矩阵，而其余的部分具有完全相同的列，行指标集，或者说这些行列具有相
12、同的稀疏性。这些的相邻行列所对应的方程被定义为U L 的超方程l l 叭。相应的数据结构是l U ，J U ，M A S T E R U ，P U 和Q u 。数值部分P A 和P U ( 或Q A 和Q U ) 的大小不能压缩，但是符号部分n ，J A 和I U ，J U 用S h e r m a n 压缩方法表示后，减少了所需的总内存空间。需要注意的是，U ，L 的超方程定义与A 的超方程定义不同。如果没有外存储策略，大型问题的一般求解方法是依靠虚拟内存系统，允许操作系统在内存和磁盘之间移动数据。这种方法的优点是不需要修改内存算法的程序，但是在大型题目中应用的经验表明它的速度很慢，所以本文
13、针对非对称矩阵的特殊性【1 0 。1 2 l 扩展了文献中针对有限元计算的外存算法。将矩阵A 和分解因子上、下三角阵U 与L 划分为相似大小的块。根据消元的要求，按照至少有两块可以同时安置于内存来决定块的大小。而在本文解法的P C 机版本实现中，块的大小要满足5 块以上可以同时放在内存中。3 稀琉L D U 分解当矩阵A 是各阶主子式大于零时，求解器一般可分为几个独立的阶段：符号刚度矩阵的组装、填充元优化、符号分解、数值刚度矩阵的组装、数值分解以及消元与回代。本文解法和有限元线性方程组解法在L D U 分解过程中的相应阶段【1o 】有很大的不同。其中，使用混合行主元和列主元的存储方法，使符号分
14、解只需要在矩阵上、下三角阵中的个进行。对称的非零元位置使我们可以使用为对称方程组开发的优化算法，如A M D 和M E T I S 。所以本文的非对称解法在优化和符号分解的阶段是和对称情形基本相同的。本文解法在对称情形下做尽量少的更改，以期提供一个广泛适用的从对称求解器过渡；至U t I E 对称求解器的解决方案。其中，符号组装过程要考虑到不同算法输出格式与有限元单刚矩阵的不同；数值组装过程要考虑矩阵上、下两部分的数值数组P A 和Q A ，所以数值部分所需的内存空间是对称情形的近两倍；数值分解过程则要考虑上、下三角相互分解；消元与回代过程则分别使用上、下三角，而对称情形时只用下三角。3 1
15、基本的L D U 分解由于使用混合行主元和列主元的存储方法，本文的L D U 分解为l e f t - l o o k i n g 的形式( 或者称为b a c k w a r d - r e f e r e n c e ) z 3 I 。算法1 描述了矩阵A 的l e f t l o o k i n g 形式的L D U 分解算法。对于下三角部分，子过程R o w C o l u m n T a s k ( k j ) 在L 和U 中对哆异消元，即目标第，列行用上下三角阵U L 的第k1 2 9行列的倍数来消元。子过程R o w C o l u m n T a s k j ) 把目标歹行列非
16、对角线元素除以该行列的对角线元素。和对称的情形类似，非对称情形的数值分解也需要跟踪对当前分解有贡献的行，列号。文献1 0 1 提出了对称情形的直接更新法( d i r e c tu p d a t i n gm e t h o d ) 来替代传统的分散一收集( s c a t t e r - g a t h e r ) 操作。在算法2 中，将这种方法应用到非对称的情形，它用前面的行列来更新第，行列的数值 P C U ( I ) ，Q C U ( I ) II C U ( j 一1 ) ，I C U ( j ) 。在这种情况下，U P D 口：n e q ) 作为从射行歹0 的列号绗号到它们的P C U Q C U 位置的映射。这种方法的另一个优点是，当分解过程是以L 和U 两部分对 称行进的时候，U P D 对于L 和U 是相同的。3 2 带有循环展开技术的L D U 分解优异的数值方法能在计算机上大大地增加求解效率，但它不是影响计算效率的唯一因
