- Added problem size constructor to decompositions that did not have one. It preallocates member data structures.

- Updated unit tests to check above constructor.
- In the compute() method of decompositions: Made temporary matrices/vectors class members to avoid heap allocations during compute() (when dynamic matrices are used, of course).

These  changes can speed up decomposition computation time when a solver instance is used to solve multiple same-sized problems. An added benefit is that the compute() method can now be invoked in contexts were heap allocations are forbidden, such as in real-time control loops.

CAVEAT: Not all of the decompositions in the Eigenvalues module have a heap-allocation-free compute() method. A future patch may address this issue, but some required API changes need to be incorporated first.

29 files changed, 396 insertions, 121 deletions

diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index f92a72c7b..206ccef4d 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -66,6 +66,7 @@ template<typename _MatrixType> class LDLT
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
     typedef typename ei_plain_col_type<MatrixType, int>::type IntColVectorType;
+    typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> TmpMatrixType;
 
     /** \brief Default Constructor.
       *
@@ -80,12 +81,17 @@ template<typename _MatrixType> class LDLT
       * according to the specified problem \a size.
       * \sa LDLT()
       */
-    LDLT(int size) : m_matrix(size,size), m_p(size), m_transpositions(size), m_isInitialized(false) {}
+    LDLT(int size) : m_matrix(size, size),
+                     m_p(size),
+                     m_transpositions(size),
+                     m_temporary(size),
+                     m_isInitialized(false) {}
 
     LDLT(const MatrixType& matrix)
       : m_matrix(matrix.rows(), matrix.cols()),
         m_p(matrix.rows()),
         m_transpositions(matrix.rows()),
+        m_temporary(matrix.rows()),
         m_isInitialized(false)
     {
       compute(matrix);
@@ -175,6 +181,7 @@ template<typename _MatrixType> class LDLT
     MatrixType m_matrix;
     IntColVectorType m_p;
     IntColVectorType m_transpositions; // FIXME do we really need to store permanently the transpositions?
+    TmpMatrixType m_temporary;
     int m_sign;
     bool m_isInitialized;
 };
@@ -206,7 +213,7 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
   // By using a temorary, packet-aligned products are guarenteed. In the LLT
   // case this is unnecessary because the diagonal is included and will always
   // have optimal alignment.
