significantly extend the tutorial of sparse matrices

8 files changed, 303 insertions, 71 deletions

diff --git a/doc/C09_TutorialSparse.dox b/doc/C09_TutorialSparse.dox
index 86a6d1430..03ef2949e 100644
--- a/doc/C09_TutorialSparse.dox
+++ b/doc/C09_TutorialSparse.dox
@@ -8,9 +8,17 @@ namespace Eigen {
 
 \b Table \b of \b contents \n
   - \ref TutorialSparseIntro
+  - \ref TutorialSparseExample "Example"
+  - \ref TutorialSparseSparseMatrix
   - \ref TutorialSparseFilling
-  - \ref TutorialSparseFeatureSet
   - \ref TutorialSparseDirectSolvers
+  - \ref TutorialSparseFeatureSet
+    - \ref TutorialSparse_BasicOps
+    - \ref TutorialSparse_Products
+    - \ref TutorialSparse_TriangularSelfadjoint
+    - \ref TutorialSparse_Submat
+
+
 <hr>
 
 Manipulating and solving sparse problems involves various modules which are summarized below:
@@ -27,13 +35,14 @@ Manipulating and solving sparse problems involves various modules which are summ
 
 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
+\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 OuterStarts: stores for each colmun (resp. row) the index of the first non zero in the previous two arrays.
  - \c InnerNNZs: 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.
@@ -53,7 +62,7 @@ and one of its possible sparse, \b column \b major representation:
 <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>OuterStarts:</td><td>0</td><td>3</td><td>5</td><td>8</td><td>10</td><td>\em 12 </td></tr>
 <tr><td>InnerNNZs:</td>    <td>2</td><td>2</td><td>1</td><td>1</td><td> 2</td><td></td></tr>
 </table>
 
@@ -65,12 +74,12 @@ On the other hand, inserting elements with increasing inner indices in a given i
 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].
