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

2 files changed, 76 insertions, 34 deletions

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;
     }
   }