1 namespace Eigen { 2 3 /** \page TutorialSparse Tutorial page 9 - Sparse Matrix 4 \ingroup Tutorial 5 6 \li \b Previous: \ref TutorialGeometry 7 \li \b Next: \ref TutorialMapClass 8 9 \b Table \b of \b contents \n 10 - \ref TutorialSparseIntro 11 - \ref TutorialSparseExample "Example" 12 - \ref TutorialSparseSparseMatrix 13 - \ref TutorialSparseFilling 14 - \ref TutorialSparseDirectSolvers 15 - \ref TutorialSparseFeatureSet 16 - \ref TutorialSparse_BasicOps 17 - \ref TutorialSparse_Products 18 - \ref TutorialSparse_TriangularSelfadjoint 19 - \ref TutorialSparse_Submat 20 21 22 <hr> 23 24 Manipulating and solving sparse problems involves various modules which are summarized below: 25 26 <table class="manual"> 27 <tr><th>Module</th><th>Header file</th><th>Contents</th></tr> 28 <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> 29 <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> 30 <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> 31 <tr><td></td><td>\code#include <Eigen/Sparse>\endcode</td><td>Includes all the above modules</td></tr> 32 </table> 33 34 \section TutorialSparseIntro Sparse matrix representation 35 36 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. 37 38 \b The \b %SparseMatrix \b class 39 40 The class SparseMatrix is the main sparse matrix representation of Eigen's sparse module; it offers high performance and low memory usage. 41 It implements a more versatile variant of the widely-used Compressed Column (or Row) Storage scheme. 42 It consists of four compact arrays: 43 - \c Values: stores the coefficient values of the non-zeros. 44 - \c InnerIndices: stores the row (resp. column) indices of the non-zeros. 45 - \c OuterStarts: stores for each column (resp. row) the index of the first non-zero in the previous two arrays. 46 - \c InnerNNZs: stores the number of non-zeros of each column (resp. row). 47 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. 48 The word \c outer refers to the other direction. 49 50 This storage scheme is better explained on an example. The following matrix 51 <table class="manual"> 52 <tr><td> 0</td><td>3</td><td> 0</td><td>0</td><td> 0</td></tr> 53 <tr><td>22</td><td>0</td><td> 0</td><td>0</td><td>17</td></tr> 54 <tr><td> 7</td><td>5</td><td> 0</td><td>1</td><td> 0</td></tr> 55 <tr><td> 0</td><td>0</td><td> 0</td><td>0</td><td> 0</td></tr> 56 <tr><td> 0</td><td>0</td><td>14</td><td>0</td><td> 8</td></tr> 57 </table> 58 59 and one of its possible sparse, \b column \b major representation: 60 <table class="manual"> 61 <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> 62 <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> 63 </table> 64 <table class="manual"> 65 <tr><td>OuterStarts:</td><td>0</td><td>3</td><td>5</td><td>8</td><td>10</td><td>\em 12 </td></tr> 66 <tr><td>InnerNNZs:</td> <td>2</td><td>2</td><td>1</td><td>1</td><td> 2</td><td></td></tr> 67 </table> 68 69 Currently the elements of a given inner vector are guaranteed to be always sorted by increasing inner indices. 70 The \c "_" indicates available free space to quickly insert new elements. 71 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. 72 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. 73 74 The case where no empty space is available is a special case, and is refered as the \em compressed mode. 75 It corresponds to the widely used Compressed Column (or Row) Storage schemes (CCS or CRS). 76 Any SparseMatrix can be turned to this form by calling the SparseMatrix::makeCompressed() function. 77 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]. 78 Therefore, in practice a call to SparseMatrix::makeCompressed() frees this buffer. 79 80 It is worth noting that most of our wrappers to external libraries requires compressed matrices as inputs. 81 82 The results of %Eigen's operations always produces \b compressed sparse matrices. 83 On the other hand, the insertion of a new element into a SparseMatrix converts this later to the \b uncompressed mode. 84 85 Here is the previous matrix represented in compressed mode: 86 <table class="manual"> 87 <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> 88 <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> 89 </table> 90 <table class="manual"> 91 <tr><td>OuterStarts:</td><td>0</td><td>2</td><td>4</td><td>5</td><td>6</td><td>\em 8 </td></tr> 92 </table> 93 94 A SparseVector is a special case of a SparseMatrix where only the \c Values and \c InnerIndices arrays are stored. 95 There is no notion of compressed/uncompressed mode for a SparseVector. 96 97 98 \section TutorialSparseExample First example 99 100 Before describing each individual class, let's start with the following typical example: solving the Lapace equation \f$ \nabla u = 0 \f$ on a regular 2D grid using a finite difference scheme and Dirichlet boundary conditions. 101 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 unknowns (in our case, the values of the pixels), \f$ b \f$ is the right hand side vector resulting from the boundary conditions, and \f$ A \f$ is an \f$ m \times m \f$ matrix containing only a few non-zero elements resulting from the discretization of the Laplacian operator. 102 103 <table class="manual"> 104 <tr><td> 105 \include Tutorial_sparse_example.cpp 106 </td> 107 <td> 108 \image html Tutorial_sparse_example.jpeg 109 </td></tr></table> 110 111 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. 112 113 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. 114 The raw and flat list of non-zero entries is then converted to a true SparseMatrix object \c A. 115 Note that the elements of the list do not have to be sorted, and possible duplicate entries will be summed up. 116 117 The last step consists of effectively solving the assembled problem. 118 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 counterpart for dense objects. 119 120 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. 121 122 Describing 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. 123 124 125 126 127 \section TutorialSparseSparseMatrix The SparseMatrix class 128 129 \b %Matrix \b and \b vector \b properties \n 130 131 The SparseMatrix and SparseVector classes take three template arguments: 132 * the scalar type (e.g., double) 133 * the storage order (ColMajor or RowMajor, the default is RowMajor) 134 * the inner index type (default is \c int). 135 136 As for dense Matrix objects, constructors takes the size of the object. 137 Here are some examples: 138 139 \code 140 SparseMatrix<std::complex<float> > mat(1000,2000); // declares a 1000x2000 column-major compressed sparse matrix of complex<float> 141 SparseMatrix<double,RowMajor> mat(1000,2000); // declares a 1000x2000 row-major compressed sparse matrix of double 142 SparseVector<std::complex<float> > vec(1000); // declares a column sparse vector of complex<float> of size 1000 143 SparseVector<double,RowMajor> vec(1000); // declares a row sparse vector of double of size 1000 144 \endcode 145 146 In the rest of the tutorial, \c mat and \c vec represent any sparse-matrix and sparse-vector objects, respectively. 147 148 The dimensions of a matrix can be queried using the following functions: 149 <table class="manual"> 150 <tr><td>Standard \n dimensions</td><td>\code 151 mat.rows() 152 mat.cols()\endcode</td> 153 <td>\code 154 vec.size() \endcode</td> 155 </tr> 156 <tr><td>Sizes along the \n inner/outer dimensions</td><td>\code 157 mat.innerSize() 158 mat.outerSize()\endcode</td> 159 <td></td> 160 </tr> 161 <tr><td>Number of non \n zero coefficients</td><td>\code 162 mat.nonZeros() \endcode</td> 163 <td>\code 164 vec.nonZeros() \endcode</td></tr> 165 </table> 166 167 168 \b Iterating \b over \b the \b nonzero \b coefficients \n 169 170 Random access to the elements of a sparse object can be done through the \c coeffRef(i,j) function. 171 However, this function involves a quite expensive binary search. 172 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. 173 Here is an example: 174 <table class="manual"> 175 <tr><td> 176 \code 177 SparseMatrix<double> mat(rows,cols); 178 for (int k=0; k<mat.outerSize(); ++k) 179 for (SparseMatrix<double>::InnerIterator it(mat,k); it; ++it) 180 { 181 it.value(); 182 it.row(); // row index 183 it.col(); // col index (here it is equal to k) 184 it.index(); // inner index, here it is equal to it.row() 185 } 186 \endcode 187 </td><td> 188 \code 189 SparseVector<double> vec(size); 190 for (SparseVector<double>::InnerIterator it(vec); it; ++it) 191 { 192 it.value(); // == vec[ it.index() ] 193 it.index(); 194 } 195 \endcode 196 </td></tr> 197 </table> 198 For a writable expression, the referenced value can be modified using the valueRef() function. 199 If the type of the sparse matrix or vector depends on a template parameter, then the \c typename keyword is 200 required to indicate that \c InnerIterator denotes a type; see \ref TopicTemplateKeyword for details. 201 202 203 \section TutorialSparseFilling Filling a sparse matrix 204 205 Because of the special storage scheme of a SparseMatrix, special care has to be taken when adding new nonzero entries. 206 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. 207 208 The simplest way to create a sparse matrix while guaranteeing good performance is thus to first build a list of so-called \em triplets, and then convert it to a SparseMatrix. 209 210 Here is a typical usage example: 211 \code 212 typedef Eigen::Triplet<double> T; 213 std::vector<T> tripletList; 214 triplets.reserve(estimation_of_entries); 215 for(...) 216 { 217 // ... 218 tripletList.push_back(T(i,j,v_ij)); 219 } 220 SparseMatrixType mat(rows,cols); 221 mat.setFromTriplets(tripletList.begin(), tripletList.end()); 222 // mat is ready to go! 223 \endcode 224 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(). 225 See the SparseMatrix::setFromTriplets() function and class Triplet for more details. 226 227 228 In some cases, however, slightly higher performance, and lower memory consumption can be reached by directly inserting the non-zeros into the destination matrix. 229 A typical scenario of this approach is illustrated bellow: 230 \code 231 1: SparseMatrix<double> mat(rows,cols); // default is column major 232 2: mat.reserve(VectorXi::Constant(cols,6)); 233 3: for each i,j such that v_ij != 0 234 4: mat.insert(i,j) = v_ij; // alternative: mat.coeffRef(i,j) += v_ij; 235 5: mat.makeCompressed(); // optional 236 \endcode 237 238 - 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 \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. 239 - 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. 240 - 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. 241 - The line 5 suppresses the remaining empty space and transforms the matrix into a compressed column storage. 242 243 244 \section TutorialSparseDirectSolvers Solving linear problems 245 246 %Eigen currently provides a limited set of built-in solvers, as well as wrappers to external solver libraries. 247 They are summarized in the following table: 248 249 <table class="manual"> 250 <tr><th>Class</th><th>Module</th><th>Solver kind</th><th>Matrix kind</th><th>Features related to performance</th> 251 <th>Dependencies,License</th><th class="width20em"><p>Notes</p></th></tr> 252 <tr><td>SimplicialLLT </td><td>\link SparseCholesky_Module SparseCholesky \endlink</td><td>Direct LLt factorization</td><td>SPD</td><td>Fill-in reducing</td> 253 <td>built-in, LGPL</td> 254 <td>SimplicialLDLT is often preferable</td></tr> 255 <tr><td>SimplicialLDLT </td><td>\link SparseCholesky_Module SparseCholesky \endlink</td><td>Direct LDLt factorization</td><td>SPD</td><td>Fill-in reducing</td> 256 <td>built-in, LGPL</td> 257 <td>Recommended for very sparse and not too large problems (e.g., 2D Poisson eq.)</td></tr> 258 <tr><td>ConjugateGradient</td><td>\link IterativeLinearSolvers_Module IterativeLinearSolvers \endlink</td><td>Classic iterative CG</td><td>SPD</td><td>Preconditionning</td> 259 <td>built-in, LGPL</td> 260 <td>Recommended for large symmetric problems (e.g., 3D Poisson eq.)</td></tr> 261 <tr><td>BiCGSTAB</td><td>\link IterativeLinearSolvers_Module IterativeLinearSolvers \endlink</td><td>Iterative stabilized bi-conjugate gradient</td><td>Square</td><td>Preconditionning</td> 262 <td>built-in, LGPL</td> 263 <td>Might not always converge</td></tr> 264 265 266 <tr><td>PastixLLT \n PastixLDLT \n PastixLU</td><td>\link PaStiXSupport_Module PaStiXSupport \endlink</td><td>Direct LLt, LDLt, LU factorizations</td><td>SPD \n SPD \n Square</td><td>Fill-in reducing, Leverage fast dense algebra, Multithreading</td> 267 <td>Requires the <a href="http://pastix.gforge.inria.fr">PaStiX</a> package, \b CeCILL-C </td> 268 <td>optimized for tough problems and symmetric patterns</td></tr> 269 <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> 270 <td>Requires the <a href="http://www.cise.ufl.edu/research/sparse/SuiteSparse/">SuiteSparse</a> package, \b GPL </td> 271 <td></td></tr> 272 <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> 273 <td>Requires the <a href="http://www.cise.ufl.edu/research/sparse/SuiteSparse/">SuiteSparse</a> package, \b GPL </td> 274 <td></td></tr> 275 <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> 276 <td>Requires the <a href="http://crd-legacy.lbl.gov/~xiaoye/SuperLU/">SuperLU</a> library, (BSD-like)</td> 277 <td></td></tr> 278 </table> 279 280 Here \c SPD means symmetric positive definite. 281 282 All these solvers follow the same general concept. 283 Here is a typical and general example: 284 \code 285 #include <Eigen/RequiredModuleName> 286 // ... 287 SparseMatrix<double> A; 288 // fill A 289 VectorXd b, x; 290 // fill b 291 // solve Ax = b 292 SolverClassName<SparseMatrix<double> > solver; 293 solver.compute(A); 294 if(solver.info()!=Succeeded) { 295 // decomposition failed 296 return; 297 } 298 x = solver.solve(b); 299 if(solver.info()!=Succeeded) { 300 // solving failed 301 return; 302 } 303 // solve for another right hand side: 304 x1 = solver.solve(b1); 305 \endcode 306 307 For \c SPD solvers, a second optional template argument allows to specify which triangular part have to be used, e.g.: 308 309 \code 310 #include <Eigen/IterativeLinearSolvers> 311 312 ConjugateGradient<SparseMatrix<double>, Eigen::Upper> solver; 313 x = solver.compute(A).solve(b); 314 \endcode 315 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. 316 317 In the case where multiple problems with the same sparcity pattern have to be solved, then the "compute" step can be decomposed as follow: 318 \code 319 SolverClassName<SparseMatrix<double> > solver; 320 solver.analyzePattern(A); // for this step the numerical values of A are not used 321 solver.factorize(A); 322 x1 = solver.solve(b1); 323 x2 = solver.solve(b2); 324 ... 325 A = ...; // modify the values of the nonzeros of A, the nonzeros pattern must stay unchanged 326 solver.factorize(A); 327 x1 = solver.solve(b1); 328 x2 = solver.solve(b2); 329 ... 330 \endcode 331 The compute() method is equivalent to calling both analyzePattern() and factorize(). 332 333 Finally, each solver provides some specific features, such as determinant, access to the factors, controls of the iterations, and so on. 334 More details are availble in the documentations of the respective classes. 335 336 337 \section TutorialSparseFeatureSet Supported operators and functions 338 339 Because of their special storage format, sparse matrices cannot offer the same level of flexbility than dense matrices. 340 In Eigen's sparse module we chose to expose only the subset of the dense matrix API which can be efficiently implemented. 341 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. 342 343 \subsection TutorialSparse_BasicOps Basic operations 344 345 %Sparse expressions support most of the unary and binary coefficient wise operations: 346 \code 347 sm1.real() sm1.imag() -sm1 0.5*sm1 348 sm1+sm2 sm1-sm2 sm1.cwiseProduct(sm2) 349 \endcode 350 However, a strong restriction is that the storage orders must match. For instance, in the following example: 351 \code 352 sm4 = sm1 + sm2 + sm3; 353 \endcode 354 sm1, sm2, and sm3 must all be row-major or all column major. 355 On the other hand, there is no restriction on the target matrix sm4. 356 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: 357 \code 358 SparseMatrix<double> A, B; 359 B = SparseMatrix<double>(A.transpose()) + A; 360 \endcode 361 362 Binary coefficient wise operators can also mix sparse and dense expressions: 363 \code 364 sm2 = sm1.cwiseProduct(dm1); 365 dm2 = sm1 + dm1; 366 \endcode 367 368 369 %Sparse expressions also support transposition: 370 \code 371 sm1 = sm2.transpose(); 372 sm1 = sm2.adjoint(); 373 \endcode 374 However, there is no transposeInPlace() method. 375 376 377 \subsection TutorialSparse_Products Matrix products 378 379 %Eigen supports various kind of sparse matrix products which are summarize below: 380 - \b sparse-dense: 381 \code 382 dv2 = sm1 * dv1; 383 dm2 = dm1 * sm1.adjoint(); 384 dm2 = 2. * sm1 * dm1; 385 \endcode 386 - \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(): 387 \code 388 dm2 = sm1.selfadjointView<>() * dm1; // if all coefficients of A are stored 389 dm2 = A.selfadjointView<Upper>() * dm1; // if only the upper part of A is stored 390 dm2 = A.selfadjointView<Lower>() * dm1; // if only the lower part of A is stored 391 \endcode 392 - \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: 393 \code 394 sm3 = sm1 * sm2; 395 sm3 = 4 * sm1.adjoint() * sm2; 396 \endcode 397 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: 398 \code 399 sm3 = (sm1 * sm2).prune(); // removes numerical zeros 400 sm3 = (sm1 * sm2).prune(ref); // removes elements much smaller than ref 401 sm3 = (sm1 * sm2).prune(ref,epsilon); // removes elements smaller than ref*epsilon 402 \endcode 403 404 - \b permutations. Finally, permutations can be applied to sparse matrices too: 405 \code 406 PermutationMatrix<Dynamic,Dynamic> P = ...; 407 sm2 = P * sm1; 408 sm2 = sm1 * P.inverse(); 409 sm2 = sm1.transpose() * P; 410 \endcode 411 412 413 \subsection TutorialSparse_TriangularSelfadjoint Triangular and selfadjoint views 414 415 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: 416 \code 417 dm2 = sm1.triangularView<Lower>(dm1); 418 dv2 = sm1.transpose().triangularView<Upper>(dv1); 419 \endcode 420 421 The selfadjointView() function permits various operations: 422 - optimized sparse-dense matrix products: 423 \code 424 dm2 = sm1.selfadjointView<>() * dm1; // if all coefficients of A are stored 425 dm2 = A.selfadjointView<Upper>() * dm1; // if only the upper part of A is stored 426 dm2 = A.selfadjointView<Lower>() * dm1; // if only the lower part of A is stored 427 \endcode 428 - copy of triangular parts: 429 \code 430 sm2 = sm1.selfadjointView<Upper>(); // makes a full selfadjoint matrix from the upper triangular part 431 sm2.selfadjointView<Lower>() = sm1.selfadjointView<Upper>(); // copies the upper triangular part to the lower triangular part 432 \endcode 433 - application of symmetric permutations: 434 \code 435 PermutationMatrix<Dynamic,Dynamic> P = ...; 436 sm2 = A.selfadjointView<Upper>().twistedBy(P); // compute P S P' from the upper triangular part of A, and make it a full matrix 437 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 438 \endcode 439 440 \subsection TutorialSparse_Submat Sub-matrices 441 442 %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: 443 \code 444 sm1.innerVector(j); // returns an expression of the j-th column (resp. row) of the matrix if sm1 is col-major (resp. row-major) 445 sm1.innerVectors(j, nb); // returns an expression of the nb columns (resp. row) starting from the j-th column (resp. row) 446 // of the matrix if sm1 is col-major (resp. row-major) 447 sm1.middleRows(j, nb); // for row major matrices only, get a range of nb rows 448 sm1.middleCols(j, nb); // for column major matrices only, get a range of nb columns 449 \endcode 450 451 \li \b Next: \ref TutorialMapClass 452 453 */ 454 455 } 456
