Home | History | Annotate | Download | only in doc 
      1 namespace Eigen {
      2 
      3 /** \eigenManualPage TutorialMatrixArithmetic Matrix and vector arithmetic
      4 
      5 This page aims to provide an overview and some details on how to perform arithmetic
      6 between matrices, vectors and scalars with Eigen.
      7 
      8 \eigenAutoToc
      9 
     10 \section TutorialArithmeticIntroduction Introduction
     11 
     12 Eigen offers matrix/vector arithmetic operations either through overloads of common C++ arithmetic operators such as +, -, *,
     13 or through special methods such as dot(), cross(), etc.
     14 For the Matrix class (matrices and vectors), operators are only overloaded to support
     15 linear-algebraic operations. For example, \c matrix1 \c * \c matrix2 means matrix-matrix product,
     16 and \c vector \c + \c scalar is just not allowed. If you want to perform all kinds of array operations,
     17 not linear algebra, see the \ref TutorialArrayClass "next page".
     18 
     19 \section TutorialArithmeticAddSub Addition and subtraction
     20 
     21 The left hand side and right hand side must, of course, have the same numbers of rows and of columns. They must
     22 also have the same \c Scalar type, as Eigen doesn't do automatic type promotion. The operators at hand here are:
     23 \li binary operator + as in \c a+b
     24 \li binary operator - as in \c a-b
     25 \li unary operator - as in \c -a
     26 \li compound operator += as in \c a+=b
     27 \li compound operator -= as in \c a-=b
     28 
     29 <table class="example">
     30 <tr><th>Example:</th><th>Output:</th></tr>
     31 <tr><td>
     32 \include tut_arithmetic_add_sub.cpp
     33 </td>
     34 <td>
     35 \verbinclude tut_arithmetic_add_sub.out
     36 </td></tr></table>
     37 
     38 \section TutorialArithmeticScalarMulDiv Scalar multiplication and division
     39 
     40 Multiplication and division by a scalar is very simple too. The operators at hand here are:
     41 \li binary operator * as in \c matrix*scalar
     42 \li binary operator * as in \c scalar*matrix
     43 \li binary operator / as in \c matrix/scalar
     44 \li compound operator *= as in \c matrix*=scalar
     45 \li compound operator /= as in \c matrix/=scalar
     46 
     47 <table class="example">
     48 <tr><th>Example:</th><th>Output:</th></tr>
     49 <tr><td>
     50 \include tut_arithmetic_scalar_mul_div.cpp
     51 </td>
     52 <td>
     53 \verbinclude tut_arithmetic_scalar_mul_div.out
     54 </td></tr></table>
     55 
     56 
     57 \section TutorialArithmeticMentionXprTemplates A note about expression templates
     58 
     59 This is an advanced topic that we explain on \ref TopicEigenExpressionTemplates "this page",
     60 but it is useful to just mention it now. In Eigen, arithmetic operators such as \c operator+ don't
     61 perform any computation by themselves, they just return an "expression object" describing the computation to be
     62 performed. The actual computation happens later, when the whole expression is evaluated, typically in \c operator=.
     63 While this might sound heavy, any modern optimizing compiler is able to optimize away that abstraction and
     64 the result is perfectly optimized code. For example, when you do:
     65 \code
     66 VectorXf a(50), b(50), c(50), d(50);
     67 ...
     68 a = 3*b + 4*c + 5*d;
     69 \endcode
     70 Eigen compiles it to just one for loop, so that the arrays are traversed only once. Simplifying (e.g. ignoring
     71 SIMD optimizations), this loop looks like this:
     72 \code
     73 for(int i = 0; i < 50; ++i)
     74   a[i] = 3*b[i] + 4*c[i] + 5*d[i];
     75 \endcode
     76 Thus, you should not be afraid of using relatively large arithmetic expressions with Eigen: it only gives Eigen
     77 more opportunities for optimization.
     78 
     79 \section TutorialArithmeticTranspose Transposition and conjugation
     80 
     81 The transpose \f$ a^T \f$, conjugate \f$ \bar{a} \f$, and adjoint (i.e., conjugate transpose) \f$ a^* \f$ of a matrix or vector \f$ a \f$ are obtained by the member functions \link DenseBase::transpose() transpose()\endlink, \link MatrixBase::conjugate() conjugate()\endlink, and \link MatrixBase::adjoint() adjoint()\endlink, respectively.
     82 
     83 <table class="example">
     84 <tr><th>Example:</th><th>Output:</th></tr>
     85 <tr><td>
     86 \include tut_arithmetic_transpose_conjugate.cpp
     87 </td>
     88 <td>
     89 \verbinclude tut_arithmetic_transpose_conjugate.out
     90 </td></tr></table>
     91 
     92 For real matrices, \c conjugate() is a no-operation, and so \c adjoint() is equivalent to \c transpose().
     93 
     94 As for basic arithmetic operators, \c transpose() and \c adjoint() simply return a proxy object without doing the actual transposition. If you do <tt>b = a.transpose()</tt>, then the transpose is evaluated at the same time as the result is written into \c b. However, there is a complication here. If you do <tt>a = a.transpose()</tt>, then Eigen starts writing the result into \c a before the evaluation of the transpose is finished. Therefore, the instruction <tt>a = a.transpose()</tt> does not replace \c a with its transpose, as one would expect:
     95 <table class="example">
     96 <tr><th>Example:</th><th>Output:</th></tr>
     97 <tr><td>
     98 \include tut_arithmetic_transpose_aliasing.cpp
     99 </td>
    100 <td>
    101 \verbinclude tut_arithmetic_transpose_aliasing.out
    102 </td></tr></table>
    103 This is the so-called \ref TopicAliasing "aliasing issue". In "debug mode", i.e., when \ref TopicAssertions "assertions" have not been disabled, such common pitfalls are automatically detected. 
    104 
    105 For \em in-place transposition, as for instance in <tt>a = a.transpose()</tt>, simply use the \link DenseBase::transposeInPlace() transposeInPlace()\endlink  function:
    106 <table class="example">
    107 <tr><th>Example:</th><th>Output:</th></tr>
    108 <tr><td>
    109 \include tut_arithmetic_transpose_inplace.cpp
    110 </td>
    111 <td>
    112 \verbinclude tut_arithmetic_transpose_inplace.out
    113 </td></tr></table>
    114 There is also the \link MatrixBase::adjointInPlace() adjointInPlace()\endlink function for complex matrices.
    115 
    116 \section TutorialArithmeticMatrixMul Matrix-matrix and matrix-vector multiplication
    117 
    118 Matrix-matrix multiplication is again done with \c operator*. Since vectors are a special
    119 case of matrices, they are implicitly handled there too, so matrix-vector product is really just a special
    120 case of matrix-matrix product, and so is vector-vector outer product. Thus, all these cases are handled by just
    121 two operators:
    122 \li binary operator * as in \c a*b
    123 \li compound operator *= as in \c a*=b (this multiplies on the right: \c a*=b is equivalent to <tt>a = a*b</tt>)
    124 
    125 <table class="example">
    126 <tr><th>Example:</th><th>Output:</th></tr>
    127 <tr><td>
    128 \include tut_arithmetic_matrix_mul.cpp
    129 </td>
    130 <td>
    131 \verbinclude tut_arithmetic_matrix_mul.out
    132 </td></tr></table>
    133 
    134 Note: if you read the above paragraph on expression templates and are worried that doing \c m=m*m might cause
    135 aliasing issues, be reassured for now: Eigen treats matrix multiplication as a special case and takes care of
    136 introducing a temporary here, so it will compile \c m=m*m as:
    137 \code
    138 tmp = m*m;
    139 m = tmp;
    140 \endcode
    141 If you know your matrix product can be safely evaluated into the destination matrix without aliasing issue, then you can use the \link MatrixBase::noalias() noalias()\endlink function to avoid the temporary, e.g.:
    142 \code
    143 c.noalias() += a * b;
    144 \endcode
    145 For more details on this topic, see the page on \ref TopicAliasing "aliasing".
    146 
    147 \b Note: for BLAS users worried about performance, expressions such as <tt>c.noalias() -= 2 * a.adjoint() * b;</tt> are fully optimized and trigger a single gemm-like function call.
    148 
    149 \section TutorialArithmeticDotAndCross Dot product and cross product
    150 
    151 For dot product and cross product, you need the \link MatrixBase::dot() dot()\endlink and \link MatrixBase::cross() cross()\endlink methods. Of course, the dot product can also be obtained as a 1x1 matrix as u.adjoint()*v.
    152 <table class="example">
    153 <tr><th>Example:</th><th>Output:</th></tr>
    154 <tr><td>
    155 \include tut_arithmetic_dot_cross.cpp
    156 </td>
    157 <td>
    158 \verbinclude tut_arithmetic_dot_cross.out
    159 </td></tr></table>
    160 
    161 Remember that cross product is only for vectors of size 3. Dot product is for vectors of any sizes.
    162 When using complex numbers, Eigen's dot product is conjugate-linear in the first variable and linear in the
    163 second variable.
    164 
    165 \section TutorialArithmeticRedux Basic arithmetic reduction operations
    166 Eigen also provides some reduction operations to reduce a given matrix or vector to a single value such as the sum (computed by \link DenseBase::sum() sum()\endlink), product (\link DenseBase::prod() prod()\endlink), or the maximum (\link DenseBase::maxCoeff() maxCoeff()\endlink) and minimum (\link DenseBase::minCoeff() minCoeff()\endlink) of all its coefficients.
    167 
    168 <table class="example">
    169 <tr><th>Example:</th><th>Output:</th></tr>
    170 <tr><td>
    171 \include tut_arithmetic_redux_basic.cpp
    172 </td>
    173 <td>
    174 \verbinclude tut_arithmetic_redux_basic.out
    175 </td></tr></table>
    176 
    177 The \em trace of a matrix, as returned by the function \link MatrixBase::trace() trace()\endlink, is the sum of the diagonal coefficients and can also be computed as efficiently using <tt>a.diagonal().sum()</tt>, as we will see later on.
    178 
    179 There also exist variants of the \c minCoeff and \c maxCoeff functions returning the coordinates of the respective coefficient via the arguments:
    180 
    181 <table class="example">
    182 <tr><th>Example:</th><th>Output:</th></tr>
    183 <tr><td>
    184 \include tut_arithmetic_redux_minmax.cpp
    185 </td>
    186 <td>
    187 \verbinclude tut_arithmetic_redux_minmax.out
    188 </td></tr></table>
    189 
    190 
    191 \section TutorialArithmeticValidity Validity of operations
    192 Eigen checks the validity of the operations that you perform. When possible,
    193 it checks them at compile time, producing compilation errors. These error messages can be long and ugly,
    194 but Eigen writes the important message in UPPERCASE_LETTERS_SO_IT_STANDS_OUT. For example:
    195 \code
    196   Matrix3f m;
    197   Vector4f v;
    198   v = m*v;      // Compile-time error: YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES
    199 \endcode
    200 
    201 Of course, in many cases, for example when checking dynamic sizes, the check cannot be performed at compile time.
    202 Eigen then uses runtime assertions. This means that the program will abort with an error message when executing an illegal operation if it is run in "debug mode", and it will probably crash if assertions are turned off.
    203 
    204 \code
    205   MatrixXf m(3,3);
    206   VectorXf v(4);
    207   v = m * v; // Run-time assertion failure here: "invalid matrix product"
    208 \endcode
    209 
    210 For more details on this topic, see \ref TopicAssertions "this page".
    211 
    212 */
    213 
    214 }
    215