JacobiSVD:

  * fix preallocating constructors, allocate U and V of the right size for computation options
  * complete documentation and internal comments
  * improve unit test, test inf/nan values

3 files changed, 185 insertions, 60 deletions

diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h
index e7156a3fb..e50c9ed15 100644
--- a/Eigen/src/SVD/JacobiSVD.h
+++ b/Eigen/src/SVD/JacobiSVD.h
@@ -25,14 +25,19 @@
 #ifndef EIGEN_JACOBISVD_H
 #define EIGEN_JACOBISVD_H
 
-// forward declarations (needed by ICC)
-// the empty bodies are required by MSVC
+// forward declaration (needed by ICC)
+// the empty body is required by MSVC
 template<typename MatrixType, int QRPreconditioner,
          bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
 struct ei_svd_precondition_2x2_block_to_be_real {};
 
 
-/*** QR preconditioners (R-SVD) ***/
+/*** QR preconditioners (R-SVD)
+ ***
+ *** Their role is to reduce the problem of computing the SVD to the case of a square matrix.
+ *** This approach, known as R-SVD, is an optimization for rectangular-enough matrices, and is a requirement for
+ *** JacobiSVD which by itself is only able to work on square matrices.
+ ***/
 
 enum { PreconditionIfMoreColsThanRows, PreconditionIfMoreRowsThanCols };
 
@@ -40,11 +45,11 @@ template<typename MatrixType, int QRPreconditioner, int Case>
 struct ei_qr_preconditioner_should_do_anything
 {
   enum { a = MatrixType::RowsAtCompileTime != Dynamic &&
-              MatrixType::ColsAtCompileTime != Dynamic &&
-              MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
+             MatrixType::ColsAtCompileTime != Dynamic &&
+             MatrixType::ColsAtCompileTime <= MatrixType::RowsAtCompileTime,
          b = MatrixType::RowsAtCompileTime != Dynamic &&
-              MatrixType::ColsAtCompileTime != Dynamic &&
-              MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
+             MatrixType::ColsAtCompileTime != Dynamic &&
+             MatrixType::RowsAtCompileTime <= MatrixType::ColsAtCompileTime,
          ret = !( (QRPreconditioner == NoQRPreconditioner) ||
                   (Case == PreconditionIfMoreColsThanRows && bool(a)) ||
                   (Case == PreconditionIfMoreRowsThanCols && bool(b)) )
@@ -64,6 +69,8 @@ struct ei_qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
   }
 };
 
+/*** preconditioner using FullPivHouseholderQR ***/
+
 template<typename MatrixType>
 struct ei_qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
 {
@@ -101,6 +108,8 @@ struct ei_qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner,
   }
 };
 
+/*** preconditioner using ColPivHouseholderQR ***/
+
 template<typename MatrixType>
 struct ei_qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
 {
@@ -146,6 +155,8 @@ struct ei_qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner,
   }
 };
 
