<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<meta name="generator" content="Doxygen 1.8.5"/>
<title>Eigen: Sparse matrix manipulations</title>
<link href="tabs.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="dynsections.js"></script>
<link href="navtree.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="resize.js"></script>
<script type="text/javascript" src="navtree.js"></script>
<script type="text/javascript">
$(document).ready(initResizable);
$(window).load(resizeHeight);
</script>
<link href="search/search.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="search/search.js"></script>
<script type="text/javascript">
$(document).ready(function() { searchBox.OnSelectItem(0); });
</script>
<link href="doxygen.css" rel="stylesheet" type="text/css" />
<link href="eigendoxy.css" rel="stylesheet" type="text/css">
<!-- -->
<script type="text/javascript" src="eigen_navtree_hacks.js"></script>
<!-- <script type="text/javascript"> -->
<!-- </script> -->
</head>
<body>
<div id="top"><!-- do not remove this div, it is closed by doxygen! -->
<!-- <a name="top"></a> -->
<div id="titlearea">
<table cellspacing="0" cellpadding="0">
<tbody>
<tr style="height: 56px;">
<td id="projectlogo"><img alt="Logo" src="Eigen_Silly_Professor_64x64.png"/></td>
<td style="padding-left: 0.5em;">
<div id="projectname"><a href="http://eigen.tuxfamily.org">Eigen</a>
 <span id="projectnumber">3.2.0</span>