-  Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> _temporary(size);
+  m_temporary.resize(size);
 
   for (int j = 0; j < size; ++j)
   {
@@ -251,11 +258,11 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
 
     int endSize = size - j - 1;
     if (endSize > 0) {
-      _temporary.tail(endSize).noalias() = m_matrix.block(j+1,0, endSize, j)
+      m_temporary.tail(endSize).noalias() = m_matrix.block(j+1,0, endSize, j)
                                 * m_matrix.col(j).head(j).conjugate();
 
       m_matrix.row(j).tail(endSize) = m_matrix.row(j).tail(endSize).conjugate()
-                                    - _temporary.tail(endSize).transpose();
+                                    - m_temporary.tail(endSize).transpose();
 
       if(ei_abs(Djj) > cutoff)
       {
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index 8b01e0ccb..22d0c91c8 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -82,6 +82,15 @@ template<typename _MatrixType, int _UpLo> class LLT
     */
     LLT() : m_matrix(), m_isInitialized(false) {}
 
+    /** \brief Default Constructor with memory preallocation
+      *
+      * Like the default constructor but with preallocation of the internal data
+      * according to the specified problem \a size.
+      * \sa LLT()
+      */
+    LLT(int size) : m_matrix(size, size),
+                    m_isInitialized(false) {}
+
     LLT(const MatrixType& matrix)
       : m_matrix(matrix.rows(), matrix.cols()),
         m_isInitialized(false)
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
index 0b2e0b0ec..8e5f1310a 100644
--- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -95,7 +95,24 @@ template<typename _MatrixType> class ComplexEigenSolver
       * The default constructor is useful in cases in which the user intends to
       * perform decompositions via compute().
       */
-    ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false)
+    ComplexEigenSolver()
+            : m_eivec(),
+              m_eivalues(),
+              m_schur(),
+              m_isInitialized(false)
+    {}
+    
+    /** \brief Default Constructor with memory preallocation
+      *
+      * Like the default constructor but with preallocation of the internal data
+      * according to the specified problem \a size.
+      * \sa ComplexEigenSolver()
+      */
+    ComplexEigenSolver(int size)
+            : m_eivec(size, size),
+              m_eivalues(size),
+              m_schur(size),
+              m_isInitialized(false)
     {}
 
     /** \brief Constructor; computes eigendecomposition of given matrix. 
@@ -107,6 +124,7 @@ template<typename _MatrixType> class ComplexEigenSolver
     ComplexEigenSolver(const MatrixType& matrix)
             : m_eivec(matrix.rows(),matrix.cols()),
               m_eivalues(matrix.cols()),
+              m_schur(matrix.rows()),
               m_isInitialized(false)
     {
       compute(matrix);
@@ -179,6 +197,7 @@ template<typename _MatrixType> class ComplexEigenSolver
   protected:
     EigenvectorType m_eivec;
     EigenvalueType m_eivalues;
+    ComplexSchur<MatrixType> m_schur;
     bool m_isInitialized;
 };
 
@@ -193,8 +212,8 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
 
   // Step 1: Do a complex Schur decomposition, A = U T U^*
   // The eigenvalues are on the diagonal of T.
-  ComplexSchur<MatrixType> schur(matrix);
-  m_eivalues = schur.matrixT().diagonal();
+  m_schur.compute(matrix);
+  m_eivalues = m_schur.matrixT().diagonal();
 
   // Step 2: Compute X such that T = X D X^(-1), where D is the diagonal of T.
   // The matrix X is unit triangular.
@@ -205,10 +224,10 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
     // Compute X(i,k) using the (i,k) entry of the equation X T = D X
     for(int i=k-1 ; i>=0 ; i--)
     {
-      X.coeffRef(i,k) = -schur.matrixT().coeff(i,k);
+      X.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
       if(k-i-1>0)
-        X.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * X.col(k).segment(i+1,k-i-1)).value();
-      ComplexScalar z = schur.matrixT().coeff(i,i) - schur.matrixT().coeff(k,k);
+        X.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * X.col(k).segment(i+1,k-i-1)).value();
+      ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
       if(z==ComplexScalar(0))
       {
 	// If the i-th and k-th eigenvalue are equal, then z equals 0. 
@@ -220,7 +239,7 @@ void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
   }
 
   // Step 3: Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
-  m_eivec = schur.matrixU() * X;
+  m_eivec = m_schur.matrixU() * X;
   // .. and normalize the eigenvectors
   for(int k=0 ; k<n ; k++)
   {
diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h
index 7c9ab03c8..8874aa257 100644
--- a/Eigen/src/Eigenvalues/ComplexSchur.h
+++ b/Eigen/src/Eigenvalues/ComplexSchur.h
@@ -96,7 +96,11 @@ template<typename _MatrixType> class ComplexSchur
       * \sa compute() for an example.
       */
     ComplexSchur(int size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
-      : m_matT(size,size), m_matU(size,size), m_isInitialized(false), m_matUisUptodate(false)
+      : m_matT(size,size),
+        m_matU(size,size),
+        m_hess(size),
+        m_isInitialized(false),
+        m_matUisUptodate(false)
     {}
 
     /** \brief Constructor; computes Schur decomposition of given matrix. 
@@ -111,6 +115,7 @@ template<typename _MatrixType> class ComplexSchur
     ComplexSchur(const MatrixType& matrix, bool skipU = false)
             : m_matT(matrix.rows(),matrix.cols()),
               m_matU(matrix.rows(),matrix.cols()),
+              m_hess(matrix.rows()),
               m_isInitialized(false),
               m_matUisUptodate(false)
     {
@@ -182,6 +187,7 @@ template<typename _MatrixType> class ComplexSchur
 
   protected:
     ComplexMatrixType m_matT, m_matU;
+    HessenbergDecomposition<MatrixType> m_hess;
     bool m_isInitialized;
     bool m_matUisUptodate;
 
@@ -300,10 +306,10 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
 
   // Reduce to Hessenberg form
   // TODO skip Q if skipU = true
-  HessenbergDecomposition<MatrixType> hess(matrix);
+  m_hess.compute(matrix);
 
-  m_matT = hess.matrixH().template cast<ComplexScalar>();
-  if(!skipU)  m_matU = hess.matrixQ().template cast<ComplexScalar>();
+  m_matT = m_hess.matrixH().template cast<ComplexScalar>();
+  if(!skipU)  m_matU = m_hess.matrixQ().template cast<ComplexScalar>();
 
   // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
 
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h
index 32dee92e9..114bfab6f 100644
--- a/Eigen/src/Eigenvalues/EigenSolver.h
+++ b/Eigen/src/Eigenvalues/EigenSolver.h
@@ -122,6 +122,17 @@ template<typename _MatrixType> class EigenSolver
       */
     EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) {}
 
+    /** \brief Default Constructor with memory preallocation
+      *
+      * Like the default constructor but with preallocation of the internal data
+      * according to the specified problem \a size.
+      * \sa EigenSolver()
+      */
+    EigenSolver(int size)
+      : m_eivec(size, size),
+        m_eivalues(size),
+        m_isInitialized(false) {}
+
     /** \brief Constructor; computes eigendecomposition of given matrix. 
       * 
       * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
index 3e8c4c64d..451604252 100644
--- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h
+++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
@@ -91,7 +91,8 @@ template<typename _MatrixType> class HessenbergDecomposition
       * \sa compute() for an example.
       */
     HessenbergDecomposition(int size = Size==Dynamic ? 2 : Size)
-      : m_matrix(size,size)
+      : m_matrix(size,size),
+        m_temp(size)
     {
       if(size>1)
         m_hCoeffs.resize(size-1);
@@ -107,12 +108,13 @@ template<typename _MatrixType> class HessenbergDecomposition
       * \sa matrixH() for an example.
       */
     HessenbergDecomposition(const MatrixType& matrix)
-      : m_matrix(matrix)
+      : m_matrix(matrix),
+        m_temp(matrix.rows())
     {
       if(matrix.rows()<2)
         return;
       m_hCoeffs.resize(matrix.rows()-1,1);
-      _compute(m_matrix, m_hCoeffs);
+      _compute(m_matrix, m_hCoeffs, m_temp);
     }
 
     /** \brief Computes Hessenberg decomposition of given matrix. 
@@ -137,7 +139,7 @@ template<typename _MatrixType> class HessenbergDecomposition
       if(matrix.rows()<2)
         return;
       m_hCoeffs.resize(matrix.rows()-1,1);
-      _compute(m_matrix, m_hCoeffs);
+      _compute(m_matrix, m_hCoeffs, m_temp);
     }
 
     /** \brief Returns the Householder coefficients.
@@ -226,13 +228,14 @@ template<typename _MatrixType> class HessenbergDecomposition
 
   private:
 
-    static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs);
     typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType;
     typedef typename NumTraits<Scalar>::Real RealScalar;
-
+    static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
+    
   protected:
     MatrixType m_matrix;
     CoeffVectorType m_hCoeffs;
+    VectorType m_temp;
 };
 
 #ifndef EIGEN_HIDE_HEAVY_CODE
@@ -250,11 +253,11 @@ template<typename _MatrixType> class HessenbergDecomposition
   * \sa packedMatrix()
   */
 template<typename MatrixType>
-void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs)
+void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
 {
   assert(matA.rows()==matA.cols());
   int n = matA.rows();
-  VectorType temp(n);
+  temp.resize(n);
   for (int i = 0; i<n-1; ++i)
   {
     // let's consider the vector v = i-th column starting at position i+1
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
index a66ab8ace..dd4a1ab7e 100644
--- a/Eigen/src/Eigenvalues/RealSchur.h
+++ b/Eigen/src/Eigenvalues/RealSchur.h
@@ -78,6 +78,7 @@ template<typename _MatrixType> class RealSchur
     typedef typename MatrixType::Scalar Scalar;
     typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
     typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
+    typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
 
     /** \brief Default constructor.
       *
@@ -92,7 +93,9 @@ template<typename _MatrixType> class RealSchur
       */
     RealSchur(int size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
             : m_matT(size, size),
-              m_matU(size, size), 
+              m_matU(size, size),
+              m_workspaceVector(size),
+              m_hess(size),
               m_isInitialized(false)
     { }
 
@@ -108,6 +111,8 @@ template<typename _MatrixType> class RealSchur
     RealSchur(const MatrixType& matrix)
             : m_matT(matrix.rows(),matrix.cols()),
               m_matU(matrix.rows(),matrix.cols()),
+              m_workspaceVector(matrix.rows()),
+              m_hess(matrix.rows()),
               m_isInitialized(false)
     {
       compute(matrix);
@@ -165,6 +170,8 @@ template<typename _MatrixType> class RealSchur
     
     MatrixType m_matT;
     MatrixType m_matU;
+    ColumnVectorType m_workspaceVector;
+    HessenbergDecomposition<MatrixType> m_hess;
     bool m_isInitialized;
 
     typedef Matrix<Scalar,3,1> Vector3s;
@@ -185,14 +192,13 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix)
 
   // Step 1. Reduce to Hessenberg form
   // TODO skip Q if skipU = true
-  HessenbergDecomposition<MatrixType> hess(matrix);
-  m_matT = hess.matrixH();
-  m_matU = hess.matrixQ();
+  m_hess.compute(matrix);
+  m_matT = m_hess.matrixH();
+  m_matU = m_hess.matrixQ();
 
   // Step 2. Reduce to real Schur form  
-  typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
-  ColumnVectorType workspaceVector(m_matU.cols());
-  Scalar* workspace = &workspaceVector.coeffRef(0);
+  m_workspaceVector.resize(m_matU.cols());
+  Scalar* workspace = &m_workspaceVector.coeffRef(0);
 
   // The matrix m_matT is divided in three parts. 
   // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index b31f77e9b..41a600290 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -58,14 +58,16 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
 
     SelfAdjointEigenSolver()
         : m_eivec(int(Size), int(Size)),
-          m_eivalues(int(Size))
+          m_eivalues(int(Size)),
+          m_subdiag(int(TridiagonalizationType::SizeMinusOne))
     {
       ei_assert(Size!=Dynamic);
     }
 
     SelfAdjointEigenSolver(int size)
         : m_eivec(size, size),
-          m_eivalues(size)
+          m_eivalues(size),
+          m_subdiag(TridiagonalizationType::SizeMinusOne)
     {}
 
     /** Constructors computing the eigenvalues of the selfadjoint matrix \a matrix,
@@ -75,8 +77,10 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       */
     SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
       : m_eivec(matrix.rows(), matrix.cols()),
-        m_eivalues(matrix.cols())
+        m_eivalues(matrix.cols()),
+        m_subdiag()
     {
+      if (matrix.rows() > 1) m_subdiag.resize(matrix.rows() - 1);
       compute(matrix, computeEigenvectors);
     }
 
@@ -89,8 +93,10 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       */
     SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
       : m_eivec(matA.rows(), matA.cols()),
-        m_eivalues(matA.cols())
+        m_eivalues(matA.cols()),
+        m_subdiag()
     {
+      if (matA.rows() > 1) m_subdiag.resize(matA.rows() - 1);
       compute(matA, matB, computeEigenvectors);
     }
 
@@ -132,6 +138,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
   protected:
     MatrixType m_eivec;
     RealVectorType m_eivalues;
+    typename TridiagonalizationType::SubDiagonalType m_subdiag;
     #ifndef NDEBUG
     bool m_eigenvectorsOk;
     #endif
@@ -187,27 +194,27 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::compute(
   // the latter avoids multiple memory allocation when the same SelfAdjointEigenSolver is used multiple times...
   // (same for diag and subdiag)
   RealVectorType& diag = m_eivalues;
-  typename TridiagonalizationType::SubDiagonalType subdiag(n-1);
-  TridiagonalizationType::decomposeInPlace(m_eivec, diag, subdiag, computeEigenvectors);
+  m_subdiag.resize(n-1);
+  TridiagonalizationType::decomposeInPlace(m_eivec, diag, m_subdiag, computeEigenvectors);
 
   int end = n-1;
   int start = 0;
   while (end>0)
   {
     for (int i = start; i<end; ++i)
-      if (ei_isMuchSmallerThan(ei_abs(subdiag[i]),(ei_abs(diag[i])+ei_abs(diag[i+1]))))
-        subdiag[i] = 0;
+      if (ei_isMuchSmallerThan(ei_abs(m_subdiag[i]),(ei_abs(diag[i])+ei_abs(diag[i+1]))))
+        m_subdiag[i] = 0;
 
     // find the largest unreduced block
-    while (end>0 && subdiag[end-1]==0)
+    while (end>0 && m_subdiag[end-1]==0)
       end--;
     if (end<=0)
       break;
     start = end - 1;
-    while (start>0 && subdiag[start-1]!=0)
+    while (start>0 && m_subdiag[start-1]!=0)
       start--;
 
-    ei_tridiagonal_qr_step(diag.data(), subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
+    ei_tridiagonal_qr_step(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
   }
 
   // Sort eigenvalues and corresponding vectors.
diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h
index a0b772888..c320242ef 100644
--- a/Eigen/src/Eigenvalues/Tridiagonalization.h
+++ b/Eigen/src/Eigenvalues/Tridiagonalization.h
@@ -210,8 +210,8 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
     // i.e., A = H A H' where H = I - h v v' and v = matA.col(i).tail(n-i-1)
     matA.col(i).coeffRef(i+1) = 1;
 
-    hCoeffs.tail(n-i-1) = (matA.corner(BottomRight,remainingSize,remainingSize).template selfadjointView<Lower>()
-                        * (ei_conj(h) * matA.col(i).tail(remainingSize)));
+    hCoeffs.tail(n-i-1).noalias() = (matA.corner(BottomRight,remainingSize,remainingSize).template selfadjointView<Lower>()
+                                  * (ei_conj(h) * matA.col(i).tail(remainingSize)));
 
     hCoeffs.tail(n-i-1) += (ei_conj(h)*Scalar(-0.5)*(hCoeffs.tail(remainingSize).dot(matA.col(i).tail(remainingSize)))) * matA.col(i).tail(n-i-1);
 
diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h
index 1fe15bfd2..241fb9d16 100644
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -81,6 +81,14 @@ template<typename _MatrixType> class FullPivLU
       */
     FullPivLU();
 
+    /** \brief Default Constructor with memory preallocation
+      *
+      * Like the default constructor but with preallocation of the internal data
+      * according to the specified problem \a size.
+      * \sa FullPivLU()
+      */
+    FullPivLU(int rows, int cols);
+
     /** Constructor.
       *
       * \param matrix the matrix of which to compute the LU decomposition.
@@ -377,6 +385,8 @@ template<typename _MatrixType> class FullPivLU
     MatrixType m_lu;
     PermutationPType m_p;
     PermutationQType m_q;
+    IntColVectorType m_rowsTranspositions;
+    IntRowVectorType m_colsTranspositions;
     int m_det_pq, m_nonzero_pivots;
     RealScalar m_maxpivot, m_prescribedThreshold;
     bool m_isInitialized, m_usePrescribedThreshold;
@@ -389,8 +399,26 @@ FullPivLU<MatrixType>::FullPivLU()
 }
 
 template<typename MatrixType>
+FullPivLU<MatrixType>::FullPivLU(int rows, int cols)
+  : m_lu(rows, cols),
+    m_p(rows),
+    m_q(cols),
+    m_rowsTranspositions(rows),
+    m_colsTranspositions(cols),
+    m_isInitialized(false),
+    m_usePrescribedThreshold(false)
+{
+}
+
+template<typename MatrixType>
 FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
-  : m_isInitialized(false), m_usePrescribedThreshold(false)
+  : m_lu(matrix.rows(), matrix.cols()),
+    m_p(matrix.rows()),
+    m_q(matrix.cols()),
+    m_rowsTranspositions(matrix.rows()),
+    m_colsTranspositions(matrix.cols()),
+    m_isInitialized(false),
+    m_usePrescribedThreshold(false)
 {
   compute(matrix);
 }
@@ -407,9 +435,9 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
 
   // will store the transpositions, before we accumulate them at the end.
   // can't accumulate on-the-fly because that will be done in reverse order for the rows.
-  IntColVectorType rows_transpositions(matrix.rows());
-  IntRowVectorType cols_transpositions(matrix.cols());
-  int number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. rows_transpositions[i]!=i
+  m_rowsTranspositions.resize(matrix.rows());
+  m_colsTranspositions.resize(matrix.cols());
+  int number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
 
   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
   m_maxpivot = RealScalar(0);
@@ -442,8 +470,8 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
       m_nonzero_pivots = k;
       for(int i = k; i < size; ++i)
       {
-        rows_transpositions.coeffRef(i) = i;
-        cols_transpositions.coeffRef(i) = i;
+        m_rowsTranspositions.coeffRef(i) = i;
+        m_colsTranspositions.coeffRef(i) = i;
       }
       break;
     }
@@ -453,8 +481,8 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
     // Now that we've found the pivot, we need to apply the row/col swaps to
     // bring it to the location (k,k).
 
-    rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
-    cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
+    m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
+    m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
     if(k != row_of_biggest_in_corner) {
       m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
       ++number_of_transpositions;
@@ -478,11 +506,11 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
 
   m_p.setIdentity(rows);
   for(int k = size-1; k >= 0; --k)
-    m_p.applyTranspositionOnTheRight(k, rows_transpositions.coeff(k));
+    m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
 
   m_q.setIdentity(cols);
   for(int k = 0; k < size; ++k)
-    m_q.applyTranspositionOnTheRight(k, cols_transpositions.coeff(k));
+    m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
 
   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
   return *this;
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index 5062fd0ce..42c74ee53 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -83,6 +83,14 @@ template<typename _MatrixType> class PartialPivLU
     */
     PartialPivLU();
 
+    /** \brief Default Constructor with memory preallocation
+      *
+      * Like the default constructor but with preallocation of the internal data
+      * according to the specified problem \a size.
+      * \sa PartialPivLU()
+      */
+    PartialPivLU(int size);
+
     /** Constructor.
       *
       * \param matrix the matrix of which to compute the LU decomposition.
@@ -176,6 +184,7 @@ template<typename _MatrixType> class PartialPivLU
   protected:
     MatrixType m_lu;
     PermutationType m_p;
+    PermutationVectorType m_rowsTranspositions;
     int m_det_p;
     bool m_isInitialized;
 };
@@ -184,6 +193,17 @@ template<typename MatrixType>
 PartialPivLU<MatrixType>::PartialPivLU()
   : m_lu(),
     m_p(),
+    m_rowsTranspositions(),
+    m_det_p(0),
+    m_isInitialized(false)
+{
+}
+
+template<typename MatrixType>
+PartialPivLU<MatrixType>::PartialPivLU(int size)
+  : m_lu(size, size),
+    m_p(size),
+    m_rowsTranspositions(size),
     m_det_p(0),
     m_isInitialized(false)
 {
@@ -191,8 +211,9 @@ PartialPivLU<MatrixType>::PartialPivLU()
 
 template<typename MatrixType>
 PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)
-  : m_lu(),
-    m_p(),
+  : m_lu(matrix.rows(), matrix.rows()),
+    m_p(matrix.rows()),
+    m_rowsTranspositions(matrix.rows()),
     m_det_p(0),
     m_isInitialized(false)
 {
@@ -384,15 +405,15 @@ PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& ma
   ei_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
   const int size = matrix.rows();
 
-  PermutationVectorType rows_transpositions(size);
+  m_rowsTranspositions.resize(size);
 
   int nb_transpositions;
-  ei_partial_lu_inplace(m_lu, rows_transpositions, nb_transpositions);
+  ei_partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
   m_det_p = (nb_transpositions%2) ? -1 : 1;
 
   m_p.setIdentity(size);
   for(int k = size-1; k >= 0; --k)
-    m_p.applyTranspositionOnTheRight(k, rows_transpositions.coeff(k));
+    m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
 
   m_isInitialized = true;
   return *this;
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
index 27f80e046..5893125cd 100644
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -70,11 +70,38 @@ template<typename _MatrixType> class ColPivHouseholderQR
     * The default constructor is useful in cases in which the user intends to
     * perform decompositions via ColPivHouseholderQR::compute(const MatrixType&).
     */
-    ColPivHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
+    ColPivHouseholderQR()
+      : m_qr(),
+        m_hCoeffs(),
+        m_colsPermutation(),
+        m_colsTranspositions(),
+        m_temp(),
+        m_colSqNorms(),
+        m_isInitialized(false) {}
+
+    /** \brief Default Constructor with memory preallocation
+      *
+      * Like the default constructor but with preallocation of the internal data
+      * according to the specified problem \a size.
+      * \sa ColPivHouseholderQR()
+      */
+    ColPivHouseholderQR(int rows, int cols)
+      : m_qr(rows, cols),
+        m_hCoeffs(std::min(rows,cols)),
+        m_colsPermutation(cols),
+        m_colsTranspositions(cols),
+        m_temp(cols),
+        m_colSqNorms(cols),
+        m_isInitialized(false),
+        m_usePrescribedThreshold(false) {}
 
     ColPivHouseholderQR(const MatrixType& matrix)
       : m_qr(matrix.rows(), matrix.cols()),
         m_hCoeffs(std::min(matrix.rows(),matrix.cols())),
+        m_colsPermutation(matrix.cols()),
+        m_colsTranspositions(matrix.cols()),
+        m_temp(matrix.cols()),
+        m_colSqNorms(matrix.cols()),
         m_isInitialized(false),
         m_usePrescribedThreshold(false)
     {
@@ -121,7 +148,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
     const PermutationType& colsPermutation() const
     {
       ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
-      return m_cols_permutation;
+      return m_colsPermutation;
     }
 
     /** \returns the absolute value of the determinant of the matrix of which
@@ -307,7 +334,10 @@ template<typename _MatrixType> class ColPivHouseholderQR
   protected:
     MatrixType m_qr;
     HCoeffsType m_hCoeffs;
-    PermutationType m_cols_permutation;
+    PermutationType m_colsPermutation;
+    IntRowVectorType m_colsTranspositions;
+    RowVectorType m_temp;
+    RealRowVectorType m_colSqNorms;
     bool m_isInitialized, m_usePrescribedThreshold;
     RealScalar m_prescribedThreshold, m_maxpivot;
     int m_nonzero_pivots;
@@ -342,35 +372,35 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
   m_qr = matrix;
   m_hCoeffs.resize(size);
 
-  RowVectorType temp(cols);
+  m_temp.resize(cols);
 
-  IntRowVectorType cols_transpositions(matrix.cols());
+  m_colsTranspositions.resize(matrix.cols());
   int number_of_transpositions = 0;
 
-  RealRowVectorType colSqNorms(cols);
+  m_colSqNorms.resize(cols);
   for(int k = 0; k < cols; ++k)
-    colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
+    m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
 
-  RealScalar threshold_helper = colSqNorms.maxCoeff() * ei_abs2(NumTraits<Scalar>::epsilon()) / rows;
+  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * ei_abs2(NumTraits<Scalar>::epsilon()) / rows;
 
   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
   m_maxpivot = RealScalar(0);
 
   for(int k = 0; k < size; ++k)
   {
-    // first, we look up in our table colSqNorms which column has the biggest squared norm
+    // first, we look up in our table m_colSqNorms which column has the biggest squared norm
     int biggest_col_index;
-    RealScalar biggest_col_sq_norm = colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
+    RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
     biggest_col_index += k;
 
-    // since our table colSqNorms accumulates imprecision at every step, we must now recompute
+    // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
     // the actual squared norm of the selected column.
     // Note that not doing so does result in solve() sometimes returning inf/nan values
     // when running the unit test with 1000 repetitions.
     biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
 
     // we store that back into our table: it can't hurt to correct our table.
-    colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
+    m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
 
     // if the current biggest column is smaller than epsilon times the initial biggest column,
     // terminate to avoid generating nan/inf values.
@@ -388,10 +418,10 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
     }
 
     // apply the transposition to the columns
-    cols_transpositions.coeffRef(k) = biggest_col_index;
+    m_colsTranspositions.coeffRef(k) = biggest_col_index;
     if(k != biggest_col_index) {
       m_qr.col(k).swap(m_qr.col(biggest_col_index));
-      std::swap(colSqNorms.coeffRef(k), colSqNorms.coeffRef(biggest_col_index));
+      std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
       ++number_of_transpositions;
     }
 
@@ -407,15 +437,15 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
 
     // apply the householder transformation
     m_qr.corner(BottomRight, rows-k, cols-k-1)
-        .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1));
+        .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
 
     // update our table of squared norms of the columns
-    colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
+    m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
   }
 
-  m_cols_permutation.setIdentity(cols);
+  m_colsPermutation.setIdentity(cols);
   for(int k = 0; k < m_nonzero_pivots; ++k)
-    m_cols_permutation.applyTranspositionOnTheRight(k, cols_transpositions.coeff(k));
+    m_colsPermutation.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
 
   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
   m_isInitialized = true;
diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h
index aeb551cc7..8df2ed1e0 100644
--- a/Eigen/src/QR/FullPivHouseholderQR.h
+++ b/Eigen/src/QR/FullPivHouseholderQR.h
@@ -69,10 +69,38 @@ template<typename _MatrixType> class FullPivHouseholderQR
       * The default constructor is useful in cases in which the user intends to
       * perform decompositions via FullPivHouseholderQR::compute(const MatrixType&).
       */
-    FullPivHouseholderQR() : m_isInitialized(false) {}
+    FullPivHouseholderQR()
+      : m_qr(),
+        m_hCoeffs(),
+        m_rows_transpositions(),
+        m_cols_transpositions(),
+        m_cols_permutation(),
+        m_temp(),
+        m_isInitialized(false) {}
+
+    /** \brief Default Constructor with memory preallocation
+      *
+      * Like the default constructor but with preallocation of the internal data
+      * according to the specified problem \a size.
+      * \sa FullPivHouseholderQR()
+      */
+    FullPivHouseholderQR(int rows, int cols)
+      : m_qr(rows, cols),
+        m_hCoeffs(std::min(rows,cols)),
+        m_rows_transpositions(rows),
+        m_cols_transpositions(cols),
+        m_cols_permutation(cols),
+        m_temp(std::min(rows,cols)),
+        m_isInitialized(false) {}
 
     FullPivHouseholderQR(const MatrixType& matrix)
-      : m_isInitialized(false)
+      : m_qr(matrix.rows(), matrix.cols()),
+        m_hCoeffs(std::min(matrix.rows(), matrix.cols())),
+        m_rows_transpositions(matrix.rows()),
+        m_cols_transpositions(matrix.cols()),
+        m_cols_permutation(matrix.cols()),
+        m_temp(std::min(matrix.rows(), matrix.cols())),
+        m_isInitialized(false)
     {
       compute(matrix);
     }
@@ -233,7 +261,9 @@ template<typename _MatrixType> class FullPivHouseholderQR
     MatrixType m_qr;
     HCoeffsType m_hCoeffs;
     IntColVectorType m_rows_transpositions;
+    IntRowVectorType m_cols_transpositions;
     PermutationType m_cols_permutation;
+    RowVectorType m_temp;
     bool m_isInitialized;
     RealScalar m_precision;
     int m_rank;
@@ -269,12 +299,12 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
   m_qr = matrix;
   m_hCoeffs.resize(size);
 
-  RowVectorType temp(cols);
+  m_temp.resize(cols);
 
   m_precision = NumTraits<Scalar>::epsilon() * size;
 
   m_rows_transpositions.resize(matrix.rows());
-  IntRowVectorType cols_transpositions(matrix.cols());
+  m_cols_transpositions.resize(matrix.cols());
   int number_of_transpositions = 0;
 
   RealScalar biggest(0);
@@ -298,14 +328,14 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
       for(int i = k; i < size; i++)
       {
         m_rows_transpositions.coeffRef(i) = i;
-        cols_transpositions.coeffRef(i) = i;
+        m_cols_transpositions.coeffRef(i) = i;
         m_hCoeffs.coeffRef(i) = Scalar(0);
       }
       break;
     }
 
     m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
-    cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
+    m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
     if(k != row_of_biggest_in_corner) {
       m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
       ++number_of_transpositions;
@@ -320,12 +350,12 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
     m_qr.coeffRef(k,k) = beta;
 
     m_qr.corner(BottomRight, rows-k, cols-k-1)
-        .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1));
+        .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
   }
 
   m_cols_permutation.setIdentity(cols);
   for(int k = 0; k < size; ++k)
-    m_cols_permutation.applyTranspositionOnTheRight(k, cols_transpositions.coeff(k));
+    m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
 
   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
   m_isInitialized = true;
diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h
index c4abc16ff..10b9fb93f 100644
--- a/Eigen/src/QR/HouseholderQR.h
+++ b/Eigen/src/QR/HouseholderQR.h
@@ -71,11 +71,24 @@ template<typename _MatrixType> class HouseholderQR
     * The default constructor is useful in cases in which the user intends to
     * perform decompositions via HouseholderQR::compute(const MatrixType&).
     */
-    HouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
+    HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
+
+    /** \brief Default Constructor with memory preallocation
+      *
+      * Like the default constructor but with preallocation of the internal data
+      * according to the specified problem \a size.
+      * \sa HouseholderQR()
+      */
+    HouseholderQR(int rows, int cols)
+      : m_qr(rows, cols),
+        m_hCoeffs(std::min(rows,cols)),
+        m_temp(cols),
+        m_isInitialized(false) {}
 
     HouseholderQR(const MatrixType& matrix)
       : m_qr(matrix.rows(), matrix.cols()),
         m_hCoeffs(std::min(matrix.rows(),matrix.cols())),
+        m_temp(matrix.cols()),
         m_isInitialized(false)
     {
       compute(matrix);
@@ -159,6 +172,7 @@ template<typename _MatrixType> class HouseholderQR
   protected:
     MatrixType m_qr;
     HCoeffsType m_hCoeffs;
+    RowVectorType m_temp;
     bool m_isInitialized;
 };
 
@@ -190,7 +204,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
   m_qr = matrix;
   m_hCoeffs.resize(size);
 
-  RowVectorType temp(cols);
+  m_temp.resize(cols);
 
   for(int k = 0; k < size; ++k)
   {
@@ -203,7 +217,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
 
     // apply H to remaining part of m_qr from the left
     m_qr.corner(BottomRight, remainingRows, remainingCols)
-        .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingRows-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1));
+        .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingRows-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
   }
   m_isInitialized = true;
   return *this;
diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h
index b7cab023b..9323c0180 100644
--- a/Eigen/src/SVD/JacobiSVD.h
+++ b/Eigen/src/SVD/JacobiSVD.h
@@ -93,10 +93,35 @@ template<typename MatrixType, unsigned int Options> class JacobiSVD
 
   public:
 
+    /** \brief Default Constructor.
+      *
+      * The default constructor is useful in cases in which the user intends to
+      * perform decompositions via JacobiSVD::compute(const MatrixType&).
+      */
     JacobiSVD() : m_isInitialized(false) {}
 
-    JacobiSVD(const MatrixType& matrix) : m_isInitialized(false)
+
+    /** \brief Default Constructor with memory preallocation
+      *
+      * Like the default constructor but with preallocation of the internal data
+      * according to the specified problem \a size.
+      * \sa JacobiSVD()
+      */
+    JacobiSVD(int rows, int cols) : m_matrixU(rows, rows),
+                                    m_matrixV(cols, cols),
+                                    m_singularValues(std::min(rows, cols)),
+                                    m_workMatrix(rows, cols),
+                                    m_isInitialized(false) {}
+
+    JacobiSVD(const MatrixType& matrix) : m_matrixU(matrix.rows(), matrix.rows()),
+                                          m_matrixV(matrix.cols(), matrix.cols()),
+                                          m_singularValues(),
+                                          m_workMatrix(),
+                                          m_isInitialized(false)
     {
+      const int minSize = std::min(matrix.rows(), matrix.cols());
+      m_singularValues.resize(minSize);
+      m_workMatrix.resize(minSize, minSize);
       compute(matrix);
     }
 
@@ -124,6 +149,7 @@ template<typename MatrixType, unsigned int Options> class JacobiSVD
     MatrixUType m_matrixU;
     MatrixVType m_matrixV;
     SingularValuesType m_singularValues;
+    WorkMatrixType m_workMatrix;
     bool m_isInitialized;
 
     template<typename _MatrixType, unsigned int _Options, bool _IsComplex>
@@ -289,17 +315,16 @@ struct ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options, true>
 template<typename MatrixType, unsigned int Options>
 JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const MatrixType& matrix)
 {
-  WorkMatrixType work_matrix;
   int rows = matrix.rows();
   int cols = matrix.cols();
   int diagSize = std::min(rows, cols);
   m_singularValues.resize(diagSize);
   const RealScalar precision = 2 * NumTraits<Scalar>::epsilon();
 
-  if(!ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options>::run(matrix, work_matrix, *this)
-  && !ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options>::run(matrix, work_matrix, *this))
+  if(!ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options>::run(matrix, m_workMatrix, *this)
+  && !ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options>::run(matrix, m_workMatrix, *this))
   {
-    work_matrix = matrix.block(0,0,diagSize,diagSize);
+    m_workMatrix = matrix.block(0,0,diagSize,diagSize);
     if(ComputeU) m_matrixU.setIdentity(rows,rows);
     if(ComputeV) m_matrixV.setIdentity(cols,cols);
   }
@@ -312,19 +337,19 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const Ma
     {
       for(int q = 0; q < p; ++q)
       {
-        if(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p)))
-            > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*precision)
+        if(std::max(ei_abs(m_workMatrix.coeff(p,q)),ei_abs(m_workMatrix.coeff(q,p)))
+            > std::max(ei_abs(m_workMatrix.coeff(p,p)),ei_abs(m_workMatrix.coeff(q,q)))*precision)
         {
           finished = false;
-          ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options>::run(work_matrix, *this, p, q);
+          ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options>::run(m_workMatrix, *this, p, q);
 
           PlanarRotation<RealScalar> j_left, j_right;
-          ei_real_2x2_jacobi_svd(work_matrix, p, q, &j_left, &j_right);
+          ei_real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
 
-          work_matrix.applyOnTheLeft(p,q,j_left);
+          m_workMatrix.applyOnTheLeft(p,q,j_left);
           if(ComputeU) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
 
-          work_matrix.applyOnTheRight(p,q,j_right);
+          m_workMatrix.applyOnTheRight(p,q,j_right);
           if(ComputeV) m_matrixV.applyOnTheRight(p,q,j_right);
         }
       }
@@ -333,9 +358,9 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const Ma
 
   for(int i = 0; i < diagSize; ++i)
   {
-    RealScalar a = ei_abs(work_matrix.coeff(i,i));
+    RealScalar a = ei_abs(m_workMatrix.coeff(i,i));
     m_singularValues.coeffRef(i) = a;
-    if(ComputeU && (a!=RealScalar(0))) m_matrixU.col(i) *= work_matrix.coeff(i,i)/a;
+    if(ComputeU && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
   }
 
   for(int i = 0; i < diagSize; i++)
diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h
index ed0e69f91..a9e22dfd4 100644
--- a/Eigen/src/SVD/SVD.h
+++ b/Eigen/src/SVD/SVD.h
@@ -73,11 +73,25 @@ template<typename _MatrixType> class SVD
     */
     SVD() : m_matU(), m_matV(), m_sigma(), m_isInitialized(false) {}
 
-    SVD(const MatrixType& matrix)
-      : m_matU(matrix.rows(), matrix.rows()),
-        m_matV(matrix.cols(),matrix.cols()),
-        m_sigma(matrix.cols()),
-        m_isInitialized(false)
+    /** \brief Default Constructor with memory preallocation
+      *
+      * Like the default constructor but with preallocation of the internal data
+      * according to the specified problem \a size.
+      * \sa JacobiSVD()
+      */
+    SVD(int rows, int cols) : m_matU(rows, rows),
+                              m_matV(cols,cols),
+                              m_sigma(std::min(rows, cols)),
+                              m_workMatrix(rows, cols),
+                              m_rv1(cols),
+                              m_isInitialized(false) {}
+
+    SVD(const MatrixType& matrix) : m_matU(matrix.rows(), matrix.rows()),
+                                    m_matV(matrix.cols(),matrix.cols()),
+                                    m_sigma(std::min(matrix.rows(), matrix.cols())),
+                                    m_workMatrix(matrix.rows(), matrix.cols()),
+                                    m_rv1(matrix.cols()),
+                                    m_isInitialized(false)
     {
       compute(matrix);
     }
@@ -165,6 +179,8 @@ template<typename _MatrixType> class SVD
     MatrixVType m_matV;
     /** \internal */
     SingularValuesType m_sigma;
+    MatrixType m_workMatrix;
+    RowVector m_rv1;
     bool m_isInitialized;
     int m_rows, m_cols;
 };
@@ -185,11 +201,12 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
   m_matU.setZero();
   m_sigma.resize(n);
   m_matV.resize(n,n);
+  m_workMatrix = matrix;
 
   int max_iters = 30;
 
   MatrixVType& V = m_matV;
-  MatrixType A = matrix;
+  MatrixType& A = m_workMatrix;
   SingularValuesType& W = m_sigma;
 
   bool flag;
@@ -198,14 +215,14 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
   bool convergence = true;
   Scalar eps = NumTraits<Scalar>::dummy_precision();
 
-  RowVector rv1(n);
+  m_rv1.resize(n);
 
   g = scale = anorm = 0;
   // Householder reduction to bidiagonal form.
   for (i=0; i<n; i++)
   {
     l = i+2;
-    rv1[i] = scale*g;
+    m_rv1[i] = scale*g;
     g = s = scale = 0.0;
     if (i < m)
     {
@@ -246,16 +263,16 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
         g = -sign(ei_sqrt(s),f);
         h = f*g - s;
         A(i,l-1) = f-g;
-        rv1.tail(n-l+1) = A.row(i).tail(n-l+1)/h;
+        m_rv1.tail(n-l+1) = A.row(i).tail(n-l+1)/h;
         for (j=l-1; j<m; j++)
         {
           s = A.row(i).tail(n-l+1).dot(A.row(j).tail(n-l+1));
-          A.row(j).tail(n-l+1) += s*rv1.tail(n-l+1).transpose();
+          A.row(j).tail(n-l+1) += s*m_rv1.tail(n-l+1).transpose();
         }
         A.row(i).tail(n-l+1) *= scale;
       }
     }
-    anorm = std::max( anorm, (ei_abs(W[i])+ei_abs(rv1[i])) );
+    anorm = std::max( anorm, (ei_abs(W[i])+ei_abs(m_rv1[i])) );
   }
   // Accumulation of right-hand transformations.
   for (i=n-1; i>=0; i--)
@@ -277,7 +294,7 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
       V.col(i).tail(n-l).setZero();
     }
     V(i, i) = 1.0;
