(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)求解稀疏矩阵
用到再说把。。。。。。