+/*** preconditioner using HouseholderQR ***/
+
 template<typename MatrixType>
 struct ei_qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
 {
@@ -198,12 +209,36 @@ struct ei_qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, Precon
   *
   * \class JacobiSVD
   *
-  * \brief Jacobi SVD decomposition of a square matrix
+  * \brief Two-sided Jacobi SVD decomposition of a square matrix
   *
   * \param MatrixType the type of the matrix of which we are computing the SVD decomposition
   * \param QRPreconditioner this optional parameter allows to specify the type of QR decomposition that will be used internally
   *                        for the R-SVD step for non-square matrices. See discussion of possible values below.
   *
+  * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product
+  *   \f[ A = U S V^* \f]
+  * where \a U is a n-by-n unitary, \a V is a p-by-p unitary, and \a S is a n-by-p real positive matrix which is zero outside of its main diagonal;
+  * the diagonal entries of S are known as the \em singular \em values of \a A and the columns of \a U and \a V are known as the left
+  * and right \em singular \em vectors of \a A respectively.
+  *
+  * Singular values are always sorted in decreasing order.
+  *
+  * This JacobiSVD decomposition computes only the singular values by default. If you want \a U or \a V, you need to ask for them explicitly.
+  *
+  * You can ask for only \em thin \a U or \a V to be computed, meaning the following. In case of a rectangular n-by-p matrix, letting \a m be the
+  * smaller value among \a n and \a p, there are only \a m singular vectors; the remaining columns of \a U and \a V do not correspond to actual
+  * singular vectors. Asking for \em thin \a U or \a V means asking for only their \a m first columns to be formed. So \a U is then a n-by-m matrix,
+  * and \a V is then a p-by-m matrix. Notice that thin \a U and \a V are all you need for (least squares) solving.
+  *
+  * Here's an example demonstrating basic usage:
+  * \include JacobiSVD_basic.cpp
+  * Output: \verbinclude JacobiSVD_basic.out
+  * 
+  * This %JacobiSVD class a two-sided Jacobi R-SVD decomposition, ensuring optimal reliability and accuracy. The downside is that it's slower than
+  * bidiagonalizing SVD algorithms for large square matrices; however its complexity is still \f$ O(n^2p) \f$ where \a n is the smaller dimension and
+  * \a p is the greater dimension, meaning that it is still of the same order of complexity as the faster bidiagonalizing R-SVD algorithms.
+  * In particular, like any R-SVD, it takes advantage of non-squareness in that its complexity is only linear in the greater dimension.
+  *
   * The possible values for QRPreconditioner are:
   * \li ColPivHouseholderQRPreconditioner is the default. In practice it's very safe. It uses column-pivoting QR.
   * \li FullPivHouseholderQRPreconditioner, is the safest and slowest. It uses full-pivoting QR.
@@ -261,14 +296,13 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
     /** \brief Default Constructor with memory preallocation
       *
       * Like the default constructor but with preallocation of the internal data
-      * according to the specified problem \a size.
+      * according to the specified problem size.
       * \sa JacobiSVD()
       */
-    JacobiSVD(Index rows, Index cols) : m_matrixU(rows, rows),
-                                    m_matrixV(cols, cols),
-                                    m_singularValues(std::min(rows, cols)),
-                                    m_workMatrix(rows, cols),
-                                    m_isInitialized(false) {}
+    JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
+    {
+      allocate(rows, cols, computationOptions);
+    }
 
     /** \brief Constructor performing the decomposition of given matrix.
      *
@@ -286,15 +320,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
      * Thin unitaries also are only available if your matrix type has a Dynamic number of columns (for example MatrixXf).
      */
     JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
-        : m_matrixU(matrix.rows(), matrix.rows()),
-          m_matrixV(matrix.cols(), matrix.cols()),
-          m_singularValues(),
-          m_workMatrix(),
-          m_isInitialized(false)
     {
-      const Index minSize = std::min(matrix.rows(), matrix.cols());
-      m_singularValues.resize(minSize);
-      m_workMatrix.resize(minSize, minSize);
       compute(matrix, computationOptions);
     }
 
@@ -315,6 +341,15 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
      */
     JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions = 0);
 
+    /** \returns the \a U matrix.
+     *
+     * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p,
+     * the U matrix is n-by-n if you asked for ComputeFullU, and is n-by-m if you asked for ComputeThinU.
+     *
+     * The \a m first columns of \a U are the left singular vectors of the matrix being decomposed.
+     *
+     * This method asserts that you asked for \a U to be computed.
+     */
     const MatrixUType& matrixU() const
     {
       ei_assert(m_isInitialized && "JacobiSVD is not initialized.");
@@ -322,27 +357,43 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
       return m_matrixU;
     }
 
