C#，数值计算（Numerical Recipes in C#），大型稀疏线性系统（Sparse Linear Systems）的数据存储结构与源代码

《Numerical Recipes in C++》原文摘要：

A system of linear equations is called sparse if only a relatively small number of its matrix elements aij are nonzero. It is wasteful to use general methods of linear algebra on such problems, because most of the O.N 3/ arithmetic operations devoted to solving the set of equations or inverting the matrix involve zero operands. Furthermore, you might wish to work problems so large as to tax your available memory space, and it is wasteful to reserve storage for unfruitful zero elements. Note that there are two distinct (and not always compatible) goals for any sparse matrix method: saving time and/or saving space.

We considered one archetypal sparse form in §2.4, the band-diagonal matrix. In the tridiagonal case, e.g., we saw that it was possible to save both time (order N instead of N 3) and space (order N instead of N 2). The method of solution was not different in principle from the general method of LU decomposition; it was just applied cleverly, and with due attention to the bookkeeping of zero elements. Many practical schemes for dealing with sparse problems have this same character. They are fundamentally decomposition schemes, or else elimination schemes akin to Gauss-Jordan, but carefully optimized so as to minimize the number of so-called

fill-ins, initially zero elements that must become nonzero during the solution process, and for which storage must be reserved.

Direct methods for solving sparse equations, then, depend crucially on the pre- cise pattern of sparsity of the matrix. Patterns that occur frequently, or that are useful as way stations in the reduction of more general forms, already have special names and special methods of solution. We do not have space here for any detailed review of these. References listed at the end of this section will furnish you with an “in” to the specialized literature, and the following list of buzz words (and Figure 2.7.1) will at least let you hold your own at cocktail parties:

tridiagonal
band-diagonal (or banded) with bandwidth M
band triangular
block diagonal
block tridiagonal
block triangular
cyclic banded
singly (or doubly) bordered block diagonal
singly (or doubly) bordered block triangular
singly (or doubly) bordered band-diagonal
singly (or doubly) bordered band triangular
other (!)

You should also be aware of some of the special sparse forms that occur in the solu- tion of partial differential equations in two or more dimensions. See Chapter 20.

If your particular pattern of sparsity is not a simple one, then you may wish to try an analyze/factorize/operate package, which automates the procedure of figuring out how fill-ins are to be minimized. The analyze stage is done once only for each pattern of sparsity. The factorize stage is done once for each particular matrix that fits the pattern. The operate stage is performed once for each right-hand side to be used with the particular matrix. Consult [2,3] for references on this. The NAG library [4] has an analyze/factorize/operate capability. A substantial collection of routines for sparse matrix calculation is also available from IMSL [5] as the Yale Sparse Matrix Package [6].

You should be aware that the special order of interchanges and eliminations, prescribed by a sparse matrix method so as to minimize fill-ins and arithmetic op- erations, generally acts to decrease the method’s numerical stability as compared to, e.g., regular LU decomposition with pivoting. Scaling your problem so as to make its nonzero matrix elements have comparable magnitudes (if you can do it) will sometimes ameliorate this problem.

In the remainder of this section, we present some concepts that are applicable to some general classes of sparse matrices, and which do not necessarily depend on details of the pattern of sparsity.

凑字数的狗屁不通的译文：

如果线性方程组中只有相对较少的矩阵元素aij非零，则称其为稀疏。在此类问题上使用线性代数的一般方法是浪费的，因为大多数用于求解方程组或矩阵求逆的O.N 3/算术运算都涉及零操作数。此外，您可能希望处理如此大的问题，以至于占用可用内存空间，而为无效的零元素保留存储空间是浪费。请注意，任何稀疏矩阵方法都有两个不同（但并不总是兼容）的目标：节省时间和/或节省空间。

我们考虑了§2.4中的一种原型稀疏形式，带对角矩阵。例如，在三对角情况下，我们发现可以节省时间（N阶而不是N 3）和空间（N阶而不是N 2）。求解方法原则上与LU分解的一般方法没有区别；它只是巧妙地应用，并适当注意零元素的簿记。许多处理稀疏问题的实用方案都具有相同的特征。它们基本上是分解方案，或者类似于高斯-乔丹的消除方案，但经过仔细优化，以最小化所谓的

