doc/C09_TutorialSparse.dox

namespace Eigen {

/** \page TutorialSparse Tutorial page 9 - Sparse Matrix
    \ingroup Tutorial

\li \b Previous: \ref TutorialGeometry
\li \b Next: TODO

\b Table \b of \b contents \n
  - \ref TutorialSparseIntro
  - \ref TutorialSparseFilling
  - \ref TutorialSparseFeatureSet
  - \ref TutorialSparseDirectSolvers
<hr>

Manipulating and solving sparse problems involves various modules which are summarized below:

<table class="manual">
<tr><th>Module</th><th>Header file</th><th>Contents</th></tr>
<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>
</table>

\section TutorialSparseIntro Representing sparse matrices

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.

\b The \b SparseMatrix \b class
The class SparseMatrix is the main sparse matrix representation of Eigen's sparse module; it offers high performance and low memory usage.
It implements a more versatile variant of the widely-used Compressed Column (or Row) Storage scheme.
It consists of four compact arrays:
 - \c Values: stores the coefficient values of the non-zeros.
 - \c InnerIndices: stores the row (resp. column) indices of the non-zeros.
 - \c OuterIndexPtrs: stores for each colmun (resp. row) the index of the first non zero in the previous arrays.
 - \c InnerSizes: stores the number of non-zeros of each column (resp. row).
The word \c inner refers to an \em inner \em vector that is a column for a column-major matrix, or a row for a row-major matrix.
The word \c outer refers to the other direction.

This storage scheme is better explained on an example. The following matrix
<table class="manual">
<tr><td> 0</td><td>3</td><td> 0</td><td>0</td><td> 0</td></tr>
<tr><td>22</td><td>0</td><td> 0</td><td>0</td><td>17</td></tr>
<tr><td> 7</td><td>5</td><td> 0</td><td>1</td><td> 0</td></tr>
<tr><td> 0</td><td>0</td><td> 0</td><td>0</td><td> 0</td></tr>
<tr><td> 0</td><td>0</td><td>14</td><td>0</td><td> 8</td></tr>
</table>

and its sparse, \b column \b major representation:
<table class="manual">
<tr><td>Values:</td>        <td>22</td><td>7</td><td>_</td><td>3</td><td>5</td><td>14</td><td>_</td><td>_</td><td>1</td><td>_</td><td>17</td><td>8</td></tr>
<tr><td>InnerIndices:</td>  <td> 1</td><td>2</td><td>_</td><td>0</td><td>2</td><td> 4</td><td>_</td><td>_</td><td>2</td><td>_</td><td> 1</td><td>4</td></tr>
</table>
<table class="manual">
<tr><td>OuterIndexPtrs:</td><td>0</td><td>3</td><td>5</td><td>8</td><td>10</td><td>\em 12 </td></tr>
<tr><td>InnerSizes:</td>    <td>2</td><td>2</td><td>1</td><td>1</td><td> 2</td><td></td></tr>
</table>

By default the elements of a given inner vector are 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 of coordinate is therefore in O(nnz) where nnz 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 incresae the respective \c InnerSizes 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 and it corresponds to the widely used Compressed Column (or Row) Storage scheme.
Any SparseMatrix can be turned in this form by calling the SparseMatrix::makeCompressed() function.
In this case one can remark that the \c InnerSizes array is superfluous because \c InnerSizes[j] is simply equals to \c OuterIndexPtrs[j+1]-\c OuterIndexPtrs[j].
Therefore, a call to SparseMatrix::makeCompressed() frees this buffer.

Here is the previous matrix represented in compressed mode:
<table class="manual">
<tr><td>Values:</td>        <td>22</td><td>7</td><td>3</td><td>5</td><td>14</td><td>1</td><td>17</td><td>8</td></tr>
<tr><td>InnerIndices:</td>  <td> 1</td><td>2</td><td>0</td><td>2</td><td> 4</td><td>2</td><td> 1</td><td>4</td></tr>
</table>
<table class="manual">
<tr><td>OuterIndexPtrs:</td><td>0</td><td>2</td><td>4</td><td>5</td><td>6</td><td>\em 8 </td></tr>
</table>
In this mode any insertion of new non zero elements would be extremely costly.



\b Matrix \b and \b vector \b properties \n

Here mat and vec represent any sparse-matrix and sparse-vector type, respectively.

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

Although a sparse matrix could also be used to represent a sparse vector, for that purpose it is better to use the specialized SparseVector class:
\code

\endcode

<table class="manual">
<tr><td>Standard \n dimensions</td><td>\code
mat.rows()
mat.cols()\endcode</td>
<td>\code
vec.size() \endcode</td>
</tr>
<tr><td>Sizes along the \n inner/outer dimensions</td><td>\code
mat.innerSize()
mat.outerSize()\endcode</td>
<td></td>
</tr>
<tr><td>Number of non \n zero coefficients</td><td>\code
mat.nonZeros() \endcode</td>
<td>\code
vec.nonZeros() \endcode</td></tr>
</table>


\b Iterating \b over \b the \b nonzero \b coefficients \n

Iterating over the coefficients of a sparse matrix can be done only in the same order as the storage order. Here is an example:
<table class="manual">
<tr><td>
\code
SparseMatrixType mat(rows,cols);
for (int k=0; k<mat.outerSize(); ++k)
  for (SparseMatrixType::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()
  }
\endcode
</td><td>
\code
SparseVector<double> vec(size);
for (SparseVector<double>::InnerIterator it(vec); it; ++it)
{
  it.value(); // == vec[ it.index() ]
  it.index();
}
\endcode
</td></tr>
</table>

If the type of the sparse matrix or vector depends on a template parameter, then the \c typename keyword is
required to indicate that \c InnerIterator denotes a type; see \ref TopicTemplateKeyword for details.

\section TutorialSparseFilling Filling a sparse matrix

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.

A typical scenario to insert nonzeros is illustrated bellow:
\code
1: SparseMatrix<double> m(rows,cols);        // default is column major
2: aux.reserve(VectorXi::Constant(rows,6));
3: for each i,j such that v_ij != 0
4:   mat.insert(i,j) = v_ij;
5: mat.makeCompressed();                     // optional
\endcode

- 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 provifing a VectorXi or std::vector compatible object. 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.
- The line 4 performs a sorted insertion. In this example, the ideal case is when the j-th column is not full and contains non-zeros whose inner-indices are smaller than i. In this case, this operation boils down to trivial O(1) operation.
- The line 5 suppresses the left empty room and transform the matrix to a compressed column storage.


\section TutorialSparseFeatureSet Supported operators and functions

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:
\code
SparseMatrixType sm1, sm2, sm3;
sm3 = sm1.transpose() + sm2;                    // invalid, because transpose() changes the storage order
sm3 = SparseMatrixType(sm1.transpose()) + sm2;  // correct, because evaluation reformats as column-major
\endcode

Here are some examples of supported operations:
\code
sm1 *= 0.5;
sm4 = sm1 + sm2 + sm3;          // only if sm1, sm2 and sm3 have the same storage order
sm3 = sm1 * sm2;
dv3 = sm1 * dv2;
dm3 = sm1 * dm2;
dm3 = dm2 * sm1;
sm3 = sm1.cwiseProduct(sm2);    // only if sm1 and sm2 have the same storage order
dv2 = sm1.triangularView<Upper>().solve(dv2);
\endcode

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:
\code
res = A.selfadjointView<>() * d;        // if all coefficients of A are stored
res = A.selfadjointView<Upper>() * d;   // if only the upper part of A is stored
res = A.selfadjointView<Lower>() * d;   // if only the lower part of A is stored
\endcode


\section TutorialSparseDirectSolvers Solving linear problems

Eigen currently provides a limited set of built-in solvers as well as wrappers to external solver libraries.
They are summarized in the following table:

<table class="manual">
<tr><td>Class    </td><td>Module</td><td>Solver kind</td><td>Matrix kind</td><td>Features related to performance</td>
    <td>Dependencies,License</td>
    <td>Notes</td></tr>

<tr><td>SimplicialLLt    </td><td>\link SparseCholesky_Module SparseCholesky \endlink</td><td>Direct LLt factorization</td><td>SPD</td><td>Fill-in reducing</td>
    <td>built-in, LGPL</td>
    <td>SimplicialLDLt is often preferable</td></tr>
<tr><td>SimplicialLDLt   </td><td>\link SparseCholesky_Module SparseCholesky \endlink</td><td>Direct LDLt factorization</td><td>SPD</td><td>Fill-in reducing</td>
    <td>built-in, LGPL</td>
    <td>Recommended for very sparse and not too large problems (e.g., 2D Poisson eq.)</td></tr>
<tr><td>ConjugateGradient</td><td>\link IterativeLinearSolvers_Module IterativeLinearSolvers \endlink</td><td>Classic iterative CG</td><td>SPD</td><td>Preconditionning</td>
    <td>built-in, LGPL</td>
    <td>Recommended for large symmetric problems (e.g., 3D Poisson eq.)</td></tr>
<tr><td>BiCGSTAB</td><td>\link IterativeLinearSolvers_Module IterativeLinearSolvers \endlink</td><td>Iterative stabilized bi-conjugate gradient</td><td>Square</td><td>Preconditionning</td>
    <td>built-in, LGPL</td>
    <td>Might not always converge</td></tr>



<tr><td>CholmodDecomposition</td><td>\link CholmodSupport_Module CholmodSupport \endlink</td><td>Direct LLT factorization</td><td>SPD</td><td>Fill-in reducing, Leverage fast dense algebra</td>
    <td>Requires the <a href="http://www.cise.ufl.edu/research/sparse/SuiteSparse/">SuiteSparse</a> package, \b GPL </td>
    <td></td></tr>
<tr><td>UmfPackLU</td><td>\link UmfPackSupport_Module UmfPackSupport \endlink</td><td>Direct LU factorization</td><td>Square</td><td>Fill-in reducing, Leverage fast dense algebra</td>
    <td>Requires the <a href="http://www.cise.ufl.edu/research/sparse/SuiteSparse/">SuiteSparse</a> package, \b GPL </td>
    <td></td></tr>
<tr><td>SuperLU</td><td>\link SuperLUSupport_Module SuperLUSupport \endlink</td><td>Direct LU factorization</td><td>Square</td><td>Fill-in reducing, Leverage fast dense algebra</td>
    <td>Requires the <a href="http://crd-legacy.lbl.gov/~xiaoye/SuperLU/">SuperLU</a> library, custom (BSD-like)</td>
    <td></td></tr>
</table>

Here \c SPD means symmetric positive definite.

All these solvers follow the same general concept.
Here is a typical and general example:
\code
#include <Eigen/RequiredModuleName>
// ...
SparseMatrix<double> A;
// fill A
VectorXd b, x;
// fill b
// solve Ax = b
SolverClassName<SparseMatrix<double> > solver;
solver.compute(A);
if(solver.info()!=Succeeded) {
  // decomposition failed
  return;
}
x = solver.solve(b);
if(solver.info()!=Succeeded) {
  // solving failed
  return;
}
// solve for another right hand side:
x1 = solver.solve(b1);
\endcode

For \c SPD solvers, a second optional template argument allows to specify which triangular part have to be used, e.g.:

\code
#include <Eigen/IterativeLinearSolvers>

ConjugateGradient<SparseMatrix<double>, Eigen::Upper> solver;
x = solver.compute(A).solve(b);
\endcode
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
solver.factorize(A);
x1 = solver.solve(b1);
x2 = solver.solve(b2);
...
\endcode
The compute() methode is equivalent to calling both analyzePattern() and factorize().

Finally, each solver provides some specific features, such as determinant, access to the factors, controls of the iterations, etc.
More details are availble on the documentation of the respective classes.

\li \b Next: TODO

*/

}