-    g = rv1[i];
+    g = m_rv1[i];
     l = i;
   }
   // Accumulation of left-hand transformations.
@@ -318,7 +335,7 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
         nm = l-1;
         // Note that rv1[1] is always zero.
         //if ((double)(ei_abs(rv1[l])+anorm) == anorm)
-        if (l==0 || ei_abs(rv1[l]) <= eps*anorm)
+        if (l==0 || ei_abs(m_rv1[l]) <= eps*anorm)
         {
           flag = false;
           break;
@@ -333,8 +350,8 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
         s = 1.0;
         for (i=l ;i<k+1; i++)
         {
-          f = s*rv1[i];
-          rv1[i] = c*rv1[i];
+          f = s*m_rv1[i];
+          m_rv1[i] = c*m_rv1[i];
           //if ((double)(ei_abs(f)+anorm) == anorm)
           if (ei_abs(f) <= eps*anorm)
             break;
@@ -363,8 +380,8 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
       x = W[l]; // Shift from bottom 2-by-2 minor.
       nm = k-1;
       y = W[nm];
-      g = rv1[nm];
-      h = rv1[k];
+      g = m_rv1[nm];
+      h = m_rv1[k];
       f = ((y-z)*(y+z) + (g-h)*(g+h))/(Scalar(2.0)*h*y);
       g = pythag(f,1.0);
       f = ((x-z)*(x+z) + h*((y/(f+sign(g,f)))-h))/x;
@@ -373,13 +390,13 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
       for (j=l; j<=nm; j++)
       {
         i = j+1;
-        g = rv1[i];
+        g = m_rv1[i];
         y = W[i];
         h = s*g;
         g = c*g;
 
         z = pythag(f,h);
-        rv1[j] = z;
+        m_rv1[j] = z;
         c = f/z;
         s = h/z;
         f = x*c + g*s;
@@ -401,8 +418,8 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix)
         x = c*y - s*g;
         A.applyOnTheRight(i,j,PlanarRotation<Scalar>(c,s));
       }
-      rv1[l] = 0.0;
-      rv1[k] = f;
+      m_rv1[l] = 0.0;
+      m_rv1[k] = f;
       W[k]   = x;
     }
   }
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index a446f5d73..0ae26c7d5 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -163,4 +163,8 @@ void test_cholesky()
   CALL_SUBTEST_7( cholesky_verify_assert<Matrix3d>() );
   CALL_SUBTEST_8( cholesky_verify_assert<MatrixXf>() );
   CALL_SUBTEST_2( cholesky_verify_assert<MatrixXd>() );
