blob: c2654aed26aac18197500f4a99242cf62b075cc2 [file] [log] [blame]
Narayan Kamathc981c482012-11-02 10:59:05 +00001namespace Eigen {
2
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -07003/** \eigenManualPage TopicAliasing Aliasing
Narayan Kamathc981c482012-11-02 10:59:05 +00004
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -07005In %Eigen, aliasing refers to assignment statement in which the same matrix (or array or vector) appears on the
Narayan Kamathc981c482012-11-02 10:59:05 +00006left and on the right of the assignment operators. Statements like <tt>mat = 2 * mat;</tt> or <tt>mat =
7mat.transpose();</tt> exhibit aliasing. The aliasing in the first example is harmless, but the aliasing in the
8second example leads to unexpected results. This page explains what aliasing is, when it is harmful, and what
9to do about it.
10
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -070011\eigenAutoToc
Narayan Kamathc981c482012-11-02 10:59:05 +000012
13
14\section TopicAliasingExamples Examples
15
16Here is a simple example exhibiting aliasing:
17
18<table class="example">
19<tr><th>Example</th><th>Output</th></tr>
20<tr><td>
21\include TopicAliasing_block.cpp
22</td>
23<td>
24\verbinclude TopicAliasing_block.out
25</td></tr></table>
26
27The output is not what one would expect. The problem is the assignment
28\code
29mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);
30\endcode
31This assignment exhibits aliasing: the coefficient \c mat(1,1) appears both in the block
32<tt>mat.bottomRightCorner(2,2)</tt> on the left-hand side of the assignment and the block
33<tt>mat.topLeftCorner(2,2)</tt> on the right-hand side. After the assignment, the (2,2) entry in the bottom
34right corner should have the value of \c mat(1,1) before the assignment, which is 5. However, the output shows
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -070035that \c mat(2,2) is actually 1. The problem is that %Eigen uses lazy evaluation (see
Narayan Kamathc981c482012-11-02 10:59:05 +000036\ref TopicEigenExpressionTemplates) for <tt>mat.topLeftCorner(2,2)</tt>. The result is similar to
37\code
38mat(1,1) = mat(0,0);
39mat(1,2) = mat(0,1);
40mat(2,1) = mat(1,0);
41mat(2,2) = mat(1,1);
42\endcode
43Thus, \c mat(2,2) is assigned the \e new value of \c mat(1,1) instead of the old value. The next section
44explains how to solve this problem by calling \link DenseBase::eval() eval()\endlink.
45
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -070046Aliasing occurs more naturally when trying to shrink a matrix. For example, the expressions <tt>vec =
47vec.head(n)</tt> and <tt>mat = mat.block(i,j,r,c)</tt> exhibit aliasing.
48
49In general, aliasing cannot be detected at compile time: if \c mat in the first example were a bit bigger,
50then the blocks would not overlap, and there would be no aliasing problem. However, %Eigen does detect some
51instances of aliasing, albeit at run time. The following example exhibiting aliasing was mentioned in \ref
52TutorialMatrixArithmetic :
Narayan Kamathc981c482012-11-02 10:59:05 +000053
54<table class="example">
55<tr><th>Example</th><th>Output</th></tr>
56<tr><td>
57\include tut_arithmetic_transpose_aliasing.cpp
58</td>
59<td>
60\verbinclude tut_arithmetic_transpose_aliasing.out
61</td></tr></table>
62
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -070063Again, the output shows the aliasing issue. However, by default %Eigen uses a run-time assertion to detect this
Narayan Kamathc981c482012-11-02 10:59:05 +000064and exits with a message like
65
66\verbatim
67void Eigen::DenseBase<Derived>::checkTransposeAliasing(const OtherDerived&) const
68[with OtherDerived = Eigen::Transpose<Eigen::Matrix<int, 2, 2, 0, 2, 2> >, Derived = Eigen::Matrix<int, 2, 2, 0, 2, 2>]:
69Assertion `(!internal::check_transpose_aliasing_selector<Scalar,internal::blas_traits<Derived>::IsTransposed,OtherDerived>::run(internal::extract_data(derived()), other))
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -070070&& "aliasing detected during transposition, use transposeInPlace() or evaluate the rhs into a temporary using .eval()"' failed.
Narayan Kamathc981c482012-11-02 10:59:05 +000071\endverbatim
72
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -070073The user can turn %Eigen's run-time assertions like the one to detect this aliasing problem off by defining the
Narayan Kamathc981c482012-11-02 10:59:05 +000074EIGEN_NO_DEBUG macro, and the above program was compiled with this macro turned off in order to illustrate the
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -070075aliasing problem. See \ref TopicAssertions for more information about %Eigen's run-time assertions.
Narayan Kamathc981c482012-11-02 10:59:05 +000076
77
78\section TopicAliasingSolution Resolving aliasing issues
79
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -070080If you understand the cause of the aliasing issue, then it is obvious what must happen to solve it: %Eigen has
Narayan Kamathc981c482012-11-02 10:59:05 +000081to evaluate the right-hand side fully into a temporary matrix/array and then assign it to the left-hand
82side. The function \link DenseBase::eval() eval() \endlink does precisely that.
83
84For example, here is the corrected version of the first example above:
85
86<table class="example">
87<tr><th>Example</th><th>Output</th></tr>
88<tr><td>
89\include TopicAliasing_block_correct.cpp
90</td>
91<td>
92\verbinclude TopicAliasing_block_correct.out
93</td></tr></table>
94
95Now, \c mat(2,2) equals 5 after the assignment, as it should be.
96
97The same solution also works for the second example, with the transpose: simply replace the line
98<tt>a = a.transpose();</tt> with <tt>a = a.transpose().eval();</tt>. However, in this common case there is a
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -070099better solution. %Eigen provides the special-purpose function
Narayan Kamathc981c482012-11-02 10:59:05 +0000100\link DenseBase::transposeInPlace() transposeInPlace() \endlink which replaces a matrix by its transpose.
101This is shown below:
102
103<table class="example">
104<tr><th>Example</th><th>Output</th></tr>
105<tr><td>
106\include tut_arithmetic_transpose_inplace.cpp
107</td>
108<td>
109\verbinclude tut_arithmetic_transpose_inplace.out
110</td></tr></table>
111
112If an xxxInPlace() function is available, then it is best to use it, because it indicates more clearly what you
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -0700113are doing. This may also allow %Eigen to optimize more aggressively. These are some of the xxxInPlace()
Narayan Kamathc981c482012-11-02 10:59:05 +0000114functions provided:
115
116<table class="manual">
117<tr><th>Original function</th><th>In-place function</th></tr>
118<tr> <td> MatrixBase::adjoint() </td> <td> MatrixBase::adjointInPlace() </td> </tr>
119<tr class="alt"> <td> DenseBase::reverse() </td> <td> DenseBase::reverseInPlace() </td> </tr>
120<tr> <td> LDLT::solve() </td> <td> LDLT::solveInPlace() </td> </tr>
121<tr class="alt"> <td> LLT::solve() </td> <td> LLT::solveInPlace() </td> </tr>
122<tr> <td> TriangularView::solve() </td> <td> TriangularView::solveInPlace() </td> </tr>
123<tr class="alt"> <td> DenseBase::transpose() </td> <td> DenseBase::transposeInPlace() </td> </tr>
124</table>
125
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -0700126In the special case where a matrix or vector is shrunk using an expression like <tt>vec = vec.head(n)</tt>,
127you can use \link PlainObjectBase::conservativeResize() conservativeResize() \endlink.
128
Narayan Kamathc981c482012-11-02 10:59:05 +0000129
130\section TopicAliasingCwise Aliasing and component-wise operations
131
132As explained above, it may be dangerous if the same matrix or array occurs on both the left-hand side and the
133right-hand side of an assignment operator, and it is then often necessary to evaluate the right-hand side
134explicitly. However, applying component-wise operations (such as matrix addition, scalar multiplication and
135array multiplication) is safe.
136
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -0700137The following example has only component-wise operations. Thus, there is no need for \link DenseBase::eval()
138eval() \endlink even though the same matrix appears on both sides of the assignments.
Narayan Kamathc981c482012-11-02 10:59:05 +0000139
140<table class="example">
141<tr><th>Example</th><th>Output</th></tr>
142<tr><td>
143\include TopicAliasing_cwise.cpp
144</td>
145<td>
146\verbinclude TopicAliasing_cwise.out
147</td></tr></table>
148
149In general, an assignment is safe if the (i,j) entry of the expression on the right-hand side depends only on
150the (i,j) entry of the matrix or array on the left-hand side and not on any other entries. In that case it is
151not necessary to evaluate the right-hand side explicitly.
152
153
154\section TopicAliasingMatrixMult Aliasing and matrix multiplication
155
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -0700156Matrix multiplication is the only operation in %Eigen that assumes aliasing by default. Thus, if \c matA is a
157matrix, then the statement <tt>matA = matA * matA;</tt> is safe. All other operations in %Eigen assume that
Narayan Kamathc981c482012-11-02 10:59:05 +0000158there are no aliasing problems, either because the result is assigned to a different matrix or because it is a
159component-wise operation.
160
161<table class="example">
162<tr><th>Example</th><th>Output</th></tr>
163<tr><td>
164\include TopicAliasing_mult1.cpp
165</td>
166<td>
167\verbinclude TopicAliasing_mult1.out
168</td></tr></table>
169
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -0700170However, this comes at a price. When executing the expression <tt>matA = matA * matA</tt>, %Eigen evaluates the
171product in a temporary matrix which is assigned to \c matA after the computation. This is fine. But %Eigen does
Narayan Kamathc981c482012-11-02 10:59:05 +0000172the same when the product is assigned to a different matrix (e.g., <tt>matB = matA * matA</tt>). In that case,
173it is more efficient to evaluate the product directly into \c matB instead of evaluating it first into a
174temporary matrix and copying that matrix to \c matB.
175
176The user can indicate with the \link MatrixBase::noalias() noalias()\endlink function that there is no
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -0700177aliasing, as follows: <tt>matB.noalias() = matA * matA</tt>. This allows %Eigen to evaluate the matrix product
Narayan Kamathc981c482012-11-02 10:59:05 +0000178<tt>matA * matA</tt> directly into \c matB.
179
180<table class="example">
181<tr><th>Example</th><th>Output</th></tr>
182<tr><td>
183\include TopicAliasing_mult2.cpp
184</td>
185<td>
186\verbinclude TopicAliasing_mult2.out
187</td></tr></table>
188
189Of course, you should not use \c noalias() when there is in fact aliasing taking place. If you do, then you
190may get wrong results:
191
192<table class="example">
193<tr><th>Example</th><th>Output</th></tr>
194<tr><td>
195\include TopicAliasing_mult3.cpp
196</td>
197<td>
198\verbinclude TopicAliasing_mult3.out
199</td></tr></table>
200
201
202\section TopicAliasingSummary Summary
203
204Aliasing occurs when the same matrix or array coefficients appear both on the left- and the right-hand side of
205an assignment operator.
206 - Aliasing is harmless with coefficient-wise computations; this includes scalar multiplication and matrix or
207 array addition.
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -0700208 - When you multiply two matrices, %Eigen assumes that aliasing occurs. If you know that there is no aliasing,
Narayan Kamathc981c482012-11-02 10:59:05 +0000209 then you can use \link MatrixBase::noalias() noalias()\endlink.
Carlos Hernandez7faaa9f2014-08-05 17:53:32 -0700210 - In all other situations, %Eigen assumes that there is no aliasing issue and thus gives the wrong result if
Narayan Kamathc981c482012-11-02 10:59:05 +0000211 aliasing does in fact occur. To prevent this, you have to use \link DenseBase::eval() eval() \endlink or
212 one of the xxxInPlace() functions.
213
214*/
215}