前戏

在应用数学，工程学和IT行业的许多领域，无论是处理图形，计算飞机的空气动力学性能还是优化物流，线性方程组的数值求解都是必不可少的步骤。 没有它，现在流行的“机器”也不会有太大进步。 而且，通常，线性系统的解决方案在所有计算成本中所占的比例最大。 很明显，此关键组件需要最大速度优化。

经常与所谓的 稀疏矩阵-零值比其他值大几个数量级的矩阵。 例如，如果要处理偏微分方程或某些其他过程，其中定义的线性关系中的出现元素仅与“邻居”相连，则这是不可避免的。 这是经典物理学中一维泊松方程的稀疏矩阵的一个可能示例 $$$-\ phi ^ {''} = f$ 在一个统一的网格上（是的，虽然其中没有太多零，但是在研磨网格时它们会很健康！）：

$$

$$\开始{pmatrix} 2＆-1＆0＆0＆0 \\ -1＆2＆-1＆0＆0 \\ 0＆-1＆2＆-1＆0 \\ 0＆0＆-1 ＆2＆-1 \\ 0＆0＆0＆-1＆2 \ end {pmatrix}$$

与它们相对的矩阵-那些不注意零数并无一例外地考虑所有分量的矩阵称为稠密矩阵。

出于各种原因，稀疏矩阵通常以压缩列格式-压缩稀疏列（CSC）表示。 此描述使用两个整数数组和一个浮点数。 让矩阵拥有一切 $$$nnz_A$ 非零和 $$$N$ 列。 矩阵元素从左到右在列中列出，没有例外。 第一个数组 $$$i_A$ 长度 $$$nnz_A$ 包含非零矩阵组件的行号。 第二个数组 $$$j_A$ 长度 $$$N + 1$ 包含...号，而不是列号，因为这样矩阵将以坐标格式（coordinate）或三元（triplet）写入。 第二个数组包含列开始的那些矩阵组件的序列号，在末尾包括一个附加的伪列。 最后，第三个数组 $$$v_A$ 长度 $$$nnz_A$ 已经包含组件本身，并且顺序与第一个数组相对应。 例如，在此假设下，对于特定矩阵，矩阵的行和列的编号从零开始

$$

$$A = \开始{pmatrix} 0＆1＆0＆4＆-1 \\ 7＆0＆2.1＆0＆3 \\ 0＆0＆0＆10＆0 \ end {pmatrix}$$

数组将是 $$$i_A = \ {1，0，1，0，2，0，1 \}$ ， $$$j_A = \ {0，1，2，3，5，7 \}$ ， $$$v_A = \ {7，1，2.1，4，10，-1，3 \}$ 。

求解线性系统的方法分为两大类-直接和迭代。 直线基于将矩阵表示为两个更简单矩阵的乘积的可能性，以便随后将解决方案分为两个简单阶段。 迭代使用线性空间的独特属性，并在古代作为世界上逐步将未知数逐步逼近所需解的方法，并且在收敛过程中，通常仅使用矩阵将它们乘以矢量。 迭代方法比直接方法便宜得多，但是在条件差的矩阵上工作缓慢。 如果钢筋混凝土的可靠性很重要，请使用直线。 我在这里，我想谈一点。

假设方阵 $$我们可以查看 $$$A = LU$ 在哪里 $$$L$ 和 $$ -分别是下三角矩阵和上三角矩阵。 第一个表示它在对角线上方有一个零，第二个表示它在对角线下方。 我们究竟是如何得到这种分解的-我们现在不感兴趣。 这是一个简单的分解示例：

$$

$$\开始{pmatrix} 1＆-1＆-1 \\ 2＆-1＆-0.5 \\ 4＆-2＆-1.5 \ end {pmatrix} = \开始{pmatrix} 1＆0＆0 \\ 2 ＆1＆0 \\ 4＆2＆1 \结束{pmatrix} \开始{pmatrix} 1＆-1＆-1 \\ 0＆1＆1.5 \\ 0＆0＆-0.5 \ end {pmatrix}$$