-    const SingularValuesType& singularValues() const
+    /** \returns the \a V matrix.
+     *
+     * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p,
+     * the V matrix is p-by-p if you asked for ComputeFullV, and is p-by-m if you asked for ComputeThinV.
+     *
+     * The \a m first columns of \a V are the right singular vectors of the matrix being decomposed.
+     *
+     * This method asserts that you asked for \a V to be computed.
+     */
+    const MatrixVType& matrixV() const
     {
       ei_assert(m_isInitialized && "JacobiSVD is not initialized.");
-      return m_singularValues;
+      ei_assert(computeV() && "This JacobiSVD decomposition didn't compute V. Did you ask for it?");
+      return m_matrixV;
     }
 
-    const MatrixVType& matrixV() const
+    /** \returns the vector of singular values.
+     *
+     * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p, the
+     * returned vector has size \a m.
+     */
+    const SingularValuesType& singularValues() const
     {
       ei_assert(m_isInitialized && "JacobiSVD is not initialized.");
-      ei_assert(computeV() && "This JacobiSVD decomposition didn't compute V. Did you ask for it?");
-      return m_matrixV;
+      return m_singularValues;
     }
 
+    /** \returns true if \a U (full or thin) is asked for in this SVD decomposition */
     inline bool computeU() const { return m_computeFullU || m_computeThinU; }
+    /** \returns true if \a V (full or thin) is asked for in this SVD decomposition */
     inline bool computeV() const { return m_computeFullV || m_computeThinV; }
 
     /** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A.
       *
       * \param b the right-hand-side of the equation to solve.
       *
-      * \note Solving requires both U and V to be computed. Thin U and V are enough, there is no need for full U or V,
+      * \note Solving requires both U and V to be computed. Thin U and V are enough, there is no need for full U or V.
       * 
       * \note SVD solving is implicitly least-squares. Thus, this method serves both purposes of exact solving and least-squares solving.
       * In other words, the returned solution is guaranteed to minimize the Euclidean norm \f$ \Vert A x - b \Vert \f$.
@@ -356,6 +407,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
       return ei_solve_retval<JacobiSVD, Rhs>(*this, b.derived());
     }
 
+    /** \returns the number of singular values that are not exactly 0 */
     Index nonzeroSingularValues() const
     {
       ei_assert(m_isInitialized && "JacobiSVD is not initialized.");
@@ -365,6 +417,9 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
     inline Index rows() const { return m_rows; }
     inline Index cols() const { return m_cols; }
 
+  private:
+    void allocate(Index rows, Index cols, unsigned int computationOptions = 0);
+
   protected:
     MatrixUType m_matrixU;
     MatrixVType m_matrixV;
@@ -373,7 +428,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
     bool m_isInitialized;
     bool m_computeFullU, m_computeThinU;
     bool m_computeFullV, m_computeThinV;
-    Index m_nonzeroSingularValues, m_rows, m_cols;
+    Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
 
     template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
     friend struct ei_svd_precondition_2x2_block_to_be_real;
@@ -382,6 +437,38 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
 };
 
 template<typename MatrixType, int QRPreconditioner>
