1、第 39 卷 第 1 期2023 年 2 月北 京 建 筑 大 学 学 报Journal of Beijing University of Civil Engineering and ArchitectureVol 39No 1Feb 2023文章编号:2096 9872(2023)01 0103 06大型稀疏线性方程组的数值解法刘长河(北京建筑大学 理学院，北京100044)摘要:在许多利用经典算法求线性方程组的数值解的过程中，系数矩阵中的零元素对计算结果没有影响，也就没有存储的必要。如果是大型稀疏线性方程组，这样可以节省大量的存储空间。为此，提出一种在 MATLAB 语言环境中仅储存系数矩
2、阵中非零元素的方法:利用 3 个 1 维数组储存系数矩阵中的非零元素及其在矩阵中的位置(行号，列号)。在编程时，忽略零元素参与的运算，可使计算量大大减少。这 2 个方面的改进使得利用经典算法求解大型稀疏线性方程组成为可能。借助于 Jacobi 迭代法进行的一系列数值实验，验证了这一探索的可行性。关键词:稀疏矩阵;大型矩阵;线性方程组;数值解中图分类号:O241.6文献标志码:ADOI:10 19740/j 2096-9872 2023 01 13开放科学(资源服务)标识码(OSID):收稿日期:2022 11 25基金项目:中国博士后科学基金项目(2018M641301)。作者简介:刘长河(1
3、966)，男，副教授，理学博士，研究方向:应用数学。Numerical Methods for Solving the Large-Scale Sparse Linear EquationsLIU Changhe(School of Science，Beijing University of Civil Engineering and Architecture，Beijing 100044)Abstract:The zero elements in coefficient matrix have no effect on the calculation results in the proce
4、ssof many classical algorithms used to solved the linear simultaneous equations，so there is no need tostore them If it is a large sparse system of linear equations，it can save a lot of storage space For thispurpose，a method in which three 1-dimentional arrays are used to store the elements and their
5、 positions(row number，column number)in the matrix is proposed by storing only the non-zero elements in thecoefficient matrix in the MATLAB language environment In addition，the calculation amount can begreatly reduced by ignoring the operation of zero elements in programming These two improvementsmak
6、e it possible to solve large sparse linear equations using classical algorithms With the help of aseries of numerical experiments carried out by Jacobi iterative method，the feasibility of this explorationis verifiedKeywords:sparse matrix;large-scale matrix;system of linear equations;numerical soluti
7、on1问题的提出求线性方程组(1)数值解的算法分为 2 类:直接解法和迭代解法。Ax=b(1)当系数矩阵 A 是低阶矩阵时，这些算法都能取得理想的计算结果。北 京 建 筑 大 学 学 报第 39 卷然而，在科学、工程和社会资产统计等计算中，经常会遇到大型线性方程组。方程组中方程与未知数的个数很多，在求其数值解时不仅需要占据大量的存储空间，而且计算量巨大，几乎所有的经典算法都束手无策。不过，大型线性方程组的系数矩阵往往是稀疏的，即其绝大多数元素为零。系数矩阵是稀疏矩阵的线性方程组称为稀疏线性方程组。虽然有些经典算法也能求大型稀疏线性方程组的数值解，但数量极少，针对性强。例如，追赶法只对三对角方程
8、组有效1。长期以来，不少学者致力于对大型稀疏线性方程组 的 数 值 解 法 的 研 究，取 得 了 不 少 成 果。CUTHILL 和 MCKEE 于 1969 年提出了一种缩减矩阵宽带的有效方法 CM 算法。此算法后来被改进为更加有效的 CM 算法2，广泛应用于许多直接解法中。但是，CM 算法的基本思想是将原方程组分解成相互独立的尽可能小的方程组来求解，用到图论的知识，非常复杂。直到近年，仍有不少学者对这一问题进行深入研究，并提出了求解大型稀疏线性方程组的一些新算法3 7。但是，这些算法大多要求系数矩阵 A 具有一定的特点。比如，A 是对称正定矩阵或要求其非零元素的分布具有较为明显的规律。有
9、些算法甚至与特定的物理模型密切相关，可移植性不强。1997 年，王崧等8 利用 C 语言中指针和结构的灵活性，提出了一种链式三元表组法，能够高效地处理高阶稀疏矩阵的运算问题。然而，对于众多非计算机专业的科技工作者来讲，要熟练掌握 C 语言的编程技巧并非易事。也许正是此缘故，他们只在文章中进行了详细的理论阐述和分析，并没有进行实例验证，真实计算效果缺乏实践支撑。目前，科技工作者多采用 MATLAB 语言进行数值计算。这是一种基于C+语言的高级语言，更加简单，更加符合数学表达式的书写格式，是当前流行的计算机软件之一。笔者在 MATLAB 语言环境下，利用 3 个 1 维数组，将系数矩阵 A 中非零
10、元素存储起来，并能很方便地调用。这种方法在计算大型稀疏线性方程组时具有明显的优点:第一，只存储系数矩阵的非零元素，大大节省了数据存储空间。第二，忽略了零元素参与的大量无效计算，极大地减少了计算量。利用本文的方法进行数据存储和编程，可以将一些原本只能计算中小型线性方程组数值解的经典算法用于求解大型稀疏线性方程组。最后，笔者利用本文的方法，借助 Jacobi 迭代法，在微机上对从含有数千至数万个未知数的稀疏线性方程组进行了求数值解的反复实验，取得了令人满意的结果。2MATLAB 软件中“find()”函数介绍矩阵中零元素所占的比例可由其稀疏度表示:稀疏度 s=零元素个数矩阵中元素总数稀疏矩阵即稀疏
11、度接近于 1 的矩阵。MATLAB 软件中提供了“find()”函数，能够给出矩阵中所有非零元素的值和其在矩阵中的位置(行号，列号)。为简单计，以 4 阶矩阵为例，简要介绍此函数(本文用到)的功能。A=110131402223000330042044(2)执行 I，J，V=find(A)后的结果(V表示 V 的转置，其余类同):V=1122421323331444;I=12412314;J=12233344。可以看出，V 是把矩阵 A 的所有非零元素按照“第 1 列从上到下第 2 列从上到下”顺序排成的 1 维数组。I，J 分别是数组 V 中相应数据在矩阵 A 中所处的“行标”和“列标”所组成
12、的 1 维数组。这样，对 A 执行“find()”函数后，得到 3 个 1维数组:I，J，V。这 3 个数组将矩阵 A 中的非零元素的数值及它们在矩阵中的位置表示出来。在计算过程中，如果零元素不对计算结果产生影响，就可以忽略这些运算，从而减少计算次数。此时，没有必要占用内存去存储这些零元素。这种情况下，利用 3 个 1 维数组:I，J，V 代替矩阵 A(2 维数组)，可以大大减小存储空间和计算量。事实上，利用 1 个 2 维数组存储 n 阶矩阵 A 需要 n2个存储单元，而 3 个 1 维数组:I，J，V 只占用 3(1 s)n2个存储单元。当系数矩阵是稀疏矩阵，如稀疏度s=0.999 时，3
13、 个 1 维数组 I，J，V 的存储空间之和只是 2 维数组的 3。401第 1 期刘长河:大型稀疏线性方程组的数值解法3系数矩阵中非零元素的存储与调用在求线性方程组(1)数值解时，经常对系数矩阵进行初等行变换。首先需要确定矩阵 A 中各行的非零元素;其次，要确定各非零元素的列号。这种需求可以利用“find(A)”来实现。以式(2)的矩阵A 为例加以说明。执行 c，r，v=find(A)，则有如下结果:v=1113142223334244;r=11122344;c=13423324。可以看出，v 是把矩阵 A 的所有非零元素按照“第 1 行从左到右第 2 行从左到右”顺序排成的 1 维数组。r
14、，c 分别是数组 v 中相应元素的行标号，列标号所组成的 1 维数组。为建立数组 v 中元素在矩阵 A 中的位置(行号，列号)与其在数组 v 中的位置(序号)之间的关系，再引入以下 2 个 1 维数组。3.11 维数组 ghf0s 储存矩阵 A 中各行非零元素的数目数组 ghf0s 中的元素个数等于矩阵 A 的行数(此时，A 的行数为 4)，其各元素可通过执行下述程序得到:for i=1 4ghf0s(i)=length(find(r=i);(3)end运行结果:ghf0s=3212。矩阵 A 中各行的非零元素数目分别是3，2，1，2。3.21 维数组 hf0slj 储存矩阵 A 中前 i(i
15、=1，2，3，4)行的非零元素的数目累计数组 hf0slj 的元素个数等于矩阵 A 的行数(此时，A 的行数为 4)，可通过执行下述程序得到:for i=1 4hf0slj(i)=sum(ghf0s(1 i);(4)end执行结果:hf0slj=3568。矩阵 A 中前 1，2，3，4 行非零元素数目累计分别是 3，5，6，8。数组 hf0slj 的第 1 个元素即矩阵 A 中第 1 行非零元素的数目，最后一个元素即矩阵 A 中各行非零元素的总数。这样，对 于 式(2)矩 阵 A，执 行 c，r，v=find(A)和程序(3)、(4)后，得到的 3 个数组 v，hf0slj 和 c 见表 1。
16、表 1由式(2)矩阵 A 所得的 3 个 1 维向量组Tab 1Three 1-dimensional vector groups obtainedfrom the matrix A in formula(2)数组名数组长度数组元素v81113142223334244c813423324hf0slj43568为清晰起见，表 1 中使用了下划线。例如，“2223”表示这 2 个数同在矩阵 A 中第 2 行;其余类同。可以看出，矩阵 A 中的非零元素存储在 1 维数组 v 中，非零元素所在的列号存储在 1 维数组 c 中。具体如下:A 中第 1 行非零元素:v(1)，v(hf0slj(1);第 1 行非零元素列号:c(1)，c(hf0slj(1)。A 中第 2 行非零元素:v(hf0slj(1)+1)，v(hf0slj(2);第 2 行非零元素列号:c(hf0slj(1)+1)，c(hf0slj(2)。A 中第 3 行非零元素:v(hf0slj(2)+1)，v(hf0slj(3);第 3 行非零元素列号:c(hf0slj(2)+1)，c(hf0slj(3)。A 中第 4 行非零元素:v(h