【Eigen】Chapter3 稀疏线性代数 Sparse Linear Algebra

（1）稀疏矩阵操作

​ 1）稀疏矩阵格式

​ 在许多应用中，矩阵其只有少量非零系数，这样的矩阵称为稀疏矩阵(Sparse Matrix)。在这种情况下，采用稀疏矩阵的存储方式，即仅存储非零系数，可以减少内存消耗并提高性能。

​ 在Eigen中，SparseMatrix<>模板类是用于表示稀疏矩阵，注意稠密的用的是Matrix<>，SparseMatrix<>提供了高性能和低内存使用率。它实现了广泛使用的压缩列（或行）存储方案的更通用的变体。

​ 它由四个紧凑数组组成：

​ 1 Values: 存储非零的系数值。注意，压缩列存储按照一列一列的顺序存储非零元素，同样，压缩行按照一行一行的顺序存储元素；

​ 2 InnerIndices: 存储非零的行（或列）下标。简单的说，我们可以先按照列或者行的顺序存储values，然后对于values中的每一个元素，确定其所在的行or列的下标；

​ 3 OuterStarts: 为每列（分别为行）存储前两个数组中第一个非零的索引。注意，这个数组存储的Values和InnerIndices中的位置；

​ 4 InnerNNZs:存储每列（分别为行）的非零元素个数。

​ 注意：术语inner指的是内部向量，该向量是列主矩阵的列，或行主矩阵的行。术语outer是指另一个方向，即行or列；

​ 在一个示例中可以更好地解释此存储方案。以下矩阵

// 0	3	0	0	0
// 22	0	0	0	17
// 7	5	0	1	0
// 0	0	0	0	0
// 0	0	14	0	8

​ and one of its possible sparse, column major representation:

Values: 		22	7	_	3	5	14	_	_	1	_	17	8
InnerIndices: 	1	2	_	0	2	4	_	_	2	_	1	4 

OuterStarts:	0	3	5	8	10	12
InnerNNZs: 		2	2	1	1	2

​ 注释：OutStarts表示的每一列中的起始非零元素在Values or InnerIndices中的索引，例如,第一列非零元素22在Values中的第0索引，第二列中的3在Values中的第3索引，第三列中的14在Values中的第5索引…

​ OutStarts多了一个元素，比如上述方阵有5列，那么应该只有5个元素，但是它有6个元素，难道是存储了最后一个元素？

​ 2）第一个例子

​ 在描述每个类之前，让我们从以下典型示例开始：

​ Δ u = 0 \Delta u = 0 Δu=0使用有限差分方案和Dirichlet边界条件在规则2D网格上求解拉普拉斯方程。这样的问题可以在数学上表示为Ax = b形式的线性问题，其中x是一个m维向量，在本例中，它表示像素值, b是从边界条件中获得，并且A是一个mXm的方阵，它是从离散化的拉普拉斯算子中获得的稀疏矩阵，即，只含有少数的非0元素。

​ （说实话，有点复杂）

#include <Eigen/Sparse>
#include <vector>
#include <iostream>
 
typedef Eigen::SparseMatrix<double> SpMat; // declares a column-major sparse matrix type of double
typedef Eigen::Triplet<double> T;
 
void buildProblem(std::vector<T>& coefficients, Eigen::VectorXd& b, int n);
void saveAsBitmap(const Eigen::VectorXd& x, int n, const char* filename);
 
int main(int argc, char** argv)
{
  if(argc!=2) {
    std::cerr << "Error: expected one and only one argument.\n";
    return -1;
  }
  
  int n = 300;  // size of the image
  int m = n*n;  // number of unknowns (=number of pixels)
 
  // Assembly:
  std::vector<T> coefficients;            // list of non-zeros coefficients
  Eigen::VectorXd b(m);                   // the right hand side-vector resulting from the constraints
  buildProblem(coefficients, b, n);
 
  SpMat A(m,m);
  A.setFromTriplets(coefficients.begin(), coefficients.end());
 
  // Solving:
  Eigen::SimplicialCholesky<SpMat> chol(A);  // performs a Cholesky factorization of A
  Eigen::VectorXd x = chol.solve(b);         // use the factorization to solve for the given right hand side
 
  // Export the result to a file:
  saveAsBitmap(x, n, argv[1]);
 
  return 0;
}

​ 3）稀疏矩阵类

SparseMatrix<std::complex<float> > mat(1000,2000);         // declares a 1000x2000 column-major compressed sparse matrix of complex<float>
SparseMatrix<double,RowMajor> mat(1000,2000);              // declares a 1000x2000 row-major compressed sparse matrix of double
SparseVector<std::complex<float> > vec(1000);              // declares a column sparse vector of complex<float> of size 1000
SparseVector<double,RowMajor> vec(1000);                   // declares a row sparse vector of double of size 1000

// 维数可以使用下面的方法获得：
// Standard    dimensions，标准的维度
// mat.rows()
// mat.cols()
// vec.size()

// Sizes along the  inner/outer dimensions，内部和外部维度的大小
// mat.innerSize()
// mat.outerSize()

// Number of non zero coefficients，非零系数的个数
// mat.nonZeros()
// vec.nonZeros()

​ 迭代非零系数:

​ 随机访问的稀疏对象的元素可以通过来完成coeffRef(i,j)功能。但是，此功能涉及相当昂贵的二进制搜索。在大多数情况下，一个人只想遍历非零元素。这是通过在外部尺寸上进行标准循环，然后通过InnerIterator对当前内部向量的非零值进行迭代来实现的。因此，必须以与存储顺序相同的顺序访问非零条目。这是一个例子：

SparseMatrix<double> mat(rows,cols);
for (int k=0; k<mat.outerSize(); ++k)
{
  for (SparseMatrix<double>::InnerIterator it(mat,k); it; ++it)
  {
    it.value();
    it.row();   // row index
    it.col();   // col index (here it is equal to k)
    it.index(); // inner index, here it is equal to it.row()
  }
}

（2）求解稀疏矩阵

​ 用到再说把。。。。。。