</div>
</td>
<td> <div id="MSearchBox" class="MSearchBoxInactive">
<span class="left">
<img id="MSearchSelect" src="search/mag_sel.png"
onmouseover="return searchBox.OnSearchSelectShow()"
onmouseout="return searchBox.OnSearchSelectHide()"
alt=""/>
<input type="text" id="MSearchField" value="Search" accesskey="S"
onfocus="searchBox.OnSearchFieldFocus(true)"
onblur="searchBox.OnSearchFieldFocus(false)"
onkeyup="searchBox.OnSearchFieldChange(event)"/>
</span><span class="right">
<a id="MSearchClose" href="javascript:searchBox.CloseResultsWindow()"><img id="MSearchCloseImg" border="0" src="search/close.png" alt=""/></a>
</span>
</div>
</td>
</tr>
</tbody>
</table>
</div>
<!-- end header part -->
<!-- Generated by Doxygen 1.8.5 -->
<script type="text/javascript">
var searchBox = new SearchBox("searchBox", "search",false,'Search');
</script>
</div><!-- top -->
<div id="side-nav" class="ui-resizable side-nav-resizable">
<div id="nav-tree">
<div id="nav-tree-contents">
<div id="nav-sync" class="sync"></div>
</div>
</div>
<div id="splitbar" style="-moz-user-select:none;"
class="ui-resizable-handle">
</div>
</div>
<script type="text/javascript">
$(document).ready(function(){initNavTree('group__TutorialSparse.html','');});
</script>
<div id="doc-content">
<!-- window showing the filter options -->
<div id="MSearchSelectWindow"
onmouseover="return searchBox.OnSearchSelectShow()"
onmouseout="return searchBox.OnSearchSelectHide()"
onkeydown="return searchBox.OnSearchSelectKey(event)">
<a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(0)"><span class="SelectionMark"> </span>All</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(1)"><span class="SelectionMark"> </span>Classes</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(2)"><span class="SelectionMark"> </span>Namespaces</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(3)"><span class="SelectionMark"> </span>Functions</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(4)"><span class="SelectionMark"> </span>Variables</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(5)"><span class="SelectionMark"> </span>Typedefs</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(6)"><span class="SelectionMark"> </span>Enumerations</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(7)"><span class="SelectionMark"> </span>Enumerator</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(8)"><span class="SelectionMark"> </span>Friends</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(9)"><span class="SelectionMark"> </span>Groups</a><a class="SelectItem" href="javascript:void(0)" onclick="searchBox.OnSelectItem(10)"><span class="SelectionMark"> </span>Pages</a></div>
<!-- iframe showing the search results (closed by default) -->
<div id="MSearchResultsWindow">
<iframe src="javascript:void(0)" frameborder="0"
name="MSearchResults" id="MSearchResults">
</iframe>
</div>
<div class="header">
<div class="headertitle">
<div class="title">Sparse matrix manipulations<div class="ingroups"><a class="el" href="group__Sparse__chapter.html">Sparse linear algebra</a></div></div> </div>
</div><!--header-->
<div class="contents">
<p>Manipulating and solving sparse problems involves various modules which are summarized below:</p>
<table class="manual">
<tr>
<th>Module</th><th>Header file</th><th>Contents </th></tr>
<tr>
<td><a class="el" href="group__SparseCore__Module.html">SparseCore </a></td><td><div class="fragment"><div class="line"><span class="preprocessor">#include <Eigen/SparseCore></span></div>
</div><!-- fragment --></td><td><a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation. ">SparseMatrix</a> and <a class="el" href="classEigen_1_1SparseVector.html" title="a sparse vector class ">SparseVector</a> classes, matrix assembly, basic sparse linear algebra (including sparse triangular solvers) </td></tr>
<tr>
<td><a class="el" href="group__SparseCholesky__Module.html">SparseCholesky </a></td><td><div class="fragment"><div class="line"><span class="preprocessor">#include <Eigen/SparseCholesky></span></div>
</div><!-- fragment --></td><td>Direct sparse <a class="el" href="classEigen_1_1LLT.html" title="Standard Cholesky decomposition (LL^T) of a matrix and associated features. ">LLT</a> and <a class="el" href="classEigen_1_1LDLT.html" title="Robust Cholesky decomposition of a matrix with pivoting. ">LDLT</a> Cholesky factorization to solve sparse self-adjoint positive definite problems </td></tr>
<tr>
<td><a class="el" href="group__SparseLU__Module.html">SparseLU </a></td><td><div class="fragment"><div class="line"><span class="preprocessor">#include<Eigen/SparseLU></span> </div>
</div><!-- fragment --> </td><td>Sparse LU factorization to solve general square sparse systems </td></tr>
<tr>
<td><a class="el" href="group__SparseQR__Module.html">SparseQR </a></td><td><div class="fragment"><div class="line"><span class="preprocessor">#include<Eigen/SparseQR></span></div>
</div><!-- fragment --> </td><td>Sparse QR factorization for solving sparse linear least-squares problems </td></tr>
<tr>
<td><a class="el" href="group__IterativeLinearSolvers__Module.html">IterativeLinearSolvers </a></td><td><div class="fragment"><div class="line"><span class="preprocessor">#include <Eigen/IterativeLinearSolvers></span></div>
</div><!-- fragment --></td><td>Iterative solvers to solve large general linear square problems (including self-adjoint positive definite problems) </td></tr>
<tr>
<td><a class="el" href="group__Sparse__Module.html">Sparse </a></td><td><div class="fragment"><div class="line"><span class="preprocessor">#include <Eigen/Sparse></span></div>
</div><!-- fragment --></td><td>Includes all the above modules </td></tr>
</table>
<h1><a class="anchor" id="TutorialSparseIntro"></a>
Sparse matrix format</h1>
<p>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.</p>
<p><b>The</b> <b>SparseMatrix</b> <b>class</b> </p>
<p>The class <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation. ">SparseMatrix</a> is the main sparse matrix representation of <a class="el" href="namespaceEigen.html" title="Namespace containing all symbols from the Eigen library. ">Eigen</a>'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:</p>
<ul>
<li><code>Values:</code> stores the coefficient values of the non-zeros.</li>
<li><code>InnerIndices:</code> stores the row (resp. column) indices of the non-zeros.</li>
<li><code>OuterStarts:</code> stores for each column (resp. row) the index of the first non-zero in the previous two arrays.</li>
<li><code>InnerNNZs:</code> stores the number of non-zeros of each column (resp. row). The word <code>inner</code> refers to an <em>inner</em> <em>vector</em> that is a column for a column-major matrix, or a row for a row-major matrix. The word <code>outer</code> refers to the other direction.</li>
</ul>
<p>This storage scheme is better explained on an example. The following matrix </p>
<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>
<p>and one of its possible sparse, <b>column</b> <b>major</b> representation: </p>
<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>OuterStarts:</td><td>0</td><td>3</td><td>5</td><td>8</td><td>10</td><td><em>12</em> </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>
<p>Currently the elements of a given inner vector are guaranteed to be always sorted by increasing inner indices. The <code>"_"</code> 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 <code>InnerNNZs</code> entry that is a O(1) operation.</p>
<p>The case where no empty space is available is a special case, and is refered as the <em>compressed</em> mode. It corresponds to the widely used Compressed Column (or Row) Storage schemes (CCS or CRS). Any <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation. ">SparseMatrix</a> can be turned to this form by calling the <a class="el" href="classEigen_1_1SparseMatrix.html#a5daf3de834937aeaa7ebc7f18cea746b">SparseMatrix::makeCompressed()</a> function. In this case, one can remark that the <code>InnerNNZs</code> array is redundant with <code>OuterStarts</code> because we the equality: <code>InnerNNZs</code>[j] = <code>OuterStarts</code>[j+1]-<code>OuterStarts</code>[j]. Therefore, in practice a call to <a class="el" href="classEigen_1_1SparseMatrix.html#a5daf3de834937aeaa7ebc7f18cea746b">SparseMatrix::makeCompressed()</a> frees this buffer.</p>
<p>It is worth noting that most of our wrappers to external libraries requires compressed matrices as inputs.</p>
<p>The results of Eigen's operations always produces <b>compressed</b> sparse matrices. On the other hand, the insertion of a new element into a <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation. ">SparseMatrix</a> converts this later to the <b>uncompressed</b> mode.</p>
<p>Here is the previous matrix represented in compressed mode: </p>
<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>OuterStarts:</td><td>0</td><td>2</td><td>4</td><td>5</td><td>6</td><td><em>8</em> </td></tr>
</table>
<p>A <a class="el" href="classEigen_1_1SparseVector.html" title="a sparse vector class ">SparseVector</a> is a special case of a <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation. ">SparseMatrix</a> where only the <code>Values</code> and <code>InnerIndices</code> arrays are stored. There is no notion of compressed/uncompressed mode for a <a class="el" href="classEigen_1_1SparseVector.html" title="a sparse vector class ">SparseVector</a>.</p>
<h1><a class="anchor" id="TutorialSparseExample"></a>
First example</h1>
<p>Before describing each individual class, let's start with the following typical example: solving the Lapace equation <img class="formulaInl" alt="$ \nabla u = 0 $" src="form_198.png"/> on 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 <img class="formulaInl" alt="$ Ax=b $" src="form_199.png"/> where <img class="formulaInl" alt="$ x $" src="form_17.png"/> is the vector of <code>m</code> unknowns (in our case, the values of the pixels), <img class="formulaInl" alt="$ b $" src="form_200.png"/> is the right hand side vector resulting from the boundary conditions, and <img class="formulaInl" alt="$ A $" src="form_1.png"/> is an <img class="formulaInl" alt="$ m \times m $" src="form_201.png"/> matrix containing only a few non-zero elements resulting from the discretization of the Laplacian operator.</p>
<table class="manual">
<tr>
<td><div class="fragment"><div class="line"><span class="preprocessor">#include <Eigen/Sparse></span></div>
<div class="line"><span class="preprocessor">#include <vector></span></div>
<div class="line"></div>
<div class="line"><span class="keyword">typedef</span> <a class="code" href="classEigen_1_1SparseMatrix.html">Eigen::SparseMatrix<double></a> SpMat; <span class="comment">// declares a column-major sparse matrix type of double</span></div>
<div class="line"><span class="keyword">typedef</span> <a class="code" href="classEigen_1_1Triplet.html">Eigen::Triplet<double></a> T;</div>
<div class="line"></div>
<div class="line"><span class="keywordtype">void</span> buildProblem(std::vector<T>& coefficients, <a class="code" href="classEigen_1_1Matrix.html">Eigen::VectorXd</a>& b, <span class="keywordtype">int</span> n);</div>
<div class="line"><span class="keywordtype">void</span> saveAsBitmap(<span class="keyword">const</span> <a class="code" href="classEigen_1_1Matrix.html">Eigen::VectorXd</a>& x, <span class="keywordtype">int</span> n, <span class="keyword">const</span> <span class="keywordtype">char</span>* filename);</div>
<div class="line"></div>
<div class="line"><span class="keywordtype">int</span> main(<span class="keywordtype">int</span> argc, <span class="keywordtype">char</span>** argv)</div>
<div class="line">{</div>
<div class="line"> <span class="keywordtype">int</span> n = 300; <span class="comment">// size of the image</span></div>
<div class="line"> <span class="keywordtype">int</span> m = n*n; <span class="comment">// number of unknows (=number of pixels)</span></div>
<div class="line"></div>
<div class="line"> <span class="comment">// Assembly:</span></div>
<div class="line"> std::vector<T> coefficients; <span class="comment">// list of non-zeros coefficients</span></div>
<div class="line"> <a class="code" href="classEigen_1_1Matrix.html">Eigen::VectorXd</a> b(m); <span class="comment">// the right hand side-vector resulting from the constraints</span></div>
<div class="line"> buildProblem(coefficients, b, n);</div>
<div class="line"></div>
<div class="line"> SpMat A(m,m);</div>
<div class="line"> A.setFromTriplets(coefficients.begin(), coefficients.end());</div>
<div class="line"></div>
<div class="line"> <span class="comment">// Solving:</span></div>
<div class="line"> <a class="code" href="classEigen_1_1SimplicialCholesky.html">Eigen::SimplicialCholesky<SpMat></a> chol(A); <span class="comment">// performs a Cholesky factorization of A</span></div>
<div class="line"> <a class="code" href="classEigen_1_1Matrix.html">Eigen::VectorXd</a> x = chol.solve(b); <span class="comment">// use the factorization to solve for the given right hand side</span></div>
<div class="line"></div>
<div class="line"> <span class="comment">// Export the result to a file:</span></div>
<div class="line"> saveAsBitmap(x, n, argv[1]);</div>
<div class="line"></div>
<div class="line"> <span class="keywordflow">return</span> 0;</div>
<div class="line">}</div>
<div class="line"></div>
</div><!-- fragment --> </td><td><div class="image">
<img src="Tutorial_sparse_example.jpeg" alt="Tutorial_sparse_example.jpeg"/>
</div>
</td></tr>
</table>
<p>In this example, we start by defining a column-major sparse matrix type of double <code>SparseMatrix<double></code>, and a triplet list of the same scalar type <code>Triplet<double></code>. A triplet is a simple object representing a non-zero entry as the triplet: <code>row</code> index, <code>column</code> index, <code>value</code>.</p>
<p>In the main function, we declare a list <code>coefficients</code> of triplets (as a std vector) and the right hand side vector <img class="formulaInl" alt="$ b $" src="form_200.png"/> which are filled by the <em>buildProblem</em> function. The raw and flat list of non-zero entries is then converted to a true <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation. ">SparseMatrix</a> object <code>A</code>. Note that the elements of the list do not have to be sorted, and possible duplicate entries will be summed up.</p>
<p>The last step consists of effectively solving the assembled problem. Since the resulting matrix <code>A</code> is symmetric by construction, we can perform a direct Cholesky factorization via the <a class="el" href="classEigen_1_1SimplicialLDLT.html" title="A direct sparse LDLT Cholesky factorizations without square root. ">SimplicialLDLT</a> class which behaves like its <a class="el" href="classEigen_1_1LDLT.html" title="Robust Cholesky decomposition of a matrix with pivoting. ">LDLT</a> counterpart for dense objects.</p>
<p>The resulting vector <code>x</code> contains the pixel values as a 1D array which is saved to a jpeg file shown on the right of the code above.</p>
<p>Describing the <em>buildProblem</em> and <em>save</em> functions is out of the scope of this tutorial. They are given <a class="el" href="TutorialSparse_example_details.html">here</a> for the curious and reproducibility purpose.</p>
<h1><a class="anchor" id="TutorialSparseSparseMatrix"></a>
The SparseMatrix class</h1>
<p><b>Matrix</b> <b>and</b> <b>vector</b> <b>properties</b> <br/>
The <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation. ">SparseMatrix</a> and <a class="el" href="classEigen_1_1SparseVector.html" title="a sparse vector class ">SparseVector</a> classes take three template arguments: the scalar type (e.g., double) the storage order (ColMajor or RowMajor, the default is ColMajor) the inner index type (default is <code>int</code>).</p>
<p>As for dense <a class="el" href="classEigen_1_1Matrix.html" title="The matrix class, also used for vectors and row-vectors. ">Matrix</a> objects, constructors takes the size of the object. Here are some examples:</p>
<div class="fragment"><div class="line">SparseMatrix<std::complex<float> > mat(1000,2000); <span class="comment">// declares a 1000x2000 column-major compressed sparse matrix of complex<float></span></div>
<div class="line">SparseMatrix<double,RowMajor> mat(1000,2000); <span class="comment">// declares a 1000x2000 row-major compressed sparse matrix of double</span></div>
<div class="line">SparseVector<std::complex<float> > vec(1000); <span class="comment">// declares a column sparse vector of complex<float> of size 1000</span></div>
<div class="line">SparseVector<double,RowMajor> vec(1000); <span class="comment">// declares a row sparse vector of double of size 1000</span></div>
</div><!-- fragment --><p>In the rest of the tutorial, <code>mat</code> and <code>vec</code> represent any sparse-matrix and sparse-vector objects, respectively.</p>
<p>The dimensions of a matrix can be queried using the following functions: </p>
<table class="manual">
<tr>
<td>Standard <br/>
dimensions</td><td><div class="fragment"><div class="line">mat.rows()</div>
<div class="line">mat.cols()</div>
</div><!-- fragment --> </td><td><div class="fragment"><div class="line">vec.size() </div>
</div><!-- fragment --> </td></tr>
<tr>
<td>Sizes along the <br/>
inner/outer dimensions</td><td><div class="fragment"><div class="line">mat.<a class="code" href="classEigen_1_1DenseBase.html#acd9791e1914ad8761992e16032d49c54">innerSize</a>()</div>
<div class="line">mat.<a class="code" href="classEigen_1_1DenseBase.html#af080d3d3f82e0d4391f19af22a5eedb8">outerSize</a>()</div>
</div><!-- fragment --> </td><td></td></tr>
<tr>
<td>Number of non <br/>
zero coefficients</td><td><div class="fragment"><div class="line">mat.<a class="code" href="classEigen_1_1DenseBase.html#abe29bfb5f2d88cf2a50ffb577ec6f0a4">nonZeros</a>() </div>
</div><!-- fragment --> </td><td><div class="fragment"><div class="line">vec.nonZeros() </div>
</div><!-- fragment --> </td></tr>
</table>
<p><b>Iterating</b> <b>over</b> <b>the</b> <b>nonzero</b> <b>coefficients</b> <br/>
Random access to the elements of a sparse object can be done through the <code>coeffRef(i,j)</code> 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: </p>
<table class="manual">
<tr>
<td><div class="fragment"><div class="line">SparseMatrix<double> mat(rows,cols);</div>
<div class="line"><span class="keywordflow">for</span> (<span class="keywordtype">int</span> k=0; k<mat.<a class="code" href="classEigen_1_1DenseBase.html#af080d3d3f82e0d4391f19af22a5eedb8">outerSize</a>(); ++k)</div>
<div class="line"> <span class="keywordflow">for</span> (SparseMatrix<double>::InnerIterator it(mat,k); it; ++it)</div>
<div class="line"> {</div>
<div class="line"> it.value();</div>
<div class="line"> it.row(); <span class="comment">// row index</span></div>
<div class="line"> it.col(); <span class="comment">// col index (here it is equal to k)</span></div>
<div class="line"> it.index(); <span class="comment">// inner index, here it is equal to it.row()</span></div>
<div class="line"> }</div>
</div><!-- fragment --> </td><td><div class="fragment"><div class="line">SparseVector<double> vec(size);</div>
<div class="line"><span class="keywordflow">for</span> (SparseVector<double>::InnerIterator it(vec); it; ++it)</div>
<div class="line">{</div>
<div class="line"> it.value(); <span class="comment">// == vec[ it.index() ]</span></div>
<div class="line"> it.index();</div>
<div class="line">}</div>
</div><!-- fragment --> </td></tr>
</table>
<p>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 <code>typename</code> keyword is required to indicate that <code>InnerIterator</code> denotes a type; see <a class="el" href="TopicTemplateKeyword.html">The template and typename keywords in C++</a> for details.</p>
<h1><a class="anchor" id="TutorialSparseFilling"></a>
Filling a sparse matrix</h1>
<p>Because of the special storage scheme of a <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation. ">SparseMatrix</a>, special care has to be taken when adding new nonzero entries. For instance, the cost of a single purely random insertion into a <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation. ">SparseMatrix</a> is <code>O(nnz)</code>, where <code>nnz</code> is the current number of non-zero coefficients.</p>
<p>The simplest way to create a sparse matrix while guaranteeing good performance is thus to first build a list of so-called <em>triplets</em>, and then convert it to a <a class="el" href="classEigen_1_1SparseMatrix.html" title="A versatible sparse matrix representation. ">SparseMatrix</a>.</p>
<p>Here is a typical usage example: </p>
<div class="fragment"><div class="line"><span class="keyword">typedef</span> <a class="code" href="classEigen_1_1Triplet.html">Eigen::Triplet<double></a> T;</div>
<div class="line">std::vector<T> tripletList;</div>
<div class="line">tripletList.reserve(estimation_of_entries);</div>
<div class="line"><span class="keywordflow">for</span>(...)</div>
<div class="line">{</div>
<div class="line"> <span class="comment">// ...</span></div>
<div class="line"> tripletList.push_back(T(i,j,v_ij));</div>
<div class="line">}</div>
<div class="line">SparseMatrixType mat(rows,cols);</div>
<div class="line">mat.setFromTriplets(tripletList.begin(), tripletList.end());</div>
<div class="line"><span class="comment">// mat is ready to go!</span></div>
</div><!-- fragment --><p> The <code>std::vector</code> of triplets might contain the elements in arbitrary order, and might even contain duplicated elements that will be summed up by setFromTriplets(). See the <a class="el" href="classEigen_1_1SparseMatrix.html#a5bcf3187e372ff7cea1e8f61152ae49b">SparseMatrix::setFromTriplets()</a> function and class <a class="el" href="classEigen_1_1Triplet.html" title="A small structure to hold a non zero as a triplet (i,j,value). ">Triplet</a> for more details.</p>
<p>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: </p>
<div class="fragment"><div class="line">1: SparseMatrix<double> mat(rows,cols); <span class="comment">// default is column major</span></div>
<div class="line">2: mat.reserve(VectorXi::Constant(cols,6));</div>
<div class="line">3: <span class="keywordflow">for each</span> i,j such that v_ij != 0</div>
<div class="line">4: mat.insert(i,j) = v_ij; <span class="comment">// alternative: mat.coeffRef(i,j) += v_ij;</span></div>
<div class="line">5: mat.makeCompressed(); <span class="comment">// optional</span></div>
</div><!-- fragment --><ul>
<li>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-zeros 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 <code>j-th</code> 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.</li>
<li>The line 4 performs a sorted insertion. In this example, the ideal case is when the <code>j-th</code> column is not full and contains non-zeros whose inner-indices are smaller than <code>i</code>. In this case, this operation boils down to trivial O(1) operation.</li>
<li>When calling insert(i,j) the element <code>i</code> <code></code>,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.</li>
<li>The line 5 suppresses the remaining empty space and transforms the matrix into a compressed column storage.</li>
</ul>
<h1><a class="anchor" id="TutorialSparseFeatureSet"></a>
Supported operators and functions</h1>
<p>Because of their special storage format, sparse matrices cannot offer the same level of flexibility than dense matrices. In <a class="el" href="namespaceEigen.html" title="Namespace containing all symbols from the Eigen library. ">Eigen</a>'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</em> denotes a sparse matrix, <em>sv</em> a sparse vector, <em>dm</em> a dense matrix, and <em>dv</em> a dense vector.</p>
<h2><a class="anchor" id="TutorialSparse_BasicOps"></a>
Basic operations</h2>
<p>Sparse expressions support most of the unary and binary coefficient wise operations: </p>
<div class="fragment"><div class="line">sm1.real() sm1.imag() -sm1 0.5*sm1</div>
<div class="line">sm1+sm2 sm1-sm2 sm1.cwiseProduct(sm2)</div>
</div><!-- fragment --><p> However, a strong restriction is that the storage orders must match. For instance, in the following example: </p>
<div class="fragment"><div class="line">sm4 = sm1 + sm2 + sm3;</div>
</div><!-- fragment --><p> 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 <img class="formulaInl" alt="$ A^T + A $" src="form_202.png"/>, the matrix <img class="formulaInl" alt="$ A^T $" src="form_203.png"/> must be evaluated into a temporary matrix of compatible storage order: </p>
<div class="fragment"><div class="line">SparseMatrix<double> A, B;</div>
<div class="line">B = SparseMatrix<double>(A.transpose()) + A;</div>
</div><!-- fragment --><p>Binary coefficient wise operators can also mix sparse and dense expressions: </p>
<div class="fragment"><div class="line">sm2 = sm1.cwiseProduct(dm1);</div>
<div class="line">dm2 = sm1 + dm1;</div>
</div><!-- fragment --><p>Sparse expressions also support transposition: </p>
<div class="fragment"><div class="line">sm1 = sm2.transpose();</div>
<div class="line">sm1 = sm2.adjoint();</div>
</div><!-- fragment --><p> However, there is no transposeInPlace() method.</p>
<h2><a class="anchor" id="TutorialSparse_Products"></a>
Matrix products</h2>
<p>Eigen supports various kind of sparse matrix products which are summarize below:</p>
<ul>
<li><b>sparse-dense</b>: <div class="fragment"><div class="line">dv2 = sm1 * dv1;</div>
<div class="line">dm2 = dm1 * sm1.adjoint();</div>
<div class="line">dm2 = 2. * sm1 * dm1;</div>
</div><!-- fragment --></li>
<li><b>symmetric</b> <b>sparse-dense</b>. The product of a sparse symmetric matrix with a dense matrix (or vector) can also be optimized by specifying the symmetry with selfadjointView(): <div class="fragment"><div class="line">dm2 = sm1.selfadjointView<>() * dm1; <span class="comment">// if all coefficients of A are stored</span></div>
<div class="line">dm2 = A.selfadjointView<<a class="code" href="group__enums.html#ggab59c1bec446b10af208f977a871d910bae70afef0d3ff7aca74e17e85ff6c9f2e">Upper</a>>() * dm1; <span class="comment">// if only the upper part of A is stored</span></div>
<div class="line">dm2 = A.selfadjointView<<a class="code" href="group__enums.html#ggab59c1bec446b10af208f977a871d910baf886b397626076218462d53d50eb96bc">Lower</a>>() * dm1; <span class="comment">// if only the lower part of A is stored</span></div>
</div><!-- fragment --></li>
<li><b>sparse-sparse</b>. For sparse-sparse products, two different algorithms are available. The default one is conservative and preserve the explicit zeros that might appear: <div class="fragment"><div class="line">sm3 = sm1 * sm2;</div>
<div class="line">sm3 = 4 * sm1.adjoint() * sm2;</div>
</div><!-- fragment --> The second algorithm prunes on the fly the explicit zeros, or the values smaller than a given threshold. It is enabled and controlled through the prune() functions: <div class="fragment"><div class="line">sm3 = (sm1 * sm2).pruned(); <span class="comment">// removes numerical zeros</span></div>
<div class="line">sm3 = (sm1 * sm2).pruned(ref); <span class="comment">// removes elements much smaller than ref</span></div>
<div class="line">sm3 = (sm1 * sm2).pruned(ref,epsilon); <span class="comment">// removes elements smaller than ref*epsilon</span></div>
</div><!-- fragment --></li>
<li><b>permutations</b>. Finally, permutations can be applied to sparse matrices too: <div class="fragment"><div class="line">PermutationMatrix<Dynamic,Dynamic> P = ...;</div>
<div class="line">sm2 = P * sm1;</div>
<div class="line">sm2 = sm1 * P.inverse();</div>
<div class="line">sm2 = sm1.transpose() * P;</div>
</div><!-- fragment --></li>
</ul>
<h2><a class="anchor" id="TutorialSparse_TriangularSelfadjoint"></a>
Triangular and selfadjoint views</h2>
<p>Just as with dense matrices, the triangularView() function can be used to address a triangular part of the matrix, and perform triangular solves with a dense right hand side: </p>
<div class="fragment"><div class="line">dm2 = sm1.triangularView<<a class="code" href="group__enums.html#ggab59c1bec446b10af208f977a871d910baf886b397626076218462d53d50eb96bc">Lower</a>>(dm1);</div>
<div class="line">dv2 = sm1.transpose().triangularView<<a class="code" href="group__enums.html#ggab59c1bec446b10af208f977a871d910bae70afef0d3ff7aca74e17e85ff6c9f2e">Upper</a>>(dv1);</div>
</div><!-- fragment --><p>The selfadjointView() function permits various operations:</p>
<ul>
<li>optimized sparse-dense matrix products: <div class="fragment"><div class="line">dm2 = sm1.selfadjointView<>() * dm1; <span class="comment">// if all coefficients of A are stored</span></div>
<div class="line">dm2 = A.selfadjointView<<a class="code" href="group__enums.html#ggab59c1bec446b10af208f977a871d910bae70afef0d3ff7aca74e17e85ff6c9f2e">Upper</a>>() * dm1; <span class="comment">// if only the upper part of A is stored</span></div>
<div class="line">dm2 = A.selfadjointView<<a class="code" href="group__enums.html#ggab59c1bec446b10af208f977a871d910baf886b397626076218462d53d50eb96bc">Lower</a>>() * dm1; <span class="comment">// if only the lower part of A is stored</span></div>
</div><!-- fragment --></li>
<li>copy of triangular parts: <div class="fragment"><div class="line">sm2 = sm1.selfadjointView<<a class="code" href="group__enums.html#ggab59c1bec446b10af208f977a871d910bae70afef0d3ff7aca74e17e85ff6c9f2e">Upper</a>>(); <span class="comment">// makes a full selfadjoint matrix from the upper triangular part</span></div>
<div class="line">sm2.selfadjointView<<a class="code" href="group__enums.html#ggab59c1bec446b10af208f977a871d910baf886b397626076218462d53d50eb96bc">Lower</a>>() = sm1.selfadjointView<<a class="code" href="group__enums.html#ggab59c1bec446b10af208f977a871d910bae70afef0d3ff7aca74e17e85ff6c9f2e">Upper</a>>(); <span class="comment">// copies the upper triangular part to the lower triangular part</span></div>
</div><!-- fragment --></li>
<li>application of symmetric permutations: <div class="fragment"><div class="line">PermutationMatrix<Dynamic,Dynamic> P = ...;</div>
<div class="line">sm2 = A.selfadjointView<<a class="code" href="group__enums.html#ggab59c1bec446b10af208f977a871d910bae70afef0d3ff7aca74e17e85ff6c9f2e">Upper</a>>().twistedBy(P); <span class="comment">// compute P S P' from the upper triangular part of A, and make it a full matrix</span></div>
<div class="line">sm2.selfadjointView<<a class="code" href="group__enums.html#ggab59c1bec446b10af208f977a871d910baf886b397626076218462d53d50eb96bc">Lower</a>>() = A.selfadjointView<<a class="code" href="group__enums.html#ggab59c1bec446b10af208f977a871d910baf886b397626076218462d53d50eb96bc">Lower</a>>().twistedBy(P); <span class="comment">// compute P S P' from the lower triangular part of A, and then only compute the lower part</span></div>
</div><!-- fragment --></li>
</ul>
<p>Please, refer to the <a class="el" href="group__SparseQuickRefPage.html">Quick Reference </a> guide for the list of supported operations. The list of linear solvers available is <a class="el" href="group__TopicSparseSystems.html">here. </a> </p>
</div><!-- contents -->
</div><!-- doc-content -->
<!-- start footer part -->
<div id="nav-path" class="navpath"><!-- id is needed for treeview function! -->
<ul>
<li class="footer">Generated on Mon Oct 28 2013 11:04:27 for Eigen by
<a href="http://www.doxygen.org/index.html">
<img class="footer" src="doxygen.png" alt="doxygen"/></a> 1.8.5 </li>
</ul>
</div>
<!-- Piwik -->
<!--
<script type="text/javascript">
var pkBaseURL = (("https:" == document.location.protocol) ? "https://stats.sylphide-consulting.com/piwik/" : "http://stats.sylphide-consulting.com/piwik/");
document.write(unescape("%3Cscript src='" + pkBaseURL + "piwik.js' type='text/javascript'%3E%3C/script%3E"));
</script><script type="text/javascript">
try {
var piwikTracker = Piwik.getTracker(pkBaseURL + "piwik.php", 20);
piwikTracker.trackPageView();
piwikTracker.enableLinkTracking();
} catch( err ) {}
</script><noscript><p><img src="http://stats.sylphide-consulting.com/piwik/piwik.php?idsite=20" style="border:0" alt="" /></p></noscript>
-->
<!-- End Piwik Tracking Code -->
</body>
</html>
