From 7f7839c12f9bbd88bee1f3c9d7f5dbbe556045d4 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 4 Jul 2016 17:18:26 +0200 Subject: Add documentation and exemples for inplace decomposition. --- Eigen/src/Cholesky/LDLT.h | 4 +- Eigen/src/Cholesky/LLT.h | 4 +- Eigen/src/LU/FullPivLU.h | 4 +- Eigen/src/LU/PartialPivLU.h | 11 +-- Eigen/src/QR/ColPivHouseholderQR.h | 4 +- Eigen/src/QR/CompleteOrthogonalDecomposition.h | 4 +- Eigen/src/QR/FullPivHouseholderQR.h | 4 +- Eigen/src/QR/HouseholderQR.h | 4 +- doc/InplaceDecomposition.dox | 114 +++++++++++++++++++++++++ doc/Manual.dox | 2 + doc/examples/TutorialInplaceLU.cpp | 61 +++++++++++++ 11 files changed, 201 insertions(+), 15 deletions(-) create mode 100644 doc/InplaceDecomposition.dox create mode 100644 doc/examples/TutorialInplaceLU.cpp diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index a31b3d6aa..69e176110 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -43,6 +43,8 @@ namespace internal { * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky * decomposition to determine whether a system of equations has a solution. * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::ldlt(), SelfAdjointView::ldlt(), class LLT */ template class LDLT @@ -112,7 +114,7 @@ template class LDLT /** \brief Constructs a LDLT factorization from a given matrix * - * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. * * \sa LDLT(const EigenBase&) */ diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index ad163c749..bd966656d 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -41,6 +41,8 @@ template struct LLT_Traits; * Example: \include LLT_example.cpp * Output: \verbinclude LLT_example.out * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT */ /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH) @@ -96,7 +98,7 @@ template class LLT /** \brief Constructs a LDLT factorization from a given matrix * - * This overloaded constructor is provided for inplace solving when + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when * \c MatrixType is a Eigen::Ref. * * \sa LLT(const EigenBase&) diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 113b8c7b8..1632d3ac3 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -52,6 +52,8 @@ template struct traits > * \include class_FullPivLU.cpp * Output: \verbinclude class_FullPivLU.out * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse() */ template class FullPivLU @@ -99,7 +101,7 @@ template class FullPivLU /** \brief Constructs a LU factorization from a given matrix * - * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. * * \sa FullPivLU(const EigenBase&) */ diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index c862d9692..87ac6a281 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -68,6 +68,8 @@ struct enable_if_ref,Derived> { * * The data of the LU decomposition can be directly accessed through the methods matrixLU(), permutationP(). * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU */ template class PartialPivLU @@ -113,17 +115,10 @@ template class PartialPivLU template explicit PartialPivLU(const EigenBase& matrix); - /** Constructor for inplace decomposition + /** Constructor for \link InplaceDecomposition inplace decomposition \endlink * * \param matrix the matrix of which to compute the LU decomposition. * - * If \c MatrixType is an Eigen::Ref, then the storage of \a matrix will be shared - * between \a matrix and \c *this and the decomposition will take place in-place. - * The memory of \a matrix will be used througrough the lifetime of \c *this. In - * particular, further calls to \c this->compute(A) will still operate on the memory - * of \a matrix meaning. This also implies that the sizes of \c A must match the - * ones of \a matrix. - * * \warning The matrix should have full rank (e.g. if it's square, it should be invertible). * If you need to deal with non-full rank, use class FullPivLU instead. */ diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index db50b5675..ee6989897 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -41,6 +41,8 @@ template struct traits > * This decomposition performs column pivoting in order to be rank-revealing and improve * numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR. * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::colPivHouseholderQr() */ template class ColPivHouseholderQR @@ -135,7 +137,7 @@ template class ColPivHouseholderQR /** \brief Constructs a QR factorization from a given matrix * - * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. * * \sa ColPivHouseholderQR(const EigenBase&) */ diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h index 967eb35dd..69ffe4f42 100644 --- a/Eigen/src/QR/CompleteOrthogonalDecomposition.h +++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h @@ -39,6 +39,8 @@ struct traits > * \b Q and \b Z are unitary matrices and \b T an upper triangular matrix of * size rank-by-rank. \b A may be rank deficient. * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::completeOrthogonalDecomposition() */ template @@ -117,7 +119,7 @@ class CompleteOrthogonalDecomposition { /** \brief Constructs a complete orthogonal decomposition from a given matrix * - * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. * * \sa CompleteOrthogonalDecomposition(const EigenBase&) */ diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index 5c2f57d04..4c85a4693 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -50,6 +50,8 @@ struct traits > * This decomposition performs a very prudent full pivoting in order to be rank-revealing and achieve optimal * numerical stability. The trade-off is that it is slower than HouseholderQR and ColPivHouseholderQR. * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::fullPivHouseholderQr() */ template class FullPivHouseholderQR @@ -136,7 +138,7 @@ template class FullPivHouseholderQR /** \brief Constructs a QR factorization from a given matrix * - * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. * * \sa FullPivHouseholderQR(const EigenBase&) */ diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index f2a9cc080..2e64fa7fe 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -37,6 +37,8 @@ namespace Eigen { * This Householder QR decomposition is faster, but less numerically stable and less feature-full than * FullPivHouseholderQR or ColPivHouseholderQR. * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::householderQr() */ template class HouseholderQR @@ -104,7 +106,7 @@ template class HouseholderQR /** \brief Constructs a QR factorization from a given matrix * - * This overloaded constructor is provided for inplace solving when + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when * \c MatrixType is a Eigen::Ref. * * \sa HouseholderQR(const EigenBase&) diff --git a/doc/InplaceDecomposition.dox b/doc/InplaceDecomposition.dox new file mode 100644 index 000000000..fc89db4a1 --- /dev/null +++ b/doc/InplaceDecomposition.dox @@ -0,0 +1,114 @@ +namespace Eigen { + +/** \eigenManualPage InplaceDecomposition Inplace matrix decompositions + +Starting from %Eigen 3.3, the LU, Cholesky, and QR decompositions can operate \em inplace, that is, directly within the given input matrix. +This feature is especially useful when dealing with huae matrices, and or when the available memory is very limited (embedded systems). + +To this end, the respective decomposition class must be instantiated with a Ref<> matrix type, and the decomposition object must be constructed with the input matrix as argument. As an example, let us consider an inplace LU decomposition with partial pivoting. + +Let's start with the basic inclusions, and declaration of a 2x2 matrix \c A: + + + + + + + +
codeoutput
\snippet TutorialInplaceLU.cpp init + \snippet TutorialInplaceLU.out init +
+ +No surprise here! Then, let's declare our inplace LU object \c lu, and check the content of the matrix \c A: + + + + + + +
\snippet TutorialInplaceLU.cpp declaration + \snippet TutorialInplaceLU.out declaration +
+ +Here, the \c lu object computes and stores the \c L and \c U factors within the memory held by the matrix \c A. +The coefficients of \c A have thus been destroyed during the factorization, and replaced by the L and U factors as one can verify: + + + + + + +
\snippet TutorialInplaceLU.cpp matrixLU + \snippet TutorialInplaceLU.out matrixLU +
+ +Then, one can use the \c lu object as usual, for instance to solve the Ax=b problem: + + + + + +
\snippet TutorialInplaceLU.cpp solve + \snippet TutorialInplaceLU.out solve +
+ +Here, since the content of the original matrix \c A has been lost, we had to declared a new matrix \c A0 to verify the result. +Since the memory is shared between \c A and \c lu, modifying the matrix \c A will make \c lu invalid. +This can easily be verified by modifying the content of \c A and trying to solve the initial problem again: + + + + + + +
\snippet TutorialInplaceLU.cpp modifyA + \snippet TutorialInplaceLU.out modifyA +
+ +Note that there is no shared pointer under the hood, it is the \b responsibility \b of \b the \b user to keep the input matrix \c A in life as long as \c lu is living. + +If one wants to update the factorization with the modified A, one has to call the compute method as usual: + + + + + +
\snippet TutorialInplaceLU.cpp recompute + \snippet TutorialInplaceLU.out recompute +
+ +Note that calling compute does not change the memory which is referenced by the \c lu object. Therefore, if the compute method is called with another matrix \c A1 different than \c A, then the content of \c A1 won't be modified. This is still the content of \c A that will be used to store the L and U factors of the matrix \c A1. +This can easily be verified as follows: + + + + + +
\snippet TutorialInplaceLU.cpp recompute_bis0 + \snippet TutorialInplaceLU.out recompute_bis0 +
+The matrix \c A1 is unchanged, and one can thus solve A1*x=b, and directly check the residual without any copy of \c A1: + + + + + +
\snippet TutorialInplaceLU.cpp recompute_bis1 + \snippet TutorialInplaceLU.out recompute_bis1 +
+ + +Here is the list of matrix decomposition supported this inplace mechanism: + +- class LLT +- class LDLT +- class PartialPivLU +- class FullPivLU +- class HouseholderQR +- class ColPivHouseholderQR +- class FullPivHouseholderQR +- class CompleteOrthogonalDecomposition + +*/ + +} \ No newline at end of file diff --git a/doc/Manual.dox b/doc/Manual.dox index 8940b0070..5167cb9e5 100644 --- a/doc/Manual.dox +++ b/doc/Manual.dox @@ -106,6 +106,8 @@ namespace Eigen { \ingroup DenseLinearSolvers_chapter */ /** \addtogroup LeastSquares \ingroup DenseLinearSolvers_chapter */ +/** \addtogroup InplaceDecomposition + \ingroup DenseLinearSolvers_chapter */ /** \addtogroup DenseLinearSolvers_Reference \ingroup DenseLinearSolvers_chapter */ diff --git a/doc/examples/TutorialInplaceLU.cpp b/doc/examples/TutorialInplaceLU.cpp new file mode 100644 index 000000000..cb9c59b60 --- /dev/null +++ b/doc/examples/TutorialInplaceLU.cpp @@ -0,0 +1,61 @@ +#include +struct init { + init() { std::cout << "[" << "init" << "]" << std::endl; } +}; +init init_obj; +// [init] +#include +#include + +using namespace std; +using namespace Eigen; + +int main() +{ + MatrixXd A(2,2); + A << 2, -1, 1, 3; + cout << "Here is the input matrix A before decomposition:\n" << A << endl; +cout << "[init]" << endl; + +cout << "[declaration]" << endl; + PartialPivLU > lu(A); + cout << "Here is the input matrix A after decomposition:\n" << A << endl; +cout << "[declaration]" << endl; + +cout << "[matrixLU]" << endl; + cout << "Here is the matrix storing the L and U factors:\n" << lu.matrixLU() << endl; +cout << "[matrixLU]" << endl; + +cout << "[solve]" << endl; + MatrixXd A0(2,2); A0 << 2, -1, 1, 3; + VectorXd b(2); b << 1, 2; + VectorXd x = lu.solve(b); + cout << "Residual: " << (A0 * x - b).norm() << endl; +cout << "[solve]" << endl; + +cout << "[modifyA]" << endl; + A << 3, 4, -2, 1; + x = lu.solve(b); + cout << "Residual: " << (A0 * x - b).norm() << endl; +cout << "[modifyA]" << endl; + +cout << "[recompute]" << endl; + A0 = A; // save A + lu.compute(A); + x = lu.solve(b); + cout << "Residual: " << (A0 * x - b).norm() << endl; +cout << "[recompute]" << endl; + +cout << "[recompute_bis0]" << endl; + MatrixXd A1(2,2); + A1 << 5,-2,3,4; + lu.compute(A1); + cout << "Here is the input matrix A1 after decomposition:\n" << A1 << endl; +cout << "[recompute_bis0]" << endl; + +cout << "[recompute_bis1]" << endl; + x = lu.solve(b); + cout << "Residual: " << (A1 * x - b).norm() << endl; +cout << "[recompute_bis1]" << endl; + +} -- cgit v1.2.3