在这种情况下如何求解方程组 $$$斧= f$ 例如，在右侧 $$$f = \开始{pmatrix} 4 \\ 2 \\ 3 \结束{pmatrix}$ ？ 第一阶段是直接行动（正解=正向替代）。 我们表示 $$$y：= Ux$ 并与系统一起工作 $$$Ly = f$ 。 因为 $$$L$ -下三角，在循环中我们依次找到所有组件 $$$y$ 从上到下：

$$

$$\开始{pmatrix} 1＆0＆0 \\ 2＆1＆0 \\ 4＆2＆1 \结束{pmatrix} \开始{pmatrix} y_1 \\ y_2 \\ y_3 \ end {pmatrix} = \开始{pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix} \表示y = \开始{pmatrix} 4 \\ -6 \\ -1 \ end {pmatrix}$$

中心思想是，一旦发现 $$$我$ 向量成分 $$$y$ 它乘以具有相同矩阵号的列 $$$L$ 然后从右侧减去。 矩阵本身似乎从左向右塌陷，随着发现越来越多的矢量分量，其尺寸减小 $$$y$ 。 此过程称为列消除。

第二阶段是向后移动（向后求解=向后替换）。 找到向量 $$$y$ 我们决定 $$$Ux = y$ 。 在这里，我们已经自下而上地研究了各个组件，但是想法仍然相同： $$$我$ 第列乘以刚刚找到的分量 $$$x_i$ 并向右移动，矩阵从右向左折叠。 下图示例中提到的矩阵说明了整个过程：

$$

$$\ small \开始{pmatrix} 1＆0＆0 \\ 2＆1＆0 \\ 4＆2＆1 \结束{pmatrix} \开始{pmatrix} y_1 \\ y_2 \\ y_3 \ end {pmatrix} = \开始{pmatrix} 4 \\ 2 \\ 3 \结束{pmatrix} \暗示\开始{pmatrix} 1＆0 \\ 2＆1 \结束{pmatrix} \开始{pmatrix} y_2 \\ y_3 \ end {pmatrix } = \开始{pmatrix} -6 \\ -13 \结束{pmatrix} \暗示\开始{pmatrix} 1 \结束{pmatrix} \开始{pmatrix} y_3 \ end {pmatrix} = \开始{pmatrix} -1 \ end {pmatrix}$$

$$

$$\ small \ begin {pmatrix} 1＆-1＆-1 \\ 0＆1＆1.5 \\ 0＆0＆-0.5 \ end {pmatrix} \ begin {pmatrix} x_1 \\ x_2 \\ x_3 \ end { pmatrix} = \开始{pmatrix} 4 \\ -6 \\ -1 \结束{pmatrix} \ Rightarrow \开始{pmatrix} 1＆-1 \\ 0＆1 \结束{pmatrix} \开始{pmatrix} x_1 \ \ x_2 \ end {pmatrix} = \ begin {pmatrix} 6 \\ -9 \ end {pmatrix} \隐含\ begin {pmatrix} 1 \ end {pmatrix} \ begin {pmatrix} x_1 \ end {pmatrix} = \ begin {pmatrix} -3 \ end {pmatrix}$$

我们的决定将是 $$$x = \开始{pmatrix} -3 \\-9 \\ 2 \结束{pmatrix}$ 。

如果矩阵是密集的，也就是说，它完全以某种数组的形式表示，一维或二维，并且随着时间的流逝访问其中的特定元素 $$$O（1）$ ，那么采用现有分解方法进行类似的求解过程并不困难，而且易于编码，因此我们甚至不会在上面浪费时间。 如果矩阵稀疏怎么办？ 原则上，这里也不难。 这是解决方案前进的C ++代码 $$$x$ 它写在右侧的位置，而无需检查输入数据的正确性（CSC数组对应于下三角矩阵）：

