<tr><td>\link Sparse_Module SparseCore \endlink</td><td>\code#include <Eigen/SparseCore>\endcode</td><td>SparseMatrix and SparseVector classes, matrix assembly, basic sparse linear algebra (including sparse triangular solvers)</td></tr>
<tr><td>\link SparseCholesky_Module SparseCholesky \endlink</td><td>\code#include <Eigen/SparseCholesky>\endcode</td><td>Direct sparse LLT and LDLT Cholesky factorization to solve sparse self-adjoint positive definite problems</td></tr>
<tr><td>\link IterativeLinearSolvers_Module IterativeLinearSolvers \endlink</td><td>\code#include <Eigen/IterativeLinearSolvers>\endcode</td><td>Iterative solvers to solve large general linear square problems (including self-adjoint positive definite problems)</td></tr>
<tr><td></td><td>\code#include <Eigen/Sparse>\endcode</td><td>Includes all the above modules</td></tr>
In many applications (e.g., finite element methods) it is common to deal with very large matrices where only a few coefficients are different from zero. In such cases, memory consumption can be reduced and performance increased by using a specialized representation storing only the nonzero coefficients. Such a matrix is called a sparse matrix.
Currently the elements of a given inner vector are guaranteed to be always sorted by increasing inner indices.
The \c "_" indicates available free space to quickly insert new elements.
Assuming no reallocation is needed, the insertion of a random element is therefore in O(nnz_j) where nnz_j is the number of nonzeros of the respective inner vector.
On the other hand, inserting elements with increasing inner indices in a given inner vector is much more efficient since this only requires to increase the respective \c InnerNNZs entry that is a O(1) operation.
The case where no empty space is available is a special case, and is refered as the \em compressed mode.
It corresponds to the widely used Compressed Column (or Row) Storage schemes (CCS or CRS).
Any SparseMatrix can be turned to this form by calling the SparseMatrix::makeCompressed() function.
In this case, one can remark that the \c InnerNNZs array is redundant with \c OuterIndexPtrs because we the equality: \c InnerNNZs[j] = \c OuterIndexPtrs[j+1]-\c OuterIndexPtrs[j].
Therefore, in practice a call to SparseMatrix::makeCompressed() frees this buffer.
It is worth noting that most of our wrappers to external libraries requires compressed matrices as inputs.
The results of Eigen's operations always produces \b compressed sparse matrices.
On the other hand, the insertion of a new element into a SparseMatrix converts this later to the \b uncompressed mode.
Because of the special storage scheme of a SparseMatrix, special care has to be taken when adding new nonzero entries.
For instance, the cost of inserting nnz non zeros in a a single purely random insertion into a SparseMatrix is O(nnz), where nnz is the current number of nonzero coefficients.
The simplest way to create a sparse matrix while guarantying good performance is to first build a list of so called \em triplets, and then convert it to a SparseMatrix.
The std::vector triplets might contain the elements in arbitrary order, and might even contain duplicated elements that will be summed up by setFromTriplets().
See the SparseMatrix::setFromTriplets() function and class Triplet for more details.
In some cases, however, slightly higher performance, and lower memory consumption can be reached by directly inserting the non zeros into the destination matrix.
A typical scenario of this approach is illustrated bellow:
- The key ingredient here is the line 2 where we reserve room for 6 non zeros per column. In many cases, the number of non zero per column or row can be easily known in advance. If it varies significantly for each inner vector, then it is possible to specify a reserve size for each inner vector by providing a vector object with an operator[](int j) returning the reserve size of the \c j-th inner vector (e.g., via a VectorXi or std::vector<int>). If only a rought estimate of the number of nonzeros per inner-vector can be obtained, it is highly recommended to overestimate it rather than the opposite. If this line is omitted, then the first insertion of a new element will reserve room for 2 elements per inner vector.
- The line 4 performs a sorted insertion. In this example, the ideal case is when the \c j-th column is not full and contains non-zeros whose inner-indices are smaller than \c i. In this case, this operation boils down to trivial O(1) operation.
- When calling insert(i,j) the element \c i \c ,j must not already exists, otherwise use the coeffRef(i,j) method that will allow to, e.g., accumulate values. This method first performs a binary search and finally calls insert(i,j) if the element does not already exist. It is more flexible than insert() but also more costly.
- The line 5 suppresses the remaining empty space and transforms the matrix into a compressed column storage.
In the following \em sm denotes a sparse matrix, \em sv a sparse vector, \em dm a dense matrix, and \em dv a dense vector.
In Eigen's sparse module we chose to expose only the subset of the dense matrix API which can be efficiently implemented. Moreover, not every combination is allowed; for instance, it is not possible to add two sparse matrices having two different storage orders. On the other hand, it is perfectly fine to evaluate a sparse matrix or expression to a matrix having a different storage order:
The product of a sparse \em symmetric matrix A with a dense matrix (or vector) d can be optimized by specifying the symmetry of A using selfadjointView:
In the above example, only the upper triangular part of the input matrix A is considered for solving. The opposite triangle might either be empty or contain arbitrary values.
In the case where multiple problems with the same sparcity pattern have to be solved, then the "compute" step can be decomposed as follow:
\code
SolverClassName<SparseMatrix<double> > solver;
solver.analyzePattern(A); // for this step the numerical values of A are not used
solver.factorize(A);
x1 = solver.solve(b1);
x2 = solver.solve(b2);
...
A = ...; // modify the values of the nonzeros of A, the nonzeros pattern must stay unchanged
