From 44a527dfa50ce9c473cbf1f446b8b6f406d4bc91 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 4 Feb 2009 09:44:44 +0000 Subject: * classify and sort the doxygen "related pages" by tweaking the filename and adding 2 categories: Troubleshooting and Advanced * use the EXCLUDE_SYMBOLS to clean the class list (insteaded of a homemade bash script) * remove the broken "exemple list" * re-structure the unsupported directory as mentionned in the ML and integrate the doc as follow: - snippets of the unsupported directory are directly imported from the main snippets/CMakefile.txt (no need to duplicate code) - add a top level "Unsupported modules" group - unsupported modules have to defined their own sub group and nest it using \ingroup Unsupported_modules - then a pair of //@{ //@} will put everything in the submodule - this is just a proposal ! --- Mainpage.dox | 4 +- disabled/buildexamplelist.sh | 15 + disabled/cleanhierarchy.sh | 5 + doc/B01_Experimental.dox | 55 ++++ doc/C01_QuickStartGuide.dox | 565 ++++++++++++++++++++++++++++++++ doc/C03_TutorialGeometry.dox | 253 ++++++++++++++ doc/C05_TutorialLinearAlgebra.dox | 188 +++++++++++ doc/C07_TutorialSparse.dox | 240 ++++++++++++++ doc/CMakeLists.txt | 2 - doc/CustomizingEigen.dox | 148 --------- doc/D01_StlContainers.dox | 52 +++ doc/D07_PassingByValue.dox | 40 +++ doc/D09_StructHavingEigenMembers.dox | 149 +++++++++ doc/D11_UnalignedArrayAssert.dox | 77 +++++ doc/Doxyfile.in | 6 +- doc/Experimental.dox | 55 ---- doc/FixedSizeVectorizable.dox | 38 --- doc/I00_CustomizingEigen.dox | 148 +++++++++ doc/I01_TopicLazyEvaluation.dox | 67 ++++ doc/I03_InsideEigenExample.dox | 500 ++++++++++++++++++++++++++++ doc/I05_FixedSizeVectorizable.dox | 38 +++ doc/InsideEigenExample.dox | 500 ---------------------------- doc/PassingByValue.dox | 40 --- doc/QuickStartGuide.dox | 565 -------------------------------- doc/StlContainers.dox | 52 --- doc/StructHavingEigenMembers.dox | 149 --------- doc/TopicLazyEvaluation.dox | 67 ---- doc/TutorialGeometry.dox | 253 -------------- doc/TutorialLinearAlgebra.dox | 188 ----------- doc/TutorialSparse.dox | 240 -------------- doc/UnalignedArrayAssert.dox | 77 ----- doc/buildexamplelist.sh | 15 - doc/cleanhierarchy.sh | 5 - doc/snippets/CMakeLists.txt | 40 +-- unsupported/AdolcSupport/AdolcForward.h | 96 ------ unsupported/Eigen/AdolcForward | 120 +++++++ unsupported/doc/unsupported.dox | 15 + 37 files changed, 2554 insertions(+), 2513 deletions(-) create mode 100755 disabled/buildexamplelist.sh create mode 100755 disabled/cleanhierarchy.sh create mode 100644 doc/B01_Experimental.dox create mode 100644 doc/C01_QuickStartGuide.dox create mode 100644 doc/C03_TutorialGeometry.dox create mode 100644 doc/C05_TutorialLinearAlgebra.dox create mode 100644 doc/C07_TutorialSparse.dox delete mode 100644 doc/CustomizingEigen.dox create mode 100644 doc/D01_StlContainers.dox create mode 100644 doc/D07_PassingByValue.dox create mode 100644 doc/D09_StructHavingEigenMembers.dox create mode 100644 doc/D11_UnalignedArrayAssert.dox delete mode 100644 doc/Experimental.dox delete mode 100644 doc/FixedSizeVectorizable.dox create mode 100644 doc/I00_CustomizingEigen.dox create mode 100644 doc/I01_TopicLazyEvaluation.dox create mode 100644 doc/I03_InsideEigenExample.dox create mode 100644 doc/I05_FixedSizeVectorizable.dox delete mode 100644 doc/InsideEigenExample.dox delete mode 100644 doc/PassingByValue.dox delete mode 100644 doc/QuickStartGuide.dox delete mode 100644 doc/StlContainers.dox delete mode 100644 doc/StructHavingEigenMembers.dox delete mode 100644 doc/TopicLazyEvaluation.dox delete mode 100644 doc/TutorialGeometry.dox delete mode 100644 doc/TutorialLinearAlgebra.dox delete mode 100644 doc/TutorialSparse.dox delete mode 100644 doc/UnalignedArrayAssert.dox delete mode 100755 doc/buildexamplelist.sh delete mode 100755 doc/cleanhierarchy.sh delete mode 100644 unsupported/AdolcSupport/AdolcForward.h create mode 100644 unsupported/Eigen/AdolcForward create mode 100644 unsupported/doc/unsupported.dox diff --git a/Mainpage.dox b/Mainpage.dox index 468045dce..b5b7f30fa 100644 --- a/Mainpage.dox +++ b/Mainpage.dox @@ -54,7 +54,7 @@ // DOXYGEN_SET_WARN_IF_UNDOCUMENTED = NO // DOXYGEN_SET_WARN_NO_PARAMDOC = NO -// DOXYGEN_SET_INPUT = @topdir@/eigen2/Eigen @topdir@/eigen2/doc @topdir@/eigen2/build/doc +// DOXYGEN_SET_INPUT = @topdir@/eigen2/Eigen @topdir@/eigen2/doc @topdir@/eigen2/build/doc @topdir@/eigen2/unsupported/Eigen // DOXYGEN_SET_EXCLUDE = *.sh *.in // DOXYGEN_SET_EXAMPLE_PATH = @topdir@/eigen2/doc/snippets/ @topdir@/eigen2/doc/examples/ @topdir@/eigen2/build/doc/examples/ @topdir@/eigen2/build/doc/snippets/ @@ -63,6 +63,8 @@ // DOXYGEN_SET_RECURSIVE = NO // DOXYGEN_SET_FILTER_SOURCE_FILES = YES +// DOXYGEN_EXCLUDE_SYMBOLS = MatrixBase<* MapBase<* RotationBase<* Matrix<* + // DOXYGEN_SET_SOURCE_BROWSER = NO // DOXYGEN_SET_INLINE_SOURCES = NO // DOXYGEN_SET_STRIP_CODE_COMMENTS = YES diff --git a/disabled/buildexamplelist.sh b/disabled/buildexamplelist.sh new file mode 100755 index 000000000..7d47ae39d --- /dev/null +++ b/disabled/buildexamplelist.sh @@ -0,0 +1,15 @@ +#!/bin/sh + +echo "namespace Eigen {" +echo "/** \page ExampleList" +echo "

Selected list of examples