+void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions)
+{
+  m_rows = rows;
+  m_cols = cols;
+  m_isInitialized = false;
+  m_computeFullU = computationOptions & ComputeFullU;
+  m_computeThinU = computationOptions & ComputeThinU;
+  m_computeFullV = computationOptions & ComputeFullV;
+  m_computeThinV = computationOptions & ComputeThinV;
+  ei_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
+  ei_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
+  ei_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
+              "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
+  if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
+  {
+      ei_assert(!(m_computeThinU || m_computeThinV) &&
+              "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
+              "Use the ColPivHouseholderQR preconditioner instead.");
+  }
+  m_diagSize = std::min(m_rows, m_cols);
+  m_singularValues.resize(m_diagSize);
+  m_matrixU.resize(m_rows, m_computeFullU ? m_rows
+                          : m_computeThinU ? m_diagSize
+                          : 0);
+  m_matrixV.resize(m_cols, m_computeFullV ? m_cols
+                          : m_computeThinV ? m_diagSize
+                          : 0);
+  m_workMatrix.resize(m_diagSize, m_diagSize);
+}
+
+
+template<typename MatrixType, int QRPreconditioner>
 struct ei_svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
 {
   typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
@@ -463,53 +550,51 @@ template<typename MatrixType, int QRPreconditioner>
 JacobiSVD<MatrixType, QRPreconditioner>&
 JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions)
 {
-  m_computeFullU = computationOptions & ComputeFullU;
-  m_computeThinU = computationOptions & ComputeThinU;
-  m_computeFullV = computationOptions & ComputeFullV;
-  m_computeThinV = computationOptions & ComputeThinV;
-  ei_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U");
-  ei_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V");
-  ei_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) &&
-            "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.");
-  if (QRPreconditioner == FullPivHouseholderQRPreconditioner)
-  {
-    ei_assert(!(m_computeThinU || m_computeThinV) &&
-              "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
-              "Use the ColPivHouseholderQR preconditioner instead.");
-  }
-  m_rows = matrix.rows();
-  m_cols = matrix.cols();
-  Index diagSize = std::min(m_rows, m_cols);
-  m_singularValues.resize(diagSize);
+  allocate(matrix.rows(), matrix.cols(), computationOptions);
+
+  // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations,
+  // only worsening the precision of U and V as we accumulate more rotations
   const RealScalar precision = 2 * NumTraits<Scalar>::epsilon();
 
+  /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
+
   if(!ei_qr_preconditioner_impl<MatrixType, QRPreconditioner, PreconditionIfMoreColsThanRows>::run(*this, matrix)
   && !ei_qr_preconditioner_impl<MatrixType, QRPreconditioner, PreconditionIfMoreRowsThanCols>::run(*this, matrix))
   {
-    m_workMatrix = matrix.block(0,0,diagSize,diagSize);
+    m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize);
     if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
-    if(m_computeThinU) m_matrixU.setIdentity(m_rows,diagSize);
+    if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
     if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
-    if(m_computeThinV) m_matrixV.setIdentity(diagSize,m_cols);
+    if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
   }
 
+  /*** step 2. The main Jacobi SVD iteration. ***/
+
   bool finished = false;
   while(!finished)
   {
     finished = true;
-    for(Index p = 1; p < diagSize; ++p)
+
+    // do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
+
+    for(Index p = 1; p < m_diagSize; ++p)
     {
       for(Index q = 0; q < p; ++q)
       {
+        // if this 2x2 sub-matrix is not diagonal already...
+        // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't
+        // keep us iterating forever.
         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, QRPreconditioner>::run(m_workMatrix, *this, p, q);
 
+          // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal
+          ei_svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q);
           PlanarRotation<RealScalar> j_left, j_right;
           ei_real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
 
+          // accumulate resulting Jacobi rotations
           m_workMatrix.applyOnTheLeft(p,q,j_left);
           if(computeU()) m_matrixU.applyOnTheRight(p,q,j_left.transpose());
 
@@ -520,19 +605,22 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
     }
   }
 
