blob: 10156d575934790767f67bbe870a08cc0668b95d [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
18 - \ref TutorialArithmeticMatrixMul
19 - \ref TutorialArithmeticDotAndCross
20 - \ref TutorialArithmeticRedux
21 - \ref TutorialArithmeticValidity
Carlos Becker9d440052010-06-25 20:16:12 -040022
Benoit Jacobe078bb22010-06-26 14:00:00 -040023\section TutorialArithmeticIntroduction Introduction
Carlos Becker9d440052010-06-25 20:16:12 -040024
Benoit Jacobe078bb22010-06-26 14:00:00 -040025Eigen offers matrix/vector arithmetic operations either through overloads of common C++ arithmetic operators such as +, -, *,
26or through special methods such as dot(), cross(), etc.
27For the Matrix class (matrices and vectors), operators are only overloaded to support
28linear-algebraic operations. For example, \c matrix1 \c * \c matrix2 means matrix-matrix product,
29and \c vector \c + \c scalar is just not allowed. If you want to perform all kinds of array operations,
30not linear algebra, see \ref TutorialArrayClass "next page".
31
32\section TutorialArithmeticMentionCommaInitializer A note about comma-initialization
33
34In examples below, we'll be initializing matrices with a very convenient device known as the \em comma-initializer.
35For now, it is enough to know this example:
36\include Tutorial_commainit_01.cpp
37Output: \verbinclude Tutorial_commainit_01.out
38The comma initializer is discussed in \ref TutorialAdvancedInitialization "this page".
39
40\section TutorialArithmeticAddSub Addition and substraction
41
42The left hand side and right hand side must, of course, have the same numbers of rows and of columns. They must
43also have the same \c Scalar type, as Eigen doesn't do automatic type promotion. The operators at hand here are:
44\li binary operator + as in \c a+b
45\li binary operator - as in \c a-b
46\li unary operator - as in \c -a
47\li compound operator += as in \c a+=b
48\li compound operator -= as in \c a-=b
49
50Example: \include tut_arithmetic_add_sub.cpp
51Output: \verbinclude tut_arithmetic_add_sub.out
52
53\section TutorialArithmeticScalarMulDiv Scalar multiplication and division
54
55Multiplication and division by a scalar is very simple too. The operators at hand here are:
56\li binary operator * as in \c matrix*scalar
57\li binary operator * as in \c scalar*matrix
58\li binary operator / as in \c matrix/scalar
59\li compound operator *= as in \c matrix*=scalar
60\li compound operator /= as in \c matrix/=scalar
61
62Example: \include tut_arithmetic_scalar_mul_div.cpp
63Output: \verbinclude tut_arithmetic_scalar_mul_div.out
64
65\section TutorialArithmeticMentionXprTemplates A note about expression templates
66
67This is an advanced topic that we explain in \ref TopicEigenExpressionTemplates "this page",
68but it is useful to just mention it now. In Eigen, arithmetic operators such as \c operator+ don't
69perform any computation by themselves, they just return an "expression object" describing the computation to be
70performed. The actual computation happens later, when the whole expression is evaluated, typically in \c operator=.
71While this might sound heavy, any modern optimizing compiler is able to optimize away that abstraction and
72the result is perfectly optimized code. For example, when you do:
Carlos Becker9d440052010-06-25 20:16:12 -040073\code
Benoit Jacobe078bb22010-06-26 14:00:00 -040074VectorXf a(50), b(50), c(50), d(50);
75...
76a = 3*b + 4*c + 5*d;
Carlos Becker9d440052010-06-25 20:16:12 -040077\endcode
Benoit Jacobe078bb22010-06-26 14:00:00 -040078Eigen compiles it to just one for loop, so that the arrays are traversed only once. Simplyfying (e.g. ignoring
79SIMD optimizations), this loop looks like this:
Carlos Becker9d440052010-06-25 20:16:12 -040080\code
Benoit Jacobe078bb22010-06-26 14:00:00 -040081for(int i = 0; i < 50; ++i)
82 a[i] = 3*b[i] + 4*c[i] + 5*d[i];
Carlos Becker9d440052010-06-25 20:16:12 -040083\endcode
Benoit Jacobe078bb22010-06-26 14:00:00 -040084Thus, you should not be afraid of using relatively large arithmetic expressions with Eigen: it only gives Eigen
85more opportunities for optimization.
86
87\section TutorialArithmeticMatrixMul Matrix-matrix and matrix-vector multiplication
88
89Matrix-matrix multiplication is again done with \c operator*. Since vectors are a special
90case of matrices, they are implicitly handled there too, so matrix-vector product is really just a special
91case of matrix-matrix product, and so is vector-vector outer product. Thus, all these cases are handled by just
92two operators:
93\li binary operator * as in \c a*b
94\li compound operator *= as in \c a*=b
95
96Example: \include tut_arithmetic_matrix_mul.cpp
97Output: \verbinclude tut_arithmetic_matrix_mul.out
98
99Note: if you read the above paragraph on expression templates and are worried that doing \c m=m*m might cause
100aliasing issues, be reassured for now: Eigen treats matrix multiplication as a special case and takes care of
101introducing a temporary here, so it will compile \c m=m*m as:
Carlos Becker9d440052010-06-25 20:16:12 -0400102\code
Benoit Jacobe078bb22010-06-26 14:00:00 -0400103tmp = m*m;
104m = tmp;
Carlos Becker9d440052010-06-25 20:16:12 -0400105\endcode
Benoit Jacobe078bb22010-06-26 14:00:00 -0400106For more details on this topic, see \ref TopicEigenExpressionTemplates "this page".
Carlos Becker9d440052010-06-25 20:16:12 -0400107
Benoit Jacobe078bb22010-06-26 14:00:00 -0400108\section TutorialArithmeticDotAndCross Dot product and cross product
Carlos Becker9d440052010-06-25 20:16:12 -0400109
Benoit Jacobe078bb22010-06-26 14:00:00 -0400110The above-discussed \c operator* does not allow to compute dot and cross products. For that, you need the dot() and cross() methods.
111Example: \include tut_arithmetic_dot_cross.cpp
112Output: \verbinclude tut_arithmetic_dot_cross.out
Carlos Becker9d440052010-06-25 20:16:12 -0400113
Benoit Jacobe078bb22010-06-26 14:00:00 -0400114Remember that cross product is only for vectors of size 3. Dot product is for vectors of any sizes.
115When using complex numbers, Eigen's dot product is conjugate-linear in the first variable and linear in the
116second variable.
Carlos Becker9d440052010-06-25 20:16:12 -0400117
Benoit Jacobe078bb22010-06-26 14:00:00 -0400118\section TutorialArithmeticRedux Basic arithmetic reduction operations
119Eigen also provides some reduction operations to obtain values such as the sum or the maximum
120or minimum of all the coefficients in a given matrix or vector.
Carlos Becker9d440052010-06-25 20:16:12 -0400121
Benoit Jacobe078bb22010-06-26 14:00:00 -0400122TODO: convert this from table format to tutorial/examples format.
Carlos Becker9d440052010-06-25 20:16:12 -0400123
124<table class="tutorial_code" align="center">
125<tr><td align="center">\b Reduction \b operation</td><td align="center">\b Usage \b example</td></tr>
126<tr><td>
127Sum of all the coefficients in a matrix</td><td>\code
128MatrixXf m;
129float totalSum = m.sum();\endcode</td></tr>
130<tr><td>
131Maximum coefficient in a matrix</td><td>\code
132MatrixXf m;
133int row, col;
134
135// minimum value will be stored in minValue
136// and the row and column where it was found in row and col,
137// (these two parameters are optional)
138float minValue = m.minCoeff(&row,&col);\endcode</td></tr>
139<tr><td>
140Maximum coefficient in a matrix</td><td>\code
141MatrixXf m;
142int row, col;
143
144// maximum value will be stored in maxValue
145// and the row and column where it was found in row and col,
146// (these two parameters are optional)
147float maxValue = m.maxCoeff(&row,&col);\endcode</td></tr>
148<tr><td>
149Product between all coefficients in a matrix</td><td>\code
150MatrixXf m;
151
152float product = m.prod();\endcode</td></tr>
153<tr><td>
154Mean of coefficients in a matrix</td><td>\code
155MatrixXf m;
156
157float mean = m.mean();\endcode</td></tr>
158<tr><td>
159Matrix's trace</td><td>\code
160MatrixXf m;
161
162float trace = m.trace();\endcode</td></tr>
163</table>
164
Benoit Jacobe078bb22010-06-26 14:00:00 -0400165\subsection TutorialArithmeticValidity Validity of operations
166Eigen checks the validity of the operations that you perform. When possible,
167it checks them at compile-time, producing compilation errors. These error messages can be long and ugly,
168but Eigen writes the important message in UPPERCASE_LETTERS_SO_IT_STANDS_OUT. For example:
169\code
170 Matrix3f m;
171 Vector4f v;
172 v = m*v; // Compile-time error: YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES
173\endcode
Carlos Becker9d440052010-06-25 20:16:12 -0400174
Benoit Jacobe078bb22010-06-26 14:00:00 -0400175Of course, in many cases, for example when checking dynamic sizes, the check cannot be performed at compile time.
176Eigen then uses runtime assertions. This means that executing an illegal operation will result in a crash at runtime,
177with an error message.
Carlos Becker9d440052010-06-25 20:16:12 -0400178
Benoit Jacobe078bb22010-06-26 14:00:00 -0400179\code
180 MatrixXf m(3,3);
181 VectorXf v(4);
182 v = m * v; // Run-time assertion failure here: "invalid matrix product"
183\endcode
184
185For more details on this topic, see \ref TopicAssertions "this page".
186
187\li \b Next: \ref TutorialArrayClass
Carlos Becker9d440052010-06-25 20:16:12 -0400188
189*/
190
191}