" + +grep \\addexample $1/Eigen/src/*/*.h -R | cut -d \\ -f 2- | \ +while read example; +do +anchor=`echo "$example" | cut -d " " -f 2` +text=`echo "$example" | cut -d " " -f 4-` +echo "\\\li \\\ref $anchor \"$text\"" +done +echo "*/" +echo "}" \ No newline at end of file diff --git a/disabled/cleanhierarchy.sh b/disabled/cleanhierarchy.sh new file mode 100755 index 000000000..fe44891b9 --- /dev/null +++ b/disabled/cleanhierarchy.sh @@ -0,0 +1,5 @@ +#!/bin/sh + +sed -i 's/^.li.*MatrixBase\<.*gt.*a.$/ /g' $1 +sed -i 's/^.li.*MapBase\<.*gt.*a.$/ /g' $1 +sed -i 's/^.li.*RotationBase\<.*gt.*a.$/ /g' $1 diff --git a/doc/B01_Experimental.dox b/doc/B01_Experimental.dox new file mode 100644 index 000000000..6d8b90d5a --- /dev/null +++ b/doc/B01_Experimental.dox @@ -0,0 +1,55 @@ +namespace Eigen { + +/** \page Experimental Experimental parts of Eigen + +\b Table \b of \b contents + - \ref summary + - \ref modules + - \ref core + +\section summary Summary + +With the 2.0 release, Eigen's API is, to a large extent, stable. However, we wish to retain the freedom to make API incompatible changes. To that effect, we call many parts of Eigen "experimental" which means that they are not subject to API stability guarantee. + +Our goal is that for the 2.1 release (expected in July 2009) most of these parts become API-stable too. + +We are aware that API stability is a major concern for our users. That's why it's a priority for us to reach it, but at the same time we're being serious about not calling Eigen API-stable too early. + +Experimental features may at any time: +\li be removed; +\li be subject to an API incompatible change; +\li introduce API or ABI incompatible changes in your own code if you let them affect your API or ABI. + +\section modules Experimental modules + +The following modules are considered entirely experimental, and we make no firm API stability guarantee about them for the time being: +\li SVD +\li QR +\li Cholesky +\li Sparse +\li Geometry (this one should be mostly stable, but it's a little too early to make a formal guarantee) + +\section core Experimental parts of the Core module + +In the Core module, the only classes subject to ABI stability guarantee (meaning that you can use it for data members in your public ABI) is: +\li Matrix +\li Map + +All other classes offer no ABI guarantee, e.g. the layout of their data can be changed. + +The only classes subject to (even partial) API stability guarantee (meaning that you can safely construct and use objects) are: +\li MatrixBase : partial API stability (see below) +\li Matrix : full API stability (except for experimental stuff inherited from MatrixBase) +\li Map : full API stability (except for experimental stuff inherited from MatrixBase) + +All other classes offer no direct API guarantee, e.g. their methods can be changed; however notice that most classes inherit MatrixBase and that this is where most of their API comes from -- so in practice most of the API is stable. + +A few MatrixBase methods are considered experimental, hence not part of any API stability guarantee: +\li all methods documented as internal +\li all methods hidden in the Doxygen documentation +\li all methods marked as experimental +\li all methods defined in experimental modules + +*/ + +} diff --git a/doc/C01_QuickStartGuide.dox b/doc/C01_QuickStartGuide.dox new file mode 100644 index 000000000..d7bec3bb7 --- /dev/null +++ b/doc/C01_QuickStartGuide.dox @@ -0,0 +1,565 @@ +namespace Eigen { + +/** \page TutorialCore Tutorial 1/3 - Core features + \ingroup Tutorial + +
\ref index "Overview" + | \b Core \b features + | \ref TutorialGeometry "Geometry" + | \ref TutorialAdvancedLinearAlgebra "Advanced linear algebra" + | \ref TutorialSparse "Sparse matrix" +
+ +\b Table \b of \b contents + - \ref TutorialCoreGettingStarted + - \ref TutorialCoreSimpleExampleFixedSize + - \ref TutorialCoreSimpleExampleDynamicSize + - \ref TutorialCoreMatrixTypes + - \ref TutorialCoreCoefficients + - \ref TutorialCoreMatrixInitialization + - \ref TutorialCoreArithmeticOperators + - \ref TutorialCoreReductions + - \ref TutorialCoreMatrixBlocks + - \ref TutorialCoreDiagonalMatrices + - \ref TutorialCoreTransposeAdjoint + - \ref TutorialCoreDotNorm + - \ref TutorialCoreTriangularMatrix + - \ref TutorialCoreSpecialTopics +\n + +
+ +\section TutorialCoreGettingStarted Getting started + +In order to use Eigen, you just need to download and extract Eigen's source code. It is not necessary to use CMake or install anything. + +Here are some quick compilation instructions with GCC. To quickly test an example program, just do + +\code g++ -I /path/to/eigen2/ my_program.cpp -o my_program \endcode + +There is no library to link to. For good performance, add the \c -O2 compile-flag. Note however that this makes it impossible to debug inside Eigen code, as many functions get inlined. In some cases, performance can be further improved by disabling Eigen assertions: use \c -DEIGEN_NO_DEBUG or \c -DNDEBUG to disable them. + +On the x86 architecture, the SSE2 instruction set is not enabled by default. Use \c -msse2 to enable it, and Eigen will then automatically enable its vectorized paths. On x86-64 and AltiVec-based architectures, vectorization is enabled by default. + + +\warning \redstar In most cases it is enough to include the \c Eigen/Core header only to get started with Eigen. However, some features presented in this tutorial require the Array module to be included (\c \#include \c ). Those features are highlighted with a red star \redstar. + +\section TutorialCoreSimpleExampleFixedSize Simple example with fixed-size matrices and vectors + +By fixed-size, we mean that the number of rows and columns are fixed at compile-time. In this case, Eigen avoids dynamic memory allocation, and unroll loops when that makes sense. This is useful for very small sizes: typically up to 4x4, sometimes up to 16x16. + + +
+\include Tutorial_simple_example_fixed_size.cpp + +output: +\include Tutorial_simple_example_fixed_size.out +
+ +top\section TutorialCoreSimpleExampleDynamicSize Simple example with dynamic-size matrices and vectors + +By dynamic-size, we mean that the numbers of rows and columns are not fixed at compile-time. In this case, they are stored as runtime variables and the arrays are dynamically allocated. + + +
+\include Tutorial_simple_example_dynamic_size.cpp + +output: +\include Tutorial_simple_example_dynamic_size.out +
+ + + + + + +top\section TutorialCoreMatrixTypes Matrix and vector types + +In Eigen, all kinds of dense matrices and vectors are represented by the template class Matrix. In most cases, you can simply use one of the \ref matrixtypedefs "convenience typedefs". + +The template class Matrix takes a number of template parameters, but for now it is enough to understand the 3 first ones (and the others can then be left unspecified): + +\code Matrix \endcode + +\li \c Scalar is the scalar type, i.e. the type of the coefficients. That is, if you want a vector of floats, choose \c float here. +\li \c RowsAtCompileTime and \c ColsAtCompileTime are the number of rows and columns of the matrix as known at compile-time. + +For example, \c Vector3d is a typedef for \code Matrix \endcode + +For dynamic-size, that is in order to left the number of rows or of columns unspecified at compile-time, use the special value Eigen::Dynamic. For example, \c VectorXd is a typedef for \code Matrix \endcode + + +top\section TutorialCoreCoefficients Coefficient access + +Eigen supports the following syntaxes for read and write coefficient access: + +\code +matrix(i,j); +vector(i) +vector[i] +vector.x() // first coefficient +vector.y() // second coefficient +vector.z() // third coefficient +vector.w() // fourth coefficient +\endcode + +Notice that these coefficient access methods have assertions checking the ranges. So if you do a lot of coefficient access, these assertion can have an important cost. There are then two possibilities if you want avoid paying this cost: +\li Either you can disable assertions altogether, by defining EIGEN_NO_DEBUG or NDEBUG. Notice that some IDEs like MS Visual Studio define NDEBUG automatically in "Release Mode". +\li Or you can disable the checks on a case-by-case basis by using the coeff() and coeffRef() methods: see MatrixBase::coeff(int,int) const, MatrixBase::coeffRef(int,int), etc. + +top\section TutorialCoreMatrixInitialization Matrix and vector creation and initialization + + +\subsection TutorialPredefMat Predefined Matrices +Eigen offers several static methods to create special matrix expressions, and non-static methods to assign these expressions to existing matrices: + + + + + + + + + + + + + + +
Fixed-size matrix or vectorDynamic-size matrixDynamic-size vector
+\code +Matrix3f x; + +x = Matrix3f::Zero(); +x = Matrix3f::Ones(); +x = Matrix3f::Constant(value); +x = Matrix3f::Identity(); +x = Matrix3f::Random(); + +x.setZero(); +x.setOnes(); +x.setIdentity(); +x.setConstant(value); +x.setRandom(); +\endcode + +\code +MatrixXf x; + +x = MatrixXf::Zero(rows, cols); +x = MatrixXf::Ones(rows, cols); +x = MatrixXf::Constant(rows, cols, value); +x = MatrixXf::Identity(rows, cols); +x = MatrixXf::Random(rows, cols); + +x.setZero(rows, cols); +x.setOnes(rows, cols); +x.setConstant(rows, cols, value); +x.setIdentity(rows, cols); +x.setRandom(rows, cols); +\endcode + +\code +VectorXf x; + +x = VectorXf::Zero(size); +x = VectorXf::Ones(size); +x = VectorXf::Constant(size, value); +x = VectorXf::Identity(size); +x = VectorXf::Random(size); + +x.setZero(size); +x.setOnes(size); +x.setConstant(size, value); +x.setIdentity(size); +x.setRandom(size); +\endcode +
\redstar the Random() and setRandom() functions require the inclusion of the Array module (\c \#include \c )
Basis vectors \link MatrixBase::Unit [details]\endlink
\code +Vector3f::UnitX() // 1 0 0 +Vector3f::UnitY() // 0 1 0 +Vector3f::UnitZ() // 0 0 1 +\endcode\code +VectorXf::Unit(size,i) +VectorXf::Unit(4,1) == Vector4f(0,1,0,0) + == Vector4f::UnitY() +\endcode +
+ +Here is an usage example: + +
+\code +cout << MatrixXf::Constant(2, 3, sqrt(2)) << endl; +RowVector3i v; +v.setConstant(6); +cout << "v = " << v << endl; +\endcode + +output: +\code +1.41 1.41 1.41 +1.41 1.41 1.41 + +v = 6 6 6 +\endcode +
+ + +\subsection TutorialCasting Casting + +In Eigen, any matrices of same size and same scalar type are all naturally compatible. The scalar type can be explicitely casted to another one using the template MatrixBase::cast() function: +\code +Matrix3d md(1,2,3); +Matrix3f mf = md.cast(); +\endcode +Note that casting to the same scalar type in an expression is free. + +The destination matrix is automatically resized in any assignment: +\code +MatrixXf res(10,10); +Matrix3f a, b; +res = a+b; // OK: res is resized to size 3x3 +\endcode +Of course, fixed-size matrices can't be resized. + + +\subsection TutorialMap Map +Any memory buffer can be mapped as an Eigen expression: +
+\code +std::vector stlarray(10); +Map(&stlarray[0], stlarray.size()).setOnes(); +int data[4] = 1, 2, 3, 4; +Matrix2i mat2x2(data); +MatrixXi mat2x2 = Map(data); +MatrixXi mat2x2 = Map(data,2,2); +\endcode +
+ + + +\subsection TutorialCommaInit Comma initializer +Eigen also offers a \ref MatrixBaseCommaInitRef "comma initializer syntax" which allows you to set all the coefficients of a matrix to specific values: + +
+\include Tutorial_commainit_01.cpp + +output: +\verbinclude Tutorial_commainit_01.out +
+ +Not excited by the above example? Then look at the following one where the matrix is set by blocks: + +
+\include Tutorial_commainit_02.cpp + +output: +\verbinclude Tutorial_commainit_02.out +
+ +\b Side \b note: here \link CommaInitializer::finished() .finished() \endlink +is used to get the actual matrix object once the comma initialization +of our temporary submatrix is done. Note that despite the apparent complexity of such an expression, +Eigen's comma initializer usually compiles to very optimized code without any overhead. + + + + + + +top\section TutorialCoreArithmeticOperators Arithmetic Operators + +In short, all arithmetic operators can be used right away as in the following example. Note however that arithmetic operators are only given their usual meaning from mathematics tradition. For other operations, such as taking the coefficient-wise product of two vectors, see the discussion of \link Cwise .cwise() \endlink below. Anyway, here is an example demonstrating basic arithmetic operators: +\code +mat4 -= mat1*1.5 + mat2 * (mat3/4); +\endcode +which includes two matrix scalar products ("mat1*1.5" and "mat3/4"), a matrix-matrix product ("mat2 * (mat3/4)"), +a matrix addition ("+") and substraction with assignment ("-="). + + + + + +
+matrix/vector product\code +col2 = mat1 * col1; +row2 = row1 * mat1; row1 *= mat1; +mat3 = mat1 * mat2; mat3 *= mat1; \endcode +
+add/subtract\code +mat3 = mat1 + mat2; mat3 += mat1; +mat3 = mat1 - mat2; mat3 -= mat1;\endcode +
+scalar product\code +mat3 = mat1 * s1; mat3 = s1 * mat1; mat3 *= s1; +mat3 = mat1 / s1; mat3 /= s1;\endcode +
+ +In Eigen, only traditional mathematical operators can be used right away. +But don't worry, thanks to the \link Cwise .cwise() \endlink operator prefix, +Eigen's matrices are also very powerful as a numerical container supporting +most common coefficient-wise operators. + + +
+ + + + + + +
Coefficient wise \link Cwise::operator*() product \endlink\code mat3 = mat1.cwise() * mat2; \endcode +
+Add a scalar to all coefficients \redstar\code +mat3 = mat1.cwise() + scalar; +mat3.cwise() += scalar; +mat3.cwise() -= scalar; +\endcode +
+Coefficient wise \link Cwise::operator/() division \endlink \redstar\code +mat3 = mat1.cwise() / mat2; \endcode +
+Coefficient wise \link Cwise::inverse() reciprocal \endlink \redstar\code +mat3 = mat1.cwise().inverse(); \endcode +
+Coefficient wise comparisons \redstar \n +(support all operators)\code +mat3 = mat1.cwise() < mat2; +mat3 = mat1.cwise() <= mat2; +mat3 = mat1.cwise() > mat2; +etc. +\endcode +
+		 + + + +
+\b Trigo \redstar: \n +\link Cwise::sin sin \endlink, \link Cwise::cos cos \endlink\code +mat3 = mat1.cwise().sin(); +etc. +\endcode +
+\b Power \redstar: \n \link Cwise::pow() pow \endlink, +\link Cwise::square square \endlink, +\link Cwise::cube cube \endlink, \n +\link Cwise::sqrt sqrt \endlink, +\link Cwise::exp exp \endlink, +\link Cwise::log log \endlink \code +mat3 = mat1.cwise().square(); +mat3 = mat1.cwise().pow(5); +mat3 = mat1.cwise().log(); +etc. +\endcode +
+\link Cwise::min min \endlink, \link Cwise::max max \endlink, \n +absolute value (\link Cwise::abs() abs \endlink, \link Cwise::abs2() abs2 \endlink) +\code +mat3 = mat1.cwise().min(mat2); +mat3 = mat1.cwise().max(mat2); +mat3 = mat1.cwise().abs(); +mat3 = mat1.cwise().abs2(); +\endcode
+
+\redstar Those functions require the inclusion of the Array module (\c \#include \c ). + +\b Side \b note: If you think that the \c .cwise() syntax is too verbose for your own taste and prefer to have non-conventional mathematical operators directly available, then feel free to extend MatrixBase as described \ref ExtendingMatrixBase "here". + +So far, we saw the notation \code mat1*mat2 \endcode for matrix product, and \code mat1.cwise()*mat2 \endcode for coefficient-wise product. What about other kinds of products, which in some other libraries also use arithmetic operators? In Eigen, they are accessed as follows -- note that here we are anticipating on further sections, for convenience. + + + + +
\link MatrixBase::dot() dot product \endlink (inner product)\code +scalar = vec1.dot(vec2);\endcode +
+outer product\code +mat = vec1 * vec2.transpose();\endcode +
+\link MatrixBase::cross() cross product \endlink\code +#include +vec3 = vec1.cross(vec2);\endcode
+ + + + +top\section TutorialCoreReductions Reductions + +Eigen provides several reduction methods such as: +\link MatrixBase::minCoeff() minCoeff() \endlink, \link MatrixBase::maxCoeff() maxCoeff() \endlink, +\link MatrixBase::sum() sum() \endlink, \link MatrixBase::trace() trace() \endlink, +\link MatrixBase::norm() norm() \endlink, \link MatrixBase::squaredNorm() squaredNorm() \endlink, +\link MatrixBase::all() all() \endlink \redstar,and \link MatrixBase::any() any() \endlink \redstar. +All reduction operations can be done matrix-wise, +\link MatrixBase::colwise() column-wise \endlink \redstar or +\link MatrixBase::rowwise() row-wise \endlink \redstar. Usage example: + + + + +
\code + 5 3 1 +mat = 2 7 8 + 9 4 6 \endcode + \code mat.minCoeff(); \endcode\code 1 \endcode
\code mat.colwise().minCoeff(); \endcode\code 2 3 1 \endcode
\code mat.rowwise().minCoeff(); \endcode\code +1 +2 +4 +\endcode
+ +Also note that maxCoeff and minCoeff can takes optional arguments returning the coordinates of the respective min/max coeff: \link MatrixBase::maxCoeff(int*,int*) const maxCoeff(int* i, int* j) \endlink, \link MatrixBase::minCoeff(int*,int*) const minCoeff(int* i, int* j) \endlink. + +\b Side \b note: The all() and any() functions are especially useful in combinaison with coeff-wise comparison operators (\ref CwiseAll "example"). + + + + + +top\section TutorialCoreMatrixBlocks Matrix blocks + +Read-write access to a \link MatrixBase::col(int) column \endlink +or a \link MatrixBase::row(int) row \endlink of a matrix: +\code +mat1.row(i) = mat2.col(j); +mat1.col(j1).swap(mat1.col(j2)); +\endcode + +Read-write access to sub-vectors: + + + + + + + + + + + + + + + + + + + + +
Default versionsOptimized versions when the size \n is known at compile time
\code vec1.start(n)\endcode\code vec1.start()\endcodethe first \c n coeffs
\code vec1.end(n)\endcode\code vec1.end()\endcodethe last \c n coeffs
\code vec1.segment(pos,n)\endcode\code vec1.segment(pos)\endcodethe \c size coeffs in \n the range [\c pos : \c pos + \c n [
+ +Read-write access to sub-matrices:
\code mat1.block(i,j,rows,cols)\endcode + \link MatrixBase::block(int,int,int,int) (more) \endlink\code mat1.block(i,j)\endcode + \link MatrixBase::block(int,int) (more) \endlinkthe \c rows x \c cols sub-matrix \n starting from position (\c i,\c j)
\code + mat1.corner(TopLeft,rows,cols) + mat1.corner(TopRight,rows,cols) + mat1.corner(BottomLeft,rows,cols) + mat1.corner(BottomRight,rows,cols)\endcode + \link MatrixBase::corner(CornerType,int,int) (more) \endlink\code + mat1.corner(TopLeft) + mat1.corner(TopRight) + mat1.corner(BottomLeft) + mat1.corner(BottomRight)\endcode + \link MatrixBase::corner(CornerType) (more) \endlinkthe \c rows x \c cols sub-matrix \n taken in one of the four corners
\code +mat4x4.minor(i,j) = mat3x3; +mat3x3 = mat4x4.minor(i,j);\endcode + +\link MatrixBase::minor() minor \endlink (read-write)
+ + + +top\section TutorialCoreDiagonalMatrices Diagonal matrices + + + + + + +
+\link MatrixBase::asDiagonal() make a diagonal matrix \endlink from a vector \n +this product is automatically optimized !\code +mat3 = mat1 * vec2.asDiagonal();\endcode +
Access \link MatrixBase::diagonal() the diagonal of a matrix \endlink as a vector (read/write)\code + vec1 = mat1.diagonal(); + mat1.diagonal() = vec1; + \endcode +
+ +top\section TutorialCoreTransposeAdjoint Transpose and Adjoint operations + + + + +
+\link MatrixBase::transpose() transposition \endlink (read-write)\code +mat3 = mat1.transpose() * mat2; +mat3.transpose() = mat1 * mat2.transpose(); +\endcode +
+\link MatrixBase::adjoint() adjoint \endlink (read only)\n\code +mat3 = mat1.adjoint() * mat2; +\endcode +
+ +top\section TutorialCoreDotNorm Dot-product, vector norm, normalization + + + + + +
+\link MatrixBase::dot() Dot-product \endlink of two vectors +\code vec1.dot(vec2);\endcode +
+\link MatrixBase::norm() norm \endlink of a vector \n +\link MatrixBase::squaredNorm() squared norm \endlink of a vector +\code vec.norm(); \endcode \n \code vec.squaredNorm() \endcode +
+returns a \link MatrixBase::normalized() normalized \endlink vector \n +\link MatrixBase::normalize() normalize \endlink a vector +\code +vec3 = vec1.normalized(); +vec1.normalize();\endcode +
+ +top\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.. + + + + + + + +
+Extract triangular matrices \n from a given matrix m: +\code +m.part() +m.part() +m.part() +m.part() +m.part() +m.part()\endcode +
+Write to triangular parts \n of a matrix m: +\code +m1.part() = m2; +m1.part() = m2; +m1.part() = m2; +m1.part() = m2;\endcode +
+Special: take advantage of symmetry \n (selfadjointness) when copying \n an expression into a matrix +\code +m.part() = someSelfadjointMatrix; +m1.part() = m2 + m2.adjoint(); // m2 + m2.adjoint() is selfadjoint \endcode +
+ + +top\section TutorialCoreSpecialTopics Special Topics + +\ref TopicLazyEvaluation "Lazy Evaluation and Aliasing": Thanks to expression templates, Eigen is able to apply lazy evaluation wherever that is beneficial. + +*/ + +} diff --git a/doc/C03_TutorialGeometry.dox b/doc/C03_TutorialGeometry.dox new file mode 100644 index 000000000..b3c860d53 --- /dev/null +++ b/doc/C03_TutorialGeometry.dox @@ -0,0 +1,253 @@ +namespace Eigen { + +/** \page TutorialGeometry Tutorial 2/3 - Geometry + \ingroup Tutorial + +
\ref index "Overview" + | \ref TutorialCore "Core features" + | \b Geometry + | \ref TutorialAdvancedLinearAlgebra "Advanced linear algebra" + | \ref TutorialSparse "Sparse matrix" +
+ +In this tutorial chapter we will shortly introduce the many possibilities offered by the \ref GeometryModule "geometry module", +namely 2D and 3D rotations and affine transformations. + +\b Table \b of \b contents + - \ref TutorialGeoElementaryTransformations + - \ref TutorialGeoCommontransformationAPI + - \ref TutorialGeoTransform + - \ref TutorialGeoEulerAngles + +Eigen's Geometry module provides two different kinds of geometric transformations: + - Abstract transformations, such as rotations (represented by \ref AngleAxis "angle and axis" or by a \ref Quaternion "quaternion"), \ref Translation "translations", \ref Scaling "scalings". These transformations are NOT represented as matrices, but you can nevertheless mix them with matrices and vectors in expressions, and convert them to matrices if you wish. + - Affine transformation matrices: see the Transform class. These are really matrices. + +\note If you are working with OpenGL 4x4 matrices then Transform3f and Transform3d are what you want. Since Eigen defaults to column-major storage, you can directly use the Transform::data() method to pass your transformation matrix to OpenGL. + +You can construct a Transform from an abstract transformation, like this: +\code + Transform t(AngleAxis(angle,axis)); +\endcode +or like this: +\code + Transform t; + t = AngleAxis(angle,axis); +\endcode +But note that unfortunately, because of how C++ works, you can \b not do this: +\code + Transform t = AngleAxis(angle,axis); +\endcode +\b Explanation: In the C++ language, this would require Transform to have a non-explicit conversion constructor from AngleAxis, but we really don't want to allow implicit casting here. + + +\section TutorialGeoElementaryTransformations Transformation types + + + + + + + + + + +
Transformation typeTypical initialization code
+\ref Rotation2D "2D rotation" from an angle\code +Rotation2D rot2(angle_in_radian);\endcode
+3D rotation as an \ref AngleAxis "angle + axis"\code +AngleAxis aa(angle_in_radian, Vector3f(ax,ay,az));\endcode
+3D rotation as a \ref Quaternion "quaternion"\code +Quaternion q = AngleAxis(angle_in_radian, axis);\endcode
+N-D Scaling\code +Scaling(sx, sy) +Scaling(sx, sy, sz) +Scaling(s) +Scaling(vecN)\endcode
+N-D Translation\code +Translation(tx, ty) +Translation(tx, ty, tz) +Translation(s) +Translation(vecN)\endcode
+N-D \ref TutorialGeoTransform "Affine transformation"\code +Transform t = concatenation_of_any_transformations; +Transform t = Translation3f(p) * AngleAxisf(a,axis) * Scaling3f(s);\endcode
+N-D Linear transformations \n +(pure rotations, \n scaling, etc.)\code +Matrix t = concatenation_of_rotations_and_scalings; +Matrix t = Rotation2Df(a) * Scaling2f(s); +Matrix t = AngleAxisf(a,axis) * Scaling3f(s);\endcode
+ +Notes on rotations\n To transform more than a single vector the preferred +representations are rotation matrices, while for other usages Quaternion is the +representation of choice as they are compact, fast and stable. Finally Rotation2D and +AngleAxis are mainly convenient types to create other rotation objects. + +Notes on Translation and Scaling\n Likewise AngleAxis, these classes were +designed to simplify the creation/initialization of linear (Matrix) and affine (Transform) +transformations. Nevertheless, unlike AngleAxis which is inefficient to use, these classes +might still be interesting to write generic and efficient algorithms taking as input any +kind of transformations. + +Any of the above transformation types can be converted to any other types of the same nature, +or to a more generic type. Here are some additional examples: + + +
\code +Rotation2Df r = Matrix2f(..); // assumes a pure rotation matrix +AngleAxisf aa = Quaternionf(..); +AngleAxisf aa = Matrix3f(..); // assumes a pure rotation matrix +Matrix2f m = Rotation2Df(..); +Matrix3f m = Quaternionf(..); Matrix3f m = Scaling3f(..); +Transform3f m = AngleAxis3f(..); Transform3f m = Scaling3f(..); +Transform3f m = Translation3f(..); Transform3f m = Matrix3f(..); +\endcode
+ + +top\section TutorialGeoCommontransformationAPI Common API across transformation types + +To some extent, Eigen's \ref GeometryModule "geometry module" allows you to write +generic algorithms working on any kind of transformation representations: + + + + + +
+Concatenation of two transformations\code +gen1 * gen2;\endcode
Apply the transformation to a vector\code +vec2 = gen1 * vec1;\endcode
Get the inverse of the transformation\code +gen2 = gen1.inverse();\endcode
Spherical interpolation \n (Rotation2D and Quaternion only)\code +rot3 = rot1.slerp(alpha,rot2);\endcode
+ + + +top\section TutorialGeoTransform Affine transformations +Generic affine transformations are represented by the Transform class which internaly +is a (Dim+1)^2 matrix. In Eigen we have chosen to not distinghish between points and +vectors such that all points are actually represented by displacement vectors from the +origin ( \f$ \mathbf{p} \equiv \mathbf{p}-0 \f$ ). With that in mind, real points and +vector distinguish when the transformation is applied. + + + + + + + +
+Apply the transformation to a \b point \code +VectorNf p1, p2; +p2 = t * p1;\endcode
+Apply the transformation to a \b vector \code +VectorNf vec1, vec2; +vec2 = t.linear() * vec1;\endcode
+Apply a \em general transformation \n to a \b normal \b vector +(explanations)\code +VectorNf n1, n2; +MatrixNf normalMatrix = t.linear().inverse().transpose(); +n2 = (normalMatrix * n1).normalized();\endcode
+Apply a transformation with \em pure \em rotation \n to a \b normal \b vector +(no scaling, no shear)\code +n2 = t.linear() * n1;\endcode
+OpenGL compatibility \b 3D \code +glLoadMatrixf(t.data());\endcode
+OpenGL compatibility \b 2D \code +Transform3f aux(Transform3f::Identity); +aux.linear().corner<2,2>(TopLeft) = t.linear(); +aux.translation().start<2>() = t.translation(); +glLoadMatrixf(aux.data());\endcode
+ +\b Component \b accessors + + + + + + +
+full read-write access to the internal matrix\code +t.matrix() = matN1xN1; // N1 means N+1 +matN1xN1 = t.matrix(); +\endcode
+coefficient accessors\code +t(i,j) = scalar; <=> t.matrix()(i,j) = scalar; +scalar = t(i,j); <=> scalar = t.matrix()(i,j); +\endcode
+translation part\code +t.translation() = vecN; +vecN = t.translation(); +\endcode
+linear part\code +t.linear() = matNxN; +matNxN = t.linear(); +\endcode
+extract the rotation matrix\code +matNxN = t.extractRotation(); +\endcode
+ + +\b Transformation \b creation \n +While transformation objects can be created and updated concatenating elementary transformations, +the Transform class also features a procedural API: + + + + + + +
\b procedurale \b API \b equivalent \b natural \b API
Translation\code +t.translate(Vector_(tx,ty,..)); +t.pretranslate(Vector_(tx,ty,..)); +\endcode\code +t *= Translation_(tx,ty,..); +t = Translation_(tx,ty,..) * t; +\endcode
\b Rotation \n In 2D and for the procedural API, any_rotation can also \n be an angle in radian\code +t.rotate(any_rotation); +t.prerotate(any_rotation); +\endcode\code +t *= any_rotation; +t = any_rotation * t; +\endcode
Scaling\code +t.scale(Vector_(sx,sy,..)); +t.scale(s); +t.prescale(Vector_(sx,sy,..)); +t.prescale(s); +\endcode\code +t *= Scaling_(sx,sy,..); +t *= Scaling_(s); +t = Scaling_(sx,sy,..) * t; +t = Scaling_(s) * t; +\endcode
Shear transformation \n ( \b 2D \b only ! )\code +t.shear(sx,sy); +t.preshear(sx,sy); +\endcode
+ +Note that in both API, any many transformations can be concatenated in a single expression as shown in the two following equivalent examples: + + + +
\code +t.pretranslate(..).rotate(..).translate(..).scale(..); +\endcode
\code +t = Translation_(..) * t * RotationType(..) * Translation_(..) * Scaling_(..); +\endcode
+ + + +top\section TutorialGeoEulerAngles Euler angles + + +
+Euler angles might be convenient to create rotation objects. +On the other hand, since there exist 24 differents convension,they are pretty confusing to use. This example shows how +to create a rotation matrix according to the 2-1-2 convention.\code +Matrix3f m; +m = AngleAxisf(angle1, Vector3f::UnitZ()) +* * AngleAxisf(angle2, Vector3f::UnitY()) +* * AngleAxisf(angle3, Vector3f::UnitZ()); +\endcode
+ +*/ + +} diff --git a/doc/C05_TutorialLinearAlgebra.dox b/doc/C05_TutorialLinearAlgebra.dox new file mode 100644 index 000000000..1b584e103 --- /dev/null +++ b/doc/C05_TutorialLinearAlgebra.dox @@ -0,0 +1,188 @@ + +namespace Eigen { + +/** \page TutorialAdvancedLinearAlgebra Tutorial 3/3 - Advanced linear algebra + \ingroup Tutorial + +
\ref index "Overview" + | \ref TutorialCore "Core features" + | \ref TutorialGeometry "Geometry" + | \b Advanced \b linear \b algebra + | \ref TutorialSparse "Sparse matrix" +
+ +\b Table \b of \b contents + - \ref TutorialAdvSolvers + - \ref TutorialAdvLU + - \ref TutorialAdvCholesky + - \ref TutorialAdvQR + - \ref TutorialAdvEigenProblems + +\section TutorialAdvSolvers Solving linear problems + +This part of the tutorial focuses on solving linear problem of the form \f$ A \mathbf{x} = b \f$, +where both \f$ A \f$ and \f$ b \f$ are known, and \f$ x \f$ is the unknown. Moreover, \f$ A \f$ +assumed to be a squared matrix. Of course, the methods described here can be used whenever an expression +involve the product of an inverse matrix with a vector or another matrix: \f$ A^{-1} B \f$. +Eigen offers various algorithms to this problem, and its choice mainly depends on the nature of +the matrix \f$ A \f$, such as its shape, size and numerical properties. + +\subsection TutorialAdvSolvers_Triangular Triangular solver +If the matrix \f$ A \f$ is triangular (upper or lower) and invertible (the coefficients of the diagonal +are all not zero), then the problem can be solved directly using MatrixBase::solveTriangular(), or better, +MatrixBase::solveTriangularInPlace(). Here is an example: + +
+\include MatrixBase_marked.cpp + +output: +\include MatrixBase_marked.out +
+ +See MatrixBase::solveTriangular() for more details. + + +\subsection TutorialAdvSolvers_Inverse Direct inversion (for small matrices) +If the matrix \f$ A \f$ is small (\f$ \leq 4 \f$) and invertible, then a good approach is to directly compute +the inverse of the matrix \f$ A \f$, and then obtain the solution \f$ x \f$ by \f$ \mathbf{x} = A^{-1} b \f$. With Eigen, +this can be implemented like this: + +\code +#include +Matrix4f A = Matrix4f::Random(); +Vector4f b = Vector4f::Random(); +Vector4f x = A.inverse() * b; +\endcode + +Note that the function inverse() is defined in the LU module. +See MatrixBase::inverse() for more details. + + +\subsection TutorialAdvSolvers_Symmetric Cholesky (for positive definite matrices) +If the matrix \f$ A \f$ is \b positive \b definite, then +the best method is to use a Cholesky decomposition. +Such positive definite matrices often arise when solving overdetermined problems in a least square sense (see below). +Eigen offers two different Cholesky decompositions: a \f$ LL^T \f$ decomposition where L is a lower triangular matrix, +and a \f$ LDL^T \f$ decomposition where L is lower triangular with unit diagonal and D is a diagonal matrix. +The latter avoids square roots and is therefore slightly more stable than the former one. +\code +#include +MatrixXf D = MatrixXf::Random(8,4); +MatrixXf A = D.transpose() * D; +VectorXf b = D.transpose() * VectorXf::Random(4); +VectorXf x; +A.llt().solve(b,&x); // using a LLT factorization +A.ldlt().solve(b,&x); // using a LDLT factorization +\endcode +when the value of the right hand side \f$ b \f$ is not needed anymore, then it is faster to use +the \em in \em place API, e.g.: +\code +A.llt().solveInPlace(b); +\endcode +In that case the value of \f$ b \f$ is replaced by the result \f$ x \f$. + +If the linear problem has to solved for various right hand side \f$ b_i \f$, but same matrix \f$ A \f$, +then it is highly recommended to perform the decomposition of \$ A \$ only once, e.g.: +\code +// ... +LLT lltOfA(A); +lltOfA.solveInPlace(b0); +lltOfA.solveInPlace(b1); +// ... +\endcode + +\sa Cholesky_Module, LLT::solve(), LLT::solveInPlace(), LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT. + + +\subsection TutorialAdvSolvers_LU LU decomposition (for most cases) +If the matrix \f$ A \f$ does not fit in any of the previous categories, or if you are unsure about the numerical +stability of your problem, then you can use the LU solver based on a decomposition of the same name : see the section \ref TutorialAdvLU below. Actually, Eigen's LU module does not implement a standard LU decomposition, but rather a so-called LU decomposition +with full pivoting and rank update which has much better numerical stability. +The API of the LU solver is the same than the Cholesky one, except that there is no \em in \em place variant: +\code +#include +MatrixXf A = MatrixXf::Random(20,20); +VectorXf b = VectorXf::Random(20); +VectorXf x; +A.lu().solve(b, &x); +\endcode + +Again, the LU decomposition can be stored to be reused or to perform other kernel operations: +\code +// ... +LU luOfA(A); +luOfA.solve(b, &x); +// ... +\endcode + +See the section \ref TutorialAdvLU below. + +\sa class LU, LU::solve(), LU_Module + + +\subsection TutorialAdvSolvers_SVD SVD solver (for singular matrices and special cases) +Finally, Eigen also offer a solver based on a singular value decomposition (SVD). Again, the API is the +same than with Cholesky or LU: +\code +#include +MatrixXf A = MatrixXf::Random(20,20); +VectorXf b = VectorXf::Random(20); +VectorXf x; +A.svd().solve(b, &x); +SVD svdOfA(A); +svdOfA.solve(b, &x); +\endcode + +\sa class SVD, SVD::solve(), SVD_Module + + + + +top\section TutorialAdvLU LU + +Eigen provides a rank-revealing LU decomposition with full pivoting, which has very good numerical stability. + +You can obtain the LU decomposition of a matrix by calling \link MatrixBase::lu() lu() \endlink, which is the easiest way if you're going to use the LU decomposition only once, as in +\code +#include +MatrixXf A = MatrixXf::Random(20,20); +VectorXf b = VectorXf::Random(20); +VectorXf x; +A.lu().solve(b, &x); +\endcode + +Alternatively, you can construct a named LU decomposition, which allows you to reuse it for more than one operation: +\code +#include +MatrixXf A = MatrixXf::Random(20,20); +Eigen::LUDecomposition lu(A); +cout << "The rank of A is" << lu.rank() << endl; +if(lu.isInvertible()) { + cout << "A is invertible, its inverse is:" << endl << lu.inverse() << endl; +} +else { + cout << "Here's a matrix whose columns form a basis of the kernel a.k.a. nullspace of A:" + << endl << lu.kernel() << endl; +} +\endcode + +\sa LU_Module, LU::solve(), class LU + +top\section TutorialAdvCholesky Cholesky +todo + +\sa Cholesky_Module, LLT::solve(), LLT::solveInPlace(), LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT + +top\section TutorialAdvQR QR +todo + +\sa QR_Module, class QR + +top\section TutorialAdvEigenProblems Eigen value problems +todo + +\sa class SelfAdjointEigenSolver, class EigenSolver + +*/ + +} diff --git a/doc/C07_TutorialSparse.dox b/doc/C07_TutorialSparse.dox new file mode 100644 index 000000000..a8bfe006e --- /dev/null +++ b/doc/C07_TutorialSparse.dox @@ -0,0 +1,240 @@ +namespace Eigen { + +/** \page TutorialSparse Tutorial - Getting started with the sparse module + \ingroup Tutorial + +
\ref index "Overview" + | \ref TutorialCore "Core features" + | \ref TutorialGeometry "Geometry" + | \ref TutorialAdvancedLinearAlgebra "Advanced linear algebra" + | \b Sparse \b matrix +
+ +\b Table \b of \b contents \n + - \ref TutorialSparseIntro + - \ref TutorialSparseFilling + - \ref TutorialSparseFeatureSet + - \ref TutorialSparseDirectSolvers +
+ +\section TutorialSparseIntro Sparse matrix representations + +In many applications (e.g., finite element methods) it is common to deal with very large matrices where only a few coefficients are different than zero. Both in term of memory consumption and performance, it is fundamental to use an adequate representation storing only nonzero coefficients. Such a matrix is called a sparse matrix. + +\b Declaring \b sparse \b matrices \b and \b vectors \n +The SparseMatrix class is the main sparse matrix representation of the Eigen's sparse module which offers high performance, low memory usage, and compatibility with most of sparse linear algebra packages. Because of its limited flexibility, we also provide a DynamicSparseMatrix variante taillored for low-level sparse matrix assembly. Both of them can be either row major or column major: + +\code +#include +SparseMatrix > m1(1000,2000); // declare a 1000x2000 col-major compressed sparse matrix of complex +SparseMatrix m2(1000,2000); // declare a 1000x2000 row-major compressed sparse matrix of double +DynamicSparseMatrix > m1(1000,2000); // declare a 1000x2000 col-major dynamic sparse matrix of complex +DynamicSparseMatrix m2(1000,2000); // declare a 1000x2000 row-major dynamic sparse matrix of double +\endcode + +Although a sparse matrix could also be used to represent a sparse vector, for that purpose it is better to use the specialized SparseVector class: +\code +SparseVector > v1(1000); // declare a column sparse vector of complex of size 1000 +SparseVector v2(1000); // declare a row sparse vector of double of size 1000 +\endcode +Note that here the size of a vector denotes its dimension and not the number of nonzero coefficients which is initially zero (like sparse matrices). + + +\b Overview \b of \b the \b internal \b sparse \b storage \n +In order to get the best of the Eigen's sparse objects, it is important to have a rough idea of the way they are internally stored. The SparseMatrix class implements the common and generic Compressed Column/Row Storage scheme. It consists of three compact arrays storing the values with their respective inner coordinates, and pointer indices to the begining of each outer vector. For instance, let \c m be a column-major sparse matrix. Then its nonzero coefficients are sequentially stored in memory in a column-major order (\em values). A second array of integer stores the respective row index of each coefficient (\em inner \em indices). Finally, a third array of integer, having the same length than the number of columns, stores the index in the previous arrays of the first element of each column (\em outer \em indices). + +Here is an example, with the matrix: + + + + + + +
03000
2200017
75010
00000
001408
+ +and its internal representation using the Compressed Column Storage format: + + + +
Values: 22735141178
Inner indices: 1202 42 14
+Outer indices:
02456\em 7
+ +As you can guess, here the storage order is even more important than with dense matrix. We will therefore often make a clear difference between the \em inner and \em outer dimensions. For instance, it is easy to loop over the coefficients of an \em inner \em vector (e.g., a column of a column-major matrix), but completely inefficient to do the same for an \em outer \em vector (e.g., a row of a col-major matrix). + +The SparseVector class implements the same compressed storage scheme but, of course, without any outer index buffer. + +Since all nonzero coefficients of such a matrix are sequentially stored in memory, random insertion of new nonzeros can be extremely costly. To overcome this limitation, Eigen's sparse module provides a DynamicSparseMatrix class which is basically implemented as an array of SparseVector. In other words, a DynamicSparseMatrix is a SparseMatrix where the values and inner-indices arrays have been splitted into multiple small and resizable arrays. Assuming the number of nonzeros per inner vector is relatively low, this slight modification allow for very fast random insertion at the cost of a slight memory overhead and a lost of compatibility with other sparse libraries used by some of our highlevel solvers. Note that the major memory overhead comes from the extra memory preallocated by each inner vector to avoid an expensive memory reallocation at every insertion. + +To summarize, it is recommanded to use a SparseMatrix whenever this is possible, and reserve the use of DynamicSparseMatrix for matrix assembly purpose when a SparseMatrix is not flexible enough. The respective pro/cons of both representations are summarized in the following table: + + + + + + + + + + + + + + + +
SparseMatrixDynamicSparseMatrix
memory usage*****
sorted insertion******
random insertion \n in sorted inner vector****
sorted insertion \n in random inner vector-***
random insertion-**
coeff wise unary operators******
coeff wise binary operators******
matrix products*****(*)
transpose*****
redux*****
*= scalar*****
Compatibility with highlevel solvers \n (TAUCS, Cholmod, SuperLU, UmfPack)***-
+ + +\b Matrix \b and \b vector \b properties \n + +Here mat and vec represents any sparse-matrix and sparse-vector types respectively. + + + + + + + + + + +
Standard \n dimensions\code +mat.rows() +mat.cols()\endcode\code +vec.size() \endcode
Sizes along the \n inner/outer dimensions\code +mat.innerSize() +mat.outerSize()\endcode
Number of non \n zero coefficiens\code +mat.nonZeros() \endcode\code +vec.nonZeros() \endcode
+ + +\b Iterating \b over \b the \b nonzero \b coefficients \n + +Iterating over the coefficients of a sparse matrix can be done only in the same order than the storage order. Here is an example: + + +
+\code +SparseMatrixType mat(rows,cols); +for (int k=0; k\ +\code +SparseVector vec(size); +for (SparseVector::InnerIterator it(vec); it; ++it) +{ + it.value(); // == vec[ it.index() ] + it.index(); +} +\endcode +
+ + +\section TutorialSparseFilling Filling a sparse matrix + +A DynamicSparseMatrix object can be set and updated just like any dense matrix using the coeffRef(row,col) method. If the coefficient is not stored yet, then it will be inserted in the matrix. Here is an example: +\code +DynamicSparseMatrix aux(1000,1000); +for (...) + for each i + for each j interacting with i + aux.coeffRef(i,j) += foo(o1,o2); +SparseMatrix mat(aux); // convert the DynamicSparseMatrix to a SparseMatrix +\endcode + +Sometimes, however, we simply want to set all the coefficients of a matrix before using it through standard matrix operations (addition, product, etc.). In that case it faster to use the low-level startFill()/fill()/fillrand()/endFill() interface. Even though this interface is availabe for both sparse matrix types, their respective restrictions slightly differ from one representation to the other. In all case, a call to startFill() set the matrix to zero, and the fill*() functions will fail if the coefficient already exist. + +As a first difference, for SparseMatrix, the fill*() functions can only be called inside a startFill()/endFill() pair, and no other member functions are allowed during the filling process, i.e., until endFill() has been called. On the other hand, a DynamicSparseMatrix is always in a stable state, and the startFill()/endFill() functions are only for compatibility purpose. + +Another difference is that the fill*() functions must be called with increasing outer indices for a SparseMatrix, while they can be random for a DynamicSparseMatrix. + +Finally, the fill() function assumes the coefficient are inserted in a sorted order per inner vector, while the fillrand() variante allows random insertions (the outer indices must still be sorted for SparseMatrix). + +Some examples: + +1 - If you can set the coefficients in exactly the same order that the storage order, then the matrix can be filled directly and very efficiently. Here is an example initializing a random, row-major sparse matrix: +\code +SparseMatrix m(rows,cols); +m.startFill(rows*cols*percent_of_non_zero); // estimate of the number of nonzeros (optional) +for (int i=0; i\ m(rows,cols); +m.startFill(rows*cols*percent_of_non_zero); // estimate of the number of nonzeros (optional) +for (int i=0; i\ m(rows,cols); +{ + RandomSetter > setter(m); + for (int k=0; k\().solveTriangular(dv2); +\endcode + +The product of a sparse matrix A time a dense matrix/vector dv with A symmetric can be optimized by telling that to Eigen: +\code +res = A.marked() * dv; // if all coefficients of A are stored +res = A.marked() * dv; // if only the upper part of A is stored +res = A.marked() * dv; // if only the lower part of A is stored +\endcode + + +\section TutorialSparseDirectSolvers Using the direct solvers + +TODO + +\subsection TutorialSparseDirectSolvers_LLT LLT +Cholmod, Taucs. + +\subsection TutorialSparseDirectSolvers_LDLT LDLT + + +\subsection TutorialSparseDirectSolvers_LU LU +SuperLU, UmfPack. + +*/ + +} diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt index d90fc887d..6934462e8 100644 --- a/doc/CMakeLists.txt +++ b/doc/CMakeLists.txt @@ -38,9 +38,7 @@ add_custom_target( ${CMAKE_CURRENT_BINARY_DIR}/html/ COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/Eigen_Silly_Professor_64x64.png ${CMAKE_CURRENT_BINARY_DIR}/html/ - COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/buildexamplelist.sh ${Eigen_SOURCE_DIR} > ${CMAKE_CURRENT_BINARY_DIR}/ExampleList.dox COMMAND doxygen - COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/cleanhierarchy.sh ${CMAKE_CURRENT_BINARY_DIR}/html/hierarchy.html WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} ) diff --git a/doc/CustomizingEigen.dox b/doc/CustomizingEigen.dox deleted file mode 100644 index 8f2d3aa25..000000000 --- a/doc/CustomizingEigen.dox +++ /dev/null @@ -1,148 +0,0 @@ -namespace Eigen { - -/** \page CustomizingEigen Customizing/Extending Eigen - -Eigen2 can be extended in several ways, for instance, by defining global methods, \ref ExtendingMatrixBase "by adding custom methods to MatrixBase", adding support to \ref CustomScalarType "custom types" etc. - -\b Table \b of \b contents - - \ref ExtendingMatrixBase - - \ref CustomScalarType - - \ref PreprocessorDirectives - -\section ExtendingMatrixBase Extending MatrixBase - -In this section we will see how to add custom methods to MatrixBase. Since all expressions and matrix types inherit MatrixBase, adding a method to MatrixBase make it immediately available to all expressions ! A typical use case is, for instance, to make Eigen compatible with another API. - -You certainly know that in C++ it is not possible to add methods to an extending class. So how that's possible ? Here the trick is to include in the declaration of MatrixBase a file defined by the preprocessor token \c EIGEN_MATRIXBASE_PLUGIN: -\code -class MatrixBase { - // ... - #ifdef EIGEN_MATRIXBASE_PLUGIN - #include EIGEN_MATRIXBASE_PLUGIN - #endif -}; -\endcode -Therefore to extend MatrixBase with you own methods you just have to create a file with your method declaration and define EIGEN_MATRIXBASE_PLUGIN before you include any Eigen's header file. - -Here is an example of such an extension file: \n -\b MatrixBaseAddons.h -\code -inline Scalar at(uint i, uint j) const { return this->operator()(i,j); } -inline Scalar& at(uint i, uint j) { return this->operator()(i,j); } -inline Scalar at(uint i) const { return this->operator[](i); } -inline Scalar& at(uint i) { return this->operator[](i); } - -inline RealScalar squaredLength() const { return squaredNorm(); } -inline RealScalar length() const { return norm(); } -inline RealScalar invLength(void) const { return fast_inv_sqrt(squaredNorm()); } - -template -inline Scalar squaredDistanceTo(const MatrixBase& other) const -{ return (derived() - other.derived()).squaredNorm(); } - -template -inline RealScalar distanceTo(const MatrixBase& other) const -{ return ei_sqrt(derived().squaredDistanceTo(other)); } - -inline void scaleTo(RealScalar l) { RealScalar vl = norm(); if (vl>1e-9) derived() *= (l/vl); } - -inline Transpose transposed() {return transpose();} -inline const Transpose transposed() const {return transpose();} - -inline uint minComponentId(void) const { int i; minCoeff(&i); return i; } -inline uint maxComponentId(void) const { int i; maxCoeff(&i); return i; } - -template -void makeFloor(const MatrixBase& other) { derived() = derived().cwise().min(other.derived()); } -template -void makeCeil(const MatrixBase& other) { derived() = derived().cwise().max(other.derived()); } - -const typename Cwise::ScalarAddReturnType -operator+(const Scalar& scalar) const { return cwise() + scalar } - -friend const typename Cwise::ScalarAddReturnType -operator+(const Scalar& scalar, const MatrixBase& mat) { return mat + scalar; } -\endcode - -Then one can the following declaration in the config.h or whatever prerequisites header file of his project: -\code -#define EIGEN_MATRIXBASE_PLUGIN "MatrixBaseAddons.h" -\endcode - - - -\section CustomScalarType Using custom scalar types - -By default, Eigen currently supports the following scalar types: \c int, \c float, \c double, \c std::complex, \c std::complex, \c long \c double, \c long \c long \c int (64 bits integers), and \c bool. The \c long \c double is especially useful on x86-64 systems or when the SSE2 instruction set is enabled because it enforces the use of x87 registers with extended accuracy. - -In order to add support for a custom type \c T you need: - 1 - make sure the common operator (+,-,*,/,etc.) are supported by the type \c T - 2 - add a specialization of struct Eigen::NumTraits (see \ref NumTraits) - 3 - define a couple of math functions for your type such as: ei_sqrt, ei_abs, etc... - (see the file Eigen/src/Core/MathFunctions.h) - -Here is a concrete example adding support for the Adolc's \c adouble type. Adolc is an automatic differentiation library. The type \c adouble is basically a real value tracking the values of any number of partial derivatives. - -\code -#ifndef ADLOCSUPPORT_H -#define ADLOCSUPPORT_H - -#define ADOLC_TAPELESS -#include -#include - -namespace Eigen { - -template<> struct NumTraits -{ - typedef adtl::adouble Real; - typedef adtl::adouble FloatingPoint; - enum { - IsComplex = 0, - HasFloatingPoint = 1, - ReadCost = 1, - AddCost = 1, - MulCost = 1 - }; -}; - -} - -// the Adolc's type adouble is defined in the adtl namespace -// therefore, the following ei_* functions *must* be defined -// in the same namespace -namespace adtl { - - inline const adouble& ei_conj(const adouble& x) { return x; } - inline const adouble& ei_real(const adouble& x) { return x; } - inline adouble ei_imag(const adouble&) { return 0.; } - inline adouble ei_abs(const adouble& x) { return fabs(x); } - inline adouble ei_abs2(const adouble& x) { return x*x; } - inline adouble ei_sqrt(const adouble& x) { return sqrt(x); } - inline adouble ei_exp(const adouble& x) { return exp(x); } - inline adouble ei_log(const adouble& x) { return log(x); } - inline adouble ei_sin(const adouble& x) { return sin(x); } - inline adouble ei_cos(const adouble& x) { return cos(x); } - inline adouble ei_pow(const adouble& x, adouble y) { return pow(x, y); } - -} - -#endif // ADLOCSUPPORT_H -\endcode - - - -\section PreprocessorDirectives Preprocessor directives - -You can control some aspects of Eigen by defining the following preprocessor tokens them before including any of Eigen's headers. - - \b EIGEN_NO_DEBUG disables Eigen assertions. Like NDEBUG but only affects Eigen's assertions. - - \b EIGEN_DONT_VECTORIZE disables explicit vectorization when defined. - - \b EIGEN_UNROLLING_LIMIT defines the maximal instruction counts to enable meta unrolling of loops. Set it to zero to disable unrolling. The default is 100. - - \b EIGEN_DEFAULT_TO_ROW_MAJOR the default storage order for matrices becomes row-major instead of column-major. - - \b EIGEN_TUNE_FOR_CPU_CACHE_SIZE represents the maximal size in Bytes of L2 blocks. Since several blocks have to stay concurently in L2 cache, this value should correspond to at most 1/4 of the size of L2 cache. - - \b EIGEN_NO_STATIC_ASSERT replaces compile time static assertions by runtime assertions - - \b EIGEN_MATRIXBASE_PLUGIN see \ref ExtendingMatrixBase - -*/ - -} diff --git a/doc/D01_StlContainers.dox b/doc/D01_StlContainers.dox new file mode 100644 index 000000000..9a2586023 --- /dev/null +++ b/doc/D01_StlContainers.dox @@ -0,0 +1,52 @@ +namespace Eigen { + +/** \page StlContainers Troubleshooting - Using STL Containers with Eigen + +\b Table \b of \b contents + - \ref summary + - \ref allocator + - \ref vector + +\section summary Executive summary + +Using STL containers on \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types" requires taking the following two steps: + +\li A 16-byte-aligned allocator must be used. Eigen does provide one ready for use: aligned_allocator. +\li If you want to use the std::vector container, you need to \#include . + +These issues arise only with \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types". For other Eigen types, such as Vector3f or MatrixXd, no special care is needed when using STL containers. + +\section allocator Using an aligned allocator + +STL containers take an optional template parameter, the allocator type. When using STL containers on \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types", you need tell the container to use an allocator that will always allocate memory at 16-byte-aligned locations. Fortunately, Eigen does provide such an allocator: Eigen::aligned_allocator. + +For example, instead of +\code +std::map +\endcode +you need to use +\code +std::map, Eigen::aligned_allocator > +\endcode +Note that here, the 3rd parameter "std::less" is just the default value, we only had to specify it because we needed to specify the allocator type, that is the 4th parameter. + +\section vector The case of std::vector + +The situation with std::vector was even worse (explanation below) so we had to specialize it for Eigen types. The upside is that our specialization takes care of specifying the aligned allocator, so you don't need to worry about it. All you need to do is to \#include . + +So as soon as you have +\code +#include +\endcode +you can simply use +\code +std::vector +\endcode +without having to worry about anything. + +\b Explanation: The resize() method of std::vector takes a value_type argument (defaulting to value_type()). So with std::vector, some Eigen::Vector4f objects will be passed by value, which discards any alignment modifiers, so a Eigen::Vector4f can be created at an unaligned location. In order to avoid that, the only solution we saw was to specialize std::vector to make it work on a slight modification of, here, Eigen::Vector4f, that is able to deal properly with this situation. + + +*/ + +} diff --git a/doc/D07_PassingByValue.dox b/doc/D07_PassingByValue.dox new file mode 100644 index 000000000..d15b375a6 --- /dev/null +++ b/doc/D07_PassingByValue.dox @@ -0,0 +1,40 @@ +namespace Eigen { + +/** \page PassingByValue Troubleshooting - Passing Eigen objects by value to functions + +Passing objects by value is almost always a very bad idea in C++, as this means useless copies, and one should pass them by reference instead. + +With Eigen, this is even more important: passing \ref FixedSizeVectorizable "fixed-size vectorizable Eigen objects" by value is not only inefficient, it can be illegal or make your program crash! And the reason is that these Eigen objects have alignment modifiers that aren't respected when they are passed by value. + +So for example, a function like this, where v is passed by value: + +\code +void my_function(Eigen::Vector2d v); +\endcode + +needs to be rewritten as follows, passing v by reference: + +\code +void my_function(const Eigen::Vector2d& v); +\endcode + +Likewise if you have a class having a Eigen object as member: + +\code +struct Foo +{ + Eigen::Vector2d v; +}; +void my_function(Foo v); +\endcode + +This function also needs to be rewritten like this: +\code +void my_function(const Foo& v); +\endcode + +Note that on the other hand, there is no problem with functions that return objects by value. + +*/ + +} diff --git a/doc/D09_StructHavingEigenMembers.dox b/doc/D09_StructHavingEigenMembers.dox new file mode 100644 index 000000000..4ff2b23aa --- /dev/null +++ b/doc/D09_StructHavingEigenMembers.dox @@ -0,0 +1,149 @@ +namespace Eigen { + +/** \page StructHavingEigenMembers Troubleshooting - Structures Having Eigen Members + +\b Table \b of \b contents + - \ref summary + - \ref what + - \ref how + - \ref why + - \ref movetotop + - \ref bugineigen + - \ref conditional + +\section summary Executive Summary + +If you define a structure having members of \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types", you must overload its "operator new" so that it generates 16-bytes-aligned pointers. Fortunately, Eigen provides you with a macro EIGEN_MAKE_ALIGNED_OPERATOR_NEW that does that for you. + +\section what What kind of code needs to be changed? + +The kind of code that needs to be changed is this: + +\code +class Foo +{ + ... + Eigen::Vector2d v; + ... +}; + +... + +Foo *foo = new Foo; +\endcode + +In other words: you have a class that has as a member a \ref FixedSizeVectorizable "fixed-size vectorizable Eigen object", and then you dynamically create an object of that class. + +\section how How should such code be modified? + +Very easy, you just need to put a EIGEN_MAKE_ALIGNED_OPERATOR_NEW macro in a public part of your class, like this: + +\code +class Foo +{ + ... + Eigen::Vector2d v; + ... +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW +}; + +... + +Foo *foo = new Foo; +\endcode + +This macro makes "new Foo" always return an aligned pointer. + +\section why Why is this needed? + +OK let's say that your code looks like this: + +\code +class Foo +{ + ... + Eigen::Vector2d v; + ... +}; + +... + +Foo *foo = new Foo; +\endcode + +A Eigen::Vector2d consists of 2 doubles, which is 128 bits. Which is exactly the size of a SSE packet, which makes it possible to use SSE for all sorts of operations on this vector. But SSE instructions (at least the ones that Eigen uses, which are the fast ones) require 128-bit alignment. Otherwise you get a segmentation fault. + +For this reason, Eigen takes care by itself to require 128-bit alignment for Eigen::Vector2d, by doing two things: +\li Eigen requires 128-bit alignment for the Eigen::Vector2d's array (of 2 doubles). With GCC, this is done with a __attribute__ ((aligned(16))). +\li Eigen overloads the "operator new" of Eigen::Vector2d so it will always return 128-bit aligned pointers. + +Thus, normally, you don't have to worry about anything, Eigen handles alignment for you... + +... except in one case. When you have a class Foo like above, and you dynamically allocate a new Foo as above, then, since Foo doesn't have aligned "operator new", the returned pointer foo is not necessarily 128-bit aligned. + +The alignment attribute of the member v is then relative to the start of the class, foo. If the foo pointer wasn't aligned, then foo->v won't be aligned either! + +The solution is to let class Foo have an aligned "operator new", as we showed in the previous section. + +\section movetotop Should I then put all the members of Eigen types at the beginning of my class? + +No, that's not needed. Since Eigen takes care of declaring 128-bit alignment, all members that need it are automatically 128-bit aligned relatively to the class. So when you have code like + +\code +class Foo +{ + double x; + Eigen::Vector2d v; +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW +}; +\endcode + +it will work just fine. You do \b not need to rewrite it as + +\code +class Foo +{ + Eigen::Vector2d v; + double x; +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW +}; +\endcode + +\section dynamicsize What about dynamic-size matrices and vectors? + +Dynamic-size matrices and vectors, such as Eigen::VectorXd, allocate dynamically their own array of coefficients, so they take care of requiring absolute alignment automatically. So they don't cause this issue. The issue discussed here is only with fixed-size matrices and vectors. + +\section bugineigen So is this a bug in Eigen? + +No, it's not our bug. It's more like an inherent problem of the C++ language -- though it must be said that any other existing language probably has the same problem. The problem is that there is no way that you can specify an aligned "operator new" that would propagate to classes having you as member data. + +\section conditional What if I want to do this conditionnally (depending on template parameters) ? + +For this situation, we offer the macro EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign). It will generate aligned operators like EIGEN_MAKE_ALIGNED_OPERATOR_NEW if NeedsToAlign is true. It will generate operators with the default alignment if NeedsToAlign is false. + +Example: + +\code +template class Foo +{ + typedef Eigen::Matrix Vector; + enum { NeedsToAlign = (sizeof(Vector)%16)==0 }; + ... + Vector v; + ... +public: + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign) +}; + +... + +Foo<4> *foo4 = new Foo<4>; // foo4 is guaranteed to be 128bit-aligned +Foo<3> *foo3 = new Foo<3>; // foo3 has only the system default alignment guarantee +\endcode + +*/ + +} diff --git a/doc/D11_UnalignedArrayAssert.dox b/doc/D11_UnalignedArrayAssert.dox new file mode 100644 index 000000000..a5ad32de9 --- /dev/null +++ b/doc/D11_UnalignedArrayAssert.dox @@ -0,0 +1,77 @@ +namespace Eigen { + +/** \page UnalignedArrayAssert Troubleshooting - Explanation of the assertion on unaligned arrays + +Hello! You are seeing this webpage because your program terminated on an assertion failure like this one: +
+my_program: path/to/eigen2/Eigen/src/Core/MatrixStorage.h:44:
+Eigen::ei_matrix_array::ei_matrix_array()
+[with T = double, int Size = 2, int MatrixOptions = 2, bool Align = true]:
+Assertion `(reinterpret_cast(array) & 0xf) == 0 && "this assertion
+is explained here: http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html
+**** READ THIS WEB PAGE !!! ****"' failed.
+
+ +There are 3 known causes for this issue. Please read on to understand them and learn how to fix them. + +\b Table \b of \b contents + - \ref c1 + - \ref c2 + - \ref c3 + - \ref explanation + +\section c1 Cause 1: Structures having Eigen objects as members + +If you have code like this, + +\code +class Foo +{ + //... + Eigen::Vector2d v; + //... +}; +//... +Foo *foo = new Foo; +\endcode + +then you need to read this separate page: \ref StructHavingEigenMembers "Structures Having Eigen Members". + +Note that here, Eigen::Vector2d is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types". + +\section c2 Cause 2: STL Containers + +If you use STL Containers such as std::vector, std::map, ..., with Eigen objects, like this, + +\code +std::vector my_vector; +std::map my_map; +\endcode + +then you need to read this separate page: \ref StlContainers "Using STL Containers with Eigen". + +Note that here, Eigen::Matrix2f is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types". + +\section c3 Cause 3: Passing Eigen objects by value + +If some function in your code is getting an Eigen object passed by value, like this, + +\code +void func(Eigen::Vector4d v); +\endcode + +then you need to read this separate page: \ref PassingByValue "Passing Eigen objects by value to functions". + +Note that here, Eigen::Vector4d is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types". + +\section explanation General explanation of this assertion + +\ref FixedSizeVectorizable "fixed-size vectorizable Eigen objects" must absolutely be created at 16-byte-aligned locations, otherwise SIMD instructions adressing them will crash. + +Eigen normally takes care of these alignment issues for you, by setting an alignment attribute on them and by overloading their "operator new". + +However there are a few corner cases where these alignment settings get overridden: they are the possible causes for this assertion. + +*/ + +} diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in index 873dd24bc..04fe712e9 100644 --- a/doc/Doxyfile.in +++ b/doc/Doxyfile.in @@ -564,7 +564,9 @@ WARN_LOGFILE = INPUT = "${Eigen_SOURCE_DIR}/Eigen" \ "${Eigen_SOURCE_DIR}/doc" \ - "${Eigen_BINARY_DIR}/doc" + "${Eigen_BINARY_DIR}/doc" \ + "${Eigen_SOURCE_DIR}/unsupported/Eigen" \ + "${Eigen_SOURCE_DIR}/unsupported/doc" # This tag can be used to specify the character encoding of the source files # that doxygen parses. Internally doxygen uses the UTF-8 encoding, which is @@ -618,7 +620,7 @@ EXCLUDE_PATTERNS = CMake* \ # wildcard * is used, a substring. Examples: ANamespace, AClass, # AClass::ANamespace, ANamespace::*Test -EXCLUDE_SYMBOLS = +EXCLUDE_SYMBOLS = MatrixBase<* MapBase<* RotationBase<* Matrix<* # The EXAMPLE_PATH tag can be used to specify one or more files or # directories that contain example code fragments that are included (see diff --git a/doc/Experimental.dox b/doc/Experimental.dox deleted file mode 100644 index 6d8b90d5a..000000000 --- a/doc/Experimental.dox +++ /dev/null @@ -1,55 +0,0 @@ -namespace Eigen { - -/** \page Experimental Experimental parts of Eigen - -\b Table \b of \b contents - - \ref summary - - \ref modules - - \ref core - -\section summary Summary - -With the 2.0 release, Eigen's API is, to a large extent, stable. However, we wish to retain the freedom to make API incompatible changes. To that effect, we call many parts of Eigen "experimental" which means that they are not subject to API stability guarantee. - -Our goal is that for the 2.1 release (expected in July 2009) most of these parts become API-stable too. - -We are aware that API stability is a major concern for our users. That's why it's a priority for us to reach it, but at the same time we're being serious about not calling Eigen API-stable too early. - -Experimental features may at any time: -\li be removed; -\li be subject to an API incompatible change; -\li introduce API or ABI incompatible changes in your own code if you let them affect your API or ABI. - -\section modules Experimental modules - -The following modules are considered entirely experimental, and we make no firm API stability guarantee about them for the time being: -\li SVD -\li QR -\li Cholesky -\li Sparse -\li Geometry (this one should be mostly stable, but it's a little too early to make a formal guarantee) - -\section core Experimental parts of the Core module - -In the Core module, the only classes subject to ABI stability guarantee (meaning that you can use it for data members in your public ABI) is: -\li Matrix -\li Map - -All other classes offer no ABI guarantee, e.g. the layout of their data can be changed. - -The only classes subject to (even partial) API stability guarantee (meaning that you can safely construct and use objects) are: -\li MatrixBase : partial API stability (see below) -\li Matrix : full API stability (except for experimental stuff inherited from MatrixBase) -\li Map : full API stability (except for experimental stuff inherited from MatrixBase) - -All other classes offer no direct API guarantee, e.g. their methods can be changed; however notice that most classes inherit MatrixBase and that this is where most of their API comes from -- so in practice most of the API is stable. - -A few MatrixBase methods are considered experimental, hence not part of any API stability guarantee: -\li all methods documented as internal -\li all methods hidden in the Doxygen documentation -\li all methods marked as experimental -\li all methods defined in experimental modules - -*/ - -} diff --git a/doc/FixedSizeVectorizable.dox b/doc/FixedSizeVectorizable.dox deleted file mode 100644 index 6793bd921..000000000 --- a/doc/FixedSizeVectorizable.dox +++ /dev/null @@ -1,38 +0,0 @@ -namespace Eigen { - -/** \page FixedSizeVectorizable Fixed-size vectorizable Eigen objects - -The goal of this page is to explain what we mean by "fixed-size vectorizable". - -\section summary Executive Summary - -An Eigen object is called "fixed-size vectorizable" if it has fixed size and that size is a multiple of 16 bytes. - -Examples include: -\li Eigen::Vector2d -\li Eigen::Vector4d -\li Eigen::Vector4f -\li Eigen::Matrix2d -\li Eigen::Matrix2f -\li Eigen::Matrix4d -\li Eigen::Matrix4f -\li Eigen::Transform3d -\li Eigen::Transform3f -\li Eigen::Quaterniond -\li Eigen::Quaternionf - -\section explanation Explanation - -First, "fixed-size" should be clear: an Eigen object has fixed size if its number of rows and its number of columns are fixed at compile-time. So for example Matrix3f has fixed size, but MatrixXf doesn't (the opposite of fixed-size is dynamic-size). - -The array of coefficients of a fixed-size Eigen object is a plain "static array", it is not dynamically allocated. For example, the data behind a Matrix4f is just a "float array[16]". - -Fixed-size objects are typically very small, which means that we want to handle them with zero runtime overhead -- both in terms of memory usage and of speed. - -Now, vectorization (both SSE and AltiVec) works with 128-bit packets. Moreover, for performance reasons, these packets need to be have 128-bit alignment. - -So it turns out that the only way that fixed-size Eigen objects can be vectorized, is if their size is a multiple of 128 bits, or 16 bytes. Eigen will then request 16-byte alignment for these object, and henceforth rely on these objects being aligned so no runtime check for alignment is performed. - -*/ - -} diff --git a/doc/I00_CustomizingEigen.dox b/doc/I00_CustomizingEigen.dox new file mode 100644 index 000000000..e8e830d37 --- /dev/null +++ b/doc/I00_CustomizingEigen.dox @@ -0,0 +1,148 @@ +namespace Eigen { + +/** \page CustomizingEigen Advanced - Customizing/Extending Eigen + +Eigen2 can be extended in several ways, for instance, by defining global methods, \ref ExtendingMatrixBase "by adding custom methods to MatrixBase", adding support to \ref CustomScalarType "custom types" etc. + +\b Table \b of \b contents + - \ref ExtendingMatrixBase + - \ref CustomScalarType + - \ref PreprocessorDirectives + +\section ExtendingMatrixBase Extending MatrixBase + +In this section we will see how to add custom methods to MatrixBase. Since all expressions and matrix types inherit MatrixBase, adding a method to MatrixBase make it immediately available to all expressions ! A typical use case is, for instance, to make Eigen compatible with another API. + +You certainly know that in C++ it is not possible to add methods to an extending class. So how that's possible ? Here the trick is to include in the declaration of MatrixBase a file defined by the preprocessor token \c EIGEN_MATRIXBASE_PLUGIN: +\code +class MatrixBase { + // ... + #ifdef EIGEN_MATRIXBASE_PLUGIN + #include EIGEN_MATRIXBASE_PLUGIN + #endif +}; +\endcode +Therefore to extend MatrixBase with you own methods you just have to create a file with your method declaration and define EIGEN_MATRIXBASE_PLUGIN before you include any Eigen's header file. + +Here is an example of such an extension file: \n +\b MatrixBaseAddons.h +\code +inline Scalar at(uint i, uint j) const { return this->operator()(i,j); } +inline Scalar& at(uint i, uint j) { return this->operator()(i,j); } +inline Scalar at(uint i) const { return this->operator[](i); } +inline Scalar& at(uint i) { return this->operator[](i); } + +inline RealScalar squaredLength() const { return squaredNorm(); } +inline RealScalar length() const { return norm(); } +inline RealScalar invLength(void) const { return fast_inv_sqrt(squaredNorm()); } + +template +inline Scalar squaredDistanceTo(const MatrixBase& other) const +{ return (derived() - other.derived()).squaredNorm(); } + +template +inline RealScalar distanceTo(const MatrixBase& other) const +{ return ei_sqrt(derived().squaredDistanceTo(other)); } + +inline void scaleTo(RealScalar l) { RealScalar vl = norm(); if (vl>1e-9) derived() *= (l/vl); } + +inline Transpose transposed() {return transpose();} +inline const Transpose transposed() const {return transpose();} + +inline uint minComponentId(void) const { int i; minCoeff(&i); return i; } +inline uint maxComponentId(void) const { int i; maxCoeff(&i); return i; } + +template +void makeFloor(const MatrixBase& other) { derived() = derived().cwise().min(other.derived()); } +template +void makeCeil(const MatrixBase& other) { derived() = derived().cwise().max(other.derived()); } + +const typename Cwise::ScalarAddReturnType +operator+(const Scalar& scalar) const { return cwise() + scalar } + +friend const typename Cwise::ScalarAddReturnType +operator+(const Scalar& scalar, const MatrixBase& mat) { return mat + scalar; } +\endcode + +Then one can the following declaration in the config.h or whatever prerequisites header file of his project: +\code +#define EIGEN_MATRIXBASE_PLUGIN "MatrixBaseAddons.h" +\endcode + + + +\section CustomScalarType Using custom scalar types + +By default, Eigen currently supports the following scalar types: \c int, \c float, \c double, \c std::complex, \c std::complex, \c long \c double, \c long \c long \c int (64 bits integers), and \c bool. The \c long \c double is especially useful on x86-64 systems or when the SSE2 instruction set is enabled because it enforces the use of x87 registers with extended accuracy. + +In order to add support for a custom type \c T you need: + 1 - make sure the common operator (+,-,*,/,etc.) are supported by the type \c T + 2 - add a specialization of struct Eigen::NumTraits (see \ref NumTraits) + 3 - define a couple of math functions for your type such as: ei_sqrt, ei_abs, etc... + (see the file Eigen/src/Core/MathFunctions.h) + +Here is a concrete example adding support for the Adolc's \c adouble type. Adolc is an automatic differentiation library. The type \c adouble is basically a real value tracking the values of any number of partial derivatives. + +\code +#ifndef ADLOCSUPPORT_H +#define ADLOCSUPPORT_H + +#define ADOLC_TAPELESS +#include +#include + +namespace Eigen { + +template<> struct NumTraits +{ + typedef adtl::adouble Real; + typedef adtl::adouble FloatingPoint; + enum { + IsComplex = 0, + HasFloatingPoint = 1, + ReadCost = 1, + AddCost = 1, + MulCost = 1 + }; +}; + +} + +// the Adolc's type adouble is defined in the adtl namespace +// therefore, the following ei_* functions *must* be defined +// in the same namespace +namespace adtl { + + inline const adouble& ei_conj(const adouble& x) { return x; } + inline const adouble& ei_real(const adouble& x) { return x; } + inline adouble ei_imag(const adouble&) { return 0.; } + inline adouble ei_abs(const adouble& x) { return fabs(x); } + inline adouble ei_abs2(const adouble& x) { return x*x; } + inline adouble ei_sqrt(const adouble& x) { return sqrt(x); } + inline adouble ei_exp(const adouble& x) { return exp(x); } + inline adouble ei_log(const adouble& x) { return log(x); } + inline adouble ei_sin(const adouble& x) { return sin(x); } + inline adouble ei_cos(const adouble& x) { return cos(x); } + inline adouble ei_pow(const adouble& x, adouble y) { return pow(x, y); } + +} + +#endif // ADLOCSUPPORT_H +\endcode + + + +\section PreprocessorDirectives Preprocessor directives + +You can control some aspects of Eigen by defining the following preprocessor tokens them before including any of Eigen's headers. + - \b EIGEN_NO_DEBUG disables Eigen assertions. Like NDEBUG but only affects Eigen's assertions. + - \b EIGEN_DONT_VECTORIZE disables explicit vectorization when defined. + - \b EIGEN_UNROLLING_LIMIT defines the maximal instruction counts to enable meta unrolling of loops. Set it to zero to disable unrolling. The default is 100. + - \b EIGEN_DEFAULT_TO_ROW_MAJOR the default storage order for matrices becomes row-major instead of column-major. + - \b EIGEN_TUNE_FOR_CPU_CACHE_SIZE represents the maximal size in Bytes of L2 blocks. Since several blocks have to stay concurently in L2 cache, this value should correspond to at most 1/4 of the size of L2 cache. + - \b EIGEN_NO_STATIC_ASSERT replaces compile time static assertions by runtime assertions + - \b EIGEN_MATRIXBASE_PLUGIN see \ref ExtendingMatrixBase + +*/ + +} diff --git a/doc/I01_TopicLazyEvaluation.dox b/doc/I01_TopicLazyEvaluation.dox new file mode 100644 index 000000000..16fe79e70 --- /dev/null +++ b/doc/I01_TopicLazyEvaluation.dox @@ -0,0 +1,67 @@ +namespace Eigen { + +/** \page TopicLazyEvaluation Advanced - Lazy Evaluation and Aliasing + +Executive summary: Eigen has intelligent compile-time mechanisms to enable lazy evaluation and removing temporaries where appropriate. +It will handle aliasing automatically in most cases, for example with matrix products. The automatic behavior can be overridden +manually by using the MatrixBase::eval() and MatrixBase::lazy() methods. + +When you write a line of code involving a complex expression such as + +\code mat1 = mat2 + mat3 * (mat4 + mat5); \endcode + +Eigen determines automatically, for each sub-expression, whether to evaluate it into a temporary variable. Indeed, in certain cases it is better to evaluate immediately a sub-expression into a temporary variable, while in other cases it is better to avoid that. + +A traditional math library without expression templates always evaluates all sub-expressions into temporaries. So with this code, + +\code vec1 = vec2 + vec3; \endcode + +a traditional library would evaluate \c vec2 + vec3 into a temporary \c vec4 and then copy \c vec4 into \c vec1. This is of course inefficient: the arrays are traversed twice, so there are a lot of useless load/store operations. + +Expression-templates-based libraries can avoid evaluating sub-expressions into temporaries, which in many cases results in large speed improvements. This is called lazy evaluation as an expression is getting evaluated as late as possible, instead of immediately. However, most other expression-templates-based libraries always choose lazy evaluation. There are two problems with that: first, lazy evaluation is not always a good choice for performance; second, lazy evaluation can be very dangerous, for example with matrix products: doing matrix = matrix*matrix gives a wrong result if the matrix product is lazy-evaluated, because of the way matrix product works. + +For these reasons, Eigen has intelligent compile-time mechanisms to determine automatically when to use lazy evaluation, and when on the contrary it should evaluate immediately into a temporary variable. + +So in the basic example, + +\code matrix1 = matrix2 + matrix3; \endcode + +Eigen chooses lazy evaluation. Thus the arrays are traversed only once, producing optimized code. If you really want to force immediate evaluation, use \link MatrixBase::eval() eval()\endlink: + +\code matrix1 = (matrix2 + matrix3).eval(); \endcode + +Here is now a more involved example: + +\code matrix1 = -matrix2 + matrix3 + 5 * matrix4; \endcode + +Eigen chooses lazy evaluation at every stage in that example, which is clearly the correct choice. In fact, lazy evaluation is the "default choice" and Eigen will choose it except in a few circumstances. + +The first circumstance in which Eigen chooses immediate evaluation, is when it sees an assignment a = b; and the expression \c b has the evaluate-before-assigning \link flags flag\endlink. The most important example of such an expression is the \link Product matrix product expression\endlink. For example, when you do + +\code matrix = matrix * matrix; \endcode + +Eigen first evaluates matrix * matrix into a temporary matrix, and then copies it into the original \c matrix. This guarantees a correct result as we saw above that lazy evaluation gives wrong results with matrix products. It also doesn't cost much, as the cost of the matrix product itself is much higher. + +What if you know what you are doing and want to force lazy evaluation? Then use \link MatrixBase::lazy() .lazy()\endlink instead. Here is an example: + +\code matrix1 = (matrix2 * matrix2).lazy(); \endcode + +Here, since we know that matrix2 is not the same matrix as matrix1, we know that lazy evaluation is not dangerous, so we may force lazy evaluation. Concretely, the effect of lazy() here is to remove the evaluate-before-assigning \link flags flag\endlink and also the evaluate-before-nesting \link flags flag\endlink which we now discuss. + +The second circumstance in which Eigen chooses immediate evaluation, is when it sees a nested expression such as a + b where \c b is already an expression having the evaluate-before-nesting \link flags flag\endlink. Again, the most important example of such an expression is the \link Product matrix product expression\endlink. For example, when you do + +\code matrix1 = matrix2 + matrix3 * matrix4; \endcode + +the product matrix3 * matrix4 gets evaluated immediately into a temporary matrix. Indeed, experiments showed that it is often beneficial for performance to evaluate immediately matrix products when they are nested into bigger expressions. + +Again, \link MatrixBase::lazy() .lazy()\endlink can be used to force lazy evaluation here. + +The third circumstance in which Eigen chooses immediate evaluation, is when its cost model shows that the total cost of an operation is reduced if a sub-expression gets evaluated into a temporary. Indeed, in certain cases, an intermediate result is sufficiently costly to compute and is reused sufficiently many times, that is worth "caching". Here is an example: + +\code matrix1 = matrix2 * (matrix3 + matrix4); \endcode + +Here, provided the matrices have at least 2 rows and 2 columns, each coefficienct of the expression matrix3 + matrix4 is going to be used several times in the matrix product. Instead of computing the sum everytime, it is much better to compute it once and store it in a temporary variable. Eigen understands this and evaluates matrix3 + matrix4 into a temporary variable before evaluating the product. + +*/ + +} \ No newline at end of file diff --git a/doc/I03_InsideEigenExample.dox b/doc/I03_InsideEigenExample.dox new file mode 100644 index 000000000..d4960e79d --- /dev/null +++ b/doc/I03_InsideEigenExample.dox @@ -0,0 +1,500 @@ +namespace Eigen { + +/** \page InsideEigenExample Advanced - What happens inside Eigen, on a simple example + +\b Table \b of \b contents + - \ref WhyInteresting + - \ref ConstructingVectors + - \ref ConstructionOfSumXpr + - \ref Assignment +\n + +
+ + +Consider the following example program: + +\code +#include + +int main() +{ + int size = 50; + // VectorXf is a vector of floats, with dynamic size. + Eigen::VectorXf u(size), v(size), w(size); + u = v + w; +} +\endcode + +The goal of this page is to understand how Eigen compiles it, assuming that SSE2 vectorization is enabled (GCC option -msse2). + +\section WhyInteresting Why it's interesting + +Maybe you think, that the above example program is so simple, that compiling it shouldn't involve anything interesting. So before starting, let us explain what is nontrivial in compiling it correctly -- that is, producing optimized code -- so that the complexity of Eigen, that we'll explain here, is really useful. + +Look at the line of code +\code + u = v + w; // (*) +\endcode + +The first important thing about compiling it, is that the arrays should be traversed only once, like +\code + for(int i = 0; i < size; i++) u[i] = v[i] + w[i]; +\endcode +The problem is that if we make a naive C++ library where the VectorXf class has an operator+ returning a VectorXf, then the line of code (*) will amount to: +\code + VectorXf tmp = v + w; + VectorXf u = tmp; +\endcode +Obviously, the introduction of the temporary \a tmp here is useless. It has a very bad effect on performance, first because the creation of \a tmp requires a dynamic memory allocation in this context, and second as there are now two for loops: +\code + for(int i = 0; i < size; i++) tmp[i] = v[i] + w[i]; + for(int i = 0; i < size; i++) u[i] = tmp[i]; +\endcode +Traversing the arrays twice instead of once is terrible for performance, as it means that we do many redundant memory accesses. + +The second important thing about compiling the above program, is to make correct use of SSE2 instructions. Notice that Eigen also supports AltiVec and that all the discussion that we make here applies also to AltiVec. + +SSE2, like AltiVec, is a set of instructions allowing to perform computations on packets of 128 bits at once. Since a float is 32 bits, this means that SSE2 instructions can handle 4 floats at once. This means that, if correctly used, they can make our computation go up to 4x faster. + +However, in the above program, we have chosen size=50, so our vectors consist of 50 float's, and 50 is not a multiple of 4. This means that we cannot hope to do all of that computation using SSE2 instructions. The second best thing, to which we should aim, is to handle the 48 first coefficients with SSE2 instructions, since 48 is the biggest multiple of 4 below 50, and then handle separately, without SSE2, the 49th and 50th coefficients. Something like this: + +\code + for(int i = 0; i < 4*(size/4); i+=4) u.packet(i) = v.packet(i) + w.packet(i); + for(int i = 4*(size/4); i < size; i++) u[i] = v[i] + w[i]; +\endcode + +So let us look line by line at our example program, and let's follow Eigen as it compiles it. + +\section ConstructingVectors Constructing vectors + +Let's analyze the first line: + +\code + Eigen::VectorXf u(size), v(size), w(size); +\endcode + +First of all, VectorXf is the following typedef: +\code + typedef Matrix VectorXf; +\endcode + +The class template Matrix is declared in src/Core/util/ForwardDeclarations.h with 6 template parameters, but the last 3 are automatically determined by the first 3. So you don't need to worry about them for now. Here, Matrix\ means a matrix of floats, with a dynamic number of rows and 1 column. + +The Matrix class inherits a base class, MatrixBase. Don't worry about it, for now it suffices to say that MatrixBase is what unifies matrices/vectors and all the expressions types -- more on that below. + +When we do +\code + Eigen::VectorXf u(size); +\endcode +the constructor that is called is Matrix::Matrix(int), in src/Core/Matrix.h. Besides some assertions, all it does is to construct the \a m_storage member, which is of type ei_matrix_storage\. + +You may wonder, isn't it overengineering to have the storage in a separate class? The reason is that the Matrix class template covers all kinds of matrices and vector: both fixed-size and dynamic-size. The storage method is not the same in these two cases. For fixed-size, the matrix coefficients are stored as a plain member array. For dynamic-size, the coefficients will be stored as a pointer to a dynamically-allocated array. Because of this, we need to abstract storage away from the Matrix class. That's ei_matrix_storage. + +Let's look at this constructor, in src/Core/MatrixStorage.h. You can see that there are many partial template specializations of ei_matrix_storages here, treating separately the cases where dimensions are Dynamic or fixed at compile-time. The partial specialization that we are looking at is: +\code +template class ei_matrix_storage +\endcode + +Here, the constructor called is ei_matrix_storage::ei_matrix_storage(int size, int rows, int columns) +with size=50, rows=50, columns=1. + +Here is this constructor: +\code +inline ei_matrix_storage(int size, int rows, int) : m_data(ei_aligned_new(size)), m_rows(rows) {} +\endcode + +Here, the \a m_data member is the actual array of coefficients of the matrix. As you see, it is dynamically allocated. Rather than calling new[] or malloc(), as you can see, we have our own ei_aligned_new defined in src/Core/util/Memory.h. What it does is that if vectorization is enabled, then it uses a platform-specific call to allocate a 128-bit-aligned array, as that is very useful for vectorization with both SSE2 and AltiVec. If vectorization is disabled, it amounts to the standard new[]. + +As you can see, the constructor also sets the \a m_rows member to \a size. Notice that there is no \a m_columns member: indeed, in this partial specialization of ei_matrix_storage, we know the number of columns at compile-time, since the _Cols template parameter is different from Dynamic. Namely, in our case, _Cols is 1, which is to say that our vector is just a matrix with 1 column. Hence, there is no need to store the number of columns as a runtime variable. + +When you call VectorXf::data() to get the pointer to the array of coefficients, it returns ei_matrix_storage::data() which returns the \a m_data member. + +When you call VectorXf::size() to get the size of the vector, this is actually a method in the base class MatrixBase. It determines that the vector is a column-vector, since ColsAtCompileTime==1 (this comes from the template parameters in the typedef VectorXf). It deduces that the size is the number of rows, so it returns VectorXf::rows(), which returns ei_matrix_storage::rows(), which returns the \a m_rows member, which was set to \a size by the constructor. + +\section ConstructionOfSumXpr Construction of the sum expression + +Now that our vectors are constructed, let's move on to the next line: + +\code +u = v + w; +\endcode + +The executive summary is that operator+ returns a "sum of vectors" expression, but doesn't actually perform the computation. It is the operator=, whose call occurs thereafter, that does the computation. + +Let us now see what Eigen does when it sees this: + +\code +v + w +\endcode + +Here, v and w are of type VectorXf, which is a typedef for a specialization of Matrix (as we explained above), which is a subclass of MatrixBase. So what is being called is + +\code +MatrixBase::operator+(const MatrixBase&) +\endcode + +The return type of this operator is +\code +CwiseBinaryOp, VectorXf, VectorXf> +\endcode +The CwiseBinaryOp class is our first encounter with an expression template. As we said, the operator+ doesn't by itself perform any computation, it just returns an abstract "sum of vectors" expression. Since there are also "difference of vectors" and "coefficient-wise product of vectors" expressions, we unify them all as "coefficient-wise binary operations", which we abbreviate as "CwiseBinaryOp". "Coefficient-wise" means that the operations is performed coefficient by coefficient. "binary" means that there are two operands -- we are adding two vectors with one another. + +Now you might ask, what if we did something like + +\code +v + w + u; +\endcode + +The first v + w would return a CwiseBinaryOp as above, so in order for this to compile, we'd need to define an operator+ also in the class CwiseBinaryOp... at this point it starts looking like a nightmare: are we going to have to define all operators in each of the expression classes (as you guessed, CwiseBinaryOp is only one of many) ? This looks like a dead end! + +The solution is that CwiseBinaryOp itself, as well as Matrix and all the other expression types, is a subclass of MatrixBase. So it is enough to define once and for all the operators in class MatrixBase. + +Since MatrixBase is the common base class of different subclasses, the aspects that depend on the subclass must be abstracted from MatrixBase. This is called polymorphism. + +The classical approach to polymorphism in C++ is by means of virtual functions. This is dynamic polymorphism. Here we don't want dynamic polymorphism because the whole design of Eigen is based around the assumption that all the complexity, all the abstraction, gets resolved at compile-time. This is crucial: if the abstraction can't get resolved at compile-time, Eigen's compile-time optimization mechanisms become useless, not to mention that if that abstraction has to be resolved at runtime it'll incur an overhead by itself. + +Here, what we want is to have a single class MatrixBase as the base of many subclasses, in such a way that each MatrixBase object (be it a matrix, or vector, or any kind of expression) knows at compile-time (as opposed to run-time) of which particular subclass it is an object (i.e. whether it is a matrix, or an expression, and what kind of expression). + +The solution is the Curiously Recurring Template Pattern. Let's do the break now. Hopefully you can read this wikipedia page during the break if needed, but it won't be allowed during the exam. + +In short, MatrixBase takes a template parameter \a Derived. Whenever we define a subclass Subclass, we actually make Subclass inherit MatrixBase\. The point is that different subclasses inherit different MatrixBase types. Thanks to this, whenever we have an object of a subclass, and we call on it some MatrixBase method, we still remember even from inside the MatrixBase method which particular subclass we're talking about. + +This means that we can put almost all the methods and operators in the base class MatrixBase, and have only the bare minimum in the subclasses. If you look at the subclasses in Eigen, like for instance the CwiseBinaryOp class, they have very few methods. There are coeff() and sometimes coeffRef() methods for access to the coefficients, there are rows() and cols() methods returning the number of rows and columns, but there isn't much more than that. All the meat is in MatrixBase, so it only needs to be coded once for all kinds of expressions, matrices, and vectors. + +So let's end this digression and come back to the piece of code from our example program that we were currently analyzing, + +\code +v + w +\endcode + +Now that MatrixBase is a good friend, let's write fully the prototype of the operator+ that gets called here (this code is from src/Core/MatrixBase.h): + +\code +template +class MatrixBase +{ + // ... + + template + const CwiseBinaryOp::Scalar>, Derived, OtherDerived> + operator+(const MatrixBase &other) const; + + // ... +}; +\endcode + +Here of course, \a Derived and \a OtherDerived are VectorXf. + +As we said, CwiseBinaryOp is also used for other operations such as substration, so it takes another template parameter determining the operation that will be applied to coefficients. This template parameter is a functor, that is, a class in which we have an operator() so it behaves like a function. Here, the functor used is ei_scalar_sum_op. It is defined in src/Core/Functors.h. + +Let us now explain the ei_traits here. The ei_scalar_sum_op class takes one template parameter: the type of the numbers to handle. Here of course we want to pass the scalar type (a.k.a. numeric type) of VectorXf, which is \c float. How do we determine which is the scalar type of \a Derived ? Throughout Eigen, all matrix and expression types define a typedef \a Scalar which gives its scalar type. For example, VectorXf::Scalar is a typedef for \c float. So here, if life was easy, we could find the numeric type of \a Derived as just +\code +typename Derived::Scalar +\endcode +Unfortunately, we can't do that here, as the compiler would complain that the type Derived hasn't yet been defined. So we use a workaround: in src/Core/util/ForwardDeclarations.h, we declared (not defined!) all our subclasses, like Matrix, and we also declared the following class template: +\code +template struct ei_traits; +\endcode +In src/Core/Matrix.h, right \em before the definition of class Matrix, we define a partial specialization of ei_traits for T=Matrix\. In this specialization of ei_traits, we define the Scalar typedef. So when we actually define Matrix, it is legal to refer to "typename ei_traits\::Scalar". + +Anyway, we have declared our operator+. In our case, where \a Derived and \a OtherDerived are VectorXf, the above declaration amounts to: +\code +class MatrixBase +{ + // ... + + const CwiseBinaryOp, VectorXf, VectorXf> + operator+(const MatrixBase &other) const; + + // ... +}; +\endcode + +Let's now jump to src/Core/CwiseBinaryOp.h to see how it is defined. As you can see there, all it does is to return a CwiseBinaryOp object, and this object is just storing references to the left-hand-side and right-hand-side expressions -- here, these are the vectors \a v and \a w. Well, the CwiseBinaryOp object is also storing an instance of the (empty) functor class, but you shouldn't worry about it as that is a minor implementation detail. + +Thus, the operator+ hasn't performed any actual computation. To summarize, the operation \a v + \a w just returned an object of type CwiseBinaryOp which did nothing else than just storing references to \a v and \a w. + +\section Assignment The assignment + +At this point, the expression \a v + \a w has finished evaluating, so, in the process of compiling the line of code +\code +u = v + w; +\endcode +we now enter the operator=. + +What operator= is being called here? The vector u is an object of class VectorXf, i.e. Matrix. In src/Core/Matrix.h, inside the definition of class Matrix, we see this: +\code + template + inline Matrix& operator=(const MatrixBase& other) + { + ei_assert(m_storage.data()!=0 && "you cannot use operator= with a non initialized matrix (instead use set()"); + return Base::operator=(other.derived()); + } +\endcode +Here, Base is a typedef for MatrixBase\. So, what is being called is the operator= of MatrixBase. Let's see its prototype in src/Core/MatrixBase.h: +\code + template + Derived& operator=(const MatrixBase& other); +\endcode +Here, \a Derived is VectorXf (since u is a VectorXf) and \a OtherDerived is CwiseBinaryOp. More specifically, as explained in the previous section, \a OtherDerived is: +\code +CwiseBinaryOp, VectorXf, VectorXf> +\endcode +So the full prototype of the operator= being called is: +\code +VectorXf& MatrixBase::operator=(const MatrixBase, VectorXf, VectorXf> > & other); +\endcode +This operator= literally reads "copying a sum of two VectorXf's into another VectorXf". + +Let's now look at the implementation of this operator=. It resides in the file src/Core/Assign.h. + +What we can see there is: +\code +template +template +inline Derived& MatrixBase + ::operator=(const MatrixBase& other) +{ + return ei_assign_selector::run(derived(), other.derived()); +} +\endcode + +OK so our next task is to understand ei_assign_selector :) + +Here is its declaration (all that is still in the same file src/Core/Assign.h) +\code +template +struct ei_assign_selector; +\endcode + +So ei_assign_selector takes 4 template parameters, but the 2 last ones are automatically determined by the 2 first ones. + +EvalBeforeAssigning is here to enforce the EvalBeforeAssigningBit. As explained here, certain expressions have this flag which makes them automatically evaluate into temporaries before assigning them to another expression. This is the case of the Product expression, in order to avoid strange aliasing effects when doing "m = m * m;" However, of course here our CwiseBinaryOp expression doesn't have the EvalBeforeAssigningBit: we said since the beginning that we didn't want a temporary to be introduced here. So if you go to src/Core/CwiseBinaryOp.h, you'll see that the Flags in ei_traits\ don't include the EvalBeforeAssigningBit. The Flags member of CwiseBinaryOp is then imported from the ei_traits by the EIGEN_GENERIC_PUBLIC_INTERFACE macro. Anyway, here the template parameter EvalBeforeAssigning has the value \c false. + +NeedToTranspose is here for the case where the user wants to copy a row-vector into a column-vector. We allow this as a special exception to the general rule that in assignments we require the dimesions to match. Anyway, here both the left-hand and right-hand sides are column vectors, in the sense that ColsAtCompileTime is equal to 1. So NeedToTranspose is \c false too. + +So, here we are in the partial specialization: +\code +ei_assign_selector +\endcode + +Here's how it is defined: +\code +template +struct ei_assign_selector { + static Derived& run(Derived& dst, const OtherDerived& other) { return dst.lazyAssign(other.derived()); } +}; +\endcode + +OK so now our next job is to understand how lazyAssign works :) + +\code +template +template +inline Derived& MatrixBase + ::lazyAssign(const MatrixBase& other) +{ + EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Derived,OtherDerived) + ei_assert(rows() == other.rows() && cols() == other.cols()); + ei_assign_impl::run(derived(),other.derived()); + return derived(); +} +\endcode + +What do we see here? Some assertions, and then the only interesting line is: +\code + ei_assign_impl::run(derived(),other.derived()); +\endcode + +OK so now we want to know what is inside ei_assign_impl. + +Here is its declaration: +\code +template::Vectorization, + int Unrolling = ei_assign_traits::Unrolling> +struct ei_assign_impl; +\endcode +Again, ei_assign_selector takes 4 template parameters, but the 2 last ones are automatically determined by the 2 first ones. + +These two parameters \a Vectorization and \a Unrolling are determined by a helper class ei_assign_traits. Its job is to determine which vectorization strategy to use (that is \a Vectorization) and which unrolling strategy to use (that is \a Unrolling). + +We'll not enter into the details of how these strategies are chosen (this is in the implementation of ei_assign_traits at the top of the same file). Let's just say that here \a Vectorization has the value \a LinearVectorization, and \a Unrolling has the value \a NoUnrolling (the latter is obvious since our vectors have dynamic size so there's no way to unroll the loop at compile-time). + +So the partial specialization of ei_assign_impl that we're looking at is: +\code +ei_assign_impl +\endcode + +Here is how it's defined: +\code +template +struct ei_assign_impl +{ + static void run(Derived1 &dst, const Derived2 &src) + { + const int size = dst.size(); + const int packetSize = ei_packet_traits::size; + const int alignedStart = ei_assign_traits::DstIsAligned ? 0 + : ei_alignmentOffset(&dst.coeffRef(0), size); + const int alignedEnd = alignedStart + ((size-alignedStart)/packetSize)*packetSize; + + for(int index = 0; index < alignedStart; index++) + dst.copyCoeff(index, src); + + for(int index = alignedStart; index < alignedEnd; index += packetSize) + { + dst.template copyPacket::SrcAlignment>(index, src); + } + + for(int index = alignedEnd; index < size; index++) + dst.copyCoeff(index, src); + } +}; +\endcode + +Here's how it works. \a LinearVectorization means that the left-hand and right-hand side expression can be accessed linearly i.e. you can refer to their coefficients by one integer \a index, as opposed to having to refer to its coefficients by two integers \a row, \a column. + +As we said at the beginning, vectorization works with blocks of 4 floats. Here, \a PacketSize is 4. + +There are two potential problems that we need to deal with: +\li first, vectorization works much better if the packets are 128-bit-aligned. This is especially important for write access. So when writing to the coefficients of \a dst, we want to group these coefficients by packets of 4 such that each of these packets is 128-bit-aligned. In general, this requires to skip a few coefficients at the beginning of \a dst. This is the purpose of \a alignedStart. We then copy these first few coefficients one by one, not by packets. However, in our case, the \a dst expression is a VectorXf and remember that in the construction of the vectors we allocated aligned arrays. Thanks to \a DstIsAligned, Eigen remembers that without having to do any runtime check, so \a alignedStart is zero and this part is avoided altogether. +\li second, the number of coefficients to copy is not in general a multiple of \a packetSize. Here, there are 50 coefficients to copy and \a packetSize is 4. So we'll have to copy the last 2 coefficients one by one, not by packets. Here, \a alignedEnd is 48. + +Now come the actual loops. + +First, the vectorized part: the 48 first coefficients out of 50 will be copied by packets of 4: +\code + for(int index = alignedStart; index < alignedEnd; index += packetSize) + { + dst.template copyPacket::SrcAlignment>(index, src); + } +\endcode + +What is copyPacket? It is defined in src/Core/Coeffs.h: +\code +template +template +inline void MatrixBase::copyPacket(int index, const MatrixBase& other) +{ + ei_internal_assert(index >= 0 && index < size()); + derived().template writePacket(index, + other.derived().template packet(index)); +} +\endcode + +OK, what are writePacket() and packet() here? + +First, writePacket() here is a method on the left-hand side VectorXf. So we go to src/Core/Matrix.h to look at its definition: +\code +template +inline void writePacket(int index, const PacketScalar& x) +{ + ei_pstoret(m_storage.data() + index, x); +} +\endcode +Here, \a StoreMode is \a Aligned, indicating that we are doing a 128-bit-aligned write access, \a PacketScalar is a type representing a "SSE packet of 4 floats" and ei_pstoret is a function writing such a packet in memory. Their definitions are architecture-specific, we find them in src/Core/arch/SSE/PacketMath.h: + +The line in src/Core/arch/SSE/PacketMath.h that determines the PacketScalar type (via a typedef in Matrix.h) is: +\code +template<> struct ei_packet_traits { typedef __m128 type; enum {size=4}; }; +\endcode +Here, __m128 is a SSE-specific type. Notice that the enum \a size here is what was used to define \a packetSize above. + +And here is the implementation of ei_pstoret: +\code +template<> inline void ei_pstore(float* to, const __m128& from) { _mm_store_ps(to, from); } +\endcode +Here, __mm_store_ps is a SSE-specific intrinsic function, representing a single SSE instruction. The difference between ei_pstore and ei_pstoret is that ei_pstoret is a dispatcher handling both the aligned and unaligned cases, you find its definition in src/Core/GenericPacketMath.h: +\code +template +inline void ei_pstoret(Scalar* to, const Packet& from) +{ + if(LoadMode == Aligned) + ei_pstore(to, from); + else + ei_pstoreu(to, from); +} +\endcode + +OK, that explains how writePacket() works. Now let's look into the packet() call. Remember that we are analyzing this line of code inside copyPacket(): +\code +derived().template writePacket(index, + other.derived().template packet(index)); +\endcode + +Here, \a other is our sum expression \a v + \a w. The .derived() is just casting from MatrixBase to the subclass which here is CwiseBinaryOp. So let's go to src/Core/CwiseBinaryOp.h: +\code +class CwiseBinaryOp +{ + // ... + template + inline PacketScalar packet(int index) const + { + return m_functor.packetOp(m_lhs.template packet(index), m_rhs.template packet(index)); + } +}; +\endcode +Here, \a m_lhs is the vector \a v, and \a m_rhs is the vector \a w. So the packet() function here is Matrix::packet(). The template parameter \a LoadMode is \a Aligned. So we're looking at +\code +class Matrix +{ + // ... + template + inline PacketScalar packet(int index) const + { + return ei_ploadt(m_storage.data() + index); + } +}; +\endcode +We let you look up the definition of ei_ploadt in GenericPacketMath.h and the ei_pload in src/Core/arch/SSE/PacketMath.h. It is very similar to the above for ei_pstore. + +Let's go back to CwiseBinaryOp::packet(). Once the packets from the vectors \a v and \a w have been returned, what does this function do? It calls m_functor.packetOp() on them. What is m_functor? Here we must remember what particular template specialization of CwiseBinaryOp we're dealing with: +\code +CwiseBinaryOp, VectorXf, VectorXf> +\endcode +So m_functor is an object of the empty class ei_scalar_sum_op. As we mentioned above, don't worry about why we constructed an object of this empty class at all -- it's an implementation detail, the point is that some other functors need to store member data. + +Anyway, ei_scalar_sum_op is defined in src/Core/Functors.h: +\code +template struct ei_scalar_sum_op EIGEN_EMPTY_STRUCT { + inline const Scalar operator() (const Scalar& a, const Scalar& b) const { return a + b; } + template + inline const PacketScalar packetOp(const PacketScalar& a, const PacketScalar& b) const + { return ei_padd(a,b); } +}; +\endcode +As you can see, all what packetOp() does is to call ei_padd on the two packets. Here is the definition of ei_padd from src/Core/arch/SSE/PacketMath.h: +\code +template<> inline __m128 ei_padd(const __m128& a, const __m128& b) { return _mm_add_ps(a,b); } +\endcode +Here, _mm_add_ps is a SSE-specific intrinsic function, representing a single SSE instruction. + +To summarize, the loop +\code + for(int index = alignedStart; index < alignedEnd; index += packetSize) + { + dst.template copyPacket::SrcAlignment>(index, src); + } +\endcode +has been compiled to the following code: for \a index going from 0 to the 11 ( = 48/4 - 1), read the i-th packet (of 4 floats) from the vector v and the i-th packet from the vector w using two __mm_load_ps SSE instructions, then add them together using a __mm_add_ps instruction, then store the result using a __mm_store_ps instruction. + +There remains the second loop handling the last few (here, the last 2) coefficients: +\code + for(int index = alignedEnd; index < size; index++) + dst.copyCoeff(index, src); +\endcode +However, it works just like the one we just explained, it is just simpler because there is no SSE vectorization involved here. copyPacket() becomes copyCoeff(), packet() becomes coeff(), writePacket() becomes coeffRef(). If you followed us this far, you can probably understand this part by yourself. + +We see that all the C++ abstraction of Eigen goes away during compilation and that we indeed are precisely controlling which assembly instructions we emit. Such is the beauty of C++! Since we have such precise control over the emitted assembly instructions, but such complex logic to choose the right instructions, we can say that Eigen really behaves like an optimizing compiler. If you prefer, you could say that Eigen behaves like a script for the compiler. In a sense, C++ template metaprogramming is scripting the compiler -- and it's been shown that this scripting language is Turing-complete. See Wikipedia. + +*/ + +} diff --git a/doc/I05_FixedSizeVectorizable.dox b/doc/I05_FixedSizeVectorizable.dox new file mode 100644 index 000000000..16f0c7f50 --- /dev/null +++ b/doc/I05_FixedSizeVectorizable.dox @@ -0,0 +1,38 @@ +namespace Eigen { + +/** \page FixedSizeVectorizable Advanced - Fixed-size vectorizable Eigen objects + +The goal of this page is to explain what we mean by "fixed-size vectorizable". + +\section summary Executive Summary + +An Eigen object is called "fixed-size vectorizable" if it has fixed size and that size is a multiple of 16 bytes. + +Examples include: +\li Eigen::Vector2d +\li Eigen::Vector4d +\li Eigen::Vector4f +\li Eigen::Matrix2d +\li Eigen::Matrix2f +\li Eigen::Matrix4d +\li Eigen::Matrix4f +\li Eigen::Transform3d +\li Eigen::Transform3f +\li Eigen::Quaterniond +\li Eigen::Quaternionf + +\section explanation Explanation + +First, "fixed-size" should be clear: an Eigen object has fixed size if its number of rows and its number of columns are fixed at compile-time. So for example Matrix3f has fixed size, but MatrixXf doesn't (the opposite of fixed-size is dynamic-size). + +The array of coefficients of a fixed-size Eigen object is a plain "static array", it is not dynamically allocated. For example, the data behind a Matrix4f is just a "float array[16]". + +Fixed-size objects are typically very small, which means that we want to handle them with zero runtime overhead -- both in terms of memory usage and of speed. + +Now, vectorization (both SSE and AltiVec) works with 128-bit packets. Moreover, for performance reasons, these packets need to be have 128-bit alignment. + +So it turns out that the only way that fixed-size Eigen objects can be vectorized, is if their size is a multiple of 128 bits, or 16 bytes. Eigen will then request 16-byte alignment for these object, and henceforth rely on these objects being aligned so no runtime check for alignment is performed. + +*/ + +} diff --git a/doc/InsideEigenExample.dox b/doc/InsideEigenExample.dox deleted file mode 100644 index 3242acf4a..000000000 --- a/doc/InsideEigenExample.dox +++ /dev/null @@ -1,500 +0,0 @@ -namespace Eigen { - -/** \page InsideEigenExample What happens inside Eigen, on a simple example - -\b Table \b of \b contents - - \ref WhyInteresting - - \ref ConstructingVectors - - \ref ConstructionOfSumXpr - - \ref Assignment -\n - -
- - -Consider the following example program: - -\code -#include - -int main() -{ - int size = 50; - // VectorXf is a vector of floats, with dynamic size. - Eigen::VectorXf u(size), v(size), w(size); - u = v + w; -} -\endcode - -The goal of this page is to understand how Eigen compiles it, assuming that SSE2 vectorization is enabled (GCC option -msse2). - -\section WhyInteresting Why it's interesting - -Maybe you think, that the above example program is so simple, that compiling it shouldn't involve anything interesting. So before starting, let us explain what is nontrivial in compiling it correctly -- that is, producing optimized code -- so that the complexity of Eigen, that we'll explain here, is really useful. - -Look at the line of code -\code - u = v + w; // (*) -\endcode - -The first important thing about compiling it, is that the arrays should be traversed only once, like -\code - for(int i = 0; i < size; i++) u[i] = v[i] + w[i]; -\endcode -The problem is that if we make a naive C++ library where the VectorXf class has an operator+ returning a VectorXf, then the line of code (*) will amount to: -\code - VectorXf tmp = v + w; - VectorXf u = tmp; -\endcode -Obviously, the introduction of the temporary \a tmp here is useless. It has a very bad effect on performance, first because the creation of \a tmp requires a dynamic memory allocation in this context, and second as there are now two for loops: -\code - for(int i = 0; i < size; i++) tmp[i] = v[i] + w[i]; - for(int i = 0; i < size; i++) u[i] = tmp[i]; -\endcode -Traversing the arrays twice instead of once is terrible for performance, as it means that we do many redundant memory accesses. - -The second important thing about compiling the above program, is to make correct use of SSE2 instructions. Notice that Eigen also supports AltiVec and that all the discussion that we make here applies also to AltiVec. - -SSE2, like AltiVec, is a set of instructions allowing to perform computations on packets of 128 bits at once. Since a float is 32 bits, this means that SSE2 instructions can handle 4 floats at once. This means that, if correctly used, they can make our computation go up to 4x faster. - -However, in the above program, we have chosen size=50, so our vectors consist of 50 float's, and 50 is not a multiple of 4. This means that we cannot hope to do all of that computation using SSE2 instructions. The second best thing, to which we should aim, is to handle the 48 first coefficients with SSE2 instructions, since 48 is the biggest multiple of 4 below 50, and then handle separately, without SSE2, the 49th and 50th coefficients. Something like this: - -\code - for(int i = 0; i < 4*(size/4); i+=4) u.packet(i) = v.packet(i) + w.packet(i); - for(int i = 4*(size/4); i < size; i++) u[i] = v[i] + w[i]; -\endcode - -So let us look line by line at our example program, and let's follow Eigen as it compiles it. - -\section ConstructingVectors Constructing vectors - -Let's analyze the first line: - -\code - Eigen::VectorXf u(size), v(size), w(size); -\endcode - -First of all, VectorXf is the following typedef: -\code - typedef Matrix VectorXf; -\endcode - -The class template Matrix is declared in src/Core/util/ForwardDeclarations.h with 6 template parameters, but the last 3 are automatically determined by the first 3. So you don't need to worry about them for now. Here, Matrix\ means a matrix of floats, with a dynamic number of rows and 1 column. - -The Matrix class inherits a base class, MatrixBase. Don't worry about it, for now it suffices to say that MatrixBase is what unifies matrices/vectors and all the expressions types -- more on that below. - -When we do -\code - Eigen::VectorXf u(size); -\endcode -the constructor that is called is Matrix::Matrix(int), in src/Core/Matrix.h. Besides some assertions, all it does is to construct the \a m_storage member, which is of type ei_matrix_storage\. - -You may wonder, isn't it overengineering to have the storage in a separate class? The reason is that the Matrix class template covers all kinds of matrices and vector: both fixed-size and dynamic-size. The storage method is not the same in these two cases. For fixed-size, the matrix coefficients are stored as a plain member array. For dynamic-size, the coefficients will be stored as a pointer to a dynamically-allocated array. Because of this, we need to abstract storage away from the Matrix class. That's ei_matrix_storage. - -Let's look at this constructor, in src/Core/MatrixStorage.h. You can see that there are many partial template specializations of ei_matrix_storages here, treating separately the cases where dimensions are Dynamic or fixed at compile-time. The partial specialization that we are looking at is: -\code -template class ei_matrix_storage -\endcode - -Here, the constructor called is ei_matrix_storage::ei_matrix_storage(int size, int rows, int columns) -with size=50, rows=50, columns=1. - -Here is this constructor: -\code -inline ei_matrix_storage(int size, int rows, int) : m_data(ei_aligned_new(size)), m_rows(rows) {} -\endcode - -Here, the \a m_data member is the actual array of coefficients of the matrix. As you see, it is dynamically allocated. Rather than calling new[] or malloc(), as you can see, we have our own ei_aligned_new defined in src/Core/util/Memory.h. What it does is that if vectorization is enabled, then it uses a platform-specific call to allocate a 128-bit-aligned array, as that is very useful for vectorization with both SSE2 and AltiVec. If vectorization is disabled, it amounts to the standard new[]. - -As you can see, the constructor also sets the \a m_rows member to \a size. Notice that there is no \a m_columns member: indeed, in this partial specialization of ei_matrix_storage, we know the number of columns at compile-time, since the _Cols template parameter is different from Dynamic. Namely, in our case, _Cols is 1, which is to say that our vector is just a matrix with 1 column. Hence, there is no need to store the number of columns as a runtime variable. - -When you call VectorXf::data() to get the pointer to the array of coefficients, it returns ei_matrix_storage::data() which returns the \a m_data member. - -When you call VectorXf::size() to get the size of the vector, this is actually a method in the base class MatrixBase. It determines that the vector is a column-vector, since ColsAtCompileTime==1 (this comes from the template parameters in the typedef VectorXf). It deduces that the size is the number of rows, so it returns VectorXf::rows(), which returns ei_matrix_storage::rows(), which returns the \a m_rows member, which was set to \a size by the constructor. - -\section ConstructionOfSumXpr Construction of the sum expression - -Now that our vectors are constructed, let's move on to the next line: - -\code -u = v + w; -\endcode - -The executive summary is that operator+ returns a "sum of vectors" expression, but doesn't actually perform the computation. It is the operator=, whose call occurs thereafter, that does the computation. - -Let us now see what Eigen does when it sees this: - -\code -v + w -\endcode - -Here, v and w are of type VectorXf, which is a typedef for a specialization of Matrix (as we explained above), which is a subclass of MatrixBase. So what is being called is - -\code -MatrixBase::operator+(const MatrixBase&) -\endcode - -The return type of this operator is -\code -CwiseBinaryOp, VectorXf, VectorXf> -\endcode -The CwiseBinaryOp class is our first encounter with an expression template. As we said, the operator+ doesn't by itself perform any computation, it just returns an abstract "sum of vectors" expression. Since there are also "difference of vectors" and "coefficient-wise product of vectors" expressions, we unify them all as "coefficient-wise binary operations", which we abbreviate as "CwiseBinaryOp". "Coefficient-wise" means that the operations is performed coefficient by coefficient. "binary" means that there are two operands -- we are adding two vectors with one another. - -Now you might ask, what if we did something like - -\code -v + w + u; -\endcode - -The first v + w would return a CwiseBinaryOp as above, so in order for this to compile, we'd need to define an operator+ also in the class CwiseBinaryOp... at this point it starts looking like a nightmare: are we going to have to define all operators in each of the expression classes (as you guessed, CwiseBinaryOp is only one of many) ? This looks like a dead end! - -The solution is that CwiseBinaryOp itself, as well as Matrix and all the other expression types, is a subclass of MatrixBase. So it is enough to define once and for all the operators in class MatrixBase. - -Since MatrixBase is the common base class of different subclasses, the aspects that depend on the subclass must be abstracted from MatrixBase. This is called polymorphism. - -The classical approach to polymorphism in C++ is by means of virtual functions. This is dynamic polymorphism. Here we don't want dynamic polymorphism because the whole design of Eigen is based around the assumption that all the complexity, all the abstraction, gets resolved at compile-time. This is crucial: if the abstraction can't get resolved at compile-time, Eigen's compile-time optimization mechanisms become useless, not to mention that if that abstraction has to be resolved at runtime it'll incur an overhead by itself. - -Here, what we want is to have a single class MatrixBase as the base of many subclasses, in such a way that each MatrixBase object (be it a matrix, or vector, or any kind of expression) knows at compile-time (as opposed to run-time) of which particular subclass it is an object (i.e. whether it is a matrix, or an expression, and what kind of expression). - -The solution is the Curiously Recurring Template Pattern. Let's do the break now. Hopefully you can read this wikipedia page during the break if needed, but it won't be allowed during the exam. - -In short, MatrixBase takes a template parameter \a Derived. Whenever we define a subclass Subclass, we actually make Subclass inherit MatrixBase\. The point is that different subclasses inherit different MatrixBase types. Thanks to this, whenever we have an object of a subclass, and we call on it some MatrixBase method, we still remember even from inside the MatrixBase method which particular subclass we're talking about. - -This means that we can put almost all the methods and operators in the base class MatrixBase, and have only the bare minimum in the subclasses. If you look at the subclasses in Eigen, like for instance the CwiseBinaryOp class, they have very few methods. There are coeff() and sometimes coeffRef() methods for access to the coefficients, there are rows() and cols() methods returning the number of rows and columns, but there isn't much more than that. All the meat is in MatrixBase, so it only needs to be coded once for all kinds of expressions, matrices, and vectors. - -So let's end this digression and come back to the piece of code from our example program that we were currently analyzing, - -\code -v + w -\endcode - -Now that MatrixBase is a good friend, let's write fully the prototype of the operator+ that gets called here (this code is from src/Core/MatrixBase.h): - -\code -template -class MatrixBase -{ - // ... - - template - const CwiseBinaryOp::Scalar>, Derived, OtherDerived> - operator+(const MatrixBase &other) const; - - // ... -}; -\endcode - -Here of course, \a Derived and \a OtherDerived are VectorXf. - -As we said, CwiseBinaryOp is also used for other operations such as substration, so it takes another template parameter determining the operation that will be applied to coefficients. This template parameter is a functor, that is, a class in which we have an operator() so it behaves like a function. Here, the functor used is ei_scalar_sum_op. It is defined in src/Core/Functors.h. - -Let us now explain the ei_traits here. The ei_scalar_sum_op class takes one template parameter: the type of the numbers to handle. Here of course we want to pass the scalar type (a.k.a. numeric type) of VectorXf, which is \c float. How do we determine which is the scalar type of \a Derived ? Throughout Eigen, all matrix and expression types define a typedef \a Scalar which gives its scalar type. For example, VectorXf::Scalar is a typedef for \c float. So here, if life was easy, we could find the numeric type of \a Derived as just -\code -typename Derived::Scalar -\endcode -Unfortunately, we can't do that here, as the compiler would complain that the type Derived hasn't yet been defined. So we use a workaround: in src/Core/util/ForwardDeclarations.h, we declared (not defined!) all our subclasses, like Matrix, and we also declared the following class template: -\code -template struct ei_traits; -\endcode -In src/Core/Matrix.h, right \em before the definition of class Matrix, we define a partial specialization of ei_traits for T=Matrix\. In this specialization of ei_traits, we define the Scalar typedef. So when we actually define Matrix, it is legal to refer to "typename ei_traits\::Scalar". - -Anyway, we have declared our operator+. In our case, where \a Derived and \a OtherDerived are VectorXf, the above declaration amounts to: -\code -class MatrixBase -{ - // ... - - const CwiseBinaryOp, VectorXf, VectorXf> - operator+(const MatrixBase &other) const; - - // ... -}; -\endcode - -Let's now jump to src/Core/CwiseBinaryOp.h to see how it is defined. As you can see there, all it does is to return a CwiseBinaryOp object, and this object is just storing references to the left-hand-side and right-hand-side expressions -- here, these are the vectors \a v and \a w. Well, the CwiseBinaryOp object is also storing an instance of the (empty) functor class, but you shouldn't worry about it as that is a minor implementation detail. - -Thus, the operator+ hasn't performed any actual computation. To summarize, the operation \a v + \a w just returned an object of type CwiseBinaryOp which did nothing else than just storing references to \a v and \a w. - -\section Assignment The assignment - -At this point, the expression \a v + \a w has finished evaluating, so, in the process of compiling the line of code -\code -u = v + w; -\endcode -we now enter the operator=. - -What operator= is being called here? The vector u is an object of class VectorXf, i.e. Matrix. In src/Core/Matrix.h, inside the definition of class Matrix, we see this: -\code - template - inline Matrix& operator=(const MatrixBase& other) - { - ei_assert(m_storage.data()!=0 && "you cannot use operator= with a non initialized matrix (instead use set()"); - return Base::operator=(other.derived()); - } -\endcode -Here, Base is a typedef for MatrixBase\. So, what is being called is the operator= of MatrixBase. Let's see its prototype in src/Core/MatrixBase.h: -\code - template - Derived& operator=(const MatrixBase& other); -\endcode -Here, \a Derived is VectorXf (since u is a VectorXf) and \a OtherDerived is CwiseBinaryOp. More specifically, as explained in the previous section, \a OtherDerived is: -\code -CwiseBinaryOp, VectorXf, VectorXf> -\endcode -So the full prototype of the operator= being called is: -\code -VectorXf& MatrixBase::operator=(const MatrixBase, VectorXf, VectorXf> > & other); -\endcode -This operator= literally reads "copying a sum of two VectorXf's into another VectorXf". - -Let's now look at the implementation of this operator=. It resides in the file src/Core/Assign.h. - -What we can see there is: -\code -template -template -inline Derived& MatrixBase - ::operator=(const MatrixBase& other) -{ - return ei_assign_selector::run(derived(), other.derived()); -} -\endcode - -OK so our next task is to understand ei_assign_selector :) - -Here is its declaration (all that is still in the same file src/Core/Assign.h) -\code -template -struct ei_assign_selector; -\endcode - -So ei_assign_selector takes 4 template parameters, but the 2 last ones are automatically determined by the 2 first ones. - -EvalBeforeAssigning is here to enforce the EvalBeforeAssigningBit. As explained here, certain expressions have this flag which makes them automatically evaluate into temporaries before assigning them to another expression. This is the case of the Product expression, in order to avoid strange aliasing effects when doing "m = m * m;" However, of course here our CwiseBinaryOp expression doesn't have the EvalBeforeAssigningBit: we said since the beginning that we didn't want a temporary to be introduced here. So if you go to src/Core/CwiseBinaryOp.h, you'll see that the Flags in ei_traits\ don't include the EvalBeforeAssigningBit. The Flags member of CwiseBinaryOp is then imported from the ei_traits by the EIGEN_GENERIC_PUBLIC_INTERFACE macro. Anyway, here the template parameter EvalBeforeAssigning has the value \c false. - -NeedToTranspose is here for the case where the user wants to copy a row-vector into a column-vector. We allow this as a special exception to the general rule that in assignments we require the dimesions to match. Anyway, here both the left-hand and right-hand sides are column vectors, in the sense that ColsAtCompileTime is equal to 1. So NeedToTranspose is \c false too. - -So, here we are in the partial specialization: -\code -ei_assign_selector -\endcode - -Here's how it is defined: -\code -template -struct ei_assign_selector { - static Derived& run(Derived& dst, const OtherDerived& other) { return dst.lazyAssign(other.derived()); } -}; -\endcode - -OK so now our next job is to understand how lazyAssign works :) - -\code -template -template -inline Derived& MatrixBase - ::lazyAssign(const MatrixBase& other) -{ - EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Derived,OtherDerived) - ei_assert(rows() == other.rows() && cols() == other.cols()); - ei_assign_impl::run(derived(),other.derived()); - return derived(); -} -\endcode - -What do we see here? Some assertions, and then the only interesting line is: -\code - ei_assign_impl::run(derived(),other.derived()); -\endcode - -OK so now we want to know what is inside ei_assign_impl. - -Here is its declaration: -\code -template::Vectorization, - int Unrolling = ei_assign_traits::Unrolling> -struct ei_assign_impl; -\endcode -Again, ei_assign_selector takes 4 template parameters, but the 2 last ones are automatically determined by the 2 first ones. - -These two parameters \a Vectorization and \a Unrolling are determined by a helper class ei_assign_traits. Its job is to determine which vectorization strategy to use (that is \a Vectorization) and which unrolling strategy to use (that is \a Unrolling). - -We'll not enter into the details of how these strategies are chosen (this is in the implementation of ei_assign_traits at the top of the same file). Let's just say that here \a Vectorization has the value \a LinearVectorization, and \a Unrolling has the value \a NoUnrolling (the latter is obvious since our vectors have dynamic size so there's no way to unroll the loop at compile-time). - -So the partial specialization of ei_assign_impl that we're looking at is: -\code -ei_assign_impl -\endcode - -Here is how it's defined: -\code -template -struct ei_assign_impl -{ - static void run(Derived1 &dst, const Derived2 &src) - { - const int size = dst.size(); - const int packetSize = ei_packet_traits::size; - const int alignedStart = ei_assign_traits::DstIsAligned ? 0 - : ei_alignmentOffset(&dst.coeffRef(0), size); - const int alignedEnd = alignedStart + ((size-alignedStart)/packetSize)*packetSize; - - for(int index = 0; index < alignedStart; index++) - dst.copyCoeff(index, src); - - for(int index = alignedStart; index < alignedEnd; index += packetSize) - { - dst.template copyPacket::SrcAlignment>(index, src); - } - - for(int index = alignedEnd; index < size; index++) - dst.copyCoeff(index, src); - } -}; -\endcode - -Here's how it works. \a LinearVectorization means that the left-hand and right-hand side expression can be accessed linearly i.e. you can refer to their coefficients by one integer \a index, as opposed to having to refer to its coefficients by two integers \a row, \a column. - -As we said at the beginning, vectorization works with blocks of 4 floats. Here, \a PacketSize is 4. - -There are two potential problems that we need to deal with: -\li first, vectorization works much better if the packets are 128-bit-aligned. This is especially important for write access. So when writing to the coefficients of \a dst, we want to group these coefficients by packets of 4 such that each of these packets is 128-bit-aligned. In general, this requires to skip a few coefficients at the beginning of \a dst. This is the purpose of \a alignedStart. We then copy these first few coefficients one by one, not by packets. However, in our case, the \a dst expression is a VectorXf and remember that in the construction of the vectors we allocated aligned arrays. Thanks to \a DstIsAligned, Eigen remembers that without having to do any runtime check, so \a alignedStart is zero and this part is avoided altogether. -\li second, the number of coefficients to copy is not in general a multiple of \a packetSize. Here, there are 50 coefficients to copy and \a packetSize is 4. So we'll have to copy the last 2 coefficients one by one, not by packets. Here, \a alignedEnd is 48. - -Now come the actual loops. - -First, the vectorized part: the 48 first coefficients out of 50 will be copied by packets of 4: -\code - for(int index = alignedStart; index < alignedEnd; index += packetSize) - { - dst.template copyPacket::SrcAlignment>(index, src); - } -\endcode - -What is copyPacket? It is defined in src/Core/Coeffs.h: -\code -template -template -inline void MatrixBase::copyPacket(int index, const MatrixBase& other) -{ - ei_internal_assert(index >= 0 && index < size()); - derived().template writePacket(index, - other.derived().template packet(index)); -} -\endcode - -OK, what are writePacket() and packet() here? - -First, writePacket() here is a method on the left-hand side VectorXf. So we go to src/Core/Matrix.h to look at its definition: -\code -template -inline void writePacket(int index, const PacketScalar& x) -{ - ei_pstoret(m_storage.data() + index, x); -} -\endcode -Here, \a StoreMode is \a Aligned, indicating that we are doing a 128-bit-aligned write access, \a PacketScalar is a type representing a "SSE packet of 4 floats" and ei_pstoret is a function writing such a packet in memory. Their definitions are architecture-specific, we find them in src/Core/arch/SSE/PacketMath.h: - -The line in src/Core/arch/SSE/PacketMath.h that determines the PacketScalar type (via a typedef in Matrix.h) is: -\code -template<> struct ei_packet_traits { typedef __m128 type; enum {size=4}; }; -\endcode -Here, __m128 is a SSE-specific type. Notice that the enum \a size here is what was used to define \a packetSize above. - -And here is the implementation of ei_pstoret: -\code -template<> inline void ei_pstore(float* to, const __m128& from) { _mm_store_ps(to, from); } -\endcode -Here, __mm_store_ps is a SSE-specific intrinsic function, representing a single SSE instruction. The difference between ei_pstore and ei_pstoret is that ei_pstoret is a dispatcher handling both the aligned and unaligned cases, you find its definition in src/Core/GenericPacketMath.h: -\code -template -inline void ei_pstoret(Scalar* to, const Packet& from) -{ - if(LoadMode == Aligned) - ei_pstore(to, from); - else - ei_pstoreu(to, from); -} -\endcode - -OK, that explains how writePacket() works. Now let's look into the packet() call. Remember that we are analyzing this line of code inside copyPacket(): -\code -derived().template writePacket(index, - other.derived().template packet(index)); -\endcode - -Here, \a other is our sum expression \a v + \a w. The .derived() is just casting from MatrixBase to the subclass which here is CwiseBinaryOp. So let's go to src/Core/CwiseBinaryOp.h: -\code -class CwiseBinaryOp -{ - // ... - template - inline PacketScalar packet(int index) const - { - return m_functor.packetOp(m_lhs.template packet(index), m_rhs.template packet(index)); - } -}; -\endcode -Here, \a m_lhs is the vector \a v, and \a m_rhs is the vector \a w. So the packet() function here is Matrix::packet(). The template parameter \a LoadMode is \a Aligned. So we're looking at -\code -class Matrix -{ - // ... - template - inline PacketScalar packet(int index) const - { - return ei_ploadt(m_storage.data() + index); - } -}; -\endcode -We let you look up the definition of ei_ploadt in GenericPacketMath.h and the ei_pload in src/Core/arch/SSE/PacketMath.h. It is very similar to the above for ei_pstore. - -Let's go back to CwiseBinaryOp::packet(). Once the packets from the vectors \a v and \a w have been returned, what does this function do? It calls m_functor.packetOp() on them. What is m_functor? Here we must remember what particular template specialization of CwiseBinaryOp we're dealing with: -\code -CwiseBinaryOp, VectorXf, VectorXf> -\endcode -So m_functor is an object of the empty class ei_scalar_sum_op. As we mentioned above, don't worry about why we constructed an object of this empty class at all -- it's an implementation detail, the point is that some other functors need to store member data. - -Anyway, ei_scalar_sum_op is defined in src/Core/Functors.h: -\code -template struct ei_scalar_sum_op EIGEN_EMPTY_STRUCT { - inline const Scalar operator() (const Scalar& a, const Scalar& b) const { return a + b; } - template - inline const PacketScalar packetOp(const PacketScalar& a, const PacketScalar& b) const - { return ei_padd(a,b); } -}; -\endcode -As you can see, all what packetOp() does is to call ei_padd on the two packets. Here is the definition of ei_padd from src/Core/arch/SSE/PacketMath.h: -\code -template<> inline __m128 ei_padd(const __m128& a, const __m128& b) { return _mm_add_ps(a,b); } -\endcode -Here, _mm_add_ps is a SSE-specific intrinsic function, representing a single SSE instruction. - -To summarize, the loop -\code - for(int index = alignedStart; index < alignedEnd; index += packetSize) - { - dst.template copyPacket::SrcAlignment>(index, src); - } -\endcode -has been compiled to the following code: for \a index going from 0 to the 11 ( = 48/4 - 1), read the i-th packet (of 4 floats) from the vector v and the i-th packet from the vector w using two __mm_load_ps SSE instructions, then add them together using a __mm_add_ps instruction, then store the result using a __mm_store_ps instruction. - -There remains the second loop handling the last few (here, the last 2) coefficients: -\code - for(int index = alignedEnd; index < size; index++) - dst.copyCoeff(index, src); -\endcode -However, it works just like the one we just explained, it is just simpler because there is no SSE vectorization involved here. copyPacket() becomes copyCoeff(), packet() becomes coeff(), writePacket() becomes coeffRef(). If you followed us this far, you can probably understand this part by yourself. - -We see that all the C++ abstraction of Eigen goes away during compilation and that we indeed are precisely controlling which assembly instructions we emit. Such is the beauty of C++! Since we have such precise control over the emitted assembly instructions, but such complex logic to choose the right instructions, we can say that Eigen really behaves like an optimizing compiler. If you prefer, you could say that Eigen behaves like a script for the compiler. In a sense, C++ template metaprogramming is scripting the compiler -- and it's been shown that this scripting language is Turing-complete. See Wikipedia. - -*/ - -} diff --git a/doc/PassingByValue.dox b/doc/PassingByValue.dox deleted file mode 100644 index 6eaf9ccb9..000000000 --- a/doc/PassingByValue.dox +++ /dev/null @@ -1,40 +0,0 @@ -namespace Eigen { - -/** \page PassingByValue Passing Eigen objects by value to functions - -Passing objects by value is almost always a very bad idea in C++, as this means useless copies, and one should pass them by reference instead. - -With Eigen, this is even more important: passing \ref FixedSizeVectorizable "fixed-size vectorizable Eigen objects" by value is not only inefficient, it can be illegal or make your program crash! And the reason is that these Eigen objects have alignment modifiers that aren't respected when they are passed by value. - -So for example, a function like this, where v is passed by value: - -\code -void my_function(Eigen::Vector2d v); -\endcode - -needs to be rewritten as follows, passing v by reference: - -\code -void my_function(const Eigen::Vector2d& v); -\endcode - -Likewise if you have a class having a Eigen object as member: - -\code -struct Foo -{ - Eigen::Vector2d v; -}; -void my_function(Foo v); -\endcode - -This function also needs to be rewritten like this: -\code -void my_function(const Foo& v); -\endcode - -Note that on the other hand, there is no problem with functions that return objects by value. - -*/ - -} diff --git a/doc/QuickStartGuide.dox b/doc/QuickStartGuide.dox deleted file mode 100644 index d7bec3bb7..000000000 --- a/doc/QuickStartGuide.dox +++ /dev/null @@ -1,565 +0,0 @@ -namespace Eigen { - -/** \page TutorialCore Tutorial 1/3 - Core features - \ingroup Tutorial - -
\ref index "Overview" - | \b Core \b features - | \ref TutorialGeometry "Geometry" - | \ref TutorialAdvancedLinearAlgebra "Advanced linear algebra" - | \ref TutorialSparse "Sparse matrix" -
- -\b Table \b of \b contents - - \ref TutorialCoreGettingStarted - - \ref TutorialCoreSimpleExampleFixedSize - - \ref TutorialCoreSimpleExampleDynamicSize - - \ref TutorialCoreMatrixTypes - - \ref TutorialCoreCoefficients - - \ref TutorialCoreMatrixInitialization - - \ref TutorialCoreArithmeticOperators - - \ref TutorialCoreReductions - - \ref TutorialCoreMatrixBlocks - - \ref TutorialCoreDiagonalMatrices - - \ref TutorialCoreTransposeAdjoint - - \ref TutorialCoreDotNorm - - \ref TutorialCoreTriangularMatrix - - \ref TutorialCoreSpecialTopics -\n - -
- -\section TutorialCoreGettingStarted Getting started - -In order to use Eigen, you just need to download and extract Eigen's source code. It is not necessary to use CMake or install anything. - -Here are some quick compilation instructions with GCC. To quickly test an example program, just do - -\code g++ -I /path/to/eigen2/ my_program.cpp -o my_program \endcode - -There is no library to link to. For good performance, add the \c -O2 compile-flag. Note however that this makes it impossible to debug inside Eigen code, as many functions get inlined. In some cases, performance can be further improved by disabling Eigen assertions: use \c -DEIGEN_NO_DEBUG or \c -DNDEBUG to disable them. - -On the x86 architecture, the SSE2 instruction set is not enabled by default. Use \c -msse2 to enable it, and Eigen will then automatically enable its vectorized paths. On x86-64 and AltiVec-based architectures, vectorization is enabled by default. - - -\warning \redstar In most cases it is enough to include the \c Eigen/Core header only to get started with Eigen. However, some features presented in this tutorial require the Array module to be included (\c \#include \c ). Those features are highlighted with a red star \redstar. - -\section TutorialCoreSimpleExampleFixedSize Simple example with fixed-size matrices and vectors - -By fixed-size, we mean that the number of rows and columns are fixed at compile-time. In this case, Eigen avoids dynamic memory allocation, and unroll loops when that makes sense. This is useful for very small sizes: typically up to 4x4, sometimes up to 16x16. - - -
-\include Tutorial_simple_example_fixed_size.cpp - -output: -\include Tutorial_simple_example_fixed_size.out -
- -top\section TutorialCoreSimpleExampleDynamicSize Simple example with dynamic-size matrices and vectors - -By dynamic-size, we mean that the numbers of rows and columns are not fixed at compile-time. In this case, they are stored as runtime variables and the arrays are dynamically allocated. - - -
-\include Tutorial_simple_example_dynamic_size.cpp - -output: -\include Tutorial_simple_example_dynamic_size.out -
- - - - - - -top\section TutorialCoreMatrixTypes Matrix and vector types - -In Eigen, all kinds of dense matrices and vectors are represented by the template class Matrix. In most cases, you can simply use one of the \ref matrixtypedefs "convenience typedefs". - -The template class Matrix takes a number of template parameters, but for now it is enough to understand the 3 first ones (and the others can then be left unspecified): - -\code Matrix \endcode - -\li \c Scalar is the scalar type, i.e. the type of the coefficients. That is, if you want a vector of floats, choose \c float here. -\li \c RowsAtCompileTime and \c ColsAtCompileTime are the number of rows and columns of the matrix as known at compile-time. - -For example, \c Vector3d is a typedef for \code Matrix \endcode - -For dynamic-size, that is in order to left the number of rows or of columns unspecified at compile-time, use the special value Eigen::Dynamic. For example, \c VectorXd is a typedef for \code Matrix \endcode - - -top\section TutorialCoreCoefficients Coefficient access - -Eigen supports the following syntaxes for read and write coefficient access: - -\code -matrix(i,j); -vector(i) -vector[i] -vector.x() // first coefficient -vector.y() // second coefficient -vector.z() // third coefficient -vector.w() // fourth coefficient -\endcode - -Notice that these coefficient access methods have assertions checking the ranges. So if you do a lot of coefficient access, these assertion can have an important cost. There are then two possibilities if you want avoid paying this cost: -\li Either you can disable assertions altogether, by defining EIGEN_NO_DEBUG or NDEBUG. Notice that some IDEs like MS Visual Studio define NDEBUG automatically in "Release Mode". -\li Or you can disable the checks on a case-by-case basis by using the coeff() and coeffRef() methods: see MatrixBase::coeff(int,int) const, MatrixBase::coeffRef(int,int), etc. - -top\section TutorialCoreMatrixInitialization Matrix and vector creation and initialization - - -\subsection TutorialPredefMat Predefined Matrices -Eigen offers several static methods to create special matrix expressions, and non-static methods to assign these expressions to existing matrices: - - - - - - - - - - - - - - -
Fixed-size matrix or vectorDynamic-size matrixDynamic-size vector
-\code -Matrix3f x; - -x = Matrix3f::Zero(); -x = Matrix3f::Ones(); -x = Matrix3f::Constant(value); -x = Matrix3f::Identity(); -x = Matrix3f::Random(); - -x.setZero(); -x.setOnes(); -x.setIdentity(); -x.setConstant(value); -x.setRandom(); -\endcode - -\code -MatrixXf x; - -x = MatrixXf::Zero(rows, cols); -x = MatrixXf::Ones(rows, cols); -x = MatrixXf::Constant(rows, cols, value); -x = MatrixXf::Identity(rows, cols); -x = MatrixXf::Random(rows, cols); - -x.setZero(rows, cols); -x.setOnes(rows, cols); -x.setConstant(rows, cols, value); -x.setIdentity(rows, cols); -x.setRandom(rows, cols); -\endcode - -\code -VectorXf x; - -x = VectorXf::Zero(size); -x = VectorXf::Ones(size); -x = VectorXf::Constant(size, value); -x = VectorXf::Identity(size); -x = VectorXf::Random(size); - -x.setZero(size); -x.setOnes(size); -x.setConstant(size, value); -x.setIdentity(size); -x.setRandom(size); -\endcode -
\redstar the Random() and setRandom() functions require the inclusion of the Array module (\c \#include \c )
Basis vectors \link MatrixBase::Unit [details]\endlink
\code -Vector3f::UnitX() // 1 0 0 -Vector3f::UnitY() // 0 1 0 -Vector3f::UnitZ() // 0 0 1 -\endcode\code -VectorXf::Unit(size,i) -VectorXf::Unit(4,1) == Vector4f(0,1,0,0) - == Vector4f::UnitY() -\endcode -
- -Here is an usage example: - -
-\code -cout << MatrixXf::Constant(2, 3, sqrt(2)) << endl; -RowVector3i v; -v.setConstant(6); -cout << "v = " << v << endl; -\endcode - -output: -\code -1.41 1.41 1.41 -1.41 1.41 1.41 - -v = 6 6 6 -\endcode -
- - -\subsection TutorialCasting Casting - -In Eigen, any matrices of same size and same scalar type are all naturally compatible. The scalar type can be explicitely casted to another one using the template MatrixBase::cast() function: -\code -Matrix3d md(1,2,3); -Matrix3f mf = md.cast(); -\endcode -Note that casting to the same scalar type in an expression is free. - -The destination matrix is automatically resized in any assignment: -\code -MatrixXf res(10,10); -Matrix3f a, b; -res = a+b; // OK: res is resized to size 3x3 -\endcode -Of course, fixed-size matrices can't be resized. - - -\subsection TutorialMap Map -Any memory buffer can be mapped as an Eigen expression: -
-\code -std::vector stlarray(10); -Map(&stlarray[0], stlarray.size()).setOnes(); -int data[4] = 1, 2, 3, 4; -Matrix2i mat2x2(data); -MatrixXi mat2x2 = Map(data); -MatrixXi mat2x2 = Map(data,2,2); -\endcode -
- - - -\subsection TutorialCommaInit Comma initializer -Eigen also offers a \ref MatrixBaseCommaInitRef "comma initializer syntax" which allows you to set all the coefficients of a matrix to specific values: - -
-\include Tutorial_commainit_01.cpp - -output: -\verbinclude Tutorial_commainit_01.out -
- -Not excited by the above example? Then look at the following one where the matrix is set by blocks: - -
-\include Tutorial_commainit_02.cpp - -output: -\verbinclude Tutorial_commainit_02.out -
- -\b Side \b note: here \link CommaInitializer::finished() .finished() \endlink -is used to get the actual matrix object once the comma initialization -of our temporary submatrix is done. Note that despite the apparent complexity of such an expression, -Eigen's comma initializer usually compiles to very optimized code without any overhead. - - - - - - -top\section TutorialCoreArithmeticOperators Arithmetic Operators - -In short, all arithmetic operators can be used right away as in the following example. Note however that arithmetic operators are only given their usual meaning from mathematics tradition. For other operations, such as taking the coefficient-wise product of two vectors, see the discussion of \link Cwise .cwise() \endlink below. Anyway, here is an example demonstrating basic arithmetic operators: -\code -mat4 -= mat1*1.5 + mat2 * (mat3/4); -\endcode -which includes two matrix scalar products ("mat1*1.5" and "mat3/4"), a matrix-matrix product ("mat2 * (mat3/4)"), -a matrix addition ("+") and substraction with assignment ("-="). - - - - - -
-matrix/vector product\code -col2 = mat1 * col1; -row2 = row1 * mat1; row1 *= mat1; -mat3 = mat1 * mat2; mat3 *= mat1; \endcode -
-add/subtract\code -mat3 = mat1 + mat2; mat3 += mat1; -mat3 = mat1 - mat2; mat3 -= mat1;\endcode -
-scalar product\code -mat3 = mat1 * s1; mat3 = s1 * mat1; mat3 *= s1; -mat3 = mat1 / s1; mat3 /= s1;\endcode -
- -In Eigen, only traditional mathematical operators can be used right away. -But don't worry, thanks to the \link Cwise .cwise() \endlink operator prefix, -Eigen's matrices are also very powerful as a numerical container supporting -most common coefficient-wise operators. - - -
- - - - - - -
Coefficient wise \link Cwise::operator*() product \endlink\code mat3 = mat1.cwise() * mat2; \endcode -
-Add a scalar to all coefficients \redstar\code -mat3 = mat1.cwise() + scalar; -mat3.cwise() += scalar; -mat3.cwise() -= scalar; -\endcode -
-Coefficient wise \link Cwise::operator/() division \endlink \redstar\code -mat3 = mat1.cwise() / mat2; \endcode -
-Coefficient wise \link Cwise::inverse() reciprocal \endlink \redstar\code -mat3 = mat1.cwise().inverse(); \endcode -
-Coefficient wise comparisons \redstar \n -(support all operators)\code -mat3 = mat1.cwise() < mat2; -mat3 = mat1.cwise() <= mat2; -mat3 = mat1.cwise() > mat2; -etc. -\endcode -
-		 - - - -
-\b Trigo \redstar: \n -\link Cwise::sin sin \endlink, \link Cwise::cos cos \endlink\code -mat3 = mat1.cwise().sin(); -etc. -\endcode -
-\b Power \redstar: \n \link Cwise::pow() pow \endlink, -\link Cwise::square square \endlink, -\link Cwise::cube cube \endlink, \n -\link Cwise::sqrt sqrt \endlink, -\link Cwise::exp exp \endlink, -\link Cwise::log log \endlink \code -mat3 = mat1.cwise().square(); -mat3 = mat1.cwise().pow(5); -mat3 = mat1.cwise().log(); -etc. -\endcode -
-\link Cwise::min min \endlink, \link Cwise::max max \endlink, \n -absolute value (\link Cwise::abs() abs \endlink, \link Cwise::abs2() abs2 \endlink) -\code -mat3 = mat1.cwise().min(mat2); -mat3 = mat1.cwise().max(mat2); -mat3 = mat1.cwise().abs(); -mat3 = mat1.cwise().abs2(); -\endcode
-
-\redstar Those functions require the inclusion of the Array module (\c \#include \c ). - -\b Side \b note: If you think that the \c .cwise() syntax is too verbose for your own taste and prefer to have non-conventional mathematical operators directly available, then feel free to extend MatrixBase as described \ref ExtendingMatrixBase "here". - -So far, we saw the notation \code mat1*mat2 \endcode for matrix product, and \code mat1.cwise()*mat2 \endcode for coefficient-wise product. What about other kinds of products, which in some other libraries also use arithmetic operators? In Eigen, they are accessed as follows -- note that here we are anticipating on further sections, for convenience. - - - - -
\link MatrixBase::dot() dot product \endlink (inner product)\code -scalar = vec1.dot(vec2);\endcode -
-outer product\code -mat = vec1 * vec2.transpose();\endcode -
-\link MatrixBase::cross() cross product \endlink\code -#include -vec3 = vec1.cross(vec2);\endcode
- - - - -top\section TutorialCoreReductions Reductions - -Eigen provides several reduction methods such as: -\link MatrixBase::minCoeff() minCoeff() \endlink, \link MatrixBase::maxCoeff() maxCoeff() \endlink, -\link MatrixBase::sum() sum() \endlink, \link MatrixBase::trace() trace() \endlink, -\link MatrixBase::norm() norm() \endlink, \link MatrixBase::squaredNorm() squaredNorm() \endlink, -\link MatrixBase::all() all() \endlink \redstar,and \link MatrixBase::any() any() \endlink \redstar. -All reduction operations can be done matrix-wise, -\link MatrixBase::colwise() column-wise \endlink \redstar or -\link MatrixBase::rowwise() row-wise \endlink \redstar. Usage example: - - - - -
\code - 5 3 1 -mat = 2 7 8 - 9 4 6 \endcode - \code mat.minCoeff(); \endcode\code 1 \endcode
\code mat.colwise().minCoeff(); \endcode\code 2 3 1 \endcode
\code mat.rowwise().minCoeff(); \endcode\code -1 -2 -4 -\endcode
- -Also note that maxCoeff and minCoeff can takes optional arguments returning the coordinates of the respective min/max coeff: \link MatrixBase::maxCoeff(int*,int*) const maxCoeff(int* i, int* j) \endlink, \link MatrixBase::minCoeff(int*,int*) const minCoeff(int* i, int* j) \endlink. - -\b Side \b note: The all() and any() functions are especially useful in combinaison with coeff-wise comparison operators (\ref CwiseAll "example"). - - - - - -top\section TutorialCoreMatrixBlocks Matrix blocks - -Read-write access to a \link MatrixBase::col(int) column \endlink -or a \link MatrixBase::row(int) row \endlink of a matrix: -\code -mat1.row(i) = mat2.col(j); -mat1.col(j1).swap(mat1.col(j2)); -\endcode - -Read-write access to sub-vectors: - - - - - - - - - - - - - - - - - - - - -
Default versionsOptimized versions when the size \n is known at compile time
\code vec1.start(n)\endcode\code vec1.start()\endcodethe first \c n coeffs
\code vec1.end(n)\endcode\code vec1.end()\endcodethe last \c n coeffs
\code vec1.segment(pos,n)\endcode\code vec1.segment(pos)\endcodethe \c size coeffs in \n the range [\c pos : \c pos + \c n [
- -Read-write access to sub-matrices:
\code mat1.block(i,j,rows,cols)\endcode - \link MatrixBase::block(int,int,int,int) (more) \endlink\code mat1.block(i,j)\endcode - \link MatrixBase::block(int,int) (more) \endlinkthe \c rows x \c cols sub-matrix \n starting from position (\c i,\c j)
\code - mat1.corner(TopLeft,rows,cols) - mat1.corner(TopRight,rows,cols) - mat1.corner(BottomLeft,rows,cols) - mat1.corner(BottomRight,rows,cols)\endcode - \link MatrixBase::corner(CornerType,int,int) (more) \endlink\code - mat1.corner(TopLeft) - mat1.corner(TopRight) - mat1.corner(BottomLeft) - mat1.corner(BottomRight)\endcode - \link MatrixBase::corner(CornerType) (more) \endlinkthe \c rows x \c cols sub-matrix \n taken in one of the four corners
\code -mat4x4.minor(i,j) = mat3x3; -mat3x3 = mat4x4.minor(i,j);\endcode - -\link MatrixBase::minor() minor \endlink (read-write)
- - - -top\section TutorialCoreDiagonalMatrices Diagonal matrices - - - - - - -
-\link MatrixBase::asDiagonal() make a diagonal matrix \endlink from a vector \n -this product is automatically optimized !\code -mat3 = mat1 * vec2.asDiagonal();\endcode -
Access \link MatrixBase::diagonal() the diagonal of a matrix \endlink as a vector (read/write)\code - vec1 = mat1.diagonal(); - mat1.diagonal() = vec1; - \endcode -
- -top\section TutorialCoreTransposeAdjoint Transpose and Adjoint operations - - - - -
-\link MatrixBase::transpose() transposition \endlink (read-write)\code -mat3 = mat1.transpose() * mat2; -mat3.transpose() = mat1 * mat2.transpose(); -\endcode -
-\link MatrixBase::adjoint() adjoint \endlink (read only)\n\code -mat3 = mat1.adjoint() * mat2; -\endcode -
- -top\section TutorialCoreDotNorm Dot-product, vector norm, normalization - - - - - -
-\link MatrixBase::dot() Dot-product \endlink of two vectors -\code vec1.dot(vec2);\endcode -
-\link MatrixBase::norm() norm \endlink of a vector \n -\link MatrixBase::squaredNorm() squared norm \endlink of a vector -\code vec.norm(); \endcode \n \code vec.squaredNorm() \endcode -
-returns a \link MatrixBase::normalized() normalized \endlink vector \n -\link MatrixBase::normalize() normalize \endlink a vector -\code -vec3 = vec1.normalized(); -vec1.normalize();\endcode -
- -top\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.. - - - - - - - -
-Extract triangular matrices \n from a given matrix m: -\code -m.part() -m.part() -m.part() -m.part() -m.part() -m.part()\endcode -
-Write to triangular parts \n of a matrix m: -\code -m1.part() = m2; -m1.part() = m2; -m1.part() = m2; -m1.part() = m2;\endcode -
-Special: take advantage of symmetry \n (selfadjointness) when copying \n an expression into a matrix -\code -m.part() = someSelfadjointMatrix; -m1.part() = m2 + m2.adjoint(); // m2 + m2.adjoint() is selfadjoint \endcode -
- - -top\section TutorialCoreSpecialTopics Special Topics - -\ref TopicLazyEvaluation "Lazy Evaluation and Aliasing": Thanks to expression templates, Eigen is able to apply lazy evaluation wherever that is beneficial. - -*/ - -} diff --git a/doc/StlContainers.dox b/doc/StlContainers.dox deleted file mode 100644 index f51ed6fd2..000000000 --- a/doc/StlContainers.dox +++ /dev/null @@ -1,52 +0,0 @@ -namespace Eigen { - -/** \page StlContainers Using STL Containers with Eigen - -\b Table \b of \b contents - - \ref summary - - \ref allocator - - \ref vector - -\section summary Executive summary - -Using STL containers on \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types" requires taking the following two steps: - -\li A 16-byte-aligned allocator must be used. Eigen does provide one ready for use: aligned_allocator. -\li If you want to use the std::vector container, you need to \#include . - -These issues arise only with \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types". For other Eigen types, such as Vector3f or MatrixXd, no special care is needed when using STL containers. - -\section allocator Using an aligned allocator - -STL containers take an optional template parameter, the allocator type. When using STL containers on \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types", you need tell the container to use an allocator that will always allocate memory at 16-byte-aligned locations. Fortunately, Eigen does provide such an allocator: Eigen::aligned_allocator. - -For example, instead of -\code -std::map -\endcode -you need to use -\code -std::map, Eigen::aligned_allocator > -\endcode -Note that here, the 3rd parameter "std::less" is just the default value, we only had to specify it because we needed to specify the allocator type, that is the 4th parameter. - -\section vector The case of std::vector - -The situation with std::vector was even worse (explanation below) so we had to specialize it for Eigen types. The upside is that our specialization takes care of specifying the aligned allocator, so you don't need to worry about it. All you need to do is to \#include . - -So as soon as you have -\code -#include -\endcode -you can simply use -\code -std::vector -\endcode -without having to worry about anything. - -\b Explanation: The resize() method of std::vector takes a value_type argument (defaulting to value_type()). So with std::vector, some Eigen::Vector4f objects will be passed by value, which discards any alignment modifiers, so a Eigen::Vector4f can be created at an unaligned location. In order to avoid that, the only solution we saw was to specialize std::vector to make it work on a slight modification of, here, Eigen::Vector4f, that is able to deal properly with this situation. - - -*/ - -} diff --git a/doc/StructHavingEigenMembers.dox b/doc/StructHavingEigenMembers.dox deleted file mode 100644 index afa197900..000000000 --- a/doc/StructHavingEigenMembers.dox +++ /dev/null @@ -1,149 +0,0 @@ -namespace Eigen { - -/** \page StructHavingEigenMembers Structures Having Eigen Members - -\b Table \b of \b contents - - \ref summary - - \ref what - - \ref how - - \ref why - - \ref movetotop - - \ref bugineigen - - \ref conditional - -\section summary Executive Summary - -If you define a structure having members of \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types", you must overload its "operator new" so that it generates 16-bytes-aligned pointers. Fortunately, Eigen provides you with a macro EIGEN_MAKE_ALIGNED_OPERATOR_NEW that does that for you. - -\section what What kind of code needs to be changed? - -The kind of code that needs to be changed is this: - -\code -class Foo -{ - ... - Eigen::Vector2d v; - ... -}; - -... - -Foo *foo = new Foo; -\endcode - -In other words: you have a class that has as a member a \ref FixedSizeVectorizable "fixed-size vectorizable Eigen object", and then you dynamically create an object of that class. - -\section how How should such code be modified? - -Very easy, you just need to put a EIGEN_MAKE_ALIGNED_OPERATOR_NEW macro in a public part of your class, like this: - -\code -class Foo -{ - ... - Eigen::Vector2d v; - ... -public: - EIGEN_MAKE_ALIGNED_OPERATOR_NEW -}; - -... - -Foo *foo = new Foo; -\endcode - -This macro makes "new Foo" always return an aligned pointer. - -\section why Why is this needed? - -OK let's say that your code looks like this: - -\code -class Foo -{ - ... - Eigen::Vector2d v; - ... -}; - -... - -Foo *foo = new Foo; -\endcode - -A Eigen::Vector2d consists of 2 doubles, which is 128 bits. Which is exactly the size of a SSE packet, which makes it possible to use SSE for all sorts of operations on this vector. But SSE instructions (at least the ones that Eigen uses, which are the fast ones) require 128-bit alignment. Otherwise you get a segmentation fault. - -For this reason, Eigen takes care by itself to require 128-bit alignment for Eigen::Vector2d, by doing two things: -\li Eigen requires 128-bit alignment for the Eigen::Vector2d's array (of 2 doubles). With GCC, this is done with a __attribute__ ((aligned(16))). -\li Eigen overloads the "operator new" of Eigen::Vector2d so it will always return 128-bit aligned pointers. - -Thus, normally, you don't have to worry about anything, Eigen handles alignment for you... - -... except in one case. When you have a class Foo like above, and you dynamically allocate a new Foo as above, then, since Foo doesn't have aligned "operator new", the returned pointer foo is not necessarily 128-bit aligned. - -The alignment attribute of the member v is then relative to the start of the class, foo. If the foo pointer wasn't aligned, then foo->v won't be aligned either! - -The solution is to let class Foo have an aligned "operator new", as we showed in the previous section. - -\section movetotop Should I then put all the members of Eigen types at the beginning of my class? - -No, that's not needed. Since Eigen takes care of declaring 128-bit alignment, all members that need it are automatically 128-bit aligned relatively to the class. So when you have code like - -\code -class Foo -{ - double x; - Eigen::Vector2d v; -public: - EIGEN_MAKE_ALIGNED_OPERATOR_NEW -}; -\endcode - -it will work just fine. You do \b not need to rewrite it as - -\code -class Foo -{ - Eigen::Vector2d v; - double x; -public: - EIGEN_MAKE_ALIGNED_OPERATOR_NEW -}; -\endcode - -\section dynamicsize What about dynamic-size matrices and vectors? - -Dynamic-size matrices and vectors, such as Eigen::VectorXd, allocate dynamically their own array of coefficients, so they take care of requiring absolute alignment automatically. So they don't cause this issue. The issue discussed here is only with fixed-size matrices and vectors. - -\section bugineigen So is this a bug in Eigen? - -No, it's not our bug. It's more like an inherent problem of the C++ language -- though it must be said that any other existing language probably has the same problem. The problem is that there is no way that you can specify an aligned "operator new" that would propagate to classes having you as member data. - -\section conditional What if I want to do this conditionnally (depending on template parameters) ? - -For this situation, we offer the macro EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign). It will generate aligned operators like EIGEN_MAKE_ALIGNED_OPERATOR_NEW if NeedsToAlign is true. It will generate operators with the default alignment if NeedsToAlign is false. - -Example: - -\code -template class Foo -{ - typedef Eigen::Matrix Vector; - enum { NeedsToAlign = (sizeof(Vector)%16)==0 }; - ... - Vector v; - ... -public: - EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign) -}; - -... - -Foo<4> *foo4 = new Foo<4>; // foo4 is guaranteed to be 128bit-aligned -Foo<3> *foo3 = new Foo<3>; // foo3 has only the system default alignment guarantee -\endcode - -*/ - -} diff --git a/doc/TopicLazyEvaluation.dox b/doc/TopicLazyEvaluation.dox deleted file mode 100644 index 7df9824ba..000000000 --- a/doc/TopicLazyEvaluation.dox +++ /dev/null @@ -1,67 +0,0 @@ -namespace Eigen { - -/** \page TopicLazyEvaluation Lazy Evaluation and Aliasing - -Executive summary: Eigen has intelligent compile-time mechanisms to enable lazy evaluation and removing temporaries where appropriate. -It will handle aliasing automatically in most cases, for example with matrix products. The automatic behavior can be overridden -manually by using the MatrixBase::eval() and MatrixBase::lazy() methods. - -When you write a line of code involving a complex expression such as - -\code mat1 = mat2 + mat3 * (mat4 + mat5); \endcode - -Eigen determines automatically, for each sub-expression, whether to evaluate it into a temporary variable. Indeed, in certain cases it is better to evaluate immediately a sub-expression into a temporary variable, while in other cases it is better to avoid that. - -A traditional math library without expression templates always evaluates all sub-expressions into temporaries. So with this code, - -\code vec1 = vec2 + vec3; \endcode - -a traditional library would evaluate \c vec2 + vec3 into a temporary \c vec4 and then copy \c vec4 into \c vec1. This is of course inefficient: the arrays are traversed twice, so there are a lot of useless load/store operations. - -Expression-templates-based libraries can avoid evaluating sub-expressions into temporaries, which in many cases results in large speed improvements. This is called lazy evaluation as an expression is getting evaluated as late as possible, instead of immediately. However, most other expression-templates-based libraries always choose lazy evaluation. There are two problems with that: first, lazy evaluation is not always a good choice for performance; second, lazy evaluation can be very dangerous, for example with matrix products: doing matrix = matrix*matrix gives a wrong result if the matrix product is lazy-evaluated, because of the way matrix product works. - -For these reasons, Eigen has intelligent compile-time mechanisms to determine automatically when to use lazy evaluation, and when on the contrary it should evaluate immediately into a temporary variable. - -So in the basic example, - -\code matrix1 = matrix2 + matrix3; \endcode - -Eigen chooses lazy evaluation. Thus the arrays are traversed only once, producing optimized code. If you really want to force immediate evaluation, use \link MatrixBase::eval() eval()\endlink: - -\code matrix1 = (matrix2 + matrix3).eval(); \endcode - -Here is now a more involved example: - -\code matrix1 = -matrix2 + matrix3 + 5 * matrix4; \endcode - -Eigen chooses lazy evaluation at every stage in that example, which is clearly the correct choice. In fact, lazy evaluation is the "default choice" and Eigen will choose it except in a few circumstances. - -The first circumstance in which Eigen chooses immediate evaluation, is when it sees an assignment a = b; and the expression \c b has the evaluate-before-assigning \link flags flag\endlink. The most important example of such an expression is the \link Product matrix product expression\endlink. For example, when you do - -\code matrix = matrix * matrix; \endcode - -Eigen first evaluates matrix * matrix into a temporary matrix, and then copies it into the original \c matrix. This guarantees a correct result as we saw above that lazy evaluation gives wrong results with matrix products. It also doesn't cost much, as the cost of the matrix product itself is much higher. - -What if you know what you are doing and want to force lazy evaluation? Then use \link MatrixBase::lazy() .lazy()\endlink instead. Here is an example: - -\code matrix1 = (matrix2 * matrix2).lazy(); \endcode - -Here, since we know that matrix2 is not the same matrix as matrix1, we know that lazy evaluation is not dangerous, so we may force lazy evaluation. Concretely, the effect of lazy() here is to remove the evaluate-before-assigning \link flags flag\endlink and also the evaluate-before-nesting \link flags flag\endlink which we now discuss. - -The second circumstance in which Eigen chooses immediate evaluation, is when it sees a nested expression such as a + b where \c b is already an expression having the evaluate-before-nesting \link flags flag\endlink. Again, the most important example of such an expression is the \link Product matrix product expression\endlink. For example, when you do - -\code matrix1 = matrix2 + matrix3 * matrix4; \endcode - -the product matrix3 * matrix4 gets evaluated immediately into a temporary matrix. Indeed, experiments showed that it is often beneficial for performance to evaluate immediately matrix products when they are nested into bigger expressions. - -Again, \link MatrixBase::lazy() .lazy()\endlink can be used to force lazy evaluation here. - -The third circumstance in which Eigen chooses immediate evaluation, is when its cost model shows that the total cost of an operation is reduced if a sub-expression gets evaluated into a temporary. Indeed, in certain cases, an intermediate result is sufficiently costly to compute and is reused sufficiently many times, that is worth "caching". Here is an example: - -\code matrix1 = matrix2 * (matrix3 + matrix4); \endcode - -Here, provided the matrices have at least 2 rows and 2 columns, each coefficienct of the expression matrix3 + matrix4 is going to be used several times in the matrix product. Instead of computing the sum everytime, it is much better to compute it once and store it in a temporary variable. Eigen understands this and evaluates matrix3 + matrix4 into a temporary variable before evaluating the product. - -*/ - -} \ No newline at end of file diff --git a/doc/TutorialGeometry.dox b/doc/TutorialGeometry.dox deleted file mode 100644 index b3c860d53..000000000 --- a/doc/TutorialGeometry.dox +++ /dev/null @@ -1,253 +0,0 @@ -namespace Eigen { - -/** \page TutorialGeometry Tutorial 2/3 - Geometry - \ingroup Tutorial - -
\ref index "Overview" - | \ref TutorialCore "Core features" - | \b Geometry - | \ref TutorialAdvancedLinearAlgebra "Advanced linear algebra" - | \ref TutorialSparse "Sparse matrix" -
- -In this tutorial chapter we will shortly introduce the many possibilities offered by the \ref GeometryModule "geometry module", -namely 2D and 3D rotations and affine transformations. - -\b Table \b of \b contents - - \ref TutorialGeoElementaryTransformations - - \ref TutorialGeoCommontransformationAPI - - \ref TutorialGeoTransform - - \ref TutorialGeoEulerAngles - -Eigen's Geometry module provides two different kinds of geometric transformations: - - Abstract transformations, such as rotations (represented by \ref AngleAxis "angle and axis" or by a \ref Quaternion "quaternion"), \ref Translation "translations", \ref Scaling "scalings". These transformations are NOT represented as matrices, but you can nevertheless mix them with matrices and vectors in expressions, and convert them to matrices if you wish. - - Affine transformation matrices: see the Transform class. These are really matrices. - -\note If you are working with OpenGL 4x4 matrices then Transform3f and Transform3d are what you want. Since Eigen defaults to column-major storage, you can directly use the Transform::data() method to pass your transformation matrix to OpenGL. - -You can construct a Transform from an abstract transformation, like this: -\code - Transform t(AngleAxis(angle,axis)); -\endcode -or like this: -\code - Transform t; - t = AngleAxis(angle,axis); -\endcode -But note that unfortunately, because of how C++ works, you can \b not do this: -\code - Transform t = AngleAxis(angle,axis); -\endcode -\b Explanation: In the C++ language, this would require Transform to have a non-explicit conversion constructor from AngleAxis, but we really don't want to allow implicit casting here. - - -\section TutorialGeoElementaryTransformations Transformation types - - - - - - - - - - -
Transformation typeTypical initialization code
-\ref Rotation2D "2D rotation" from an angle\code -Rotation2D rot2(angle_in_radian);\endcode
-3D rotation as an \ref AngleAxis "angle + axis"\code -AngleAxis aa(angle_in_radian, Vector3f(ax,ay,az));\endcode
-3D rotation as a \ref Quaternion "quaternion"\code -Quaternion q = AngleAxis(angle_in_radian, axis);\endcode
-N-D Scaling\code -Scaling(sx, sy) -Scaling(sx, sy, sz) -Scaling(s) -Scaling(vecN)\endcode
-N-D Translation\code -Translation(tx, ty) -Translation(tx, ty, tz) -Translation(s) -Translation(vecN)\endcode
-N-D \ref TutorialGeoTransform "Affine transformation"\code -Transform t = concatenation_of_any_transformations; -Transform t = Translation3f(p) * AngleAxisf(a,axis) * Scaling3f(s);\endcode
-N-D Linear transformations \n -(pure rotations, \n scaling, etc.)\code -Matrix t = concatenation_of_rotations_and_scalings; -Matrix t = Rotation2Df(a) * Scaling2f(s); -Matrix t = AngleAxisf(a,axis) * Scaling3f(s);\endcode
- -Notes on rotations\n To transform more than a single vector the preferred -representations are rotation matrices, while for other usages Quaternion is the -representation of choice as they are compact, fast and stable. Finally Rotation2D and -AngleAxis are mainly convenient types to create other rotation objects. - -Notes on Translation and Scaling\n Likewise AngleAxis, these classes were -designed to simplify the creation/initialization of linear (Matrix) and affine (Transform) -transformations. Nevertheless, unlike AngleAxis which is inefficient to use, these classes -might still be interesting to write generic and efficient algorithms taking as input any -kind of transformations. - -Any of the above transformation types can be converted to any other types of the same nature, -or to a more generic type. Here are some additional examples: - - -
\code -Rotation2Df r = Matrix2f(..); // assumes a pure rotation matrix -AngleAxisf aa = Quaternionf(..); -AngleAxisf aa = Matrix3f(..); // assumes a pure rotation matrix -Matrix2f m = Rotation2Df(..); -Matrix3f m = Quaternionf(..); Matrix3f m = Scaling3f(..); -Transform3f m = AngleAxis3f(..); Transform3f m = Scaling3f(..); -Transform3f m = Translation3f(..); Transform3f m = Matrix3f(..); -\endcode
- - -top\section TutorialGeoCommontransformationAPI Common API across transformation types - -To some extent, Eigen's \ref GeometryModule "geometry module" allows you to write -generic algorithms working on any kind of transformation representations: - - - - - -
-Concatenation of two transformations\code -gen1 * gen2;\endcode
Apply the transformation to a vector\code -vec2 = gen1 * vec1;\endcode
Get the inverse of the transformation\code -gen2 = gen1.inverse();\endcode
Spherical interpolation \n (Rotation2D and Quaternion only)\code -rot3 = rot1.slerp(alpha,rot2);\endcode
- - - -top\section TutorialGeoTransform Affine transformations -Generic affine transformations are represented by the Transform class which internaly -is a (Dim+1)^2 matrix. In Eigen we have chosen to not distinghish between points and -vectors such that all points are actually represented by displacement vectors from the -origin ( \f$ \mathbf{p} \equiv \mathbf{p}-0 \f$ ). With that in mind, real points and -vector distinguish when the transformation is applied. - - - - - - - -
-Apply the transformation to a \b point \code -VectorNf p1, p2; -p2 = t * p1;\endcode
-Apply the transformation to a \b vector \code -VectorNf vec1, vec2; -vec2 = t.linear() * vec1;\endcode
-Apply a \em general transformation \n to a \b normal \b vector -(explanations)\code -VectorNf n1, n2; -MatrixNf normalMatrix = t.linear().inverse().transpose(); -n2 = (normalMatrix * n1).normalized();\endcode
-Apply a transformation with \em pure \em rotation \n to a \b normal \b vector -(no scaling, no shear)\code -n2 = t.linear() * n1;\endcode
-OpenGL compatibility \b 3D \code -glLoadMatrixf(t.data());\endcode
-OpenGL compatibility \b 2D \code -Transform3f aux(Transform3f::Identity); -aux.linear().corner<2,2>(TopLeft) = t.linear(); -aux.translation().start<2>() = t.translation(); -glLoadMatrixf(aux.data());\endcode
- -\b Component \b accessors - - - - - - -
-full read-write access to the internal matrix\code -t.matrix() = matN1xN1; // N1 means N+1 -matN1xN1 = t.matrix(); -\endcode
-coefficient accessors\code -t(i,j) = scalar; <=> t.matrix()(i,j) = scalar; -scalar = t(i,j); <=> scalar = t.matrix()(i,j); -\endcode
-translation part\code -t.translation() = vecN; -vecN = t.translation(); -\endcode
-linear part\code -t.linear() = matNxN; -matNxN = t.linear(); -\endcode
-extract the rotation matrix\code -matNxN = t.extractRotation(); -\endcode
- - -\b Transformation \b creation \n -While transformation objects can be created and updated concatenating elementary transformations, -the Transform class also features a procedural API: - - - - - - -
\b procedurale \b API \b equivalent \b natural \b API
Translation\code -t.translate(Vector_(tx,ty,..)); -t.pretranslate(Vector_(tx,ty,..)); -\endcode\code -t *= Translation_(tx,ty,..); -t = Translation_(tx,ty,..) * t; -\endcode
\b Rotation \n In 2D and for the procedural API, any_rotation can also \n be an angle in radian\code -t.rotate(any_rotation); -t.prerotate(any_rotation); -\endcode\code -t *= any_rotation; -t = any_rotation * t; -\endcode
Scaling\code -t.scale(Vector_(sx,sy,..)); -t.scale(s); -t.prescale(Vector_(sx,sy,..)); -t.prescale(s); -\endcode\code -t *= Scaling_(sx,sy,..); -t *= Scaling_(s); -t = Scaling_(sx,sy,..) * t; -t = Scaling_(s) * t; -\endcode
Shear transformation \n ( \b 2D \b only ! )\code -t.shear(sx,sy); -t.preshear(sx,sy); -\endcode
- -Note that in both API, any many transformations can be concatenated in a single expression as shown in the two following equivalent examples: - - - -
\code -t.pretranslate(..).rotate(..).translate(..).scale(..); -\endcode
\code -t = Translation_(..) * t * RotationType(..) * Translation_(..) * Scaling_(..); -\endcode
- - - -top\section TutorialGeoEulerAngles Euler angles - - -
-Euler angles might be convenient to create rotation objects. -On the other hand, since there exist 24 differents convension,they are pretty confusing to use. This example shows how -to create a rotation matrix according to the 2-1-2 convention.\code -Matrix3f m; -m = AngleAxisf(angle1, Vector3f::UnitZ()) -* * AngleAxisf(angle2, Vector3f::UnitY()) -* * AngleAxisf(angle3, Vector3f::UnitZ()); -\endcode
- -*/ - -} diff --git a/doc/TutorialLinearAlgebra.dox b/doc/TutorialLinearAlgebra.dox deleted file mode 100644 index 1b584e103..000000000 --- a/doc/TutorialLinearAlgebra.dox +++ /dev/null @@ -1,188 +0,0 @@ - -namespace Eigen { - -/** \page TutorialAdvancedLinearAlgebra Tutorial 3/3 - Advanced linear algebra - \ingroup Tutorial - -
\ref index "Overview" - | \ref TutorialCore "Core features" - | \ref TutorialGeometry "Geometry" - | \b Advanced \b linear \b algebra - | \ref TutorialSparse "Sparse matrix" -
- -\b Table \b of \b contents - - \ref TutorialAdvSolvers - - \ref TutorialAdvLU - - \ref TutorialAdvCholesky - - \ref TutorialAdvQR - - \ref TutorialAdvEigenProblems - -\section TutorialAdvSolvers Solving linear problems - -This part of the tutorial focuses on solving linear problem of the form \f$ A \mathbf{x} = b \f$, -where both \f$ A \f$ and \f$ b \f$ are known, and \f$ x \f$ is the unknown. Moreover, \f$ A \f$ -assumed to be a squared matrix. Of course, the methods described here can be used whenever an expression -involve the product of an inverse matrix with a vector or another matrix: \f$ A^{-1} B \f$. -Eigen offers various algorithms to this problem, and its choice mainly depends on the nature of -the matrix \f$ A \f$, such as its shape, size and numerical properties. - -\subsection TutorialAdvSolvers_Triangular Triangular solver -If the matrix \f$ A \f$ is triangular (upper or lower) and invertible (the coefficients of the diagonal -are all not zero), then the problem can be solved directly using MatrixBase::solveTriangular(), or better, -MatrixBase::solveTriangularInPlace(). Here is an example: - -
-\include MatrixBase_marked.cpp - -output: -\include MatrixBase_marked.out -
- -See MatrixBase::solveTriangular() for more details. - - -\subsection TutorialAdvSolvers_Inverse Direct inversion (for small matrices) -If the matrix \f$ A \f$ is small (\f$ \leq 4 \f$) and invertible, then a good approach is to directly compute -the inverse of the matrix \f$ A \f$, and then obtain the solution \f$ x \f$ by \f$ \mathbf{x} = A^{-1} b \f$. With Eigen, -this can be implemented like this: - -\code -#include -Matrix4f A = Matrix4f::Random(); -Vector4f b = Vector4f::Random(); -Vector4f x = A.inverse() * b; -\endcode - -Note that the function inverse() is defined in the LU module. -See MatrixBase::inverse() for more details. - - -\subsection TutorialAdvSolvers_Symmetric Cholesky (for positive definite matrices) -If the matrix \f$ A \f$ is \b positive \b definite, then -the best method is to use a Cholesky decomposition. -Such positive definite matrices often arise when solving overdetermined problems in a least square sense (see below). -Eigen offers two different Cholesky decompositions: a \f$ LL^T \f$ decomposition where L is a lower triangular matrix, -and a \f$ LDL^T \f$ decomposition where L is lower triangular with unit diagonal and D is a diagonal matrix. -The latter avoids square roots and is therefore slightly more stable than the former one. -\code -#include -MatrixXf D = MatrixXf::Random(8,4); -MatrixXf A = D.transpose() * D; -VectorXf b = D.transpose() * VectorXf::Random(4); -VectorXf x; -A.llt().solve(b,&x); // using a LLT factorization -A.ldlt().solve(b,&x); // using a LDLT factorization -\endcode -when the value of the right hand side \f$ b \f$ is not needed anymore, then it is faster to use -the \em in \em place API, e.g.: -\code -A.llt().solveInPlace(b); -\endcode -In that case the value of \f$ b \f$ is replaced by the result \f$ x \f$. - -If the linear problem has to solved for various right hand side \f$ b_i \f$, but same matrix \f$ A \f$, -then it is highly recommended to perform the decomposition of \$ A \$ only once, e.g.: -\code -// ... -LLT lltOfA(A); -lltOfA.solveInPlace(b0); -lltOfA.solveInPlace(b1); -// ... -\endcode - -\sa Cholesky_Module, LLT::solve(), LLT::solveInPlace(), LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT. - - -\subsection TutorialAdvSolvers_LU LU decomposition (for most cases) -If the matrix \f$ A \f$ does not fit in any of the previous categories, or if you are unsure about the numerical -stability of your problem, then you can use the LU solver based on a decomposition of the same name : see the section \ref TutorialAdvLU below. Actually, Eigen's LU module does not implement a standard LU decomposition, but rather a so-called LU decomposition -with full pivoting and rank update which has much better numerical stability. -The API of the LU solver is the same than the Cholesky one, except that there is no \em in \em place variant: -\code -#include -MatrixXf A = MatrixXf::Random(20,20); -VectorXf b = VectorXf::Random(20); -VectorXf x; -A.lu().solve(b, &x); -\endcode - -Again, the LU decomposition can be stored to be reused or to perform other kernel operations: -\code -// ... -LU luOfA(A); -luOfA.solve(b, &x); -// ... -\endcode - -See the section \ref TutorialAdvLU below. - -\sa class LU, LU::solve(), LU_Module - - -\subsection TutorialAdvSolvers_SVD SVD solver (for singular matrices and special cases) -Finally, Eigen also offer a solver based on a singular value decomposition (SVD). Again, the API is the -same than with Cholesky or LU: -\code -#include -MatrixXf A = MatrixXf::Random(20,20); -VectorXf b = VectorXf::Random(20); -VectorXf x; -A.svd().solve(b, &x); -SVD svdOfA(A); -svdOfA.solve(b, &x); -\endcode - -\sa class SVD, SVD::solve(), SVD_Module - - - - -top\section TutorialAdvLU LU - -Eigen provides a rank-revealing LU decomposition with full pivoting, which has very good numerical stability. - -You can obtain the LU decomposition of a matrix by calling \link MatrixBase::lu() lu() \endlink, which is the easiest way if you're going to use the LU decomposition only once, as in -\code -#include -MatrixXf A = MatrixXf::Random(20,20); -VectorXf b = VectorXf::Random(20); -VectorXf x; -A.lu().solve(b, &x); -\endcode - -Alternatively, you can construct a named LU decomposition, which allows you to reuse it for more than one operation: -\code -#include -MatrixXf A = MatrixXf::Random(20,20); -Eigen::LUDecomposition lu(A); -cout << "The rank of A is" << lu.rank() << endl; -if(lu.isInvertible()) { - cout << "A is invertible, its inverse is:" << endl << lu.inverse() << endl; -} -else { - cout << "Here's a matrix whose columns form a basis of the kernel a.k.a. nullspace of A:" - << endl << lu.kernel() << endl; -} -\endcode - -\sa LU_Module, LU::solve(), class LU - -top\section TutorialAdvCholesky Cholesky -todo - -\sa Cholesky_Module, LLT::solve(), LLT::solveInPlace(), LDLT::solve(), LDLT::solveInPlace(), class LLT, class LDLT - -top\section TutorialAdvQR QR -todo - -\sa QR_Module, class QR - -top\section TutorialAdvEigenProblems Eigen value problems -todo - -\sa class SelfAdjointEigenSolver, class EigenSolver - -*/ - -} diff --git a/doc/TutorialSparse.dox b/doc/TutorialSparse.dox deleted file mode 100644 index a8bfe006e..000000000 --- a/doc/TutorialSparse.dox +++ /dev/null @@ -1,240 +0,0 @@ -namespace Eigen { - -/** \page TutorialSparse Tutorial - Getting started with the sparse module - \ingroup Tutorial - -
\ref index "Overview" - | \ref TutorialCore "Core features" - | \ref TutorialGeometry "Geometry" - | \ref TutorialAdvancedLinearAlgebra "Advanced linear algebra" - | \b Sparse \b matrix -
- -\b Table \b of \b contents \n - - \ref TutorialSparseIntro - - \ref TutorialSparseFilling - - \ref TutorialSparseFeatureSet - - \ref TutorialSparseDirectSolvers -
- -\section TutorialSparseIntro Sparse matrix representations - -In many applications (e.g., finite element methods) it is common to deal with very large matrices where only a few coefficients are different than zero. Both in term of memory consumption and performance, it is fundamental to use an adequate representation storing only nonzero coefficients. Such a matrix is called a sparse matrix. - -\b Declaring \b sparse \b matrices \b and \b vectors \n -The SparseMatrix class is the main sparse matrix representation of the Eigen's sparse module which offers high performance, low memory usage, and compatibility with most of sparse linear algebra packages. Because of its limited flexibility, we also provide a DynamicSparseMatrix variante taillored for low-level sparse matrix assembly. Both of them can be either row major or column major: - -\code -#include -SparseMatrix > m1(1000,2000); // declare a 1000x2000 col-major compressed sparse matrix of complex -SparseMatrix m2(1000,2000); // declare a 1000x2000 row-major compressed sparse matrix of double -DynamicSparseMatrix > m1(1000,2000); // declare a 1000x2000 col-major dynamic sparse matrix of complex -DynamicSparseMatrix m2(1000,2000); // declare a 1000x2000 row-major dynamic sparse matrix of double -\endcode - -Although a sparse matrix could also be used to represent a sparse vector, for that purpose it is better to use the specialized SparseVector class: -\code -SparseVector > v1(1000); // declare a column sparse vector of complex of size 1000 -SparseVector v2(1000); // declare a row sparse vector of double of size 1000 -\endcode -Note that here the size of a vector denotes its dimension and not the number of nonzero coefficients which is initially zero (like sparse matrices). - - -\b Overview \b of \b the \b internal \b sparse \b storage \n -In order to get the best of the Eigen's sparse objects, it is important to have a rough idea of the way they are internally stored. The SparseMatrix class implements the common and generic Compressed Column/Row Storage scheme. It consists of three compact arrays storing the values with their respective inner coordinates, and pointer indices to the begining of each outer vector. For instance, let \c m be a column-major sparse matrix. Then its nonzero coefficients are sequentially stored in memory in a column-major order (\em values). A second array of integer stores the respective row index of each coefficient (\em inner \em indices). Finally, a third array of integer, having the same length than the number of columns, stores the index in the previous arrays of the first element of each column (\em outer \em indices). - -Here is an example, with the matrix: - - - - - - -
03000
2200017
75010
00000
001408
- -and its internal representation using the Compressed Column Storage format: - - - -
Values: 22735141178
Inner indices: 1202 42 14
-Outer indices:
02456\em 7
- -As you can guess, here the storage order is even more important than with dense matrix. We will therefore often make a clear difference between the \em inner and \em outer dimensions. For instance, it is easy to loop over the coefficients of an \em inner \em vector (e.g., a column of a column-major matrix), but completely inefficient to do the same for an \em outer \em vector (e.g., a row of a col-major matrix). - -The SparseVector class implements the same compressed storage scheme but, of course, without any outer index buffer. - -Since all nonzero coefficients of such a matrix are sequentially stored in memory, random insertion of new nonzeros can be extremely costly. To overcome this limitation, Eigen's sparse module provides a DynamicSparseMatrix class which is basically implemented as an array of SparseVector. In other words, a DynamicSparseMatrix is a SparseMatrix where the values and inner-indices arrays have been splitted into multiple small and resizable arrays. Assuming the number of nonzeros per inner vector is relatively low, this slight modification allow for very fast random insertion at the cost of a slight memory overhead and a lost of compatibility with other sparse libraries used by some of our highlevel solvers. Note that the major memory overhead comes from the extra memory preallocated by each inner vector to avoid an expensive memory reallocation at every insertion. - -To summarize, it is recommanded to use a SparseMatrix whenever this is possible, and reserve the use of DynamicSparseMatrix for matrix assembly purpose when a SparseMatrix is not flexible enough. The respective pro/cons of both representations are summarized in the following table: - - - - - - - - - - - - - - - -
SparseMatrixDynamicSparseMatrix
memory usage*****
sorted insertion******
random insertion \n in sorted inner vector****
sorted insertion \n in random inner vector-***
random insertion-**
coeff wise unary operators******
coeff wise binary operators******
matrix products*****(*)
transpose*****
redux*****
*= scalar*****
Compatibility with highlevel solvers \n (TAUCS, Cholmod, SuperLU, UmfPack)***-
- - -\b Matrix \b and \b vector \b properties \n - -Here mat and vec represents any sparse-matrix and sparse-vector types respectively. - - - - - - - - - - -
Standard \n dimensions\code -mat.rows() -mat.cols()\endcode\code -vec.size() \endcode
Sizes along the \n inner/outer dimensions\code -mat.innerSize() -mat.outerSize()\endcode
Number of non \n zero coefficiens\code -mat.nonZeros() \endcode\code -vec.nonZeros() \endcode
- - -\b Iterating \b over \b the \b nonzero \b coefficients \n - -Iterating over the coefficients of a sparse matrix can be done only in the same order than the storage order. Here is an example: - - -
-\code -SparseMatrixType mat(rows,cols); -for (int k=0; k\ -\code -SparseVector vec(size); -for (SparseVector::InnerIterator it(vec); it; ++it) -{ - it.value(); // == vec[ it.index() ] - it.index(); -} -\endcode -
- - -\section TutorialSparseFilling Filling a sparse matrix - -A DynamicSparseMatrix object can be set and updated just like any dense matrix using the coeffRef(row,col) method. If the coefficient is not stored yet, then it will be inserted in the matrix. Here is an example: -\code -DynamicSparseMatrix aux(1000,1000); -for (...) - for each i - for each j interacting with i - aux.coeffRef(i,j) += foo(o1,o2); -SparseMatrix mat(aux); // convert the DynamicSparseMatrix to a SparseMatrix -\endcode - -Sometimes, however, we simply want to set all the coefficients of a matrix before using it through standard matrix operations (addition, product, etc.). In that case it faster to use the low-level startFill()/fill()/fillrand()/endFill() interface. Even though this interface is availabe for both sparse matrix types, their respective restrictions slightly differ from one representation to the other. In all case, a call to startFill() set the matrix to zero, and the fill*() functions will fail if the coefficient already exist. - -As a first difference, for SparseMatrix, the fill*() functions can only be called inside a startFill()/endFill() pair, and no other member functions are allowed during the filling process, i.e., until endFill() has been called. On the other hand, a DynamicSparseMatrix is always in a stable state, and the startFill()/endFill() functions are only for compatibility purpose. - -Another difference is that the fill*() functions must be called with increasing outer indices for a SparseMatrix, while they can be random for a DynamicSparseMatrix. - -Finally, the fill() function assumes the coefficient are inserted in a sorted order per inner vector, while the fillrand() variante allows random insertions (the outer indices must still be sorted for SparseMatrix). - -Some examples: - -1 - If you can set the coefficients in exactly the same order that the storage order, then the matrix can be filled directly and very efficiently. Here is an example initializing a random, row-major sparse matrix: -\code -SparseMatrix m(rows,cols); -m.startFill(rows*cols*percent_of_non_zero); // estimate of the number of nonzeros (optional) -for (int i=0; i\ m(rows,cols); -m.startFill(rows*cols*percent_of_non_zero); // estimate of the number of nonzeros (optional) -for (int i=0; i\ m(rows,cols); -{ - RandomSetter > setter(m); - for (int k=0; k\().solveTriangular(dv2); -\endcode - -The product of a sparse matrix A time a dense matrix/vector dv with A symmetric can be optimized by telling that to Eigen: -\code -res = A.marked() * dv; // if all coefficients of A are stored -res = A.marked() * dv; // if only the upper part of A is stored -res = A.marked() * dv; // if only the lower part of A is stored -\endcode - - -\section TutorialSparseDirectSolvers Using the direct solvers - -TODO - -\subsection TutorialSparseDirectSolvers_LLT LLT -Cholmod, Taucs. - -\subsection TutorialSparseDirectSolvers_LDLT LDLT - - -\subsection TutorialSparseDirectSolvers_LU LU -SuperLU, UmfPack. - -*/ - -} diff --git a/doc/UnalignedArrayAssert.dox b/doc/UnalignedArrayAssert.dox deleted file mode 100644 index 05c784ca9..000000000 --- a/doc/UnalignedArrayAssert.dox +++ /dev/null @@ -1,77 +0,0 @@ -namespace Eigen { - -/** \page UnalignedArrayAssert Explanation of the assertion on unaligned arrays - -Hello! You are seeing this webpage because your program terminated on an assertion failure like this one: -
-my_program: path/to/eigen2/Eigen/src/Core/MatrixStorage.h:44:
-Eigen::ei_matrix_array::ei_matrix_array()
-[with T = double, int Size = 2, int MatrixOptions = 2, bool Align = true]:
-Assertion `(reinterpret_cast(array) & 0xf) == 0 && "this assertion
-is explained here: http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html
-**** READ THIS WEB PAGE !!! ****"' failed.
-
- -There are 3 known causes for this issue. Please read on to understand them and learn how to fix them. - -\b Table \b of \b contents - - \ref c1 - - \ref c2 - - \ref c3 - - \ref explanation - -\section c1 Cause 1: Structures having Eigen objects as members - -If you have code like this, - -\code -class Foo -{ - //... - Eigen::Vector2d v; - //... -}; -//... -Foo *foo = new Foo; -\endcode - -then you need to read this separate page: \ref StructHavingEigenMembers "Structures Having Eigen Members". - -Note that here, Eigen::Vector2d is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types". - -\section c2 Cause 2: STL Containers - -If you use STL Containers such as std::vector, std::map, ..., with Eigen objects, like this, - -\code -std::vector my_vector; -std::map my_map; -\endcode - -then you need to read this separate page: \ref StlContainers "Using STL Containers with Eigen". - -Note that here, Eigen::Matrix2f is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types". - -\section c3 Cause 3: Passing Eigen objects by value - -If some function in your code is getting an Eigen object passed by value, like this, - -\code -void func(Eigen::Vector4d v); -\endcode - -then you need to read this separate page: \ref PassingByValue "Passing Eigen objects by value to functions". - -Note that here, Eigen::Vector4d is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types". - -\section explanation General explanation of this assertion - -\ref FixedSizeVectorizable "fixed-size vectorizable Eigen objects" must absolutely be created at 16-byte-aligned locations, otherwise SIMD instructions adressing them will crash. - -Eigen normally takes care of these alignment issues for you, by setting an alignment attribute on them and by overloading their "operator new". - -However there are a few corner cases where these alignment settings get overridden: they are the possible causes for this assertion. - -*/ - -} diff --git a/doc/buildexamplelist.sh b/doc/buildexamplelist.sh deleted file mode 100755 index 7d47ae39d..000000000 --- a/doc/buildexamplelist.sh +++ /dev/null @@ -1,15 +0,0 @@ -#!/bin/sh - -echo "namespace Eigen {" -echo "/** \page ExampleList" -echo "

Selected list of examples

" - -grep \\addexample $1/Eigen/src/*/*.h -R | cut -d \\ -f 2- | \ -while read example; -do -anchor=`echo "$example" | cut -d " " -f 2` -text=`echo "$example" | cut -d " " -f 4-` -echo "\\\li \\\ref $anchor \"$text\"" -done -echo "*/" -echo "}" \ No newline at end of file diff --git a/doc/cleanhierarchy.sh b/doc/cleanhierarchy.sh deleted file mode 100755 index fe44891b9..000000000 --- a/doc/cleanhierarchy.sh +++ /dev/null @@ -1,5 +0,0 @@ -#!/bin/sh - -sed -i 's/^.li.*MatrixBase\<.*gt.*a.$/ /g' $1 -sed -i 's/^.li.*MapBase\<.*gt.*a.$/ /g' $1 -sed -i 's/^.li.*RotationBase\<.*gt.*a.$/ /g' $1 diff --git a/doc/snippets/CMakeLists.txt b/doc/snippets/CMakeLists.txt index 72bd7770a..195ed9e99 100644 --- a/doc/snippets/CMakeLists.txt +++ b/doc/snippets/CMakeLists.txt @@ -1,25 +1,25 @@ -FILE(GLOB snippets_SRCS "*.cpp") +FILE(GLOB snippets_SRCS "*.cpp" "../../unsupported/snippets/*.cpp") ADD_CUSTOM_TARGET(all_snippets) FOREACH(snippet_src ${snippets_SRCS}) -GET_FILENAME_COMPONENT(snippet ${snippet_src} NAME_WE) -SET(compile_snippet_target compile_${snippet}) -SET(compile_snippet_src ${compile_snippet_target}.cpp) -FILE(READ ${snippet_src} snippet_source_code) -CONFIGURE_FILE(${CMAKE_CURRENT_SOURCE_DIR}/compile_snippet.cpp.in - ${CMAKE_CURRENT_BINARY_DIR}/${compile_snippet_src}) -ADD_EXECUTABLE(${compile_snippet_target} - ${CMAKE_CURRENT_BINARY_DIR}/${compile_snippet_src}) -GET_TARGET_PROPERTY(compile_snippet_executable - ${compile_snippet_target} LOCATION) -ADD_CUSTOM_COMMAND( - TARGET ${compile_snippet_target} - POST_BUILD - COMMAND ${compile_snippet_executable} - ARGS >${CMAKE_CURRENT_BINARY_DIR}/${snippet}.out -) -ADD_DEPENDENCIES(all_snippets ${compile_snippet_target}) -set_source_files_properties(${CMAKE_CURRENT_BINARY_DIR}/${compile_snippet_src} - PROPERTIES OBJECT_DEPENDS ${snippet_src}) + GET_FILENAME_COMPONENT(snippet ${snippet_src} NAME_WE) + SET(compile_snippet_target compile_${snippet}) + SET(compile_snippet_src ${compile_snippet_target}.cpp) + FILE(READ ${snippet_src} snippet_source_code) + CONFIGURE_FILE(${CMAKE_CURRENT_SOURCE_DIR}/compile_snippet.cpp.in + ${CMAKE_CURRENT_BINARY_DIR}/${compile_snippet_src}) + ADD_EXECUTABLE(${compile_snippet_target} + ${CMAKE_CURRENT_BINARY_DIR}/${compile_snippet_src}) + GET_TARGET_PROPERTY(compile_snippet_executable + ${compile_snippet_target} LOCATION) + ADD_CUSTOM_COMMAND( + TARGET ${compile_snippet_target} + POST_BUILD + COMMAND ${compile_snippet_executable} + ARGS >${CMAKE_CURRENT_BINARY_DIR}/${snippet}.out + ) + ADD_DEPENDENCIES(all_snippets ${compile_snippet_target}) + set_source_files_properties(${CMAKE_CURRENT_BINARY_DIR}/${compile_snippet_src} + PROPERTIES OBJECT_DEPENDS ${snippet_src}) ENDFOREACH(snippet_src) diff --git a/unsupported/AdolcSupport/AdolcForward.h b/unsupported/AdolcSupport/AdolcForward.h deleted file mode 100644 index 0a55311a8..000000000 --- a/unsupported/AdolcSupport/AdolcForward.h +++ /dev/null @@ -1,96 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. Eigen itself is part of the KDE project. -// -// Copyright (C) 2008-2009 Gael Guennebaud -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see . - -#ifndef EIGEN_ADLOC_FORWARD_H -#define EIGEN_ADLOC_FORWARD_H - -//-------------------------------------------------------------------------------- -// -// This file provides support for adolc's adouble type in forward mode. -// ADOL-C is a C++ automatic differentiation library, -// see http://www.math.tu-dresden.de/~adol-c/ for more information. -// -// Note that the maximal number of directions is controlled by -// the preprocessor token NUMBER_DIRECTIONS. The default is 2. -// -//-------------------------------------------------------------------------------- - -#define ADOLC_TAPELESS -#ifndef NUMBER_DIRECTIONS -# define NUMBER_DIRECTIONS 2 -#endif -#include - -// adolc defines some very stupid macros: -#if defined(malloc) -# undef malloc -#endif - -#if defined(calloc) -# undef calloc -#endif - -#if defined(realloc) -# undef realloc -#endif - -#include - -namespace Eigen { - -template<> struct NumTraits -{ - typedef adtl::adouble Real; - typedef adtl::adouble FloatingPoint; - enum { - IsComplex = 0, - HasFloatingPoint = 1, - ReadCost = 1, - AddCost = 1, - MulCost = 1 - }; -}; - -} - -// the Adolc's type adouble is defined in the adtl namespace -// therefore, the following ei_* functions *must* be defined -// in the same namespace -namespace adtl { - - inline const adouble& ei_conj(const adouble& x) { return x; } - inline const adouble& ei_real(const adouble& x) { return x; } - inline adouble ei_imag(const adouble&) { return 0.; } - inline adouble ei_abs(const adouble& x) { return fabs(x); } - inline adouble ei_abs2(const adouble& x) { return x*x; } - inline adouble ei_sqrt(const adouble& x) { return sqrt(x); } - inline adouble ei_exp(const adouble& x) { return exp(x); } - inline adouble ei_log(const adouble& x) { return log(x); } - inline adouble ei_sin(const adouble& x) { return sin(x); } - inline adouble ei_cos(const adouble& x) { return cos(x); } - inline adouble ei_pow(const adouble& x, adouble y) { return pow(x, y); } - -} - -#endif // EIGEN_ADLOC_FORWARD_H diff --git a/unsupported/Eigen/AdolcForward b/unsupported/Eigen/AdolcForward new file mode 100644 index 000000000..1845cb917 --- /dev/null +++ b/unsupported/Eigen/AdolcForward @@ -0,0 +1,120 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008-2009 Gael Guennebaud +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see . + +#ifndef EIGEN_ADLOC_FORWARD +#define EIGEN_ADLOC_FORWARD + +//-------------------------------------------------------------------------------- +// +// This file provides support for adolc's adouble type in forward mode. +// ADOL-C is a C++ automatic differentiation library, +// see http://www.math.tu-dresden.de/~adol-c/ for more information. +// +// Note that the maximal number of directions is controlled by +// the preprocessor token NUMBER_DIRECTIONS. The default is 2. +// +//-------------------------------------------------------------------------------- + +#define ADOLC_TAPELESS +#ifndef NUMBER_DIRECTIONS +# define NUMBER_DIRECTIONS 2 +#endif +#include + +// adolc defines some very stupid macros: +#if defined(malloc) +# undef malloc +#endif + +#if defined(calloc) +# undef calloc +#endif + +#if defined(realloc) +# undef realloc +#endif + +#include + +namespace Eigen { + +namespace unsupported { + +/** \ingroup Unsupported_modules + * \defgroup AdolcForward_Module Adolc forward module + * This module provides support for adolc's adouble type in forward mode. + * ADOL-C is a C++ automatic differentiation library, + * see http://www.math.tu-dresden.de/~adol-c/ for more information. + * It mainly consists in: + * - a struct Eigen::NumTraits specialization + * - overloads of ei_* math function for adtl::adouble type. + * + * Note that the maximal number of directions is controlled by + * the preprocessor token NUMBER_DIRECTIONS. The default is 2. + * + * \code + * #include + * \endcode + */ + //@{ + +} // namespace unsupported + +} // namespace Eigen + +// the Adolc's type adouble is defined in the adtl namespace +// therefore, the following ei_* functions *must* be defined +// in the same namespace +namespace adtl { + + inline const adouble& ei_conj(const adouble& x) { return x; } + inline const adouble& ei_real(const adouble& x) { return x; } + inline adouble ei_imag(const adouble&) { return 0.; } + inline adouble ei_abs(const adouble& x) { return fabs(x); } + inline adouble ei_abs2(const adouble& x) { return x*x; } + inline adouble ei_sqrt(const adouble& x) { return sqrt(x); } + inline adouble ei_exp(const adouble& x) { return exp(x); } + inline adouble ei_log(const adouble& x) { return log(x); } + inline adouble ei_sin(const adouble& x) { return sin(x); } + inline adouble ei_cos(const adouble& x) { return cos(x); } + inline adouble ei_pow(const adouble& x, adouble y) { return pow(x, y); } + +} + +namespace Eigen { namespace unsupported { /*@}*/ } } + +template<> struct EigenNumTraits +{ + typedef adtl::adouble Real; + typedef adtl::adouble FloatingPoint; + enum { + IsComplex = 0, + HasFloatingPoint = 1, + ReadCost = 1, + AddCost = 1, + MulCost = 1 + }; +}; + +#endif // EIGEN_ADLOC_FORWARD diff --git a/unsupported/doc/unsupported.dox b/unsupported/doc/unsupported.dox new file mode 100644 index 000000000..abef263ad --- /dev/null +++ b/unsupported/doc/unsupported.dox @@ -0,0 +1,15 @@ + +namespace Eigen { + +namespace unsupported { + +/** \defgroup Unsupported_modules Unsupported modules + * + * The unsupported modules are contributions from various users. They are + * provided "as is", without any support. Nevertheless, they are subject to be + * included in Eigen in the future. + */ + +} // namespace unsupported + +} // namespace Eigen \ No newline at end of file -- cgit v1.2.3