- 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.

3 files changed, 103 insertions, 29 deletions

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;