求解大型稀疏矩阵方程组的SPIKE算法

The SPIKE Algorithm for Solving Large Sparse Matrix Equations

摘要

不同于传统的LU分解算法和QR分解算法，本文研究了一种新的基于DS矩阵分解的递归SPIKE算法。SPIKE算法采用了一种新颖的分解方法来平衡通信和算法开销，相比其他方法在现代并行架构上有更好的延展性。首先，从系数矩阵的分块、DS分解、简化系数矩阵方程组的提取和求解四方面介绍了递归SPIKE算法的工作原理。然后，首次将其应用到具体的系数矩阵规模不同的线性方程组中，并与LU分解算法与QR分解算法进行了比较。三组数值实验分别给出了各个求解算法的结果和运行时间。实验结果表明，递归SPIKE算法不仅能够求解得到准确结果，而且求解速度更快。数值案例表明，递归SPIKE算法所需的计算时间约为LU算法的40%，约为QR分解算法的8%。

Abstract

A recursive SPIKE algorithm was presented based on DS matrix decomposition， which was very different from LU decomposition algorithm and QR decomposition algorithm. The SPIKE algorithm used a novel decomposition method to balance communication overhead with arithmetic cost to achieve better scalability than other methods on modern parallel architectures. Firstly， the basic operating principle of a recursive SPIKE algorithm was introduced by four aspects： partitioning coefficient matrix， DS decomposition， extraction and solution of simplified coefficient matrix equations. The SPIKE algorithm was applied to specific linear equations with different coefficient matrix for the first time. Three sets of numerical experiments gave the results and running time of each solution algorithm. The experimental results show that the recursive SPIKE algorithm can not only obtain accurate results， but also has a faster solution speed.Experiments indicate the computing time of the proposed algorithms is about 40 percent of LU decomposition algorithm， and about 8 percent of QR decomposition algorithm.

关键词

一般带状矩阵 / 三对角矩阵 / DS矩阵分解 / 递归SPIKE算法

Key words

general banded matrix / tridiagonal matrix / DS matrix decomposition / recursive SPIKE algorithm

0 引 言

偏微分方程在实际工程应用中起着至关重要的作用，它经常出现在计算机科学和工程等领域，如有限元分析^［1］、计算力学（流体、结构和流体结构的相互作用）^［2‐6］、计算纳米电子学^［7］、计算机视觉^［8］。然而，对于偏微分方程，很难求得其解析解，于是人们转而求解其数值近似解。通常通过离散化方法将大型稀疏矩阵化为系数矩阵的线性方程组来求解。大型稀疏矩阵大多为带宽较窄的大型带状矩阵。带状矩阵是指一个矩阵中只有对角线附近的元素可以是非零的，而其他位置的元素均为零。这种矩阵结构的特点是它拥有较少的非零元素，因此储存和计算时占用的资源较少，这使得它在实际问题的求解中具有很大优势。虽然求解带状矩阵线性方程组所需的储存和计算量相比一般的矩阵线性方程组要小，但是在解带状矩阵线性方程组时，当前主流的直接求解方法，如LU分解算法^［9］和QR分解算法^［10］，由于要处理大量的零元素，效率往往较低。因此，如何高效快速准确地求解带状矩阵线性方程组是一个重要的研究课题。

本文主要探究一种新的基于DS矩阵分解的求解大型带状矩阵线性方程组的递归SPIKE算法^{［11‐15］}的性能。2001年，基于DS矩阵分解的递归SPIKE算法问世，该算法通过对系数矩阵的分块形式进行简化分解来降低计算复杂度。2008年，递归SPIKE算法引入随机的DS矩阵分解来增加算法的随机性，从而进一步提高算法的计算速度。2020年，Braegan S.Spring突破了递归SPIKE算法求解大型带状矩阵线性方程组时系数矩阵阶数为{\mathrm{2}}^n的限制，使其能够求解更为一般的大型带状矩阵线性方程组^［11］。本文首次将递归SPIKE算法运用于具体的数值算例，并与当前主流的直接算法LU分解算法和QR分解算法进行比较，以此来检验基于DS矩阵分解的递归SPIKE算法的优越性，即递归SPIKE算法相比LU分解和QR分解在求解的准确性和速度上的优越性能。本文研究可以为国产化服务器的高性能计算提供更优的求解大型带状线性系统的工具。

1 递归SPIKE算法

递归SPIKE算法最早可以追溯到20世纪70年代，其首先应用于求解三对角矩阵线性方程组，后来扩展到处理带状矩阵线性方程组。一般带状矩阵可以经过一定的分块看成块状三对角矩阵，因此，递归SPIKE算法可以被看作一种求解块状三对角矩阵的区域分解算法。递归SPIKE算法的核心思想与传统的LU分解算法的思想不同，引入了一种新的DS分解算法，通过这样的矩阵分解方式，能够有效降低计算的复杂度。

1.1 系数矩阵的分块

