blob: 7a29128f3a1df9b50966394d719d6a9f155ab92a [file] [log] [blame]
Carlos Becker9d440052010-06-25 20:16:12 -04001namespace Eigen {
2
Benoit Jacobe078bb22010-06-26 14:00:00 -04003/** \page TutorialMatrixArithmetic Tutorial page 2 - %Matrix and vector arithmetic
Carlos Becker9d440052010-06-25 20:16:12 -04004 \ingroup Tutorial
5
Benoit Jacobe078bb22010-06-26 14:00:00 -04006\li \b Previous: \ref TutorialMatrixClass
7\li \b Next: \ref TutorialArrayClass
8
9This tutorial aims to provide an overview and some details on how to perform arithmetic
10between matrices, vectors and scalars with Eigen.
Carlos Becker9d440052010-06-25 20:16:12 -040011
12\b Table \b of \b contents
Benoit Jacobe078bb22010-06-26 14:00:00 -040013 - \ref TutorialArithmeticIntroduction
14 - \ref TutorialArithmeticMentionCommaInitializer
15 - \ref TutorialArithmeticAddSub
16 - \ref TutorialArithmeticScalarMulDiv
17 - \ref TutorialArithmeticMentionXprTemplates
Gael Guennebaudde1220a2010-06-28 00:05:11 +020018 - \ref TutorialArithmeticTranspose
Benoit Jacobe078bb22010-06-26 14:00:00 -040019 - \ref TutorialArithmeticMatrixMul
20 - \ref TutorialArithmeticDotAndCross
21 - \ref TutorialArithmeticRedux
22 - \ref TutorialArithmeticValidity
Carlos Becker9d440052010-06-25 20:16:12 -040023
Benoit Jacobe078bb22010-06-26 14:00:00 -040024\section TutorialArithmeticIntroduction Introduction
Carlos Becker9d440052010-06-25 20:16:12 -040025
Benoit Jacobe078bb22010-06-26 14:00:00 -040026Eigen offers matrix/vector arithmetic operations either through overloads of common C++ arithmetic operators such as +, -, *,
27or through special methods such as dot(), cross(), etc.
28For the Matrix class (matrices and vectors), operators are only overloaded to support
29linear-algebraic operations. For example, \c matrix1 \c * \c matrix2 means matrix-matrix product,
30and \c vector \c + \c scalar is just not allowed. If you want to perform all kinds of array operations,
31not linear algebra, see \ref TutorialArrayClass "next page".
32
33\section TutorialArithmeticMentionCommaInitializer A note about comma-initialization
34
35In examples below, we'll be initializing matrices with a very convenient device known as the \em comma-initializer.
36For now, it is enough to know this example:
37\include Tutorial_commainit_01.cpp
38Output: \verbinclude Tutorial_commainit_01.out
39The comma initializer is discussed in \ref TutorialAdvancedInitialization "this page".
40
41\section TutorialArithmeticAddSub Addition and substraction
42
43The left hand side and right hand side must, of course, have the same numbers of rows and of columns. They must
44also have the same \c Scalar type, as Eigen doesn't do automatic type promotion. The operators at hand here are:
45\li binary operator + as in \c a+b
46\li binary operator - as in \c a-b
47\li unary operator - as in \c -a
48\li compound operator += as in \c a+=b
49\li compound operator -= as in \c a-=b
50
51Example: \include tut_arithmetic_add_sub.cpp
52Output: \verbinclude tut_arithmetic_add_sub.out
53
54\section TutorialArithmeticScalarMulDiv Scalar multiplication and division
55
56Multiplication and division by a scalar is very simple too. The operators at hand here are:
57\li binary operator * as in \c matrix*scalar
58\li binary operator * as in \c scalar*matrix
59\li binary operator / as in \c matrix/scalar
60\li compound operator *= as in \c matrix*=scalar
61\li compound operator /= as in \c matrix/=scalar
62
63Example: \include tut_arithmetic_scalar_mul_div.cpp
64Output: \verbinclude tut_arithmetic_scalar_mul_div.out
65
66\section TutorialArithmeticMentionXprTemplates A note about expression templates
67
68This is an advanced topic that we explain in \ref TopicEigenExpressionTemplates "this page",
69but it is useful to just mention it now. In Eigen, arithmetic operators such as \c operator+ don't
70perform any computation by themselves, they just return an "expression object" describing the computation to be
71performed. The actual computation happens later, when the whole expression is evaluated, typically in \c operator=.
72While this might sound heavy, any modern optimizing compiler is able to optimize away that abstraction and
73the result is perfectly optimized code. For example, when you do:
Carlos Becker9d440052010-06-25 20:16:12 -040074\code
Benoit Jacobe078bb22010-06-26 14:00:00 -040075VectorXf a(50), b(50), c(50), d(50);
76...
77a = 3*b + 4*c + 5*d;
Carlos Becker9d440052010-06-25 20:16:12 -040078\endcode
Benoit Jacobe078bb22010-06-26 14:00:00 -040079Eigen compiles it to just one for loop, so that the arrays are traversed only once. Simplyfying (e.g. ignoring
80SIMD optimizations), this loop looks like this:
Carlos Becker9d440052010-06-25 20:16:12 -040081\code
Benoit Jacobe078bb22010-06-26 14:00:00 -040082for(int i = 0; i < 50; ++i)
83 a[i] = 3*b[i] + 4*c[i] + 5*d[i];
Carlos Becker9d440052010-06-25 20:16:12 -040084\endcode
Benoit Jacobe078bb22010-06-26 14:00:00 -040085Thus, you should not be afraid of using relatively large arithmetic expressions with Eigen: it only gives Eigen
86more opportunities for optimization.
87
Gael Guennebaudde1220a2010-06-28 00:05:11 +020088\section TutorialArithmeticTranspose Transposition and conjugation
89
90The \c transpose \f$ a^T \f$, \c conjugate \f$ \bar{a} \f$, and the \c adjoint (i.e., conjugate transpose) of the matrix or vector \f$ a \f$, are simply obtained by the functions of the same names.
91
92<table class="tutorial_code"><tr><td>
93Example: \include tut_arithmetic_transpose_conjugate.cpp
94</td>
95<td>
96Output: \include tut_arithmetic_transpose_conjugate.out
97</td></tr></table>
98
99For real matrices, \c conjugate() is a no-operation, and so \c adjoint() is 100% equivalent to \c transpose().
100
101As for basic arithmetic operators, \c transpose and \c adjoint simply return a proxy object without doing the actual transposition. Therefore, <tt>a=a.transpose()</tt> leads to an unexpected result:
102<table class="tutorial_code"><tr><td>
103Example: \include tut_arithmetic_transpose_aliasing.cpp
104</td>
105<td>
106Output: \include tut_arithmetic_transpose_aliasing.out
107</td></tr></table>
108In "debug mode", i.e., when assertions have not been disabled, such common pitfalls are automatically detected. For \em in-place transposition, simply use the transposeInPlace() function:
109<table class="tutorial_code"><tr><td>
110Example: \include tut_arithmetic_transpose_inplace.cpp
111</td>
112<td>
113Output: \include tut_arithmetic_transpose_inplace.out
114</td></tr></table>
115There is also the adjointInPlace() function for complex matrix.
116
Benoit Jacobe078bb22010-06-26 14:00:00 -0400117\section TutorialArithmeticMatrixMul Matrix-matrix and matrix-vector multiplication
118
119Matrix-matrix multiplication is again done with \c operator*. Since vectors are a special
120case of matrices, they are implicitly handled there too, so matrix-vector product is really just a special
121case of matrix-matrix product, and so is vector-vector outer product. Thus, all these cases are handled by just
122two operators:
123\li binary operator * as in \c a*b
124\li compound operator *= as in \c a*=b
125
126Example: \include tut_arithmetic_matrix_mul.cpp
127Output: \verbinclude tut_arithmetic_matrix_mul.out
128
129Note: if you read the above paragraph on expression templates and are worried that doing \c m=m*m might cause
130aliasing issues, be reassured for now: Eigen treats matrix multiplication as a special case and takes care of
131introducing a temporary here, so it will compile \c m=m*m as:
Carlos Becker9d440052010-06-25 20:16:12 -0400132\code
Benoit Jacobe078bb22010-06-26 14:00:00 -0400133tmp = m*m;
134m = tmp;
Carlos Becker9d440052010-06-25 20:16:12 -0400135\endcode
Benoit Jacobe078bb22010-06-26 14:00:00 -0400136For more details on this topic, see \ref TopicEigenExpressionTemplates "this page".
Carlos Becker9d440052010-06-25 20:16:12 -0400137
Benoit Jacobe078bb22010-06-26 14:00:00 -0400138\section TutorialArithmeticDotAndCross Dot product and cross product
Carlos Becker9d440052010-06-25 20:16:12 -0400139
Benoit Jacobe078bb22010-06-26 14:00:00 -0400140The above-discussed \c operator* does not allow to compute dot and cross products. For that, you need the dot() and cross() methods.
141Example: \include tut_arithmetic_dot_cross.cpp
142Output: \verbinclude tut_arithmetic_dot_cross.out
Carlos Becker9d440052010-06-25 20:16:12 -0400143
Benoit Jacobe078bb22010-06-26 14:00:00 -0400144Remember that cross product is only for vectors of size 3. Dot product is for vectors of any sizes.
145When using complex numbers, Eigen's dot product is conjugate-linear in the first variable and linear in the
146second variable.
Carlos Becker9d440052010-06-25 20:16:12 -0400147
Benoit Jacobe078bb22010-06-26 14:00:00 -0400148\section TutorialArithmeticRedux Basic arithmetic reduction operations
149Eigen also provides some reduction operations to obtain values such as the sum or the maximum
150or minimum of all the coefficients in a given matrix or vector.
Carlos Becker9d440052010-06-25 20:16:12 -0400151
Benoit Jacobe078bb22010-06-26 14:00:00 -0400152TODO: convert this from table format to tutorial/examples format.
Carlos Becker9d440052010-06-25 20:16:12 -0400153
154<table class="tutorial_code" align="center">
155<tr><td align="center">\b Reduction \b operation</td><td align="center">\b Usage \b example</td></tr>
156<tr><td>
157Sum of all the coefficients in a matrix</td><td>\code
158MatrixXf m;
159float totalSum = m.sum();\endcode</td></tr>
160<tr><td>
161Maximum coefficient in a matrix</td><td>\code
162MatrixXf m;
163int row, col;
164
165// minimum value will be stored in minValue
166// and the row and column where it was found in row and col,
167// (these two parameters are optional)
168float minValue = m.minCoeff(&row,&col);\endcode</td></tr>
169<tr><td>
170Maximum coefficient in a matrix</td><td>\code
171MatrixXf m;
172int row, col;
173
174// maximum value will be stored in maxValue
175// and the row and column where it was found in row and col,
176// (these two parameters are optional)
177float maxValue = m.maxCoeff(&row,&col);\endcode</td></tr>
178<tr><td>
179Product between all coefficients in a matrix</td><td>\code
180MatrixXf m;
181
182float product = m.prod();\endcode</td></tr>
183<tr><td>
184Mean of coefficients in a matrix</td><td>\code
185MatrixXf m;
186
187float mean = m.mean();\endcode</td></tr>
188<tr><td>
189Matrix's trace</td><td>\code
190MatrixXf m;
191
192float trace = m.trace();\endcode</td></tr>
193</table>
194
Benoit Jacobe078bb22010-06-26 14:00:00 -0400195\subsection TutorialArithmeticValidity Validity of operations
196Eigen checks the validity of the operations that you perform. When possible,
197it checks them at compile-time, producing compilation errors. These error messages can be long and ugly,
198but Eigen writes the important message in UPPERCASE_LETTERS_SO_IT_STANDS_OUT. For example:
199\code
200 Matrix3f m;
201 Vector4f v;
202 v = m*v; // Compile-time error: YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES
203\endcode
Carlos Becker9d440052010-06-25 20:16:12 -0400204
Benoit Jacobe078bb22010-06-26 14:00:00 -0400205Of course, in many cases, for example when checking dynamic sizes, the check cannot be performed at compile time.
206Eigen then uses runtime assertions. This means that executing an illegal operation will result in a crash at runtime,
207with an error message.
Carlos Becker9d440052010-06-25 20:16:12 -0400208
Benoit Jacobe078bb22010-06-26 14:00:00 -0400209\code
210 MatrixXf m(3,3);
211 VectorXf v(4);
212 v = m * v; // Run-time assertion failure here: "invalid matrix product"
213\endcode
214
215For more details on this topic, see \ref TopicAssertions "this page".
216
217\li \b Next: \ref TutorialArrayClass
Carlos Becker9d440052010-06-25 20:16:12 -0400218
219*/
220
221}