1*bf2c3715SXin Linamespace Eigen { 2*bf2c3715SXin Li 3*bf2c3715SXin Li/** \eigenManualPage InplaceDecomposition Inplace matrix decompositions 4*bf2c3715SXin Li 5*bf2c3715SXin LiStarting from %Eigen 3.3, the LU, Cholesky, and QR decompositions can operate \em inplace, that is, directly within the given input matrix. 6*bf2c3715SXin LiThis feature is especially useful when dealing with huge matrices, and or when the available memory is very limited (embedded systems). 7*bf2c3715SXin Li 8*bf2c3715SXin LiTo 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. 9*bf2c3715SXin Li 10*bf2c3715SXin LiLet's start with the basic inclusions, and declaration of a 2x2 matrix \c A: 11*bf2c3715SXin Li 12*bf2c3715SXin Li<table class="example"> 13*bf2c3715SXin Li<tr><th>code</th><th>output</th></tr> 14*bf2c3715SXin Li<tr> 15*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.cpp init 16*bf2c3715SXin Li </td> 17*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.out init 18*bf2c3715SXin Li </td> 19*bf2c3715SXin Li</tr> 20*bf2c3715SXin Li</table> 21*bf2c3715SXin Li 22*bf2c3715SXin LiNo surprise here! Then, let's declare our inplace LU object \c lu, and check the content of the matrix \c A: 23*bf2c3715SXin Li 24*bf2c3715SXin Li<table class="example"> 25*bf2c3715SXin Li<tr> 26*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.cpp declaration 27*bf2c3715SXin Li </td> 28*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.out declaration 29*bf2c3715SXin Li </td> 30*bf2c3715SXin Li</tr> 31*bf2c3715SXin Li</table> 32*bf2c3715SXin Li 33*bf2c3715SXin LiHere, the \c lu object computes and stores the \c L and \c U factors within the memory held by the matrix \c A. 34*bf2c3715SXin LiThe coefficients of \c A have thus been destroyed during the factorization, and replaced by the L and U factors as one can verify: 35*bf2c3715SXin Li 36*bf2c3715SXin Li<table class="example"> 37*bf2c3715SXin Li<tr> 38*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.cpp matrixLU 39*bf2c3715SXin Li </td> 40*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.out matrixLU 41*bf2c3715SXin Li </td> 42*bf2c3715SXin Li</tr> 43*bf2c3715SXin Li</table> 44*bf2c3715SXin Li 45*bf2c3715SXin LiThen, one can use the \c lu object as usual, for instance to solve the Ax=b problem: 46*bf2c3715SXin Li<table class="example"> 47*bf2c3715SXin Li<tr> 48*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.cpp solve 49*bf2c3715SXin Li </td> 50*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.out solve 51*bf2c3715SXin Li </td> 52*bf2c3715SXin Li</tr> 53*bf2c3715SXin Li</table> 54*bf2c3715SXin Li 55*bf2c3715SXin LiHere, 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. 56*bf2c3715SXin Li 57*bf2c3715SXin LiSince the memory is shared between \c A and \c lu, modifying the matrix \c A will make \c lu invalid. 58*bf2c3715SXin LiThis can easily be verified by modifying the content of \c A and trying to solve the initial problem again: 59*bf2c3715SXin Li 60*bf2c3715SXin Li<table class="example"> 61*bf2c3715SXin Li<tr> 62*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.cpp modifyA 63*bf2c3715SXin Li </td> 64*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.out modifyA 65*bf2c3715SXin Li </td> 66*bf2c3715SXin Li</tr> 67*bf2c3715SXin Li</table> 68*bf2c3715SXin Li 69*bf2c3715SXin LiNote 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. 70*bf2c3715SXin Li 71*bf2c3715SXin LiIf one wants to update the factorization with the modified A, one has to call the compute method as usual: 72*bf2c3715SXin Li<table class="example"> 73*bf2c3715SXin Li<tr> 74*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.cpp recompute 75*bf2c3715SXin Li </td> 76*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.out recompute 77*bf2c3715SXin Li </td> 78*bf2c3715SXin Li</tr> 79*bf2c3715SXin Li</table> 80*bf2c3715SXin Li 81*bf2c3715SXin LiNote 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. 82*bf2c3715SXin LiThis can easily be verified as follows: 83*bf2c3715SXin Li<table class="example"> 84*bf2c3715SXin Li<tr> 85*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.cpp recompute_bis0 86*bf2c3715SXin Li </td> 87*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.out recompute_bis0 88*bf2c3715SXin Li </td> 89*bf2c3715SXin Li</tr> 90*bf2c3715SXin Li</table> 91*bf2c3715SXin LiThe matrix \c A1 is unchanged, and one can thus solve A1*x=b, and directly check the residual without any copy of \c A1: 92*bf2c3715SXin Li<table class="example"> 93*bf2c3715SXin Li<tr> 94*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.cpp recompute_bis1 95*bf2c3715SXin Li </td> 96*bf2c3715SXin Li <td>\snippet TutorialInplaceLU.out recompute_bis1 97*bf2c3715SXin Li </td> 98*bf2c3715SXin Li</tr> 99*bf2c3715SXin Li</table> 100*bf2c3715SXin Li 101*bf2c3715SXin Li 102*bf2c3715SXin LiHere is the list of matrix decompositions supporting this inplace mechanism: 103*bf2c3715SXin Li 104*bf2c3715SXin Li- class LLT 105*bf2c3715SXin Li- class LDLT 106*bf2c3715SXin Li- class PartialPivLU 107*bf2c3715SXin Li- class FullPivLU 108*bf2c3715SXin Li- class HouseholderQR 109*bf2c3715SXin Li- class ColPivHouseholderQR 110*bf2c3715SXin Li- class FullPivHouseholderQR 111*bf2c3715SXin Li- class CompleteOrthogonalDecomposition 112*bf2c3715SXin Li 113*bf2c3715SXin Li*/ 114*bf2c3715SXin Li 115*bf2c3715SXin Li}