填充，最初为零元素，在求解过程中必须变为非零，并且必须为其保留存储。

因此，求解稀疏方程的直接方法在很大程度上取决于矩阵的稀疏性的精确模式。频繁出现的模式，或在简化更一般形式时用作中途站的模式，已经有了特殊的名称和特殊的求解方法。我们这里没有空间对这些进行任何详细审查。本节末尾列出的参考文献将为您提供专业文献的“in”，以下流行语列表（和图2.7.1）至少可以让您在鸡尾酒会上保持自己的风格：

●对角

●带宽为M的带对角（或带状）

●带状三角形

●块对角线

●块三对角

●块三角形

●环状带状

●单（或双）边界块对角

●单（或双）边界块三角形

●单（或双）边带对角

●单（或双）边带三角形

●其他（！）

您还应该了解在二维或多维偏微分方程求解中出现的一些特殊稀疏形式。见第20章。

如果您的特定稀疏模式不是一个简单的模式，那么您可能希望尝试一个分析/因子分解/操作包，该包自动计算如何最小化填充的过程。分析阶段仅对每个稀疏模式执行一次。分解阶段对符合模式的每个特定矩阵执行一次。对于要与特定矩阵一起使用的每个右侧，执行一次操作阶段。参考[2,3]。NAG库[4]具有分析/因子分解/操作能力。IMSL[5]还提供了大量稀疏矩阵计算例程，作为耶鲁稀疏矩阵包[6]。

您应该知道，稀疏矩阵方法规定的特殊交换和消去顺序（以最小化填充和算术运算）通常会降低该方法的数值稳定性，例如，与旋转的规则LU分解相比。缩放问题，使其非零矩阵元素具有可比的大小（如果可以的话），有时会改善此问题。

在本节的其余部分中，我们提出了一些适用于某些一般稀疏矩阵类的概念，这些概念不一定取决于稀疏模式的细节。

源程序（POWER BY 315SOFT.COM）

/02 Solution of Linear Algebraic Equations/NRsparseCol.cs

using System;namespace Legalsoft.Truffer
{/// <summary>/// Sparse vector data structure./// </summary>public class NRsparseCol{/// <summary>/// Number of rows./// </summary>public int nrows { get; set; } = 0;/// <summary>/// Maximum number of nonzeros./// </summary>public int nvals { get; set; } = 0;/// <summary>/// Row indices of nonzeros./// </summary>public int[] row_ind { get; set; }/// <summary>/// Array of nonzero values./// </summary>public double[] val { get; set; }public NRsparseCol(){}public NRsparseCol(int m, int nnvals){this.nrows = m;this.nvals = nnvals;this.row_ind = new int[nnvals];this.val = new double[nnvals];}public void resize(int m, int nnvals){nrows = m;nvals = nnvals;row_ind = new int[nnvals];val = new double[nnvals];}}
}

/02 Solution of Linear Algebraic Equations/NRsparseMat.cs

using System;namespace Legalsoft.Truffer
{/// <summary>/// Sparse matrix data structure for compressed column storage./// </summary>public class NRsparseMat{/// <summary>/// Number of rows./// </summary>public int nrows { get; set; } = 0;/// <summary>/// Number of columns./// </summary>public int ncols { get; set; } = 0;/// <summary>/// Maximum number of nonzeros./// </summary>public int nvals { get; set; } = 0;/// <summary>/// Pointers to start of columns. Length is ncols+1./// </summary>public int[] col_ptr { get; set; }/// <summary>/// Row indices of nonzeros/// </summary>public int[] row_ind { get; set; }/// <summary>/// Array of nonzero values./// </summary>public double[] val { get; set; }public NRsparseMat(){}public NRsparseMat(int m, int n, int nnvals){this.nrows = m;this.ncols = n;this.nvals = nnvals;this.col_ptr = new int[n + 1];this.row_ind = new int[nnvals];this.val = new double[nnvals];}public double[] ax(double[] x){double[] y = new double[nrows];for (int j = 0; j < ncols; j++){for (int i = col_ptr[j]; i < col_ptr[j + 1]; i++){y[row_ind[i]] += val[i] * x[j];}}return y;}public double[] atx(double[] x){double[] y = new double[ncols];for (int i = 0; i < ncols; i++){y[i] = 0.0;for (int j = col_ptr[i]; j < col_ptr[i + 1]; j++){y[i] += val[j] * x[row_ind[j]];}}return y;}public NRsparseMat transpose(){int m = nrows;int n = ncols;NRsparseMat at = new NRsparseMat(n, m, nvals);int[] count = new int[m];for (int i = 0; i < n; i++){for (int j = col_ptr[i]; j < col_ptr[i + 1]; j++){int k = row_ind[j];count[k]++;}}for (int j = 0; j < m; j++){at.col_ptr[j + 1] = at.col_ptr[j] + count[j];}for (int j = 0; j < m; j++){count[j] = 0;}for (int i = 0; i < n; i++){for (int j = col_ptr[i]; j < col_ptr[i + 1]; j++){int k = row_ind[j];int index = at.col_ptr[k] + count[k];at.row_ind[index] = i;at.val[index] = val[j];count[k]++;}}return at;}}
}