+
+  // Test problem size constructors
+  CALL_SUBTEST_9( LLT<MatrixXf>(10) );
+  CALL_SUBTEST_9( LDLT<MatrixXf>(10) );
 }
diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp
index 4e973c877..b3d9ac24b 100644
--- a/test/eigensolver_complex.cpp
+++ b/test/eigensolver_complex.cpp
@@ -63,4 +63,7 @@ void test_eigensolver_complex()
     CALL_SUBTEST_3( eigensolver(Matrix<std::complex<float>, 1, 1>()) );
     CALL_SUBTEST_4( eigensolver(Matrix3f()) );
   }
+
+  // Test problem size constructors
+  CALL_SUBTEST_5(ComplexEigenSolver<MatrixXf>(10));
 }
diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp
index 5dc18f141..f24c3b4ed 100644
--- a/test/eigensolver_generic.cpp
+++ b/test/eigensolver_generic.cpp
@@ -89,4 +89,7 @@ void test_eigensolver_generic()
   CALL_SUBTEST_2( eigensolver_verify_assert<MatrixXd>() );
   CALL_SUBTEST_4( eigensolver_verify_assert<Matrix2d>() );
   CALL_SUBTEST_5( eigensolver_verify_assert<MatrixXf>() );
+
+  // Test problem size constructors
+  CALL_SUBTEST_6(EigenSolver<MatrixXf>(10));
 }
diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp
index f2146cb88..70b3e6791 100644
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -129,5 +129,9 @@ void test_eigensolver_selfadjoint()
     CALL_SUBTEST_6( selfadjointeigensolver(Matrix<double,1,1>()) );
     CALL_SUBTEST_7( selfadjointeigensolver(Matrix<double,2,2>()) );
   }
+
+  // Test problem size constructors
+  CALL_SUBTEST_8(SelfAdjointEigenSolver<MatrixXf>(10));
+  CALL_SUBTEST_8(Tridiagonalization<MatrixXf>(10));
 }
 
diff --git a/test/hessenberg.cpp b/test/hessenberg.cpp
index cba9c8fda..ec1148bfc 100644
--- a/test/hessenberg.cpp
+++ b/test/hessenberg.cpp
@@ -61,4 +61,7 @@ void test_hessenberg()
   CALL_SUBTEST_3(( hessenberg<std::complex<float>,4>() ));
   CALL_SUBTEST_4(( hessenberg<float,Dynamic>(ei_random<int>(1,320)) ));
   CALL_SUBTEST_5(( hessenberg<std::complex<double>,Dynamic>(ei_random<int>(1,320)) ));
+
+  // Test problem size constructors
+  CALL_SUBTEST_6(HessenbergDecomposition<MatrixXf>(10));
 }
diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp
index c8a1495d2..bc9d93754 100644
--- a/test/jacobisvd.cpp
+++ b/test/jacobisvd.cpp
@@ -106,4 +106,7 @@ void test_jacobisvd()
   CALL_SUBTEST_3(( svd_verify_assert<Matrix3d>() ));
   CALL_SUBTEST_9(( svd_verify_assert<MatrixXf>() ));
   CALL_SUBTEST_11(( svd_verify_assert<MatrixXd>() ));
+
+  // Test problem size constructors
+  CALL_SUBTEST_12( JacobiSVD<MatrixXf>(10, 20) );
 }
diff --git a/test/lu.cpp b/test/lu.cpp
index 37e2990d2..22dca76d2 100644
--- a/test/lu.cpp
+++ b/test/lu.cpp
@@ -212,5 +212,9 @@ void test_lu()
     CALL_SUBTEST_6( lu_verify_assert<MatrixXcd>() );
 
     CALL_SUBTEST_7(( lu_non_invertible<Matrix<float,Dynamic,16> >() ));
+
+    // Test problem size constructors
+    CALL_SUBTEST_9( PartialPivLU<MatrixXf>(10) );
+    CALL_SUBTEST_9( FullPivLU<MatrixXf>(10, 20); );
   }
 }
diff --git a/test/qr.cpp b/test/qr.cpp
index 1ce1f46e6..90209639f 100644
--- a/test/qr.cpp
+++ b/test/qr.cpp
@@ -133,4 +133,7 @@ void test_qr()
   CALL_SUBTEST_6(qr_verify_assert<MatrixXd>());
   CALL_SUBTEST_7(qr_verify_assert<MatrixXcf>());
   CALL_SUBTEST_8(qr_verify_assert<MatrixXcd>());
