InplaceDecomposition.dox 3.7 KB

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