Move D&C SVD to official SVD module.

Eigen/SVD

@@ -12,20 +12,25 @@
   *
   *
   * This module provides SVD decomposition for matrices (both real and complex).
-  * This decomposition is accessible via the following MatrixBase method:
+  * Two decomposition algorithms are provided:
+  *  - JacobiSVD implementing two-sided Jacobi iterations is numerically very accurate, fast for small matrices, but very slow for larger ones.
+  *  - BDCSVD implementing a recursive divide & conquer strategy on top of an upper-bidiagonalization which remains fast for large problems.
+  * These decompositions are accessible via the respective classes and following MatrixBase methods:
   *  - MatrixBase::jacobiSvd()
+  *  - MatrixBase::bdcSvd()
   *
   * \code
   * #include <Eigen/SVD>
   * \endcode
   */
 
+#include "src/SVD/UpperBidiagonalization.h"
 #include "src/SVD/SVDBase.h"
 #include "src/SVD/JacobiSVD.h"
+#include "src/SVD/BDCSVD.h"
 #if defined(EIGEN_USE_LAPACKE) && !defined(EIGEN_USE_LAPACKE_STRICT)
 #include "src/SVD/JacobiSVD_MKL.h"
 #endif
-#include "src/SVD/UpperBidiagonalization.h"
 
 #include "src/Core/util/ReenableStupidWarnings.h"
 

Eigen/src/Core/MatrixBase.h

@@ -390,6 +390,7 @@ template<typename Derived> class MatrixBase
 /////////// SVD module ///////////
 
     JacobiSVD<PlainObject> jacobiSvd(unsigned int computationOptions = 0) const;
+    BDCSVD<PlainObject>    bdcSvd(unsigned int computationOptions = 0) const;
 
 /////////// Geometry module ///////////
 

Eigen/src/Core/util/ForwardDeclarations.h

@@ -251,6 +251,7 @@ template<typename MatrixType> class HouseholderQR;
 template<typename MatrixType> class ColPivHouseholderQR;
 template<typename MatrixType> class FullPivHouseholderQR;
 template<typename MatrixType, int QRPreconditioner = ColPivHouseholderQRPreconditioner> class JacobiSVD;
+template<typename MatrixType> class BDCSVD;
 template<typename MatrixType, int UpLo = Lower> class LLT;
 template<typename MatrixType, int UpLo = Lower> class LDLT;
 template<typename VectorsType, typename CoeffsType, int Side=OnTheLeft> class HouseholderSequence;

Eigen/src/SVD/BDCSVD.h

@@ -514,6 +514,7 @@ void BDCSVD<MatrixType>::divide (Index firstCol, Index lastCol, Index firstRowW,
 //
 // TODO Opportunities for optimization: better root finding algo, better stopping criterion, better
 // handling of round-off errors, be consistent in ordering
+// For instance, to solve the secular equation using FMM, see http://www.stat.uchicago.edu/~lekheng/courses/302/classics/greengard-rokhlin.pdf
 template <typename MatrixType>
 void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V)
 {
@@ -528,7 +529,10 @@ void BDCSVD<MatrixType>::computeSVDofM(Index firstCol, Index n, MatrixXr& U, Vec
   U.resize(n+1, n+1);
   if (m_compV) V.resize(n, n);
 
-  if (col0.hasNaN() || diag.hasNaN()) { std::cout << "\n\nHAS NAN\n\n"; return; }
+#ifdef EIGEN_BDCSVD_DEBUG_VERBOSE
+  if (col0.hasNaN() || diag.hasNaN())
+    std::cout << "\n\nHAS NAN\n\n";
+#endif
   
   // Many singular values might have been deflated, the zero ones have been moved to the end,
   // but others are interleaved and we must ignore them at this stage.
@@ -1149,21 +1153,18 @@ void BDCSVD<MatrixType>::deflation(Index firstCol, Index lastCol, Index k, Index
 }//end deflation
 
 
-  /** \svd_module
+/** \svd_module
   *
-   * \return the singular value decomposition of \c *this computed by 
+  * \return the singular value decomposition of \c *this computed by Divide & Conquer algorithm
-   *  BDC Algorithm
   *
   * \sa class BDCSVD
   */
-/*
 template<typename Derived>
 BDCSVD<typename MatrixBase<Derived>::PlainObject>
 MatrixBase<Derived>::bdcSvd(unsigned int computationOptions) const
 {
   return BDCSVD<PlainObject>(*this, computationOptions);
 }
-*/
 
 } // end namespace Eigen
 

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)

test/bdcsvd.cpp

@@ -15,7 +15,7 @@
 #define EIGEN_RUNTIME_NO_MALLOC
 
 #include "main.h"
-#include <unsupported/Eigen/BDCSVD>
+#include <Eigen/SVD>
 #include <iostream>
 #include <Eigen/LU>
 
@@ -35,18 +35,18 @@ void bdcsvd(const MatrixType& a = MatrixType(), bool pickrandom = true)
   CALL_SUBTEST(( svd_test_all_computation_options<BDCSVD<MatrixType> >(m, false)  ));
 }
 
