bug #843: fix jacobisvd for complexes and extend respective unit test to chack with random tricky matrices

1 files changed, 67 insertions, 21 deletions

diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp
index d441a6eca..36721b496 100644
--- a/test/jacobisvd.cpp
+++ b/test/jacobisvd.cpp
@@ -67,6 +67,7 @@ template<typename MatrixType, int QRPreconditioner>
 void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions)
 {
   typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
   typedef typename MatrixType::Index Index;
   Index rows = m.rows();
   Index cols = m.cols();
@@ -81,9 +82,37 @@ void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions)
 
   RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols));
   JacobiSVD<MatrixType, QRPreconditioner> svd(m, computationOptions);
+
+       if(internal::is_same<RealScalar,double>::value) svd.setThreshold(1e-8);
+  else if(internal::is_same<RealScalar,float>::value)  svd.setThreshold(1e-4);
+  
   SolutionType x = svd.solve(rhs);
+  
+  RealScalar residual = (m*x-rhs).norm();
+  // Check that there is no significantly better solution in the neighborhood of x
+  if(!test_isMuchSmallerThan(residual,rhs.norm()))
+  {
+    // If the residual is very small, then we have an exact solution, so we are already good.
+    for(int k=0;k<x.rows();++k)
+    {
+      SolutionType y(x);
+      y.row(k).array() += 2*NumTraits<RealScalar>::epsilon();
+      RealScalar residual_y = (m*y-rhs).norm();
+      VERIFY( test_isApprox(residual_y,residual) || residual < residual_y );
+      
+      y.row(k) = x.row(k).array() - 2*NumTraits<RealScalar>::epsilon();
+      residual_y = (m*y-rhs).norm();
+      VERIFY( test_isApprox(residual_y,residual) || residual < residual_y );
+    }
+  }
+  
   // evaluate normal equation which works also for least-squares solutions
-  VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs);
+  if(internal::is_same<RealScalar,double>::value)
+  {
+    // This test is not stable with single precision.
+    // This is probably because squaring m signicantly affects the precision.
+    VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs);
+  }
   
   // check minimal norm solutions
   {
@@ -139,10 +168,9 @@ void jacobisvd_test_all_computation_options(const MatrixType& m)
   if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols())
     return;
   JacobiSVD<MatrixType, QRPreconditioner> fullSvd(m, ComputeFullU|ComputeFullV);
-
-  jacobisvd_check_full(m, fullSvd);
-  jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeFullV);
-
+  CALL_SUBTEST(( jacobisvd_check_full(m, fullSvd) ));
+  CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeFullV) ));
+  
   #if defined __INTEL_COMPILER
   // remark #111: statement is unreachable
   #pragma warning disable 111
@@ -150,20 +178,20 @@ void jacobisvd_test_all_computation_options(const MatrixType& m)
   if(QRPreconditioner == FullPivHouseholderQRPreconditioner)
     return;
 
-  jacobisvd_compare_to_full(m, ComputeFullU, fullSvd);
-  jacobisvd_compare_to_full(m, ComputeFullV, fullSvd);
-  jacobisvd_compare_to_full(m, 0, fullSvd);
+  CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU, fullSvd) ));
+  CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullV, fullSvd) ));
+  CALL_SUBTEST(( jacobisvd_compare_to_full(m, 0, fullSvd) ));
 
   if (MatrixType::ColsAtCompileTime == Dynamic) {
     // thin U/V are only available with dynamic number of columns
-    jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd);
-    jacobisvd_compare_to_full(m,              ComputeThinV, fullSvd);
-    jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd);
-    jacobisvd_compare_to_full(m, ComputeThinU             , fullSvd);
-    jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd);
-    jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeThinV);
-    jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeFullV);
-    jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeThinV);
+    CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd) ));
+    CALL_SUBTEST(( jacobisvd_compare_to_full(m,              ComputeThinV, fullSvd) ));
+    CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd) ));
+    CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU             , fullSvd) ));
+    CALL_SUBTEST(( jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd) ));
+    CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeThinV) ));
+    CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeFullV) ));
+    CALL_SUBTEST(( jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeThinV) ));
 
     // test reconstruction
     typedef typename MatrixType::Index Index;
@@ -176,12 +204,29 @@ void jacobisvd_test_all_computation_options(const MatrixType& m)
 template<typename MatrixType>
 void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true)
 {
-  MatrixType m = pickrandom ? MatrixType::Random(a.rows(), a.cols()) : a;
+  MatrixType m = a;
+  if(pickrandom)
+  {
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    typedef typename MatrixType::Index Index;
+    Index diagSize = (std::min)(a.rows(), a.cols());
+    RealScalar s = std::numeric_limits<RealScalar>::max_exponent10/4;
+    s = internal::random<RealScalar>(1,s);
+    Matrix<RealScalar,Dynamic,1> d =  Matrix<RealScalar,Dynamic,1>::Random(diagSize);
+    for(Index k=0; k<diagSize; ++k)
+      d(k) = d(k)*std::pow(RealScalar(10),internal::random<RealScalar>(-s,s));
+    m = Matrix<Scalar,Dynamic,Dynamic>::Random(a.rows(),diagSize) * d.asDiagonal() * Matrix<Scalar,Dynamic,Dynamic>::Random(diagSize,a.cols());
+    // cancel some coeffs
+    Index n  = internal::random<Index>(0,m.size()-1);
+    for(Index i=0; i<n; ++i)
+      m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = Scalar(0);
+  }
 
-  jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m);
-  jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m);
-  jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m);
-  jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m);
+  CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m) ));
+  CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m) ));
+  CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m) ));
+  CALL_SUBTEST(( jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m) ));
 }
 
 template<typename MatrixType> void jacobisvd_verify_assert(const MatrixType& m)
@@ -384,6 +429,7 @@ void test_jacobisvd()
     TEST_SET_BUT_UNUSED_VARIABLE(r)
     TEST_SET_BUT_UNUSED_VARIABLE(c)
     
+    CALL_SUBTEST_10(( jacobisvd<MatrixXd>(MatrixXd(r,c)) ));
     CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(r,c)) ));
     CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(r,c)) ));
     (void) r;