+In this case, one can remark that the \c InnerNNZs array is redundant with \c OuterStarts because we the equality: \c InnerNNZs[j] = \c OuterStarts[j+1]-\c OuterStarts[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.
+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.
 
 Here is the previous matrix represented in compressed mode:
@@ -79,25 +88,64 @@ Here is the previous matrix represented in compressed mode:
 <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>
+<tr><td>OuterStarts:</td><td>0</td><td>2</td><td>4</td><td>5</td><td>6</td><td>\em 8 </td></tr>
 </table>
 
 A SparseVector is a special case of a SparseMatrix where only the \c Values and \c InnerIndices arrays are stored.
 There is no notion of compressed/uncompressed mode for a SparseVector.
 
 
-\b Matrix \b and \b vector \b properties \n
+\section TutorialSparseExample First example
+
+Before describing each individual classes, lets start with the following typical example solving for the Lapace equation \f$ \nabla u = 0 \f$ onto a regular 2D grid using a finite difference scheme and Dirichlet boundary conditions.
+Such problem can be mathematically expressed as a linear problem of the form \f$ Ax=b \f$ where \f$ x \f$ is the vector of \c m unknows (in our case, the values of the pixels), \f$ b \f$ is the right hand side vector resuting from the boundary conditions, and \f$ A \f$ is a \f$ m \times m \f$ matrix containing only a few non-zeros elements resulting from the discretization of the Laplacian operator.
+
+<table class="manual">
+<tr><td>
+\include Tutorial_sparse_example.cpp
+</td>
+<td>
+\image html Tutorial_sparse_example.jpeg
+</td></tr></table>
+
+In this example, we start by defining a column-major sparse matrix type of double \c SparseMatrix<double>, and a triplet list of the same scalar type \c  Triplet<double>. A triplet is a simple object representing a non-zero entry as the triplet: \c row index, \c column index, \c value.
+
+In the main function, we declare a list \c coefficients of triplets (as a std vector) and the right hand side vector \f$ b \f$ which are filled by the \a buildProblem function.
+The raw and flat list of non-zero entries is then converted to a true SparseMatrix object \c A.
+Note that the elements of the list do not have to be sorted, and possible duplicate entries will be summed up.
+
+The last step consists in effectively solving the such assembled problem.
+Since the resulting matrix \c A is symmetric by construction, we can perform a direct Cholesky factorization via the SimplicialLDLT class which behaves like its LDLT counter part for dense objects.
+
+The resulting vector \c x contains the pixel values as a 1D array which is saved to a jpeg file shown on the right of the code above.
+
+Commenting the \a buildProblem and \a save functions is out of the scope of this tutorial. They are given \ref TutorialSparse_example_details "here" for the curious and reproducibility purpose.
+
+
+
 
-Here mat and vec represent any sparse-matrix and sparse-vector type, respectively.
+\section TutorialSparseSparseMatrix The SparseMatrix class
+
+\b %Matrix \b and \b vector \b properties \n
+
+The SparseMatrix and SparseVector classes take three template arguments:
+ * the scalar type (e.g., double)
+ * the storage order (ColMajor or RowMajor, the default is RowMajor)
+ * the inner index type (default is \c int).
+
+As for dense Matrix objects, constructors takes the size of the object.
+Here are some examples:
 
-Declarations:
 \code
-SparseMatrix<std::complex<float> > mat(1000,2000);         // declares a 1000x2000 col-major compressed sparse matrix of complex<float>
+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
 \endcode
 
+In the rest of the tutorial, \c mat and \c vec represent any sparse-matrix and sparse-vector objects, respectively.
+
+The dimensions of a matrix can be queried using the following functions:
 <table class="manual">
 <tr><td>Standard \n dimensions</td><td>\code
 mat.rows()
@@ -119,13 +167,16 @@ vec.nonZeros() \endcode</td></tr>
 
 \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:
+Random access to the elements of a sparse object can be done through the \c coeffRef(i,j) function.
+However, this function involves a quite expensive binary search.
+In most cases, one only wants to iterate over the non-zeros elements. This is achieved by a standard loop over the outer dimension, and then by iterating over the non-zeros of the current inner vector via an InnerIterator. Thus, the non-zero entries have to be visited in the same order than the storage order.
+Here is an example:
 <table class="manual">
 <tr><td>
 \code
-SparseMatrixType mat(rows,cols);
+SparseMatrix<double> mat(rows,cols);
 for (int k=0; k<mat.outerSize(); ++k)
-  for (SparseMatrixType::InnerIterator it(mat,k); it; ++it)
+  for (SparseMatrix<double>::InnerIterator it(mat,k); it; ++it)
   {
     it.value();
     it.row();   // row index
@@ -144,22 +195,21 @@ for (SparseVector<double>::InnerIterator it(vec); it; ++it)
 \endcode
 </td></tr>
 </table>
-
+For a writable expression, the referenced value can be modified using the valueRef() function.
 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.
+For instance, the cost of a single purely random insertion into a SparseMatrix is \c O(nnz), where \c nnz is the current number of non-zero 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 simplest way to create a sparse matrix while guarantying good performance is thus to first build a list of so called \em triplets, and then convert it to a SparseMatrix.
 
 Here is a typical usage example:
 \code
-typedef Triplet<double> T;
+typedef Eigen::Triplet<double> T;
 std::vector<T> tripletList;
 triplets.reserve(estimation_of_entries);
 for(...)
@@ -167,11 +217,11 @@ for(...)
   // ...
   tripletList.push_back(T(i,j,v_ij));
 }
-SparseMatrixType m(rows,cols);
-m.setFromTriplets(tripletList.begin(), tripletList.end());
-// m is ready to go!
+SparseMatrixType mat(rows,cols);
+mat.setFromTriplets(tripletList.begin(), tripletList.end());
+// mat is ready to go!
 \endcode
-The std::vector triplets might contain the elements in arbitrary order, and might even contain duplicated elements that will be summed up by setFromTriplets().
+The \c std::vector of 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.
 
 
@@ -185,61 +235,24 @@ A typical scenario of this approach is illustrated bellow:
 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 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 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 easily be 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.
 
 
-\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;
-sm1 = sm2 * 0.5;
-sm1 = sm2.transpose();
-sm1 = sm2.adjoint();
-sm4 = sm1 + sm2 + sm3;                          // only if sm1, sm2 and sm3 have the same storage order
-sm3 = sm1 * sm2;                                // conservative sparse * sparse product preserving numerical zeros
-sm3 = (sm1 * sm2).pruned();                     // sparse * sparse product that removes numerical zeros (triggers a different algorithm)
-sm3 = (sm1 * sm2).pruned(ref);                  // sparse * sparse product that removes elements much smaller than ref
-sm3 = (sm1 * sm2).pruned(ref,epsilon);          // sparse * sparse product that removes elements smaller than ref*epsilon
-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.
+%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 class="width20em"><p>Notes</p></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>
+<tr><th>Class</th><th>Module</th><th>Solver kind</th><th>Matrix kind</th><th>Features related to performance</th>
+    <th>Dependencies,License</th><th class="width20em"><p>Notes</p></th></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>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>
@@ -251,7 +264,7 @@ They are summarized in the following table:
 
 
 
-<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>
+<tr><td>CholmodSupernodalLLT</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>
@@ -318,6 +331,121 @@ The compute() method is equivalent to calling both analyzePattern() and factoriz
 Finally, each solver provides some specific features, such as determinant, access to the factors, controls of the iterations, and so on.
 More details are availble in the documentations of the respective classes.
 
+
+\section TutorialSparseFeatureSet Supported operators and functions
+
+Because of their special storage format, sparse matrices cannot offer the same level of flexbility than dense matrices.
+In Eigen's sparse module we chose to expose only the subset of the dense matrix API which can be efficiently implemented.
+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.
+
+\subsection TutorialSparse_BasicOps Basic operations
+
+%Sparse expressions support most of the unary and binary coefficient wise operations:
+\code
+sm1.real()   sm1.imag()   -sm1                    0.5*sm1
+sm1+sm2      sm1-sm2      sm1.cwiseProduct(sm2)
+\endcode
+However, a strong restriction is that the storage orders must match. For instance, in the following example:
+\code
+sm4 = sm1 + sm2 + sm3;
+\endcode
+sm1, sm2, and sm3 must all be row-major or all column major.
+On the other hand, there is no restriction on the target matrix sm4.
+For instance, this means that for computing \f$ A^T + A \f$, the matrix \f$ A^T \f$ must be evaluated into a temporary matrix of compatible storage order:
+\code
+SparseMatrix<double> A, B;
+B = SparseMatrix<double>(A.transpose()) + A;
+\endcode
+
+Binary coefficient wise operators can also mix sparse and dense expressions:
+\code
+sm2 = sm1.cwiseProduct(dm1);
+dm2 = sm1 + dm1;
+\endcode
+
+
+%Sparse expressions also support transposition:
+\code
+sm1 = sm2.transpose();
+sm1 = sm2.adjoint();
+\endcode
+However, there is no transposeInPlace() method.
+
+
+\subsection TutorialSparse_Products Matrix products
+
+%Eigen supports various kind of sparse matrix products which are summarize below:
+  - \b sparse-dense:
+    \code
+dv2 = sm1 * dv1;
+dm2 = dm1 * sm1.adjoint();
+dm2 = 2. * sm1 * dm1;
+    \endcode
+  - \b symmetric \b sparse-dense. The product of a sparse symmetric matrix with a dense matrix (or vector) can also be optimized by specifying the symmetry with selfadjointView():
+    \code
+dm2 = sm1.selfadjointView<>() * dm1;        // if all coefficients of A are stored
+dm2 = A.selfadjointView<Upper>() * dm1;     // if only the upper part of A is stored
+dm2 = A.selfadjointView<Lower>() * dm1;     // if only the lower part of A is stored
+    \endcode
+  - \b sparse-sparse. For sparse-sparse products, two different algorithms are available. The default one is conservative and preserve the explicit zeros that might appear:
+    \code
+sm3 = sm1 * sm2;
+sm3 = 4 * sm1.adjoint() * sm2;
+    \endcode
+    The second algorithm punes on the fly the explicit zeros, or the values smaller than a given threshold. It is enabled and controlled through the prune() functions:
+    \code
+sm3 = (sm1 * sm2).prune();                  // removes numerical zeros
+sm3 = (sm1 * sm2).prune(ref);               // removes elements much smaller than ref
+sm3 = (sm1 * sm2).prune(ref,epsilon);       // removes elements much smaller than ref*epsilon
+    \endcode
+
+  - \b permutations. Finally, permutations can be applied to sparse matrices too:
+    \code
+PermutationMatrix<Dynamic,Dynamic> P = ...;
+sm2 = P * sm1;
+sm2 = sm1 * P.inverse();
+sm2 = sm1.transpose() * P;
+    \endcode
+
+
+\subsection TutorialSparse_TriangularSelfadjoint Triangular and selfadjoint views
+
+Just as dense matrices, the triangularView() function can be used to address a triangular part of the matrix, and perform triangular solves with a dense right and side:
+\code
+dm2 = sm1.triangularView<Lower>(dm1);
+dv2 = sm1.transpose().triangularView<Upper>(dv1);
+\endcode
+
+The selfadjointView() function permits various operations:
+ - optimized sparse-dense matrix products:
+    \code
+dm2 = sm1.selfadjointView<>() * dm1;        // if all coefficients of A are stored
+dm2 = A.selfadjointView<Upper>() * dm1;     // if only the upper part of A is stored
+dm2 = A.selfadjointView<Lower>() * dm1;     // if only the lower part of A is stored
+    \endcode
+ - copy of triangular parts:
+    \code
+sm2 = sm1.selfadjointView<Upper>();                               // makes a full selfadjoint matrix from the upper triangular part
+sm2.selfadjointView<Lower>() = sm1.selfadjointView<Upper>();      // copie the upper triangular part to the lower triangular part
+    \endcode
+ - application of symmetric permutations:
+ \code
+PermutationMatrix<Dynamic,Dynamic> P = ...;
+sm2 = A.selfadjointView<Upper>().twistedBy(P);                                // compute P S P' from the upper triangular part of A, and make it a full matrix
+sm2.selfadjointView<Lower>() = A.selfadjointView<Lower>().twistedBy(P);       // compute P S P' from the lower triangular part of A, and then only compute the lower part
+ \endcode
+
+\subsection TutorialSparse_Submat Sub-matrices
+
+%Sparse matrices does not support yet the addressing of arbitrary sub matrices. Currently, one can only reference a set of contiguous \em inner vectors, i.e., a set of contiguous rows for a row-major matrix, or a set of contiguous columns for a column major matrix:
+\code
+  sm1.innerVector(j);       // returns an expression of the j-th column (resp. row) of the matrix if sm1 is col-major (resp. row-major)
+  sm1.innerVectors(j, nb);  // returns an expression of the nb columns (resp. row) starting from the j-th column (resp. row)
+                            // of the matrix if sm1 is col-major (resp. row-major)
+  sm1.middleRows(j, nb);    // for row major matrices only, get a range of nb rows
+  sm1.middleCols(j, nb);    // for column major matrices only, get a range of nb columns
+\endcode
+
 \li \b Next: \ref TutorialMapClass
 
 */
diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt
index 50ce7ee0c..96bff41bf 100644
--- a/doc/CMakeLists.txt
+++ b/doc/CMakeLists.txt
@@ -36,6 +36,7 @@ set(snippets_targets "")
 add_definitions("-DEIGEN_MAKING_DOCS")
 
 add_subdirectory(examples)
+add_subdirectory(special_examples)
 add_subdirectory(snippets)
 
 add_custom_target(
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index e87c7c2ba..e9e89d486 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -592,6 +592,7 @@ RECURSIVE              = YES
 EXCLUDE                = "${Eigen_SOURCE_DIR}/Eigen/Eigen2Support" \
                          "${Eigen_SOURCE_DIR}/Eigen/src/Eigen2Support" \
                          "${Eigen_SOURCE_DIR}/doc/examples" \
+                         "${Eigen_SOURCE_DIR}/doc/special_examples" \
                          "${Eigen_SOURCE_DIR}/doc/snippets" 
 
 # The EXCLUDE_SYMLINKS tag can be used select whether or not files or
@@ -638,7 +639,9 @@ EXCLUDE_SYMBOLS        = internal::* Flagged* *InnerIterator* DenseStorage<*
 EXAMPLE_PATH           = "${Eigen_SOURCE_DIR}/doc/snippets" \
                          "${Eigen_BINARY_DIR}/doc/snippets" \
                          "${Eigen_SOURCE_DIR}/doc/examples" \
-                         "${Eigen_BINARY_DIR}/doc/examples"
+                         "${Eigen_BINARY_DIR}/doc/examples" \
+                         "${Eigen_SOURCE_DIR}/doc/special_examples" \
+                         "${Eigen_BINARY_DIR}/doc/special_examples"
 
 # If the value of the EXAMPLE_PATH tag contains directories, you can use the
 # EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp
diff --git a/doc/SparseQuickReference.dox b/doc/SparseQuickReference.dox
index 7b86e50d3..7d6eb0fa9 100644
--- a/doc/SparseQuickReference.dox
+++ b/doc/SparseQuickReference.dox
@@ -18,7 +18,7 @@ In this page, we give a quick summary of the main operations available for spars
 i.e either row major or column major. The default is column major. Most arithmetic operations on sparse matrices will assert that they have the same storage order. Moreover, when interacting with external libraries that are not yet supported by Eigen, it is important to know how to send the required matrix pointers. 
 
 \section Constructors Constructors and assignments
-SparseMatrix is the core class to build and manipulate sparse matrices in Eigen. It takes as template parameters the Scalar type and the storage order, either RowMajor or ColumnMajor. The default is ColumnMajor. ??? It is possible to modify the default storage order at compile-time with the cmake variable \b EIGEN_DEFAULT_ROW_MAJOR ???
+SparseMatrix is the core class to build and manipulate sparse matrices in Eigen. It takes as template parameters the Scalar type and the storage order, either RowMajor or ColumnMajor. The default is ColumnMajor.
 
 \code
   SparseMatrix<double> sm1(1000,1000);              // 1000x1000 compressed sparse matrix of double. 
@@ -195,4 +195,4 @@ The following functions are useful to extract a block of rows (resp. columns) fr
 
 
 */
-}
\ No newline at end of file
+}
diff --git a/doc/TutorialSparse_example_details.dox b/doc/TutorialSparse_example_details.dox
new file mode 100644
index 000000000..0438da8bb
--- /dev/null
+++ b/doc/TutorialSparse_example_details.dox
@@ -0,0 +1,4 @@
+/**
+\page TutorialSparse_example_details
+\include Tutorial_sparse_example_details.cpp
+*/
diff --git a/doc/special_examples/CMakeLists.txt b/doc/special_examples/CMakeLists.txt
new file mode 100644
index 000000000..eeeae1d2a
--- /dev/null
+++ b/doc/special_examples/CMakeLists.txt
@@ -0,0 +1,20 @@
+
+if(NOT EIGEN_TEST_NOQT)
+  find_package(Qt4)
+  if(QT4_FOUND)
+    include(${QT_USE_FILE})
+  endif()
+endif(NOT EIGEN_TEST_NOQT)
+
+
+if(QT4_FOUND)
+  add_executable(Tutorial_sparse_example Tutorial_sparse_example.cpp Tutorial_sparse_example_details.cpp)
+  target_link_libraries(Tutorial_sparse_example ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO} ${QT_QTCORE_LIBRARY} ${QT_QTGUI_LIBRARY})
+
+  add_custom_command(
+      TARGET Tutorial_sparse_example
+      POST_BUILD
+      COMMAND Tutorial_sparse_example
+      ARGS ${CMAKE_CURRENT_BINARY_DIR}/../html/Tutorial_sparse_example.jpeg
+  )
+endif(QT4_FOUND)
diff --git a/doc/special_examples/Tutorial_sparse_example.cpp b/doc/special_examples/Tutorial_sparse_example.cpp
new file mode 100644
index 000000000..002f19f01
--- /dev/null
+++ b/doc/special_examples/Tutorial_sparse_example.cpp
@@ -0,0 +1,32 @@
+#include <Eigen/Sparse>
+#include <vector>
+
+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)
+{
+  int n = 300;  // size of the image
+  int m = n*n;  // number of unknows (=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;
+}
+
diff --git a/doc/special_examples/Tutorial_sparse_example_details.cpp b/doc/special_examples/Tutorial_sparse_example_details.cpp
new file mode 100644
index 000000000..8c3020b63
--- /dev/null
+++ b/doc/special_examples/Tutorial_sparse_example_details.cpp
@@ -0,0 +1,44 @@
+#include <Eigen/Sparse>
+#include <vector>
+#include <QImage>
+
+typedef Eigen::SparseMatrix<double> SpMat; // declares a column-major sparse matrix type of double
+typedef Eigen::Triplet<double> T;
+
+void insertCoefficient(int id, int i, int j, double w, std::vector<T>& coeffs,
+                       Eigen::VectorXd& b, const Eigen::VectorXd& boundary)
+{
+  int n = boundary.size();
+  int id1 = i+j*n;
+
+        if(i==-1 || i==n) b(id) -= w * boundary(j); // constrained coeffcieint
+  else  if(j==-1 || j==n) b(id) -= w * boundary(i); // constrained coeffcieint
+  else  coeffs.push_back(T(id,id1,w));              // unknown coefficient
+}
+
+void buildProblem(std::vector<T>& coefficients, Eigen::VectorXd& b, int n)
+{
+  b.setZero();
+  Eigen::ArrayXd boundary = Eigen::ArrayXd::LinSpaced(n, 0,M_PI).sin().pow(2);
+  for(int j=0; j<n; ++j)
+  {
+    for(int i=0; i<n; ++i)
+    {
+      int id = i+j*n;
+      insertCoefficient(id, i-1,j, -1, coefficients, b, boundary);
+      insertCoefficient(id, i+1,j, -1, coefficients, b, boundary);
+      insertCoefficient(id, i,j-1, -1, coefficients, b, boundary);
+      insertCoefficient(id, i,j+1, -1, coefficients, b, boundary);
+      insertCoefficient(id, i,j,    4, coefficients, b, boundary);
+    }
+  }
+}
+
+void saveAsBitmap(const Eigen::VectorXd& x, int n, const char* filename)
+{
+  Eigen::Array<unsigned char,Eigen::Dynamic,Eigen::Dynamic> bits = (x*255).cast<unsigned char>();
+  QImage img(bits.data(), n,n,QImage::Format_Indexed8);
+  img.setColorCount(256);
+  for(int i=0;i<256;i++) img.setColor(i,qRgb(i,i,i));
+  img.save(filename);
+}