-// template<typename MatrixType>
+template<typename MatrixType>
-// void bdcsvd_method()
+void bdcsvd_method()
-// {
+{
-//   enum { Size = MatrixType::RowsAtCompileTime };
+  enum { Size = MatrixType::RowsAtCompileTime };
-//   typedef typename MatrixType::RealScalar RealScalar;
+  typedef typename MatrixType::RealScalar RealScalar;
-//   typedef Matrix<RealScalar, Size, 1> RealVecType;
+  typedef Matrix<RealScalar, Size, 1> RealVecType;
-//   MatrixType m = MatrixType::Identity();
+  MatrixType m = MatrixType::Identity();
-//   VERIFY_IS_APPROX(m.bdcSvd().singularValues(), RealVecType::Ones());
+  VERIFY_IS_APPROX(m.bdcSvd().singularValues(), RealVecType::Ones());
-//   VERIFY_RAISES_ASSERT(m.bdcSvd().matrixU());
+  VERIFY_RAISES_ASSERT(m.bdcSvd().matrixU());
-//   VERIFY_RAISES_ASSERT(m.bdcSvd().matrixV());
+  VERIFY_RAISES_ASSERT(m.bdcSvd().matrixV());
-//   VERIFY_IS_APPROX(m.bdcSvd(ComputeFullU|ComputeFullV).solve(m), m);
+  VERIFY_IS_APPROX(m.bdcSvd(ComputeFullU|ComputeFullV).solve(m), m);
-// }
+}
 
 // compare the Singular values returned with Jacobi and Bdc
 template<typename MatrixType> 
@@ -97,8 +97,8 @@ void test_bdcsvd()
   }
 
   // test matrixbase method
-//   CALL_SUBTEST_1(( bdcsvd_method<Matrix2cd>() ));
+  CALL_SUBTEST_1(( bdcsvd_method<Matrix2cd>() ));
-//   CALL_SUBTEST_3(( bdcsvd_method<Matrix3f>() ));
+  CALL_SUBTEST_3(( bdcsvd_method<Matrix3f>() ));
 
   // Test problem size constructors
   CALL_SUBTEST_7( BDCSVD<MatrixXf>(10,10) );

unsupported/Eigen/BDCSVD

@@ -1,26 +0,0 @@
-#ifndef EIGEN_BDCSVD_MODULE_H
-#define EIGEN_BDCSVD_MODULE_H
-
-#include <Eigen/SVD>
-
-#include "../../Eigen/src/Core/util/DisableStupidWarnings.h"
-
-/** \defgroup BDCSVD_Module BDCSVD module
-  *
-  *
-  *
-  * This module provides Divide & Conquer SVD decomposition for matrices (both real and complex).
-  * This decomposition is accessible via the following MatrixBase method:
-  *  - MatrixBase::bdcSvd()
-  *
-  * \code
-  * #include <Eigen/BDCSVD>
-  * \endcode
-  */
-
-#include "src/BDCSVD/BDCSVD.h"
-
-#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
-
-#endif // EIGEN_BDCSVD_MODULE_H
-/* vim: set filetype=cpp et sw=2 ts=2 ai: */

unsupported/Eigen/src/BDCSVD/CMakeLists.txt

@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_BDCSVD_SRCS "*.h")
-
-INSTALL(FILES
-  ${Eigen_BDCSVD_SRCS}
-  DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/BDCSVD COMPONENT Devel
-  )

unsupported/Eigen/src/BDCSVD/TODOBdcsvd.txt

@@ -1,13 +0,0 @@
-TO DO LIST
-
-- check more carefully single precision
-- check with duplicated singularvalues
-- no-malloc mode
-
-(optional optimization)
-- do all the allocations in the allocate part 
-- support static matrices
-- return a error at compilation time when using integer matrices (int, long, std::complex<int>, ...)
-- To solve the secular equation using FMM:
-http://www.stat.uchicago.edu/~lekheng/courses/302/classics/greengard-rokhlin.pdf
-

unsupported/Eigen/src/BDCSVD/doneInBDCSVD.txt

@@ -1,8 +0,0 @@
-This unsupported package is about a divide and conquer algorithm to compute SVD.
-
-The implementation follows as closely as possible the following reference paper : 
-http://www.cs.yale.edu/publications/techreports/tr933.pdf
-
-To solve the secular equation using FMM:
-http://www.stat.uchicago.edu/~lekheng/courses/302/classics/greengard-rokhlin.pdf
-

unsupported/Eigen/src/CMakeLists.txt

@@ -12,4 +12,3 @@ ADD_SUBDIRECTORY(Skyline)
 ADD_SUBDIRECTORY(SparseExtra)
 ADD_SUBDIRECTORY(KroneckerProduct)
 ADD_SUBDIRECTORY(Splines)
-ADD_SUBDIRECTORY(BDCSVD)

unsupported/test/CMakeLists.txt

@@ -92,7 +92,6 @@ ei_add_test(splines)
 ei_add_test(gmres)
 ei_add_test(minres)
 ei_add_test(levenberg_marquardt)
-ei_add_test(bdcsvd)
 ei_add_test(kronecker_product)
 
 option(EIGEN_TEST_CXX11 "Enable testing of C++11 features (e.g. Tensor module)." OFF)