算法1：

void forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x) { size_t j, p, n = x.size(); for (j = 0; j < n; j++) //    { x[j] /= vA[jA[j]]; for (p = jA[j]+1; p < jA[j+1]; p++) x[iA[p]] -= vA[p] * x[j] ; } } 

所有进一步的讨论将只涉及解决下三角系统的正向行程 $$$Lx = f$ 。

领带

但是如果右边是 等号右边的向量 $$$Lx = f$ 它本身是否有大量的零？ 这样就可以跳过与零位相关的计算。 在这种情况下，代码中的更改很小：

算法2：

 void fake_sparse_forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x) { size_t j, p, n = x.size(); for (j = 0; j < n; j++) //    { if(x[j]) { x[j] /= vA[jA [j]]; for (p = jA[j]+1; p < jA[j+1]; p++) x[iA[p]] -= vA[p] * x[j]; } } } 

我们添加的唯一内容是if ，其目的是将算术运算的数量减少到实际数量。 例如，如果整个右侧都由零组成，那么您无需计数：解决方案将在右侧。

一切看起来都很不错，当然可以运行，但是经过长时间的前戏之后，问题终于显现了出来-对于大型系统，该求解器的渐近性能低下。 都是因为有一个for循环的事实。 怎么了 即使if内的条件极少为真，也无法摆脱循环本身，这造成了算法的复杂性 $$$O（N）$ 在哪里 $$$N$ 是矩阵的大小。 不管现代编译器如何优化循环，这种复杂性都会使自己感觉到很大 $$$N$ 。 当整个矢量都特别令人反感 $$$f$ 完全由零组成，因为正如我们所说，在这种情况下什么也不做！ 没有一个算术运算！ 什么鬼 $$$O（N）$ 行动吗？

好吧，让我们说。 即便如此，您为什么不能忍受for run idle，因为带有实数的实际计算（即 那些落在下面的人还会少吗？ 但是事实是，这种右手边本身很少的直流程序经常在外部循环中使用，并构成了Cholesky分解的基础 $$$A = LL ^ T$ 和左眼LU分解 。 是的，是的，其中的一种扩展无法解决线性系统中的所有这些正向和反向运动，从而失去了所有实际意义。

定理 如果矩阵是对称正定（SPD），则可以表示为 $$$A = LL ^ T$ 唯一的方法 $$$L$ -在对角线上带有正元素的下三角矩阵。

对于高度稀疏的SPD矩阵，使用前瞻性的Cholesky分解。 以矩阵块形式示意性表示分解

$$

$$\开始{pmatrix} \ mathbf {L} _ {11}＆\ mathbf {0} \\ \ mathbf {l} _ {12} ^ T＆l_ {22} \ end {pmatrix} \ begin {pmatrix} \ mathbf {L} _ {11} ^ T＆\ mathbf {l} _ {12} \\ \ mathbf {0}＆l_ {22} \ end {pmatrix} = \开始{pmatrix} \ mathbf {A} _ { 11}＆\ mathbf {a} _ {12} \\ \ mathbf {a} _ {12} ^ T＆a_ {22} \ end {pmatrix}，$$

整个分解过程在逻辑上可以分为三个步骤。

算法3：

1. 较小维的Cholesky分解的上部方法 $$$\ mathbf {L} _ {11} \ mathbf {L} _ {11} ^ T = \ mathbf {A} _ {11}$ （递归！）
2. 下三角系统，右侧稀疏 $$$\ mathbf {L} _ {11} \ mathbf {l} _ {12} = \ mathbf {a} _ {12}$ （ 在这里！ ）
3. 计算 $$$l_ {22} = \ sqrt {a_ {22}-\ mathbf {l} _ {12} ^ T \ mathbf {l} _ {12}}$ （琐碎的操作）

实际上，这是以这样的方式实现的，即以一个大周期并以此顺序执行步骤3和2。 因此，矩阵从上到下对角延伸，从而增加了矩阵 $$$L$ 在循环的每次迭代中逐行。

如果在这种算法中，步骤2中的前向复杂度将为 $$$O（N）$ 在哪里 $$$N$ 下三角矩阵的大小是 $$$\ mathbf {L} _ {11}$ 在大循环的任意迭代中，整个分解的复杂度至少为 $$$O（N ^ 2）$ ！ 哦，我不会那样的！

最先进的

许多算法都以模仿人类解决问题的行动为基础。 如果您给一个人带右手的下三角线性系统，其中只有1-2是非零的，那么他当然首先会用眼睛从上到下运行右手的向量（那该死的复杂性循环） $$$O（N）$ ！ ）找到这些非零值。 然后他将只使用它们，而不会在零分量上浪费时间，因为后者不会影响决策：将零划分为矩阵的对角线分量以及将乘以零的列向右移动都没有意义。 这就是上面的算法2。没有奇迹。 但是，如果一个人立即从其他来源获得非零分量的数目怎么办？ 例如，如果右侧是其他矩阵的一列（例如Cholesky分解的情况），那么如果按顺序请求，我们可以立即访问其非零分量：

 //     j,     : for (size_t p = jA[j]; p < jA[j+1]; p++) //    vA[p]     iA[p]   j 

这种访问的复杂性是 $$$O（nnz_j）$ 在哪里 $$$nnz_j$ 是j列中非零分量的数量。 感谢上帝提供的CSC格式！ 如您所见，它不仅用于节省内存。

在这种情况下，我们可以掌握直接过程方法求解时发生的本质 $$$Lx = f$ 用于稀疏矩阵 $$$L$ 和右边 $$$f$ 。 屏住呼吸： 我们采取非零成分 $$$f_i$ 在右侧，我们找到相应的变量 $$$x_i$ 除以 $$$L_ {ii}$ 然后，将整个第i列乘以该找到的变量， 我们通过从右侧减去该列，在右侧引入其他非零值 ！ 用...图形的语言很好地描述了此过程。 而且，定向且非循环。

对于在对角线上没有零的较低三角形矩阵，我们定义一个连通图。 我们假设行和列的编号从零开始。
定义 下三角矩阵的连通图 $$$L$ 大小 $$$N$ 对角线上没有零的称为节点集 $$$V = \ {0，...，N-1 \}$ 和定向的边缘 $$$E = \ {（j，i）| L_ {ij} \ ne 0 \}$ 。

例如，在这里是下三角矩阵的连接图。

其中数字仅对应于对角线元素的序数，而点表示对角线以下的非零分量：



定义 可达性图 $$$G$ 在多个指数上 $$$W \子集V$ 叫那么多指数 $$$R_G（W）\子集V$ 在任何索引中 $$$z \在R_G（W）$$ 您可以通过从某些索引中遍历图来获得 $$$w（z）\以W$ 。

示例：来自图片的图形 $$$R_G（\ {0，4 \}）= \ {0，4，5，6 \}$ 。 很明显，它总是成立 $$$W \子集R_G（W）$ 。

如果我们将图的每个节点表示为生成它的矩阵的列号，则其边沿指向的相邻节点对应于该列中非零分量的行号。

让 $$$nnz（x）= \ {j | x_j \ ne 0 \}$ 表示对应于向量x中非零位置的一组索引。

希尔伯特定理 （否，不是用空格来命名） $$$nnz（x）$ 在哪里 $$$x$ 有一个稀疏的下三角系统的解向量 $$$Lx = f$ 右侧稀疏 $$$f$ 符合 $$$R_G（nnz（f））$ 。

我们自己加法：在定理中，我们没有考虑不太可能的可能性，即销毁列时右侧的非零数字可以减小为零，例如3-3 =0。这种效应称为数值抵消。 考虑到这些自发出现的零是浪费时间，并且将它们视为处于非零位置的所有其他数字。

假设我们可以通过索引直接访问其非零分量，则在给定的稀疏右手侧进行直接移动的有效方法如下。

1. 我们使用“深度优先搜索”方法遍历图形，从索引开始 $$$i \ in nnz（f）$ 右侧的每个非零分量。 将找到的节点写入数组 $$$R_G（nnz（f））$ 同时按照我们返回图表的顺序进行操作！ 类似于入侵者的军队：我们没有残酷地占领了这个国家，但是当他们开始把我们赶回来时，我们愤怒地摧毁了前进道路上的一切。
   值得注意的是，索引列表绝对没有必要 $$$nnz（f）$ 当它被输入到“深度优先搜索”算法的输入时，按升序排序。 您可以按任意顺序开始拍摄 $$$nnz（f）$ 。 属于集合的不同序列 $$$nnz（f）$ 索引不会影响最终的解决方案，如示例所示。
   
   此步骤不需要任何实际数组的知识。 $$$v_A$ ！ 欢迎使用直接稀疏求解器进行符号分析的世界！
   
2. 我们着手解决方案本身，拥有一个阵列 $$$R_G（nnz（f））$ 从上一步开始。 我们以与该数组的记录相反的顺序销毁这些列。 瞧！

例子

考虑一个示例，其中详细介绍了两个步骤。 假设我们具有以下形式的12 x 12矩阵：

相应的连接图具有以下形式：

令右侧的非零值仅位于位置3和5，即 $$$nnz（f）= \ {3，5 \}$ 。 让我们遍历图表，从这些索引中按书面顺序开始。 然后，“深度搜索”方法将如下所示。 首先，我们将按顺序访问索引 $$ ，不要忘记将索引标记为已访问。 在下图中，它们以蓝色阴影显示：

返回时，我们将索引输入到非零解分量的索引数组中 $$$nnz（x）：= \ {\颜色{blue} {11}，\颜色{blue} 8，\颜色{blue} 3 \}$ 。 接下来，尝试运行 $$$5 \ to8 \至...$ ，但是遇到了一个已经标记为8的节点，因此我们不会碰这条路线并转到分支 $$$5 \ to9 \至...$ 。 运行的结果将是 $$$5 \ to9 \ to10$ 。 我们无法访问节点11，因为它已被标记。 因此，返回并补充数组 $$$nnz（x）$ 绿色标记的新渔获物： $$$nnz（x）：= \ {\颜色{蓝色} {11}，\颜色{蓝色} 8，\颜色{蓝色} 3，\颜色{绿色} {10}，\颜色{绿色} 9，\颜色{green} 5 \}$ 。 这是图片：

我们想在第二次运行中访问那些浑浊的绿色节点8和11，但是不能，因为 已经在第一时间访问过。

所以数组 $$$nnz（x）= R_G（nnz（f））$ 形成。 我们进入第二步。 一个简单的步骤：我们遍历数组 $$$nnz（x）$ 以相反的顺序（从右到左），找到解矢量的相应分量 $$$x$ 除以矩阵的对角线分量 $$$L$ 并将列向右移动。 其他组件 $$$x$ 由于它们为零，因此保持不变。 示意图如下所示，底部箭头指示销毁列的顺序：

注意：按照破坏列的顺序，编号3将比编号5、9和10更晚！ 列不是按升序销毁的，这对于密集的列是错误的，即 未完成的矩阵。 但是对于稀疏矩阵，这是司空见惯的。 从本例中的非零矩阵结构可以看出，对第3列使用第5、9和10列不会使响应中的分量失真 $$$x_5$ ， $$$x_9$ ， $$$x_ {10}$ 根本没有，他们有c $$$x_3$ 只是“没有路口”。 我们的方法考虑了这一点！ 同样，在第9列之后使用第8列也不会破坏组件 $$$x_9$ 。 很棒的算法，不是吗？

如果我们先从索引5开始遍历图，然后再从索引3开始遍历，那么我们的数组将 $$$nnz（x）：= \ {\颜色{blue} {11}，\颜色{blue} 8，\颜色{blue} {10}，\颜色{blue} {9}，\颜色{blue} 5， \ color {green} 3 \}$ 。 以与该数组相反的顺序销毁这些列将提供与第一种情况完全相同的解决方案！

我们所有运算的复杂性都取决于实际算术运算的数量和右侧非零分量的数量，而不是矩阵的大小！ 我们已经实现了目标。

批评

但是！ 关键的读者会注意到，一开始我们就可以“直接访问索引右侧的非零分量”这一假设已经意味着，在我们从上到下运行右侧以找到这些索引并组织它们之前在数组中 $$$nnz（f）$ ，也就是说，已经花了 $$$O（N）$ 行动！ 此外，图运行本身要求我们首先分配最大可能长度的内存（在其他地方，您需要写下通过深度搜索找到的索引！），以免随着数组的增长而遭受进一步的重新分配 $$$nnz（x）$ ，这也需要 $$$O（N）$ 操作！ 他们为什么这么大惊小怪？

但是实际上，对于最初定义为密集向量的，具有稀疏右手侧的稀疏下三角系统的一次性解决方案，在上述所有算法上浪费开发人员时间是没有意义的。 它们甚至可能比上面算法2提出的额头方法还要慢。 但是，正如已经提到的那样，此设备对于Cholesky的因式分解是必不可少的，因此您不应向我扔西红柿。 实际上，在运行算法3之前，所有必要的最大长度的内存都会立即分配并需要 $$$O（N）$ 时间。 在进一步的色谱柱循环中 $$所有数据仅以固定长度的数组覆盖 $$$N$ ，并且仅在与重写有关的位置上，这要归功于直接访问非零元素。 正是因为如此，效率才得以提高！

C ++实现

图形本身作为代码中的数据结构不是必需的。 直接使用矩阵时，隐式使用它就足够了。 该算法将考虑所有必需的连接性。 深度搜索当然可以使用平凡递归方便地实现。 例如，此处看起来像是基于单个起始索引的递归深度搜索：

 /* j -   iA, jA -    ,    CSC top -     nnz(x);     0 result -   N,    nnz(x)   0  top-1  marked -    N;      */ void DepthFirstSearch(size_t j, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t& top, std::vector<size_t>& result, std::vector<int>& marked) { size_t p; marked[j] = 1; //   j   for (p = jA[j]; p < jA[j+1]; p++) //       j { if (!marked[iA[p]]) //  iA[p]   { DepthFirstSearch(iA[p], iA, jA, top, result, marked); //      iA[p] } } result[top++] = j; //  j   nnz(x) } 

如果在第一次调用DepthFirstSearch将top变量传递为零，则在完成整个递归过程后， top变量将等于在result数组中找到的索引数。 作业：编写另一个函数，该函数在右侧采用一组非零成分的索引，并将它们顺序传递给DepthFirstSearch函数。 没有这个，算法就不完整。 请注意：实数数组 $$$vA$ 我们根本不会将其传递给函数，因为 在符号分析过程中不需要它。

尽管它很简单，但实现上的缺陷却很明显：对于大型系统，堆栈溢出指日可待。 那么，有一个基于循环而不是递归的选项。 很难理解，但已经适合所有场合：

 /* j -   N -   iA, jA -    ,    CSC top -     nnz(x) result -   N,     nnz(x)   top  ret-1 ,  ret -    marked -    N work_stack -     N */ size_t DepthFirstSearch(size_t j, size_t N, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t top, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack) { size_t p, head = N - 1; int done; result[N - 1] = j; //    while (head < N) { j = result[head]; //  j     if (!marked[j]) { marked[j] = 1; //   j   work_stack[head] = jA[j]; } done = 1; //    j      for (p = work_stack[head] + 1; p < jA[j+1]; p++) //    j { //   c   iA[p] if (marked[iA[p]]) continue; //  iA[p]  ,   work_stack[head] = p; //       j result[--head] = iA[p]; //       iA[p] done = 0; //   j    break; } if (done) //      j  { head++; //  j    result[top++] = j; //  j    } } return (top); } 

这就是非零解矢量结构的生成器已经看起来像的样子 $$$nnz（x）$ ：

 /* iA, jA -   ,    CSC iF -        f nnzf -      f result -   N,    nnz(x)   0  ret-1 ,  ret -    marked -    N,    work_stack -     N */ size_t NnzX(const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<size_t>& iF, size_t nnzf, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack) { size_t p, N, top; N = jA.size() - 1; top = 0; for (p = 0; p < nnzf; p++) if (!marked[iF[p]]) //        top = DepthFirstSearch(iF[p], N, iA, jA, vA, top, result, marked, work_stack); for (p = 0; p < top; p++) marked[result[p]] = 0; //   return (top); } 

最后，将所有内容合而为一，我们针对稀有右手边的情况编写了下三角求解器：

 /* iA, jA, vA -  ,    CSC iF -        f nnzf -      f vF -       f result -   N,    nnz(x)   0  ret-1 ,  ret -    marked -    N,    work_stack -     N x -    N,    ;     */ size_t lower_triangular_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, const std::vector<size_t>& iF, size_t nnzf, const std::vector<double>& vF, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack, std::vector<double>& x) { size_t top, p, j; ptrdiff_t px; top = NnzX(iA, jA, iF, nnzf, result, marked, work_stack); for (p = 0; p < nnzf; p++) x[iF[p]] = vF[p]; //    for (px = top; px > -1; px--) //     { j = result[px]; // x(j)   x [j] /= vA[jA[j]]; //   x(j) for (p = jA[j]+1; p < jA[j+1]; p++) { x[iA[p]] -= vA[p]*x[j]; //  j-  } } return (top) ; } 

我们看到我们的循环仅遍历数组的索引 $$$nnz（x）$ ，而不是整套 $$$0，1，...，N-1$ 。 做完了！

有一种实现不使用标签数组marked来节省内存。相反，使用了另外一组索引。$$$V_1$ 不与布景相交 $$$V = \{0, ..., N - 1\}$通过简单的代数运算将其中的一对一映射作为节点标记过程。但是，在我们廉价的内存时代，将其保存在单个长度的数组中$$$N$ 似乎完全多余。

总结

通常，通过直接方法求解线性方程组的稀疏系统的过程分为三个阶段：

1. 人物分析
2. 基于sivol分析数据的数值分解
3. 所获得的右手边密集的三角系统的解

第二步，数值分解，是最消耗资源的部分，消耗了估计时间的大部分（> 90％）。第一步的目的是减少第二步的高成本。这篇文章中提供了一个符号分析的示例。但是，这是第一步，需要开发人员最长的开发时间和最大的精神成本。好的符号分析需要图形和树木的理论知识以及“算法本能​​”的掌握。第二步易于实施。
好吧，第三步，无论是在实施上还是在估计的时间内，在大多数情况下都是最朴实的。

德克萨斯A＆M大学计算机系教授Tim Davis很好地介绍了稀疏系统的直接方法，题目是““ 稀疏线性系统的直接方法 ”。此处描述了上述算法。虽然我本人曾经在同一所大学（尽管在不同的系中工作），但我对Davis并不熟悉。如果我没记错的话，Davis亲自参与了求解器的开发（！）。在Matlab使用。此外，他参与甚至顺便说在谷歌街景（OLS）发电机的画面，在比自己同教员上市莫属了Bjarne Stroustrup的- .. C ++的创造者

也许有一天我将在其中写有关稀疏矩阵或数值方法的其他内容 一般。对所有人都好！