+
+  // Test problem size constructors
+  CALL_SUBTEST_12(HouseholderQR<MatrixXf>(10, 20));
 }
diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp
index 0ca897b52..a34feedd9 100644
--- a/test/qr_colpivoting.cpp
+++ b/test/qr_colpivoting.cpp
@@ -157,4 +157,7 @@ void test_qr_colpivoting()
   CALL_SUBTEST_2(qr_verify_assert<MatrixXd>());
   CALL_SUBTEST_6(qr_verify_assert<MatrixXcf>());
   CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>());
+
+  // Test problem size constructors
+  CALL_SUBTEST_9(ColPivHouseholderQR<MatrixXf>(10, 20));
 }
diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp
index 7ad3af1fe..82c42c759 100644
--- a/test/qr_fullpivoting.cpp
+++ b/test/qr_fullpivoting.cpp
@@ -139,4 +139,7 @@ void test_qr_fullpivoting()
   CALL_SUBTEST_2(qr_verify_assert<MatrixXd>());
   CALL_SUBTEST_4(qr_verify_assert<MatrixXcf>());
   CALL_SUBTEST_3(qr_verify_assert<MatrixXcd>());
+
+  // Test problem size constructors
+  CALL_SUBTEST_7(FullPivHouseholderQR<MatrixXf>(10, 20));
 }