求解\mathit{A}\mathit{x}=\mathit{b}，其中矩阵\mathit{A}是一个n={\mathrm{2}}^m(m\in {\mathbf{N}}^{+})阶的带状矩阵。一些带宽较窄的n阶带状矩阵都可以划分为块状三对角矩阵。假定划分对角线上对角块的数量为p，那么每个带状对角块{\mathit{A}}_j(j=\mathrm{1},\cdots ,p)的阶数为n_j（或者大致为n/p）。对于已经划定好的分区，{\mathit{B}}_j(j=\mathrm{1},\cdots ,p-\mathrm{1})和{\mathit{C}}_j(j=\mathrm{2},\cdots ,p)是与对角线上对角块相耦合的次对角线上的块状矩阵，具体为

\mathit{A}=\left[\begin{array}{cccc}{\mathit{A}}_{\mathrm{1}}& {\mathit{B}}_{\mathrm{1}}& & \\ {\mathit{C}}_{\mathrm{2}}& {\mathit{A}}_{\mathrm{2}}& \ddots & \\ & \ddots & \ddots & {\mathit{B}}_{p-\mathrm{1}}\\ & & {\mathit{C}}_p& {\mathit{A}}_p\end{array}\right],（1）

式中：{\mathit{B}}_j和{\mathit{C}}_j的阶数均为n_j\times k，k=\frac{\mathrm{矩阵}\mathrm{的带}\mathrm{宽数}-\mathrm{1}}{\mathrm{2}}，且矩阵{\mathit{B}}_j和{\mathit{C}}_j包含着大量的零元素，具体为

{\mathit{B}}_j=\left[\begin{array}{cc}\mathrm{0}& \mathrm{0}\\ {\overline{\mathit{B}}}_j& \mathrm{0}\end{array}\right],{\mathit{C}}_j=\left[\begin{array}{cc}\mathrm{0}& {\overline{\mathit{C}}}_j\\ \mathrm{0}& \mathrm{0}\end{array}\right],（2）

式中：{\overline{\mathit{B}}}_j和{\overline{\mathit{C}}}_j是k阶矩阵。

1.2 系数矩阵的DS分解

基于以上的分块形式，接下来我们将对系数矩阵 进行DS分解，且具体分解形式为

\begin{array}{c}\mathit{A}=\mathit{D}\mathit{S}=\\ \left[\begin{array}{cccc}{\mathit{D}}_{\mathrm{1}}& & & \\ & {\mathit{D}}_{\mathrm{2}}& & \\ & & \ddots & \\ & & & {\mathit{D}}_p\end{array}\right]\left[\begin{array}{cccc}{\mathit{I}}_{\mathrm{1}}& {\mathit{V}}_{\mathrm{1}}& & \\ {\mathit{W}}_{\mathrm{2}}& {\mathit{I}}_{\mathrm{2}}& \ddots & \\ & \ddots & \ddots & {\mathit{V}}_{p-\mathrm{1}}\\ & & {\mathit{W}}_p& {\mathit{I}}_p\end{array}\right],\end{array}（3）

式中：{\mathit{I}}_j表示一个矩阵阶数为n_j的单位矩阵且{\mathit{D}}_j\equiv {\mathit{A}}_j，而{\mathit{V}}_j和{\mathit{W}}_j的非零元素所在的列形成了大小为n_j\times k（又称尖峰）的高而窄的子矩阵。由上述DS分解的矩阵乘法可以得到{\mathit{V}}_j和{\mathit{W}}_j的表达式

{\mathit{V}}_j=({\mathit{A}}_j{)}^{-\mathrm{1}}{\mathit{B}}_j,{\mathit{W}}_j=({\mathit{A}}_j{)}^{-\mathrm{1}}{\mathit{C}}_j\mathrm{。}

{\mathit{V}}_j和{\mathit{W}}_j可以通过求解矩阵方程（4）得到。

{\mathit{A}}_j\left[\begin{array}{cc}{\mathit{V}}_j& {\mathit{W}}_j\end{array}\right]=\left[\begin{array}{cc}\mathrm{0}& {\dot{\mathit{C}}}_j\\ ⋮& \mathrm{0}\\ \mathrm{0}& ⋮\\ {\dot{\mathit{B}}}_j& \mathrm{0}\end{array}\right]\mathrm{。}（4）

1.3 简化系数矩阵方程组的提取

复杂矩阵线性方程组的求解经过前面的分解后缩减到两个步骤，即\mathit{D}\mathit{G}=\mathit{b}和\mathit{S}\mathit{x}=\mathit{G}。