C#，数值计算（Numerical Recipes in C#），大型稀疏线性系统（Sparse Linear Systems）的数据存储结构与源代码相关推荐

大型稀疏线性方程组求解技术——工业仿真的底层核心
背景在工业仿真领域,对各种现实世界的问题进行数值模拟时,如流体动力学分析.电磁场仿真.结构力学应力应变分析等,其控制方程通常是偏微分方程组,在经过不同方法的隐式离散之后最终都可转化为大型稀疏线性 ...
c++稀疏表sparse table的实现算法(附完整源码)
C++稀疏表sparse table的实现算法 C++稀疏表sparse table的实现算法完整源码(定义,实现,main函数测试) C++稀疏表sparse table的实现算法完整源码(定义,实 ...
08_clickhouse主键/索引的工作机制(MergeTree的稀疏索引、索引的生成过程、索引的查询)，数据标记的工作机制（数据存储、数据标记、数据查询、数据查询示例）（学习笔记）
5.主键/索引的工作机制 5.1.MergeTree的稀疏索引 5.2.索引的生成过程 5.3.索引的查询 6.数据标记的工作机制 6.1.数据存储 6.2.数据标记 6.3.数据查询 6.4.数据查 ...
cuSPARSE库：（十四）求解稀疏三角形线性系统（solution of sparse triangular linear systems）
在cuSPARSE中,求解稀疏三角形线性系统(solution of sparse triangular linear systems)分两步实现: (1)分析阶段:调用函数csrsv_analysi ...
稀疏性(sparse)知识点
稀疏性(sparse) 定义:Sparse表示为模型内的参数中,只用很少的几个非零元素或只有很少的几个远大于零的元素. WHY: 为什么模型中要包含稀疏性的特征呢? 例子:考研学霸有10000的词汇量 ...
Eigen求解大型稀疏对称矩阵（Cholesky分解）
参考自Eigen文档代码如下: #include <Eigen/Sparse>typedef Eigen::SparseMatrix<double> SpMat; typed ...
稀疏编码(Sparse Coding)（二）
为了更进一步的清晰理解大脑皮层对信号编码的工作机制(策略),需要把他们转成数学语言,因为数学语言作为一种严谨的语言,可以利用它推导出期望和要寻找的程式.本节就使用概率推理(bayes views)的方 ...
live messenger与稀疏文件—Sparse File Bit
今天进行磁盘整理,发现一个奇怪的文件SimilarityTable_1:下面是我的C盘整理后的结果卷 (C:) 卷的大小 ...
python多维数据存储_在Python中存储和重新加载大型多维数据集
我将运行大量的模拟,产生大量的数据,这些数据需要在以后存储和访问.我的模拟程序的输出数据被写入文本文件(每个模拟一个).我计划编写一个Python程序来读取这些文本文件,然后以更便于以后分析的格式存储 ...

C#，数值计算（Numerical Recipes in C#），大型稀疏线性系统（Sparse Linear Systems）的数据存储结构与源代码

C#，数值计算（Numerical Recipes in C#），大型稀疏线性系统（Sparse Linear Systems）的数据存储结构与源代码相关推荐

最新文章

热门文章