diff options
| -rw-r--r-- | doc/C01_QuickStartGuide.dox | 86 | ||||
| -rw-r--r-- | doc/I02_HiPerformance.dox | 97 |
2 files changed, 165 insertions, 18 deletions
diff --git a/doc/C01_QuickStartGuide.dox b/doc/C01_QuickStartGuide.dox index 2f8e9f5c8..3d98e14f5 100644 --- a/doc/C01_QuickStartGuide.dox +++ b/doc/C01_QuickStartGuide.dox @@ -24,6 +24,7 @@ namespace Eigen { - \ref TutorialCoreTransposeAdjoint - \ref TutorialCoreDotNorm - \ref TutorialCoreTriangularMatrix + - \ref TutorialCoreSelfadjointMatrix - \ref TutorialCoreSpecialTopics \n @@ -577,35 +578,88 @@ vec1.normalize();\endcode <a href="#" class="top">top</a>\section TutorialCoreTriangularMatrix Dealing with triangular matrices -Read/write access to special parts of a matrix can be achieved. See \link MatrixBase::part() const this \endlink for read access and \link MatrixBase::part() this \endlink for write access.. +Currently, Eigen does not provide any explcit triangular matrix, with storage class. Instead, we +can reference a triangular part of a square matrix or expression to perform special treatment on it. +This is achieved by the class TriangularView and the MatrixBase::triangularView template function. +Note that the opposite triangular part of the matrix is never referenced, and so it can, e.g., store +a second triangular matrix. <table class="tutorial_code"> <tr><td> -Extract triangular matrices \n from a given matrix m: +Reference a read/write triangular part of a given \n +matrix (or expression) m with optional unit diagonal: </td><td>\code -m.part<Eigen::UpperTriangular>() -m.part<Eigen::StrictlyUpperTriangular>() -m.part<Eigen::UnitUpperTriangular>() -m.part<Eigen::LowerTriangular>() -m.part<Eigen::StrictlyLowerTriangular>() -m.part<Eigen::UnitLowerTriangular>()\endcode +m.triangularView<Eigen::UpperTriangular>() +m.triangularView<Eigen::UnitUpperTriangular>() +m.triangularView<Eigen::LowerTriangular>() +m.triangularView<Eigen::UnitLowerTriangular>()\endcode </td></tr> <tr><td> -Write to triangular parts \n of a matrix m: +Writting to a specific triangular part:\n (only the referenced triangular part is evaluated) </td><td>\code -m1.part<Eigen::UpperTriangular>() = m2; -m1.part<Eigen::StrictlyUpperTriangular>() = m2; -m1.part<Eigen::LowerTriangular>() = m2; -m1.part<Eigen::StrictlyLowerTriangular>() = m2;\endcode +m1.triangularView<Eigen::LowerTriangular>() = m2 + m3 \endcode </td></tr> <tr><td> -Special: take advantage of symmetry \n (selfadjointness) when copying \n an expression into a matrix +Convertion to a dense matrix setting the opposite triangular part to zero: </td><td>\code -m.part<Eigen::SelfAdjoint>() = someSelfadjointMatrix; -m1.part<Eigen::SelfAdjoint>() = m2 + m2.adjoint(); // m2 + m2.adjoint() is selfadjoint \endcode +m2 = m1.triangularView<Eigen::UnitUpperTriangular>()\endcode </td></tr> +<tr><td> +Products: +</td><td>\code +m3 += s1 * m1.adjoint().triangularView<Eigen::UnitUpperTriangular>() * m2 +m3 -= s1 * m2.conjugate() * m1.adjoint().triangularView<Eigen::LowerTriangular>() \endcode +</td></tr> +<tr><td> +Solving linear equations:\n(\f$ m_2 := m_1^{-1} m_2 \f$) +</td><td>\code +m1.triangularView<Eigen::UnitLowerTriangular>().solveInPlace(m2) +m1.adjoint().triangularView<Eigen::UpperTriangular>().solveInPlace(m2)\endcode +</td></tr> +</table> +<a href="#" class="top">top</a>\section TutorialCoreSelfadjointMatrix Dealing with symmetric/selfadjoint matrices + +Just as for triangular matrix, you can reference any triangular part of a square matrix to see it a selfadjoint +matrix to perform special and optimized operations. Again the opposite triangular is never referenced and can be +used to store other information. + +<table class="tutorial_code"> +<tr><td> +Conversion to a dense matrix: +</td><td>\code +m2 = m.selfadjointView<Eigen::LowerTriangular>();\endcode +</td></tr> +<tr><td> +Product with another general matrix or vector: +</td><td>\code +m3 = s1 * m1.conjugate().selfadjointView<Eigen::UpperTriangular>() * m3; +m3 -= s1 * m3.adjoint() * m1.selfadjointView<Eigen::UpperTriangular>();\endcode +</td></tr> +<tr><td> +Rank 1 and rank K update: +</td><td>\code +// fast version of m1 += s1 * m2 * m2.adjoint(): +m1.selfadjointView<Eigen::UpperTriangular>().rankUpdate(m2,s1); +// fast version of m1 -= m2.adjoint() * m2: +m1.selfadjointView<Eigen::LowerTriangular>().rankUpdate(m2.adjoint(),-1); \endcode +</td></tr> +<tr><td> +Rank 2 update: (\f$ m += s u v^* + s v u^* \f$) +</td><td>\code +m.selfadjointView<Eigen::UpperTriangular>().rankUpdate(u,v,s); +\endcode +</td></tr> +<tr><td> +Solving linear equations:\n(\f$ m_2 := m_1^{-1} m_2 \f$) +</td><td>\code +// via a standard Cholesky factorization +m1.selfadjointView<Eigen::UpperTriangular>().llt().solveInPlace(m2); +// via a Cholesky factorization with pivoting +m1.selfadjointView<Eigen::UpperTriangular>().ldlt().solveInPlace(m2); +\endcode +</td></tr> </table> diff --git a/doc/I02_HiPerformance.dox b/doc/I02_HiPerformance.dox index 012e7d71b..8f23b5d19 100644 --- a/doc/I02_HiPerformance.dox +++ b/doc/I02_HiPerformance.dox @@ -3,6 +3,99 @@ namespace Eigen { /** \page HiPerformance Advanced - Using Eigen with high performance +In general achieving good performance with Eigen does no require any special effort: +simply write your expressions in the most high level way. This is especially true +for small fixed size matrices. For large matrices, however, it might useful to +take some care when writing your expressions in order to minimize useless evaluations +and optimize the performance. +In this page we will give a brief overview of the Eigen's internal mechanism to simplify +and evaluate complex expressions, and discuss the current limitations. +In particular we will focus on expressions matching level 2 and 3 BLAS routines, i.e, +all kind of matrix products and triangular solvers. + +Indeed, in Eigen we have implemented a set of highly optimized routines which are very similar +to BLAS's ones. Unlike BLAS, those routines are made available to user via a high level and +natural API. Each of these routines can perform in a single evaluation a wide variety of expressions. +Given an expression, the challenge is then to map it to a minimal set of primitives. +As explained latter, this mechanism has some limitations, and knowing them will allow +you to write faster code by making your expressions more Eigen friendly. + +\section GEMM General Matrix-Matrix product (GEMM) + +Let's start with the most common primitive: the matrix product of general dense matrices. +In the BLAS world this corresponds to the GEMM routine. Our equivalent primitive can +perform the following operation: +\f$ C += \alpha op1(A) * op2(B) \f$ +where A, B, and C are column and/or row major matrices (or sub-matrices), +alpha is a scalar value, and op1, op2 can be transpose, adjoint, conjugate, or the identity. +When Eigen detects a matrix product, it analyzes both sides of the product to extract a +unique scalar factor alpha, and for each side its effective storage (order and shape) and conjugate state. +More precisely each side is simplified by iteratively removing trivial expressions such as scalar multiple, +negate and conjugate. Transpose and Block expressions are not evaluated and only modify the storage order +and shape. All other expressions are immediately evaluated. +For instance, the following expression: +\code m1 -= (s1 * m2.adjoint() * (-(s3*m3).conjugate()*s2)).lazy() \endcode +is automatically simplified to: +\code m1 += (s1*s2*conj(s3)) * m2.adjoint() * m3.conjugate() \endcode +which exactly matches our GEMM routine. + +\subsection GEMM_Limitations Limitations +Unfortunately, this simplification mechanism is not perfect yet and not all expressions which could be +handled by a single GEMM-like call are correctly detected. +<table class="tutorial_code"> +<tr> +<td>Not optimal expression</td> +<td>Evaluated as</td> +<td>Optimal version (single evaluation)</td> +<td>Comments</td> +</tr> +<tr> +<td>\code m1 += m2 * m3; \endcode</td> +<td>\code temp = m2 * m3; m1 += temp; \endcode</td> +<td>\code m1 += (m2 * m3).lazy(); \endcode</td> +<td>Use .lazy() to tell Eigen the result and right-hand-sides do not alias.</td> +</tr> +<tr> +<td>\code m1 += (s1 * (m2 * m3)).lazy(); \endcode</td> +<td>\code temp = (m2 * m3).lazy(); m1 += s1 * temp; \endcode</td> +<td>\code m1 += (s1 * m2 * m3).lazy(); \endcode</td> +<td>This is because m2 * m3 is immediately evaluated by the scalar product. <br> + Make sure the matrix product is the top most expression.</td> +</tr> +<tr> +<td>\code m1 = m1 + m2 * m3; \endcode</td> +<td>\code temp = (m2 * m3).lazy(); m1 = m1 + temp; \endcode</td> +<td>\code m1 += (m2 * m3).lazy(); \endcode</td> +<td>Here there is no way to detect at compile time that the two m1 are the same, + and so the matrix product will be immediately evaluated.</td> +</tr> +<tr> +<td>\code m1 += ((s1 * m2).transpose() * m3).lazy(); \endcode</td> +<td>\code temp = (s1*m2).transpose(); m1 = (temp * m3).lazy(); \endcode</td> +<td>\code m1 += (s1 * m2.transpose() * m3).lazy(); \endcode</td> +<td>This is because our expression analyzer stops at the first expression which cannot + be converted to a scalar multiple of a conjugate and therefore the nested scalar + multiple cannot be properly extracted.</td> +</tr> +<tr> +<td>\code m1 += (m2.conjugate().transpose() * m3).lazy(); \endcode</td> +<td>\code temp = m2.conjugate().transpose(); m1 += (temp * m3).lazy(); \endcode</td> +<td>\code m1 += (m2.adjoint() * m3).lazy(); \endcode</td> +<td>Same reason. Use .adjoint() or .transpose().conjugate()</td> +</tr> +<tr> +<td>\code m1 += ((s1*m2).block(....) * m3).lazy(); \endcode</td> +<td>\code temp = (s1*m2).block(....); m1 += (temp * m3).lazy(); \endcode</td> +<td>\code m1 += (s1 * m2.block(....) * m3).lazy(); \endcode</td> +<td>Same reason.</td> +</tr> +</table> + +Of course all these remarks hold for all other kind of products that we will describe in the following paragraphs. + + + + <table class="tutorial_code"> <tr> <td>BLAS equivalent routine</td> @@ -20,7 +113,7 @@ namespace Eigen { <td>GEMM</td> <td>m1 += s1 * m2.adjoint() * m3</td> <td>m1 += (s1 * m2).adjoint() * m3</td> -<td>This is because our expression analyser stops at the first transpose expression and cannot extract the nested scalar multiple.</td> +<td>This is because our expression analyzer stops at the first transpose expression and cannot extract the nested scalar multiple.</td> </tr> <tr> <td>GEMM</td> @@ -36,7 +129,7 @@ namespace Eigen { </tr> <tr> <td>SYR</td> -<td>m.sefadjointView<LowerTriangular>().rankUpdate(v,s)</td> +<td>m.seductive<LowerTriangular>().rankUpdate(v,s)</td> <td></td> <td>Computes m += s * v * v.adjoint()</td> </tr> |