求解\mathit{D}\mathit{G}=\mathit{b}中的线性方程组所得向量\mathit{G}将作为求解\mathit{S}\mathit{x}=\mathit{G}的线性方程组的右端向量。在求解\mathit{D}\mathit{G}=\mathit{b}的线性方程组时，利用矩阵\mathit{D}的特殊对角块形式，将线性方程组\mathit{D}\mathit{G}=\mathit{b}化为p个阶数较小的系数矩阵线性方程组来求解，这将在很大程度上降低计算的复杂度。如果将DS分解阶段和简化线性方程组求解阶段分离，则求解\mathit{D}\mathit{G}=\mathit{b}的过程可以与分解过程中的尖峰矩阵{\mathit{V}}_j和{\mathit{W}}_j的生成相结合。

将求解线性方程组\mathit{S}\mathit{x}=\mathit{G}进一步化简为求解一个系数矩阵阶数更小的线性方程组

\widehat{\mathit{S}}\widehat{\mathit{x}}=\widehat{\mathit{G}},（5）

式中：\widehat{\mathit{S}}由矩阵\mathit{S}每个分块的正上方和正下方的k行组成。实际上，尖峰矩阵{\mathit{V}}_j和{\mathit{W}}_j也可以划分为

{\mathit{V}}_j=\left[\begin{array}{c}{\mathit{V}}_j^{(t)}\\ {\mathit{V}}_j^{\text{'}}\\ {\mathit{V}}_j^{(b)}\end{array}\right],{\mathit{W}}_j=\left[\begin{array}{c}{\mathit{W}}_j^{(t)}\\ {\mathit{W}}_j^{\mathrm{\text{'}}}\\ {\mathit{W}}_j^{(b)}\end{array}\right],（6）

式中：{\mathit{V}}_j^{(t)}、{\mathit{V}}_j^{\mathrm{\text{'}}}、{\mathit{V}}_j^{(b)}和{\mathit{W}}_j^{(t)}、{\mathit{W}}_j^{\mathrm{\text{'}}}、{\mathit{W}}_j^{(b)}是{\mathit{V}}_j和{\mathit{W}}_j的顶部k行、中间n_j-\mathrm{2}k和底部k行。这里，

{\mathit{V}}_j^{(b)}=\left[\begin{array}{cc}\mathrm{0}& {\mathit{I}}_k\end{array}\right]{\mathit{V}}_j,{\mathit{W}}_j^{(t)}=\left[\begin{array}{cc}{\mathit{I}}_k& \mathrm{0}\end{array}\right]{\mathit{W}}_j,

和 {\mathit{V}}_j^{(t)}=\left[\begin{array}{cc}{\mathit{I}}_k& \mathrm{0}\end{array}\right]{\mathit{V}}_j,{\mathit{W}}_j^{(b)}=\left[\begin{array}{cc}\mathrm{0}& {\mathit{I}}_k\end{array}\right]{\mathit{W}}_j,

类似地，如果{\mathit{x}}_j和{\mathit{G}}_j是列向量\mathit{x}和\mathit{G}的第j个分区，有

{\mathit{x}}_j=\left[\begin{array}{c}{\mathit{x}}_j^{(t)}\\ {\mathit{x}}_j^{\text{'}}\\ {\mathit{x}}_j^{(b)}\end{array}\right],{\mathit{G}}_j=\left[\begin{array}{c}{\mathit{G}}_j^{(t)}\\ {\mathit{G}}_j^{\text{'}}\\ {\mathit{G}}_j^{(b)}\end{array}\right]\mathrm{。}（7）

从\mathit{S}\mathit{x}=\mathit{G}中提取出简化的线性方程组（5），该线性方程组仅包含{\mathit{V}}_j、{\mathit{W}}_j、{\mathit{x}}_j和{\mathit{G}}_j顶部和底部k行。这个简化的线性方程组的系数矩阵\widehat{\mathit{S}}包含p-\mathrm{1}个对角线对角块，其第i个对角块表示为

\left[\begin{array}{cc}{\mathit{I}}_k& {\mathit{V}}_i^{(b)}\\ {\mathit{W}}_{i+\mathrm{1}}^{(t)}& {\mathit{I}}_k\end{array}\right],（8）

相应的次对角块表示为

\left[\begin{array}{cc}{\mathit{W}}_i^{(b)}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}\end{array}\right],\left[\begin{array}{cc}\mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\mathit{V}}_{i+\mathrm{1}}^{(t)}\end{array}\right],（9）

该缩减线性方程组中的解向量块和右端向量块表示为

\left[\begin{array}{c}{\mathit{x}}_i^{(b)}\\ {\mathit{x}}_{i+\mathrm{1}}^{(t)}\end{array}\right],\left[\begin{array}{c}{\mathit{G}}_i^{(b)}\\ {\mathit{G}}_{i+\mathrm{1}}^{(t)}\end{array}\right]\mathrm{。}（10）

如果已经求解出缩减线性方程组（5），得到原线性方程组的部分解，就可以快速检索获得原线性方程组的全局解，具体求解方程为

