Move D&C SVD to official SVD module.

2 files changed, 112 insertions, 0 deletions

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 63ed2c7a4..f57d8ce36 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -211,6 +211,7 @@ ei_add_test(real_qz)
 ei_add_test(eigensolver_generalized_real)
 ei_add_test(jacobi)
 ei_add_test(jacobisvd)
+ei_add_test(bdcsvd)
 ei_add_test(householder)
 ei_add_test(geo_orthomethods)
 ei_add_test(geo_quaternion)
diff --git a/test/bdcsvd.cpp b/test/bdcsvd.cpp
new file mode 100644
index 000000000..d58d259d6
--- /dev/null
+++ b/test/bdcsvd.cpp
@@ -0,0 +1,111 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com>
+// Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr>
+// Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr>
+// Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/
+
+// discard stack allocation as that too bypasses malloc
+#define EIGEN_STACK_ALLOCATION_LIMIT 0
+#define EIGEN_RUNTIME_NO_MALLOC
+
+#include "main.h"
+#include <Eigen/SVD>
+#include <iostream>
+#include <Eigen/LU>
+
+
+#define SVD_DEFAULT(M) BDCSVD<M>
+#define SVD_FOR_MIN_NORM(M) BDCSVD<M>
+#include "../../test/svd_common.h"
+
+// Check all variants of JacobiSVD
+template<typename MatrixType>
+void bdcsvd(const MatrixType& a = MatrixType(), bool pickrandom = true)
+{
+  MatrixType m = a;
+  if(pickrandom)
+    svd_fill_random(m);
+
+  CALL_SUBTEST(( svd_test_all_computation_options<BDCSVD<MatrixType> >(m, false)  ));
+}
+
+template<typename MatrixType>
+void bdcsvd_method()
+{
+  enum { Size = MatrixType::RowsAtCompileTime };
+  typedef typename MatrixType::RealScalar RealScalar;
+  typedef Matrix<RealScalar, Size, 1> RealVecType;
+  MatrixType m = MatrixType::Identity();
+  VERIFY_IS_APPROX(m.bdcSvd().singularValues(), RealVecType::Ones());
+  VERIFY_RAISES_ASSERT(m.bdcSvd().matrixU());
+  VERIFY_RAISES_ASSERT(m.bdcSvd().matrixV());
+  VERIFY_IS_APPROX(m.bdcSvd(ComputeFullU|ComputeFullV).solve(m), m);
+}
+
+// compare the Singular values returned with Jacobi and Bdc
+template<typename MatrixType> 
+void compare_bdc_jacobi(const MatrixType& a = MatrixType(), unsigned int computationOptions = 0)
+{
+  MatrixType m = MatrixType::Random(a.rows(), a.cols());
+  BDCSVD<MatrixType> bdc_svd(m);
+  JacobiSVD<MatrixType> jacobi_svd(m);
+  VERIFY_IS_APPROX(bdc_svd.singularValues(), jacobi_svd.singularValues());
+  if(computationOptions & ComputeFullU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU());
+  if(computationOptions & ComputeThinU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU());
+  if(computationOptions & ComputeFullV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV());
+  if(computationOptions & ComputeThinV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV());
+}
+
+void test_bdcsvd()
+{
+  CALL_SUBTEST_3(( svd_verify_assert<BDCSVD<Matrix3f>  >(Matrix3f()) ));
+  CALL_SUBTEST_4(( svd_verify_assert<BDCSVD<Matrix4d>  >(Matrix4d()) ));
+  CALL_SUBTEST_7(( svd_verify_assert<BDCSVD<MatrixXf>  >(MatrixXf(10,12)) ));
+  CALL_SUBTEST_8(( svd_verify_assert<BDCSVD<MatrixXcd> >(MatrixXcd(7,5)) ));
+  
+  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> >() ));
+
+    int r = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2),
+        c = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2);
+    
+    TEST_SET_BUT_UNUSED_VARIABLE(r)
+    TEST_SET_BUT_UNUSED_VARIABLE(c)
+    
+    CALL_SUBTEST_6((  bdcsvd(Matrix<double,Dynamic,2>(r,2)) ));
+    CALL_SUBTEST_7((  bdcsvd(MatrixXf(r,c)) ));
+    CALL_SUBTEST_7((  compare_bdc_jacobi(MatrixXf(r,c)) ));
+    CALL_SUBTEST_10(( bdcsvd(MatrixXd(r,c)) ));
+    CALL_SUBTEST_10(( compare_bdc_jacobi(MatrixXd(r,c)) ));
+    CALL_SUBTEST_8((  bdcsvd(MatrixXcd(r,c)) ));
+    CALL_SUBTEST_8((  compare_bdc_jacobi(MatrixXcd(r,c)) ));
+
+    // Test on inf/nan matrix
+    CALL_SUBTEST_7(  (svd_inf_nan<BDCSVD<MatrixXf>, MatrixXf>()) );
+    CALL_SUBTEST_10( (svd_inf_nan<BDCSVD<MatrixXd>, MatrixXd>()) );
+  }
+
+  // test matrixbase method
+  CALL_SUBTEST_1(( bdcsvd_method<Matrix2cd>() ));
+  CALL_SUBTEST_3(( bdcsvd_method<Matrix3f>() ));
+
+  // Test problem size constructors
+  CALL_SUBTEST_7( BDCSVD<MatrixXf>(10,10) );
+
+  // Check that preallocation avoids subsequent mallocs
+  CALL_SUBTEST_9( svd_preallocate() );
+
+  CALL_SUBTEST_2( svd_underoverflow() );
+}
+