D&C SVD: add scaling to avoid overflow, fix handling of fixed size matrices

test/svd_common.h

diff --git a/test/svd_common.h b/test/svd_common.h
index 4631939..b880efc 100644
--- a/test/svd_common.h
+++ b/test/svd_common.h
@@ -301,7 +301,7 @@
                  0, 5.60844e-313;
   SVD_DEFAULT(Matrix2d) svd;
   svd.compute(M,ComputeFullU|ComputeFullV);
-  svd_check_full(M,svd);
+  CALL_SUBTEST( svd_check_full(M,svd) );
   
   // Check all 2x2 matrices made with the following coefficients:
   VectorXd value_set(9);
@@ -312,7 +312,7 @@
   {
     M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3));
     svd.compute(M,ComputeFullU|ComputeFullV);
-    svd_check_full(M,svd);
+    CALL_SUBTEST( svd_check_full(M,svd) );
     
     id(k)++;
     if(id(k)>=value_set.size())
@@ -336,7 +336,7 @@
 
   SVD_DEFAULT(Matrix3d) svd3;
   svd3.compute(M3,ComputeFullU|ComputeFullV); // just check we don't loop indefinitely
-  svd_check_full(M3,svd3);
+  CALL_SUBTEST( svd_check_full(M3,svd3) );
 }
 
 // void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true)

unsupported/Eigen/src/BDCSVD/BDCSVD.h

diff --git a/unsupported/Eigen/src/BDCSVD/BDCSVD.h b/unsupported/Eigen/src/BDCSVD/BDCSVD.h
index b13ee41..5bf9b0a 100644
--- a/unsupported/Eigen/src/BDCSVD/BDCSVD.h
+++ b/unsupported/Eigen/src/BDCSVD/BDCSVD.h
@@ -236,26 +236,29 @@
 {
   allocate(matrix.rows(), matrix.cols(), computationOptions);
   using std::abs;
-
-  //**** step 1 Bidiagonalization
-  MatrixType copy;
-  if (m_isTranspose) copy = matrix.adjoint();
-  else               copy = matrix;
   
+  //**** step 0 - Copy the input matrix and apply scaling to reduce over/under-flows
+  RealScalar scale = matrix.cwiseAbs().maxCoeff();
+  if(scale==RealScalar(0)) scale = RealScalar(1);
+  MatrixX copy;
+  if (m_isTranspose) copy = matrix.adjoint()/scale;
+  else               copy = matrix/scale;
+  
+  //**** step 1 - Bidiagonalization
   internal::UpperBidiagonalization<MatrixX> bid(copy);
 
-  //**** step 2 Divide
+  //**** step 2 - Divide & Conquer
   m_naiveU.setZero();
   m_naiveV.setZero();
   m_computed.topRows(m_diagSize) = bid.bidiagonal().toDenseMatrix().transpose();
   m_computed.template bottomRows<1>().setZero();
   divide(0, m_diagSize - 1, 0, 0, 0);
 
-  //**** step 3 copy
+  //**** step 3 - Copy singular values and vectors
   for (int i=0; i<m_diagSize; i++)
   {
     RealScalar a = abs(m_computed.coeff(i, i));
-    m_singularValues.coeffRef(i) = a;
+    m_singularValues.coeffRef(i) = a * scale;
     if (a == 0)
     {
       m_nonzeroSingularValues = i;
@@ -268,11 +271,13 @@
       break;
     }
   }
-//   std::cout << "m_naiveU\n" << m_naiveU << "\n\n";
-//   std::cout << "m_naiveV\n" << m_naiveV << "\n\n";
+#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
+  std::cout << "m_naiveU\n" << m_naiveU << "\n\n";
+  std::cout << "m_naiveV\n" << m_naiveV << "\n\n";
+#endif
   if(m_isTranspose) copyUV(bid.householderV(), bid.householderU(), m_naiveV, m_naiveU);
   else              copyUV(bid.householderU(), bid.householderV(), m_naiveU, m_naiveV);
-//   std::cout << "DONE\n";
+
   m_isInitialized = true;
   return *this;
 }// end compute
@@ -569,13 +574,13 @@
   computeSingVecs(zhat, diag, singVals, shifts, mus, U, V);
   
 #ifdef  EIGEN_BDCSVD_DEBUG_VERBOSE
-  std::cout << "U^T U: " << (U.transpose() * U - MatrixType(MatrixType::Identity(U.cols(),U.cols()))).norm() << "\n";
-  std::cout << "V^T V: " << (V.transpose() * V - MatrixType(MatrixType::Identity(V.cols(),V.cols()))).norm() << "\n";
+  std::cout << "U^T U: " << (U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() << "\n";
+  std::cout << "V^T V: " << (V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() << "\n";
 #endif
   
 #ifdef EIGEN_BDCSVD_SANITY_CHECKS
-  assert((U.transpose() * U - MatrixType(MatrixType::Identity(U.cols(),U.cols()))).norm() < 1e-14 * n);
-  assert((V.transpose() * V - MatrixType(MatrixType::Identity(V.cols(),V.cols()))).norm() < 1e-14 * n);
+  assert((U.transpose() * U - MatrixXr(MatrixXr::Identity(U.cols(),U.cols()))).norm() < 1e-14 * n);
+  assert((V.transpose() * V - MatrixXr(MatrixXr::Identity(V.cols(),V.cols()))).norm() < 1e-14 * n);
   assert(m_naiveU.allFinite());
   assert(m_naiveV.allFinite());
   assert(m_computed.allFinite());

unsupported/test/bdcsvd.cpp

diff --git a/unsupported/test/bdcsvd.cpp b/unsupported/test/bdcsvd.cpp
index 4ad9915..95cdb6a 100644
--- a/unsupported/test/bdcsvd.cpp
+++ b/unsupported/test/bdcsvd.cpp
@@ -70,13 +70,13 @@
   CALL_SUBTEST_7(( svd_verify_assert<BDCSVD<MatrixXf>  >(MatrixXf(10,12)) ));
   CALL_SUBTEST_8(( svd_verify_assert<BDCSVD<MatrixXcd> >(MatrixXcd(7,5)) ));
   
-//   svd_all_trivial_2x2(bdcsvd<Matrix2cd>);
-//   svd_all_trivial_2x2(bdcsvd<Matrix2d>);
+  CALL_SUBTEST_1(( svd_all_trivial_2x2(bdcsvd<Matrix2cd>) ));
+  CALL_SUBTEST_1(( svd_all_trivial_2x2(bdcsvd<Matrix2d>) ));
 
   for(int i = 0; i < g_repeat; i++) {
-//     CALL_SUBTEST_3(( bdcsvd<Matrix3f>() ));
-//     CALL_SUBTEST_4(( bdcsvd<Matrix4d>() ));
-//     CALL_SUBTEST_5(( bdcsvd<Matrix<float,3,5> >() ));
+    CALL_SUBTEST_3(( bdcsvd<Matrix3f>() ));
+    CALL_SUBTEST_4(( bdcsvd<Matrix4d>() ));
+    CALL_SUBTEST_5(( bdcsvd<Matrix<float,3,5> >() ));
 
     int r = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2),
         c = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2);