-  for(Index i = 0; i < diagSize; ++i)
+  /*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/
+
+  for(Index i = 0; i < m_diagSize; ++i)
   {
     RealScalar a = ei_abs(m_workMatrix.coeff(i,i));
     m_singularValues.coeffRef(i) = a;
     if(computeU() && (a!=RealScalar(0))) m_matrixU.col(i) *= m_workMatrix.coeff(i,i)/a;
   }
 
-  m_nonzeroSingularValues = diagSize;
+  /*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
 
-  for(Index i = 0; i < diagSize; i++)
+  m_nonzeroSingularValues = m_diagSize;
+  for(Index i = 0; i < m_diagSize; i++)
   {
     Index pos;
-    RealScalar maxRemainingSingularValue = m_singularValues.tail(diagSize-i).maxCoeff(&pos);
+    RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
     if(maxRemainingSingularValue == RealScalar(0))
     {
       m_nonzeroSingularValues = i;
diff --git a/doc/snippets/JacobiSVD_basic.cpp b/doc/snippets/JacobiSVD_basic.cpp
new file mode 100644
index 000000000..ab24b9bca
--- /dev/null
+++ b/doc/snippets/JacobiSVD_basic.cpp
@@ -0,0 +1,9 @@
+MatrixXf m = MatrixXf::Random(3,2);
+cout << "Here is the matrix m:" << endl << m << endl;
+JacobiSVD<MatrixXf> svd(m, ComputeThinU | ComputeThinV);
+cout << "Its singular values are:" << endl << svd.singularValues() << endl;
+cout << "Its left singular vectors are the columns of the thin U matrix:" << endl << svd.matrixU() << endl;
+cout << "Its right singular vectors are the columns of the thin V matrix:" << endl << svd.matrixV() << endl;
+Vector3f rhs(1, 0, 0);
+cout << "Now consider this rhs vector:" << endl << rhs << endl;
+cout << "A least-squares solution of m*x = rhs is:" << endl << svd.solve(rhs) << endl;
diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp
index a6dbcf2e8..7d8c36308 100644
--- a/test/jacobisvd.cpp
+++ b/test/jacobisvd.cpp
@@ -196,6 +196,23 @@ template<typename MatrixType> void jacobisvd_verify_assert(const MatrixType& m)
   }
 }
 
+template<typename MatrixType>
+void jacobisvd_inf_nan()
+{
+  JacobiSVD<MatrixType> svd;
+  typedef typename MatrixType::Scalar Scalar;
+  Scalar some_inf = Scalar(1) / Scalar(0);
+  svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV);
+  Scalar some_nan = Scalar(0) / Scalar(0);
+  svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV);
+  MatrixType m = MatrixType::Zero(10,10);
+  m(ei_random<int>(0,9), ei_random<int>(0,9)) = some_inf;
+  svd.compute(m, ComputeFullU | ComputeFullV);
+  m = MatrixType::Zero(10,10);
+  m(ei_random<int>(0,9), ei_random<int>(0,9)) = some_nan;
+  svd.compute(m, ComputeFullU | ComputeFullV);
+}
+
 void test_jacobisvd()
 {
   CALL_SUBTEST_3(( jacobisvd_verify_assert(Matrix3f()) ));
@@ -211,10 +228,15 @@ void test_jacobisvd()
     m << 1, 0,
          1, 0;
     CALL_SUBTEST_1(( jacobisvd(m, false) ));
+
     Matrix2d n;
-    n << 1, 1,
-         1, -1;
+    n << 0, 0,
+         0, 0;
+    CALL_SUBTEST_2(( jacobisvd(n, false) ));
+    n << 0, 0,
+         0, 1;
     CALL_SUBTEST_2(( jacobisvd(n, false) ));
+    
     CALL_SUBTEST_3(( jacobisvd<Matrix3f>() ));
     CALL_SUBTEST_4(( jacobisvd<Matrix4d>() ));
     CALL_SUBTEST_5(( jacobisvd<Matrix<float,3,5> >() ));
@@ -226,8 +248,14 @@ void test_jacobisvd()
     CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(r,c)) ));
     (void) r;
     (void) c;
+
+    // Test on inf/nan matrix
+    CALL_SUBTEST_7( jacobisvd_inf_nan<MatrixXf>() );
   }
 
-  CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(ei_random<int>(200, 300), ei_random<int>(200, 300))) ));
-  CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(ei_random<int>(100, 120), ei_random<int>(100, 120))) ));
+  CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(ei_random<int>(100, 150), ei_random<int>(100, 150))) ));
+  CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(ei_random<int>(80, 100), ei_random<int>(80, 100))) ));
+
+  // Test problem size constructors
+  CALL_SUBTEST_7( JacobiSVD<MatrixXf>(10,10) );
 }