diff --git a/test/schur_complex.cpp b/test/schur_complex.cpp
index 3659a074c..b33411cf2 100644
--- a/test/schur_complex.cpp
+++ b/test/schur_complex.cpp
@@ -64,4 +64,7 @@ void test_schur_complex()
   CALL_SUBTEST_2(( schur<MatrixXcf>(ei_random<int>(1,50)) ));
   CALL_SUBTEST_3(( schur<Matrix<std::complex<float>, 1, 1> >() ));
   CALL_SUBTEST_4(( schur<Matrix<float, 3, 3, Eigen::RowMajor> >() ));
+
+  // Test problem size constructors
+  CALL_SUBTEST_5(ComplexSchur<MatrixXf>(10));
 }
diff --git a/test/schur_real.cpp b/test/schur_real.cpp
index c1747e2f5..d0aca4308 100644
--- a/test/schur_real.cpp
+++ b/test/schur_real.cpp
@@ -82,4 +82,7 @@ void test_schur_real()
   CALL_SUBTEST_2(( schur<MatrixXd>(ei_random<int>(1,50)) ));
   CALL_SUBTEST_3(( schur<Matrix<float, 1, 1> >() ));
   CALL_SUBTEST_4(( schur<Matrix<double, 3, 3, Eigen::RowMajor> >() ));
+
+  // Test problem size constructors
+  CALL_SUBTEST_5(RealSchur<MatrixXf>(10));
 }
diff --git a/test/svd.cpp b/test/svd.cpp
index 1b41d7032..9f3072d3b 100644
--- a/test/svd.cpp
+++ b/test/svd.cpp
@@ -113,4 +113,7 @@ void test_svd()
   CALL_SUBTEST_2( svd_verify_assert<Matrix4d>() );
   CALL_SUBTEST_3( svd_verify_assert<MatrixXf>() );
   CALL_SUBTEST_4( svd_verify_assert<MatrixXd>() );
+
+  // Test problem size constructors
+  CALL_SUBTEST_9( SVD<MatrixXf>(10, 20) );
 }