\left\{\begin{array}{l}{\mathit{x}}_{\mathrm{1}}^{\text{'}}={\mathit{G}}_{\mathrm{1}}^{\text{'}}-{\mathit{V}}_{\mathrm{1}}^{\text{'}}{x}_{\mathrm{2}}^{(t)},\\ {\mathit{x}}_j^{\text{'}}={\mathit{G}}_j^{\text{'}}-{\mathit{V}}_j^{\text{'}}{x}_{j+\mathrm{1}}^{(t)}-{\mathit{W}}_j^{\text{'}}{\mathit{x}}_{j-\mathrm{1}}^{(b)},j=\mathrm{2},\cdots ,p-\mathrm{1},\\ {\mathit{x}}_p^{\text{'}}={\mathit{G}}_p^{\text{'}}-{\mathit{W}}_p^{\text{'}}{x}_{p-\mathrm{1}}^{\left(b\right)}\mathrm{。}\end{array}\right.（11）

1.4 简化系数矩阵方程组的求解

并行求解简化系数矩阵方程组（5）的一种自然方法是\mathrm{K}\mathrm{r}\mathrm{y}\mathrm{l}\mathrm{o}\mathrm{v}子空间迭代法^［16］。然而对于系数矩阵为非对角占优矩阵的线性方程组，该处理方法对大量的分块矩阵无效。因此，为了实现较小的相对残差，需要进行大量的外部迭代，这反过来又会导致计算复杂度的增加。此外，如果简化的线性方程组的系数矩阵阶数较大，那么这种方案可能因内存限制而无法实现。为此，一种新的求解简化线性方程组的直接递归方法被提出，这种递归方案涉及SPIKE算法的连续迭代，能够有效降低算法的计算复杂度。

如果简化线性方程组（5）的系数矩阵只有两个分块，即p=\mathrm{2}，那么，只由一个对角块组成的系数矩阵简化线性方程组就可以直接求解，其线性方程组表示为

\left[\begin{array}{cc}{\mathit{I}}_k& {\mathit{V}}_{\mathrm{1}}^{(b)}\\ {\mathit{W}}_{\mathrm{2}}^{(t)}& {\mathit{I}}_k\end{array}\right]\left[\begin{array}{c}{\mathit{x}}_{\mathrm{1}}^{(b)}\\ {\mathit{x}}_{\mathrm{2}}^{(t)}\end{array}\right]=\left[\begin{array}{c}{\mathit{G}}_{\mathrm{1}}^{(b)}\\ {\mathit{G}}_{\mathrm{2}}^{(t)}\end{array}\right]。（12）

具体求解步骤：首先计算\mathit{E}={\mathit{I}}_k-{\mathit{W}}_{\mathrm{2}}^{(t)}{\mathit{V}}_{\mathrm{1}}^{(b)},然后通过\mathit{E}{\mathit{x}}_{\mathrm{2}}^{(t)}={\mathit{G}}_{\mathrm{2}}^{(t)}-{\mathit{W}}_{\mathrm{2}}^{\left(t\right)}{\mathit{G}}_{\mathrm{1}}^{\left(b\right)}来获取{\mathit{x}}_{\mathrm{2}}^{(t)}，最后计算{\mathit{x}}_{\mathrm{1}}^{(b)}={\mathit{G}}_{\mathrm{1}}^{(b)}-{\mathit{V}}_{\mathrm{1}}^{\left(b\right)}{\mathit{x}}_{\mathrm{2}}^{\left(t\right)}。

{\mathit{x}}_{\mathrm{1}}和{\mathit{x}}_{\mathrm{2}}其余部分的解可以通过式（11）来获得。

下面讨论系数矩阵的多分块情况，假设分块数量p={\mathrm{2}}^d，经过分解后形成SPIKE矩阵\mathit{S}，新缩减线性方程组的分区数为2的倍数，并且可以利用另一个级别的SPIKE算法。这个过程被递归地重复执行，直到最新的矩阵\mathit{S}只有两个分区，由此得到的简化线性方程组的形式如式（12）。

在实际应用中，递归方案与整体矩阵\mathit{S}无关，而是与简化线性方程组（5）中的系数矩阵\widehat{\mathit{S}}有关。这样能够简化执行SPIKE算法，减少内存占有量，同时保存所有不同的级别的尖峰矩阵（{\mathit{V}}_j和{\mathit{W}}_j）。值得注意的是，在简化线性方程组（5）中，系数矩阵\widehat{\mathit{S}}是块状三对角矩阵，对角块由式（8）给出，次对角块由式（9）给出。

如果将矩阵\mathit{S}中的第一个分区和最后一个分区的顶部k行和底部k行提取到简化线性方程组系数矩阵\widehat{\mathit{S}}中，那么，矩阵\widehat{\mathit{S}}的阶数为\mathrm{2}kp，而不是\mathrm{2}k(p-\mathrm{1})，并且阶数为\mathrm{2}kp的矩阵\widehat{\mathit{S}}的结构仍然保持块状三对角的形式。在这种情况下，每个对角线上的对角块都是阶数为\mathrm{2}k的单位矩阵，与第i个对角块相关联的次对角块为

\left[\begin{array}{cc}\mathrm{0}& {\mathit{W}}_i^{(t)}\\ \mathrm{0}& {\mathit{W}}_i^{(b)}\end{array}\right](i=\mathrm{2},\cdots ,p),

\left[\begin{array}{cc}{\mathit{V}}_i^{(t)}& \mathrm{0}\\ {\mathit{V}}_i^{(b)}& \mathrm{0}\end{array}\right]\left(i=\mathrm{1},\cdots ,p-\mathrm{1}\right)\mathrm{。}（13）

将{\mathit{V}}_i^{[\mathrm{1}]}和{\mathit{W}}_i^{[\mathrm{1}]}作为执行第一次递归的简化线性方程组的尖峰矩阵，其中，

{\mathit{V}}_i^{[\mathrm{1}]}=\left[\begin{array}{c}{\mathit{V}}_i^{(t)}\\ {\mathit{V}}_i^{(b)}\end{array}\right],{\mathit{W}}_i^{[\mathrm{1}]}=\left[\begin{array}{c}{\mathit{W}}_i^{(t)}\\ {\mathit{W}}_i^{(b)}\end{array}\right]\mathrm{。}（14）

当p=\mathrm{4}时，新的简化线性方程组的系数矩阵{\tilde{\mathit{S}}}_{\mathrm{1}}表示为

{\tilde{\mathit{S}}}_{\mathrm{1}}=\left[\begin{array}{cccc}\mathit{I}& {\mathit{V}}_{\mathrm{1}}^{[\mathrm{1}]}& & \\ {\mathit{W}}_{\mathrm{2}}^{[\mathrm{1}]}& \mathit{I}& {\mathit{V}}_{\mathrm{2}}^{[\mathrm{1}]}& \\ & {\mathit{W}}_{\mathrm{3}}^{[\mathrm{1}]}& \mathit{I}& {\mathit{V}}_{\mathrm{3}}^{[\mathrm{1}]}\\ & & {\mathit{W}}_{\mathrm{4}}^{[\mathrm{1}]}& \mathit{I}\end{array}\right]\mathrm{。}（15）

在进行SPIKE算法的第二次递归时，先将{\tilde{\mathit{S}}}_{\mathrm{1}}进行分块，每个分块的大小为\mathrm{4}k，一共有p/\mathrm{2}个对角块。然后对矩阵{\tilde{\mathit{S}}}_{\mathrm{1}}进行DS分解

{\tilde{\mathit{S}}}_{\mathrm{1}}={\mathit{D}}_{\mathrm{1}}{\tilde{\mathit{S}}}_{\mathrm{2}},（16）

式中：{\mathit{D}}_{\mathrm{1}}由矩阵{\tilde{\mathit{S}}}_{\mathrm{1}}的对角线上对角块组成，每个块的阶数为\mathrm{4}k。因此，{\tilde{\mathit{S}}}_{\mathrm{2}}是由尖峰矩阵{\mathit{V}}_i^{[\mathrm{2}]}和{\mathit{W}}_i^{[\mathrm{2}]}组成。当p=\mathrm{4}时，这些矩阵表示为

{\mathit{D}}_{\mathrm{1}}=\left[\begin{array}{cc}\left[\begin{array}{cc}{\mathit{I}}_{\mathrm{2}k}& {\mathit{V}}_{\mathrm{1}}^{[\mathrm{1}]}\\ {\mathit{W}}_{\mathrm{2}}^{[\mathrm{1}]}& {\mathit{I}}_{\mathrm{2}k}\end{array}\right]& \\ & \left[\begin{array}{cc}{\mathit{I}}_{\mathrm{2}k}& {\mathit{V}}_{\mathrm{3}}^{[\mathrm{1}]}\\ {\mathit{W}}_{\mathrm{4}}^{[\mathrm{1}]}& {\mathit{I}}_{\mathrm{2}k}\end{array}\right]\end{array}\right]

和 {\tilde{\mathit{S}}}_{\mathrm{2}}=\left[\begin{array}{cc}{\mathit{I}}_{\mathrm{4}k}& {\mathit{V}}_{\mathrm{1}}^{[\mathrm{2}]}\\ {\mathit{W}}_{\mathrm{2}}^{[\mathrm{2}]}& {\mathit{I}}_{\mathrm{4}k}\end{array}\right]\mathrm{。}

总的来说，在第j次递归的过程中，尖峰矩阵{\mathit{V}}_i^{[j]}和{\mathit{W}}_i^{[j]}（i=1~p/{\mathrm{2}}^j）的阶数为{\mathrm{2}}^{j}k\times k。因此，如果原始矩阵的分区数为p={\mathrm{2}}^d，则需要进行递归操作的次数为d-\mathrm{1}，且矩阵{\tilde{\mathit{S}}}_{\mathrm{1}}的分解形式为

{\tilde{\mathit{S}}}_{\mathrm{1}}={\mathit{D}}_{\mathrm{1}}{\mathit{D}}_{\mathrm{2}}\cdots {\mathit{D}}_{d-\mathrm{1}}{\tilde{\mathit{S}}}_d,（17）

式中：矩阵{\tilde{\mathit{S}}}_d仅含有两个尖峰矩阵{\mathit{V}}_{\mathrm{1}}^{[d]}和{\mathit{W}}_{\mathrm{2}}^{[d]}。

简化的线性矩阵方程组可以写成

{\tilde{\mathit{S}}}_d\tilde{\mathit{x}}=\mathit{B},（18）

式中：\mathit{B}为修正的右端向量，具体表示为

\mathit{B}={\mathit{D}}_{d-\mathrm{1}}^{-\mathrm{1}}\cdots {\mathit{D}}_{\mathrm{2}}^{-\mathrm{1}}{\mathit{D}}_{\mathrm{1}}^{-\mathrm{1}}\tilde{\mathit{G}}\mathrm{。}

如果假设矩阵{\tilde{\mathit{S}}}_j的尖峰块状矩阵{\mathit{V}}_i^{[j]}和{\mathit{W}}_i^{[j]} （\forall i），在给定j的情况下，可以计算第j+\mathrm{1}次递归的尖峰矩阵{\mathit{V}}_i^{[j+\mathrm{1}]}和{\mathit{W}}_i^{[j+\mathrm{1}]}。第j次递归的尖峰矩阵的顶部和底部k行分别表示为

{\mathit{V}}_{i\mathrm{\text{'}}}^{[j](b)}=\left[\begin{array}{cc}\mathrm{0}& {\mathit{I}}_k\end{array}\right]{\mathit{V}}_i^{[j]},(19)

{\mathit{W}}_{i\mathrm{\text{'}}}^{[j](t)}=\left[\begin{array}{cc}{\mathit{I}}_k& \mathrm{0}\end{array}\right]{\mathit{W}}_i^{[j]}\mathrm{。}(20)

第j+\mathrm{1}次递归的尖峰矩阵的中间\mathrm{2}k行可表示为

\left[\begin{array}{c}{\dot{\mathit{V}}}_i^{[j+\mathrm{1}]}\\ {\ddot{\mathit{V}}}_i^{[j+\mathrm{1}]}\end{array}\right]=\left[\begin{array}{ccc}\mathrm{0}& {\mathit{I}}_{\mathrm{2}k}& \mathrm{0}\end{array}\right]{\mathit{V}}_i^{[j+\mathrm{1}]},（21）

\left[\begin{array}{c}{\dot{\mathit{W}}}_i^{[j+\mathrm{1}]}\\ {\ddot{\mathit{W}}}_i^{[j+\mathrm{1}]}\end{array}\right]=\left[\begin{array}{ccc}\mathrm{0}& {\mathit{I}}_{\mathrm{2}k}& \mathrm{0}\end{array}\right]{\mathit{W}}_i^{[j+\mathrm{1}]},（22）

则可以形成简化的矩阵方程

\begin{array}{l}\left[\begin{array}{cc}{\mathit{I}}_k& {\mathit{V}}_{\mathrm{2}i-\mathrm{1}}^{[j](b)}\\ {\mathit{W}}_{\mathrm{2}i}^{[j](t)}& {\mathit{I}}_k\end{array}\right]\left[\begin{array}{c}{\dot{\mathit{V}}}_k^{[j+\mathrm{1}]}\\ {\ddot{\mathit{V}}}_k^{[j+\mathrm{1}]}\end{array}\right]=\left[\begin{array}{c}\mathrm{0}\\ {\mathit{V}}_{\mathrm{2}i}^{[j](t)}\end{array}\right],\\ i=\mathrm{1,2},\cdots ,\frac{p}{{\mathrm{2}}^{i-\mathrm{1}}}-\mathrm{1}\end{array}

和 \left[\begin{array}{cc}{\mathit{I}}_k& {\mathit{V}}_{\mathrm{2}i-\mathrm{1}}^{[j](b)}\\ {\mathit{W}}_{\mathrm{2}i}^{[j](t)}& {\mathit{I}}_k\end{array}\right]\left[\begin{array}{c}{\dot{\mathit{W}}}_i^{[j+\mathrm{1}]}\\ {\ddot{\mathit{W}}}_i^{[j+\mathrm{1}]}\end{array}\right]=\left[\begin{array}{c}{\mathit{W}}_{\mathrm{2}i-\mathrm{1}}^{[j](b)}\\ \mathrm{0}\end{array}\right],

i=\mathrm{2,3},\cdots ,\frac{p}{{\mathrm{2}}^{i-\mathrm{1}}}\mathrm{。}

与线性方程组式（12）类似，可以利用求解式（12）的方法来获得j+\mathrm{1}次递归的尖峰矩阵的中间\mathrm{2}k行。然后检索获得第i+\mathrm{1}次递归的全部尖峰矩阵

\begin{array}{l}\left[\begin{array}{cc}{\mathit{I}}_{{\mathrm{2}}^{j}k}& \mathrm{0}\end{array}\right]{\mathit{V}}_i^{[j+\mathrm{1}]}=-{\mathit{V}}_{\mathrm{2}i-\mathrm{1}}^{[j]}{\ddot{\mathit{V}}}_i^{[j+\mathrm{1}]},\\ \left[\begin{array}{cc}\mathrm{0}& {\mathit{I}}_{{\mathrm{2}}^{j}k}\end{array}\right]{\mathit{V}}_i^{[j+\mathrm{1}]}={\mathit{V}}_{\mathrm{2}i}^{[j]}-{\mathit{W}}_{\mathrm{2}i}^{[j]}{\dot{\mathit{V}}}_i^{[j+\mathrm{1}]},\end{array}

以及

\left[\begin{array}{cc}{\mathit{I}}_{{\mathrm{2}}^{j}k}& \mathrm{0}\end{array}\right]{\mathit{W}}_i^{[j+\mathrm{1}]}={\mathit{W}}_{\mathrm{2}i-\mathrm{1}}^{[j]}-{\mathit{V}}_{\mathrm{2}i-\mathrm{1}}^{[j]}{\ddot{\mathit{W}}}_i^{[j+\mathrm{1}]},\left[\begin{array}{cc}\mathrm{0}& {\mathit{I}}_{{\mathrm{2}}^{j}k}\end{array}\right]{\mathit{W}}_i^{[j+\mathrm{1}]}=-{\mathit{W}}_{\mathrm{2}i}^{\left[j\right]}{\dot{\mathit{W}}}_i^{\left[i+\mathrm{1}\right]}\mathrm{。}

接着，通过求解式（12）的方法来求解式（18）的解向量\tilde{\mathit{x}}。最后，通过式（11）检索获得整个解向量\mathit{x}。

2 数值实验

对两组带宽相同阶数不同的带状矩阵线性方程组进行数值实验，以证明递归SPIKE算法不仅能够得到预期的结果，而且相比传统的LU分解算法和QR分解算法，能够在求解速度上更具优越性。

数值实验环境配置：设备为LAPTOP-R6HU69BR；处理器为AMD Ryzen 7 6800H with Radeon Graphics，3.20 GHz；机带RAM为16.0 GB；系统为基于x64处理器的64位操作系统；平台为matlab R2022a。

LU分解算法和QR分解算法使用matlab内置函数来求解，以便将递归SPIKE算法与最优化的LU分解算法和QR分解算法进行比较。

例 1^［17］带宽为3的带状矩阵线性方程组，对应系数矩阵和右端向量分别为

\mathit{A}={\left[\begin{array}{ccccc}\mathrm{4}& \mathrm{2}& & & \\ \mathrm{1}& \mathrm{4}& \mathrm{1}& & \\ & \ddots & \ddots & \ddots & \\ & & \mathrm{1}& \mathrm{4}& \mathrm{1}\\ & & & \mathrm{2}& \mathrm{4}\end{array}\right]}_{n\times n},\mathit{b}={\left[\begin{array}{c}\mathrm{6}\\ \mathrm{6}\\ ⋮\\ \mathrm{6}\\ \mathrm{6}\end{array}\right]}_{n\times \mathrm{1}}\mathrm{。}

利用递归SPIKE算法、LU分解算法和QR分解算法分别求解系数矩阵阶数n=\mathrm{4}\mathrm{096},\mathrm{8}\mathrm{192},\mathrm{16}\mathrm{384}时的带状矩阵线性方程组，且3种算法都能求得准确的结果\mathit{x}={(\mathrm{1,1},\cdots ,\mathrm{1})}_{\mathrm{1}\times n}^{\mathrm{T}}。3种算法求解运行的时间如表 1 所示。

例2 带宽为5的带状矩阵线性方程组，对应系数矩阵和右端向量分别为

\mathit{A}={\left[\begin{array}{ccccc}\mathrm{35}& \mathrm{5}& \mathrm{1}& & \\ \mathrm{5}& \ddots & \ddots & \ddots & \\ \mathrm{1}& \ddots & \ddots & \ddots & \mathrm{1}\\ & \ddots & \ddots & \ddots & \mathrm{5}\\ & & \mathrm{1}& \mathrm{5}& \mathrm{35}\end{array}\right]}_{n\times n},\mathit{b}={\left[\begin{array}{c}\mathrm{41}\\ \mathrm{46}\\ \mathrm{47}\\ ⋮\\ \mathrm{47}\\ \mathrm{46}\\ \mathrm{41}\end{array}\right]}_{n\times \mathrm{1}}\mathrm{。}

利用递归SPIKE算法、LU分解算法和QR分解算法分别求解系数矩阵阶数n=\mathrm{4}\mathrm{096,8}\mathrm{192,16}\mathrm{384}时的带状矩阵线性方程组，且3种算法都能求得准确的结果\mathit{x}={(\mathrm{1,1},\cdots ,\mathrm{1})}_{\mathrm{1}\times n}^{\mathrm{T}}。3种算法求解运行的时间如表 2 所示。

从上面两组实验的结果可以看出，对带宽分别为3和5，矩阵阶数为4 096的带状矩阵线性方程组求解时，递归SPIKE算法的求解速度相比LU分解算法并没有很大优势，但是明显优于QR分解算法。当带状矩阵阶数变大时，即矩阵阶数为8 192和16 384时，SPIKE算法的求解速度远快于LU分解算法和QR分解算法。因此，对于求解实际工程问题中经常碰到的大型带状矩阵线性方程组，递归SPIKE算法优于传统的LU分解算法和QR分解算法。

3 结束语

与传统的矩阵分解方式不同，本文提出的递归SPIKE算法是一种基于新的DS矩阵分解方式的算法，该算法通过提取规模更小的简化线性方程组来进行求解，降低了算法的计算复杂度。数值实验结果表明，对于大型带状矩阵线性方程组的求解，递归SPIKE算法不仅能够取得良好的结果，而且求解速度更快。在后续工作中，将在SPIKE的处理细节方面进行进一步优化，以期得到更好的效果。

参考文献

[1] FREYTAG M， SHAPIRO V， TSUKANOV I.Finite element analysis in situ［J］.Finite Elements in Analysis and Design，2011，47（9）：957-972.

[2] OISHI A， YAGAWA G.Computational mechanics enhanced by deep learning［J］.Computer Methods in Applied Mechanics and Engineering，2017，327： 327-351.

[3] HERFF S， NIEMÖLLER A， MEINKE M， et al.LES of a turbulent swirl flame using a mesh adaptive level-set method with dynamic load balancing［J］.Computers & Fluids，2021，221：104900.

[4] SHAKEEL M R， MOKHEIMER E M A.Swirl flow in annular geometry with varying cross-section［J］.Engineering Applications of Computational Fluid Mechanics，2022，16（1）：1154-1172.

[5] ZEIDAN D， BÄHR P， FARBER P， et al.Numerical investigation of a mixture two-phase flow model in two-dimensional space［J］.Computers & Fluids，2019，181：90-106.

[6] 董建伟，娄光谱.量子流体动力学等温模型的拟中性极限［J］.中北大学学报（自然科学版），2012，33（2）：102-106.

[7] DONG Jianwei， LOU Guangpu.Quasineutral limit of isothermalk quantum hydrodynamic model［J］.Journal of North University of China（Natural Science Edition），2012，33（2）：102-106.（in Chinese）

[8] POLIZZI E， ABDALLAH N B.Subband decomposition approach for the simulation of quantum electron transport in nanostructures［J］.Journal of Computational Physics，2004，202：150-180.

[9] 杨彦利，苗长云，亢伉，等.输送带跑偏故障的机器视觉检测技术［J］.中北大学学报（自然科学版），2012，33（6）：667-671.

[10] YANG Yanli， MIAO Changyun， KANG Kang，et al.Machine vision inspection technique for conveyor belt deviation［J］.Journal of North University of China（Natural Science Edition），2012，33（6）：667-671.（in Chinese）

[11] BARTELS R H， GOLUB G H.The simplex method of linear programming using LU decomposition［J］.Communications of the ACM，1969，12（5）： 266-268.

[12] SHARMA A， PALIWAL K K， IMOTO S，et al.Principal component analysis using QR decomposition［J］.International Journal of Machine Learning and Cybernetics，2013，4（6）：679-683.

[13] SPRING B S， POLIZZI E， SAMEH A H.A feature-complete spike dense banded solver［J］.ACM Transactions on Mathematical Software，2020，46（4）：1-35.

[14] SPRING B S.Enhanced capabilities of the spike algorithm and a new spike-openMP solver［D］.Amherst： University of Massachusetts Amherst，2014.

[15] SAMEH A H， POLIZZI E.A parallel hybrid banded system solver：the SPIKE algorithm［J］.Parallel Computing，2006，32（2）：177-194.

[16] SPRING B S， POLIZZI E， SAMEH A H.A feature complete spike banded algorithm and solver［DB/OL］.（2018-11-08）［2023-07-06］.

[17] POLIZZI E， SAMEH A H.SPIKE：A parallel environment for solving banded linear systems［J］.Computers & Fluids，2007，36（1）：113-120.

[18] STYKEL T， SIMONCINI V.Krylov subspace methods for projected Lyapunov equations［J］.Applied Numerical Mathematics，2012，62（1）：35-50.

[19] 冉瑞生.一些矩阵计算问题及其在图像识别中的应用研究［D］.成都：电子科技大学，2006.

基金资助

国家自然科学基金资助项目(11801281)

江苏省博士后科研资助项目(2020Z380)