mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-04-22 17:49:36 +08:00
update doc
This commit is contained in:
Gael Guennebaud
parent
7579360672
commit
6713c75fac
@ -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>
|
||||
|
||||
|
||||
|
||||
@ -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>
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user