带有双重循环展开的非对称稀疏线性方程组快速直接解法.pdf

带有双重循环展开的非对称稀疏线性方程组快 速直接解法· 2 苑维然1，陈璞1，刘凯欣k (1北京大学工学院北京；2北京大学工程研究院北京) 警要：本文介绍了一种新的直接解法来求解科学与工程计算中生成的大型非对称稀疏线性方程组。该 解法从现有的对称解法中演变出来，其分解过程在矩阵的上、下三角阵中对称行进。该解法中的LDU 分解算法利用了双重循环展开技术，并且由于其对称行进的求解方式，可以通过修改已有的对称矩阵分 解算法的代码来实现，这提供了从对称解法到非对称解法的快捷转换。在数值测试中求解了若干无网格 局部Pc订ov—Galerkin法(MLPG)法生成的矩阵，结果表明本文的方法可以大幅度提高了大型非对称稀疏线 性方程组的求解速度。 关键词：稀疏矩阵线性方程组无网格法高性能计算 1引言 在上个世纪，随着有限元法的发展中涌现了大量有效的对称稀疏矩阵的解法，它们适用于在不同 计算机环境下进行有限元方程组的求解【l圳。在计算力学领域，近年来发展了无网格法等新方法，与 以下特点：(1)没有单元刚度矩阵的概念；(2)系数矩阵是非对称的，但它的非零元分布几乎是对称 的：(3)各阶顺序主子式都是正的。本文利用有限元解法的研究成果来发展这一类非对称稀疏矩阵的 高效求解器。 在科学与工程计算中，MLPG法产生的线性方程组的未知量数目随着研究的深入不断变大。文献 中解中小例题的方法，例如：满矩阵或带宽矩阵的三角分解已不能适应求解大型应用问题的需求。这 一因难主要反映在求解时间和系数矩阵的大小上。快速解法是指非并行计算机上比传统的带宽解法快 数倍至数十倍的线性代数方程组的解法。尽管计算机的速度、内存、外存容量等在提高，但随着工程 实际问题复杂程度的增加和分析要求的提高，计算机性能的提高并不能完全满足大规模计算的需要。 另一方面，更快更节省存贮空间的算法也一直是计算科学所追求的目标之一。快速解法是计算力学的 核心技术之一，也是计算力学和计算数学领域中的一个非常重要的课题。 本文研究并实现了一个直接非对称矩阵求解器，它能够快速求解大型非对称稀疏线性方程组。在 求解的过程中，系数矩阵的三角分解是求解稀疏线性方程组最重要的一步，本文将着重讨论三角分解 的策略。在二十世纪70年代，不少有限元软件采用了变带宽存储的解法。最近十多年，总体刚度矩 阵的存贮逐渐变为稀疏存贮法。国际上一些著名的有限元分析软件已经将稀疏存储的三角分解作为默 认求解方案。与传统的方案比较，稀疏存储解法仅仅需要为刚度矩阵及其因子的非零元分配存储空间， 矩阵的运算都在索引形式下进行，只有非零元才被处理。稀疏存储的解法比传统的变带宽存储更快速， 需要更少的磁盘空间。当然，稀疏解法的编程远比传统的解法困难。 考虑非对称的线性方程 Ax=b (1) 其中，矩阵A=(珥，)，ai，∈R是一个大型稀疏矩阵，向量x和b分别是待求的未知向量和已知的 ’国家自然科学基金资助项目10572002) 127 输入向量。 在我们的问题中，系数矩阵A在数值上是非对称的，但是A的上三角和下三角部分在非零元的 分布上具有几乎对称的特征。对这样的方程组，可以在矩阵A中少量结构非对称的位置填入零，使 矩阵A在结构上对称，然后利用对称矩阵已有的成熟方法进行计算。 MLPG法生成的系数矩阵A是各阶顺序主子式全为正的，这种特性使矩阵A存在唯一的三角分 解如下： A=LDU (2) 其中，L和U分别是单位下三角和单位上三角矩阵，D是对角矩阵。在完成三角分解之后，式 (1)的求解可以分为三个步骤：Ly=b，Dz=y，Ux=z嘲。当A为对称矩阵，U=o。 没有任何疑问，L和U的计算是求解线性方程组(1)中计算量最高的步骤。如果矩阵A是稀疏的， 则一般L和U也是稀疏的，但是呸，／≠0中的(f，歹)集是非零元结构