doc/InplaceDecomposition.dox

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; the Hessenberg, tridiagonal, Schur, and QZ reductions and the dense eigensolvers followed.
This feature is especially useful when dealing with huge 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:

<table class="example">
<tr><th>code</th><th>output</th></tr>
<tr>
  <td>\snippet TutorialInplaceLU.cpp init
  </td>
  <td>\snippet TutorialInplaceLU.out init
  </td>
</tr>
</table>

No surprise here! Then, let's declare our inplace LU object \c lu, and check the content of the matrix \c A:

<table class="example">
<tr>
  <td>\snippet TutorialInplaceLU.cpp declaration
  </td>
  <td>\snippet TutorialInplaceLU.out declaration
  </td>
</tr>
</table>

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:

<table class="example">
<tr>
  <td>\snippet TutorialInplaceLU.cpp matrixLU
  </td>
  <td>\snippet TutorialInplaceLU.out matrixLU
  </td>
</tr>
</table>

Then, one can use the \c lu object as usual, for instance to solve the Ax=b problem:
<table class="example">
<tr>
  <td>\snippet TutorialInplaceLU.cpp solve
  </td>
  <td>\snippet TutorialInplaceLU.out solve
  </td>
</tr>
</table>

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:

<table class="example">
<tr>
  <td>\snippet TutorialInplaceLU.cpp modifyA
  </td>
  <td>\snippet TutorialInplaceLU.out modifyA
  </td>
</tr>
</table>

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:
<table class="example">
<tr>
  <td>\snippet TutorialInplaceLU.cpp recompute
  </td>
  <td>\snippet TutorialInplaceLU.out recompute
  </td>
</tr>
</table>

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:
<table class="example">
<tr>
  <td>\snippet TutorialInplaceLU.cpp recompute_bis0
 </td>
  <td>\snippet TutorialInplaceLU.out recompute_bis0
 </td>
</tr>
</table>
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:
<table class="example">
<tr>
  <td>\snippet TutorialInplaceLU.cpp recompute_bis1
  </td>
  <td>\snippet TutorialInplaceLU.out recompute_bis1
 </td>
</tr>
</table>


Here is the list of matrix decompositions supporting this inplace mechanism:

- class LLT
- class LDLT
- class BunchKaufman
- class PartialPivLU
- class FullPivLU
- class HouseholderQR
- class ColPivHouseholderQR
- class RandColPivHouseholderQR
- class FullPivHouseholderQR
- class CompleteOrthogonalDecomposition
- class RandCompleteOrthogonalDecomposition
- class HessenbergDecomposition
- class Tridiagonalization
- class RealSchur
- class ComplexSchur (complex \c MatrixType only)
- class RealQZ
- class ComplexQZ
- class SelfAdjointEigenSolver
- class GeneralizedSelfAdjointEigenSolver
- class EigenSolver
- class ComplexEigenSolver (complex \c MatrixType only)
- class GeneralizedEigenSolver

The factorizations store their packed factors in the input matrix, as in the LU example above. For the reductions
and eigensolvers, the input matrix is used as the workspace of the algorithm; what it holds afterwards depends on
the class:

- HessenbergDecomposition and Tridiagonalization: the packed representation returned by packedMatrix().
- RealSchur, ComplexSchur, RealQZ, and ComplexQZ: the (quasi-)triangular factors T, or S and T; the unitary
  factors are stored in the decomposition object.
- SelfAdjointEigenSolver and GeneralizedSelfAdjointEigenSolver: the eigenvectors when they are requested, and
  eigenvectors() then refers to the input matrix; otherwise, unspecified values. The generalized solver also
  factorizes \c B in place.
- EigenSolver, ComplexEigenSolver, and GeneralizedEigenSolver: unspecified values; the results are stored in the
  decomposition object.

RealSchur may use a separate workspace even when instantiated with Ref<>. For an n-by-n matrix with n >= 128 and
an outer stride in bytes divisible by 1024, it allocates a temporary padded matrix of n*(n+1) coefficients, performs
the Hessenberg reduction and Schur iteration there, and copies T back to the bound input matrix. This also applies
to compute(); computeFromHessenberg() uses the same padding policy for the Schur iteration. Zero inputs still
incur the temporary allocation at these sizes and strides.

To avoid this temporary, bind a Ref<> with an outer stride in bytes that is not divisible by 1024. For example,
a 512-by-512 double matrix can use a 513-by-512 column-major backing matrix and bind its top 512 rows; the outer
stride is then 513*sizeof(double) bytes. Defining EIGEN_NO_MALLOC, or disabling allocation or deallocation through
EIGEN_RUNTIME_NO_MALLOC, also suppresses padding and uses the bound storage, which may be slower. These controls
do not eliminate other allocations: construction and resizing of solver storage, forming Schur vectors, or
evaluating input expressions may still require memory. Inplace decomposition does not guarantee allocation-free
construction or computation.

Since the decomposition object then refers to the input matrix, the input matrix must outlive it, as in the LU
example. The QZ and generalized eigensolver classes take two matrices, and both must be passed as non-const
lvalues to select the inplace constructor.

*/

}