diff options
author | 2014-09-19 15:25:48 +0200 | |
---|---|---|
committer | 2014-09-19 15:25:48 +0200 | |
commit | 03dd4dd91a5d8963f56eebe3b9d2eb924bc06e02 (patch) | |
tree | 74293b3b60724f704022aa2e13c7612c6eb6df9b /unsupported | |
parent | 0a18eecab332d0dd87154b9eef7ff993a4bb625c (diff) |
Unify unit test for BDC and Jacobi SVD. This reveals some numerical issues in BDCSVD.
Diffstat (limited to 'unsupported')
-rw-r--r-- | unsupported/test/bdcsvd.cpp | 233 | ||||
-rw-r--r-- | unsupported/test/jacobisvd.cpp | 198 | ||||
-rw-r--r-- | unsupported/test/svd_common.h | 261 |
3 files changed, 67 insertions, 625 deletions
diff --git a/unsupported/test/bdcsvd.cpp b/unsupported/test/bdcsvd.cpp index 115a649b0..4ad991522 100644 --- a/unsupported/test/bdcsvd.cpp +++ b/unsupported/test/bdcsvd.cpp @@ -10,204 +10,105 @@ // 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/ -#include "svd_common.h" +// discard stack allocation as that too bypasses malloc +#define EIGEN_STACK_ALLOCATION_LIMIT 0 +#define EIGEN_RUNTIME_NO_MALLOC + +#include "main.h" +#include <unsupported/Eigen/BDCSVD> #include <iostream> #include <Eigen/LU> -// check if "svd" is the good image of "m" -template<typename MatrixType> -void bdcsvd_check_full(const MatrixType& m, const BDCSVD<MatrixType>& svd) -{ - svd_check_full< MatrixType, BDCSVD< MatrixType > >(m, svd); -} - -// Compare to a reference value -template<typename MatrixType> -void bdcsvd_compare_to_full(const MatrixType& m, - unsigned int computationOptions, - const BDCSVD<MatrixType>& referenceSvd) -{ - svd_compare_to_full< MatrixType, BDCSVD< MatrixType > >(m, computationOptions, referenceSvd); -} // end bdcsvd_compare_to_full +#define SVD_DEFAULT(M) BDCSVD<M> +// #define SVD_FOR_MIN_NORM(M) BDCSVD<M> +#define SVD_FOR_MIN_NORM(M) JacobiSVD<M,ColPivHouseholderQRPreconditioner> +#include "../../test/svd_common.h" -template<typename MatrixType> -void bdcsvd_solve(const MatrixType& m, unsigned int computationOptions) -{ - svd_solve< MatrixType, BDCSVD< MatrixType > >(m, computationOptions); -} // end template bdcsvd_solve - - -// test the computations options -template<typename MatrixType> -void bdcsvd_test_all_computation_options(const MatrixType& m) -{ - BDCSVD<MatrixType> fullSvd(m, ComputeFullU|ComputeFullV); - svd_test_computation_options_1< MatrixType, BDCSVD< MatrixType > >(m, fullSvd); - svd_test_computation_options_2< MatrixType, BDCSVD< MatrixType > >(m, fullSvd); -} // end bdcsvd_test_all_computation_options - - -// Call a test with all the computations options +// Check all variants of JacobiSVD template<typename MatrixType> void bdcsvd(const MatrixType& a = MatrixType(), bool pickrandom = true) { - MatrixType m = pickrandom ? MatrixType::Random(a.rows(), a.cols()) : a; - bdcsvd_test_all_computation_options<MatrixType>(m); -} // end template bdcsvd - - -// verify assert -template<typename MatrixType> -void bdcsvd_verify_assert(const MatrixType& m) -{ - svd_verify_assert< MatrixType, BDCSVD< MatrixType > >(m); -}// end template bdcsvd_verify_assert + MatrixType m = a; + if(pickrandom) + svd_fill_random(m); + CALL_SUBTEST(( svd_test_all_computation_options<BDCSVD<MatrixType> >(m, false) )); +} -// test weird values -template<typename MatrixType> -void bdcsvd_inf_nan() -{ - svd_inf_nan< MatrixType, BDCSVD< MatrixType > >(); -}// end template bdcsvd_inf_nan - - - -void bdcsvd_preallocate() -{ - svd_preallocate< BDCSVD< MatrixXf > >(); -} // end bdcsvd_preallocate - +// 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) { - std::cout << "debut compare" << std::endl; 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()); - std::cout << "fin compare" << std::endl; -} // end template compare_bdc_jacobi - - -// call the tests + 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() { - // test of Dynamic defined Matrix (42, 42) of float - CALL_SUBTEST_11(( bdcsvd_verify_assert<Matrix<float,Dynamic,Dynamic> > - (Matrix<float,Dynamic,Dynamic>(42,42)) )); - CALL_SUBTEST_11(( compare_bdc_jacobi<Matrix<float,Dynamic,Dynamic> > - (Matrix<float,Dynamic,Dynamic>(42,42), 0) )); - CALL_SUBTEST_11(( bdcsvd<Matrix<float,Dynamic,Dynamic> > - (Matrix<float,Dynamic,Dynamic>(42,42)) )); - - // test of Dynamic defined Matrix (50, 50) of double - CALL_SUBTEST_13(( bdcsvd_verify_assert<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(50,50)) )); - CALL_SUBTEST_13(( compare_bdc_jacobi<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(50,50), 0) )); - CALL_SUBTEST_13(( bdcsvd<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(50, 50)) )); - - // test of Dynamic defined Matrix (22, 22) of complex double - CALL_SUBTEST_14(( bdcsvd_verify_assert<Matrix<std::complex<double>,Dynamic,Dynamic> > - (Matrix<std::complex<double>,Dynamic,Dynamic>(22,22)) )); - CALL_SUBTEST_14(( compare_bdc_jacobi<Matrix<std::complex<double>,Dynamic,Dynamic> > - (Matrix<std::complex<double>, Dynamic, Dynamic> (22,22), 0) )); - CALL_SUBTEST_14(( bdcsvd<Matrix<std::complex<double>,Dynamic,Dynamic> > - (Matrix<std::complex<double>,Dynamic,Dynamic>(22, 22)) )); - - // test of Dynamic defined Matrix (10, 10) of int - //CALL_SUBTEST_15(( bdcsvd_verify_assert<Matrix<int,Dynamic,Dynamic> > - // (Matrix<int,Dynamic,Dynamic>(10,10)) )); - //CALL_SUBTEST_15(( compare_bdc_jacobi<Matrix<int,Dynamic,Dynamic> > - // (Matrix<int,Dynamic,Dynamic>(10,10), 0) )); - //CALL_SUBTEST_15(( bdcsvd<Matrix<int,Dynamic,Dynamic> > - // (Matrix<int,Dynamic,Dynamic>(10, 10)) )); + 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)) )); +// svd_all_trivial_2x2(bdcsvd<Matrix2cd>); +// svd_all_trivial_2x2(bdcsvd<Matrix2d>); - // test of Dynamic defined Matrix (8, 6) of double - - CALL_SUBTEST_16(( bdcsvd_verify_assert<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(8,6)) )); - CALL_SUBTEST_16(( compare_bdc_jacobi<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(8, 6), 0) )); - CALL_SUBTEST_16(( bdcsvd<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(8, 6)) )); - - - - // test of Dynamic defined Matrix (36, 12) of float - CALL_SUBTEST_17(( compare_bdc_jacobi<Matrix<float,Dynamic,Dynamic> > - (Matrix<float,Dynamic,Dynamic>(36, 12), 0) )); - CALL_SUBTEST_17(( bdcsvd<Matrix<float,Dynamic,Dynamic> > - (Matrix<float,Dynamic,Dynamic>(36, 12)) )); - - // test of Dynamic defined Matrix (5, 8) of double - CALL_SUBTEST_18(( compare_bdc_jacobi<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(5, 8), 0) )); - CALL_SUBTEST_18(( bdcsvd<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(5, 8)) )); - - - // non regression tests - CALL_SUBTEST_3(( bdcsvd_verify_assert(Matrix3f()) )); - CALL_SUBTEST_4(( bdcsvd_verify_assert(Matrix4d()) )); - CALL_SUBTEST_7(( bdcsvd_verify_assert(MatrixXf(10,12)) )); - CALL_SUBTEST_8(( bdcsvd_verify_assert(MatrixXcd(7,5)) )); - - // SUBTESTS 1 and 2 on specifics matrix for(int i = 0; i < g_repeat; i++) { - Matrix2cd m; - m << 0, 1, - 0, 1; - CALL_SUBTEST_1(( bdcsvd(m, false) )); - m << 1, 0, - 1, 0; - CALL_SUBTEST_1(( bdcsvd(m, false) )); - - Matrix2d n; - n << 0, 0, - 0, 0; - CALL_SUBTEST_2(( bdcsvd(n, false) )); - n << 0, 0, - 0, 1; - CALL_SUBTEST_2(( bdcsvd(n, false) )); +// 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); - // Statics matrix don't work with BDSVD yet - // bdc algo on a random 3x3 float matrix - // CALL_SUBTEST_3(( bdcsvd<Matrix3f>() )); - // bdc algo on a random 4x4 double matrix - // CALL_SUBTEST_4(( bdcsvd<Matrix4d>() )); - // bdc algo on a random 3x5 float matrix - // CALL_SUBTEST_5(( bdcsvd<Matrix<float,3,5> >() )); - - int r = internal::random<int>(1, 30), - c = internal::random<int>(1, 30); - CALL_SUBTEST_7(( bdcsvd<MatrixXf>(MatrixXf(r,c)) )); - CALL_SUBTEST_8(( bdcsvd<MatrixXcd>(MatrixXcd(r,c)) )); + 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)) )); (void) r; (void) c; // Test on inf/nan matrix - CALL_SUBTEST_7( bdcsvd_inf_nan<MatrixXf>() ); + CALL_SUBTEST_7( (svd_inf_nan<BDCSVD<MatrixXf>, MatrixXf>()) ); + CALL_SUBTEST_10( (svd_inf_nan<BDCSVD<MatrixXd>, MatrixXd>()) ); } - CALL_SUBTEST_7(( bdcsvd<MatrixXf>(MatrixXf(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) )); - CALL_SUBTEST_8(( bdcsvd<MatrixXcd>(MatrixXcd(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3))) )); + // 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) ); -} // end test_bdcsvd + // Check that preallocation avoids subsequent mallocs + CALL_SUBTEST_9( svd_preallocate() ); + + CALL_SUBTEST_2( svd_underoverflow() ); +} + diff --git a/unsupported/test/jacobisvd.cpp b/unsupported/test/jacobisvd.cpp deleted file mode 100644 index b4e884eee..000000000 --- a/unsupported/test/jacobisvd.cpp +++ /dev/null @@ -1,198 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> -// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> -// -// 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/. - -#include "svd_common.h" - -template<typename MatrixType, int QRPreconditioner> -void jacobisvd_check_full(const MatrixType& m, const JacobiSVD<MatrixType, QRPreconditioner>& svd) -{ - svd_check_full<MatrixType, JacobiSVD<MatrixType, QRPreconditioner > >(m, svd); -} - -template<typename MatrixType, int QRPreconditioner> -void jacobisvd_compare_to_full(const MatrixType& m, - unsigned int computationOptions, - const JacobiSVD<MatrixType, QRPreconditioner>& referenceSvd) -{ - svd_compare_to_full<MatrixType, JacobiSVD<MatrixType, QRPreconditioner> >(m, computationOptions, referenceSvd); -} - - -template<typename MatrixType, int QRPreconditioner> -void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions) -{ - svd_solve< MatrixType, JacobiSVD< MatrixType, QRPreconditioner > >(m, computationOptions); -} - - - -template<typename MatrixType, int QRPreconditioner> -void jacobisvd_test_all_computation_options(const MatrixType& m) -{ - - if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols()) - return; - - JacobiSVD< MatrixType, QRPreconditioner > fullSvd(m, ComputeFullU|ComputeFullV); - svd_test_computation_options_1< MatrixType, JacobiSVD< MatrixType, QRPreconditioner > >(m, fullSvd); - - if(QRPreconditioner == FullPivHouseholderQRPreconditioner) - return; - svd_test_computation_options_2< MatrixType, JacobiSVD< MatrixType, QRPreconditioner > >(m, fullSvd); - -} - -template<typename MatrixType> -void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true) -{ - MatrixType m = pickrandom ? MatrixType::Random(a.rows(), a.cols()) : a; - - 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); -} - - -template<typename MatrixType> -void jacobisvd_verify_assert(const MatrixType& m) -{ - - svd_verify_assert<MatrixType, JacobiSVD< MatrixType > >(m); - - typedef typename MatrixType::Index Index; - Index rows = m.rows(); - Index cols = m.cols(); - - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime - }; - - MatrixType a = MatrixType::Zero(rows, cols); - a.setZero(); - - if (ColsAtCompileTime == Dynamic) - { - JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner> svd_fullqr; - VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeFullU|ComputeThinV)) - VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeThinV)) - VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeFullV)) - } -} - -template<typename MatrixType> -void jacobisvd_method() -{ - enum { Size = MatrixType::RowsAtCompileTime }; - typedef typename MatrixType::RealScalar RealScalar; - typedef Matrix<RealScalar, Size, 1> RealVecType; - MatrixType m = MatrixType::Identity(); - VERIFY_IS_APPROX(m.jacobiSvd().singularValues(), RealVecType::Ones()); - VERIFY_RAISES_ASSERT(m.jacobiSvd().matrixU()); - VERIFY_RAISES_ASSERT(m.jacobiSvd().matrixV()); - VERIFY_IS_APPROX(m.jacobiSvd(ComputeFullU|ComputeFullV).solve(m), m); -} - - - -template<typename MatrixType> -void jacobisvd_inf_nan() -{ - svd_inf_nan<MatrixType, JacobiSVD< MatrixType > >(); -} - - -// Regression test for bug 286: JacobiSVD loops indefinitely with some -// matrices containing denormal numbers. -void jacobisvd_bug286() -{ -#if defined __INTEL_COMPILER -// shut up warning #239: floating point underflow -#pragma warning push -#pragma warning disable 239 -#endif - Matrix2d M; - M << -7.90884e-313, -4.94e-324, - 0, 5.60844e-313; -#if defined __INTEL_COMPILER -#pragma warning pop -#endif - JacobiSVD<Matrix2d> svd; - svd.compute(M); // just check we don't loop indefinitely -} - - -void jacobisvd_preallocate() -{ - svd_preallocate< JacobiSVD <MatrixXf> >(); -} - -void test_jacobisvd() -{ - CALL_SUBTEST_11(( jacobisvd<Matrix<double,Dynamic,Dynamic> > - (Matrix<double,Dynamic,Dynamic>(16, 6)) )); - - CALL_SUBTEST_3(( jacobisvd_verify_assert(Matrix3f()) )); - CALL_SUBTEST_4(( jacobisvd_verify_assert(Matrix4d()) )); - CALL_SUBTEST_7(( jacobisvd_verify_assert(MatrixXf(10,12)) )); - CALL_SUBTEST_8(( jacobisvd_verify_assert(MatrixXcd(7,5)) )); - - for(int i = 0; i < g_repeat; i++) { - Matrix2cd m; - m << 0, 1, - 0, 1; - CALL_SUBTEST_1(( jacobisvd(m, false) )); - m << 1, 0, - 1, 0; - CALL_SUBTEST_1(( jacobisvd(m, false) )); - - Matrix2d n; - 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> >() )); - CALL_SUBTEST_6(( jacobisvd<Matrix<double,Dynamic,2> >(Matrix<double,Dynamic,2>(10,2)) )); - - int r = internal::random<int>(1, 30), - c = internal::random<int>(1, 30); - CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(r,c)) )); - 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(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) )); - CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3))) )); - - - // test matrixbase method - CALL_SUBTEST_1(( jacobisvd_method<Matrix2cd>() )); - CALL_SUBTEST_3(( jacobisvd_method<Matrix3f>() )); - - - // Test problem size constructors - CALL_SUBTEST_7( JacobiSVD<MatrixXf>(10,10) ); - - // Check that preallocation avoids subsequent mallocs - CALL_SUBTEST_9( jacobisvd_preallocate() ); - - // Regression check for bug 286 - CALL_SUBTEST_2( jacobisvd_bug286() ); -} diff --git a/unsupported/test/svd_common.h b/unsupported/test/svd_common.h deleted file mode 100644 index 6a96e7c74..000000000 --- a/unsupported/test/svd_common.h +++ /dev/null @@ -1,261 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> -// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> -// -// 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 <unsupported/Eigen/BDCSVD> -#include <Eigen/LU> - - -// check if "svd" is the good image of "m" -template<typename MatrixType, typename SVD> -void svd_check_full(const MatrixType& m, const SVD& svd) -{ - typedef typename MatrixType::Index Index; - Index rows = m.rows(); - Index cols = m.cols(); - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime - }; - - typedef typename MatrixType::Scalar Scalar; - typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType; - typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType; - - - MatrixType sigma = MatrixType::Zero(rows, cols); - sigma.diagonal() = svd.singularValues().template cast<Scalar>(); - MatrixUType u = svd.matrixU(); - MatrixVType v = svd.matrixV(); - VERIFY_IS_APPROX(m, u * sigma * v.adjoint()); - VERIFY_IS_UNITARY(u); - VERIFY_IS_UNITARY(v); -} // end svd_check_full - - - -// Compare to a reference value -template<typename MatrixType, typename SVD> -void svd_compare_to_full(const MatrixType& m, - unsigned int computationOptions, - const SVD& referenceSvd) -{ - typedef typename MatrixType::Index Index; - Index rows = m.rows(); - Index cols = m.cols(); - Index diagSize = (std::min)(rows, cols); - - SVD svd(m, computationOptions); - - VERIFY_IS_APPROX(svd.singularValues(), referenceSvd.singularValues()); - if(computationOptions & ComputeFullU) - VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU()); - if(computationOptions & ComputeThinU) - VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize)); - if(computationOptions & ComputeFullV) - VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV()); - if(computationOptions & ComputeThinV) - VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize)); -} // end svd_compare_to_full - - - -template<typename MatrixType, typename SVD> -void svd_solve(const MatrixType& m, unsigned int computationOptions) -{ - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::Index Index; - Index rows = m.rows(); - Index cols = m.cols(); - - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime - }; - - typedef Matrix<Scalar, RowsAtCompileTime, Dynamic> RhsType; - typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType; - - RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols)); - SVD svd(m, computationOptions); - SolutionType x = svd.solve(rhs); - // evaluate normal equation which works also for least-squares solutions - VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs); -} // end svd_solve - - -// test computations options -// 2 functions because Jacobisvd can return before the second function -template<typename MatrixType, typename SVD> -void svd_test_computation_options_1(const MatrixType& m, const SVD& fullSvd) -{ - svd_check_full< MatrixType, SVD >(m, fullSvd); - svd_solve< MatrixType, SVD >(m, ComputeFullU | ComputeFullV); -} - - -template<typename MatrixType, typename SVD> -void svd_test_computation_options_2(const MatrixType& m, const SVD& fullSvd) -{ - svd_compare_to_full< MatrixType, SVD >(m, ComputeFullU, fullSvd); - svd_compare_to_full< MatrixType, SVD >(m, ComputeFullV, fullSvd); - svd_compare_to_full< MatrixType, SVD >(m, 0, fullSvd); - - if (MatrixType::ColsAtCompileTime == Dynamic) { - // thin U/V are only available with dynamic number of columns - - svd_compare_to_full< MatrixType, SVD >(m, ComputeFullU|ComputeThinV, fullSvd); - svd_compare_to_full< MatrixType, SVD >(m, ComputeThinV, fullSvd); - svd_compare_to_full< MatrixType, SVD >(m, ComputeThinU|ComputeFullV, fullSvd); - svd_compare_to_full< MatrixType, SVD >(m, ComputeThinU , fullSvd); - svd_compare_to_full< MatrixType, SVD >(m, ComputeThinU|ComputeThinV, fullSvd); - svd_solve<MatrixType, SVD>(m, ComputeFullU | ComputeThinV); - svd_solve<MatrixType, SVD>(m, ComputeThinU | ComputeFullV); - svd_solve<MatrixType, SVD>(m, ComputeThinU | ComputeThinV); - - typedef typename MatrixType::Index Index; - Index diagSize = (std::min)(m.rows(), m.cols()); - SVD svd(m, ComputeThinU | ComputeThinV); - VERIFY_IS_APPROX(m, svd.matrixU().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint()); - } -} - -template<typename MatrixType, typename SVD> -void svd_verify_assert(const MatrixType& m) -{ - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::Index Index; - Index rows = m.rows(); - Index cols = m.cols(); - - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime - }; - - typedef Matrix<Scalar, RowsAtCompileTime, 1> RhsType; - RhsType rhs(rows); - SVD svd; - VERIFY_RAISES_ASSERT(svd.matrixU()) - VERIFY_RAISES_ASSERT(svd.singularValues()) - VERIFY_RAISES_ASSERT(svd.matrixV()) - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - MatrixType a = MatrixType::Zero(rows, cols); - a.setZero(); - svd.compute(a, 0); - VERIFY_RAISES_ASSERT(svd.matrixU()) - VERIFY_RAISES_ASSERT(svd.matrixV()) - svd.singularValues(); - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - - if (ColsAtCompileTime == Dynamic) - { - svd.compute(a, ComputeThinU); - svd.matrixU(); - VERIFY_RAISES_ASSERT(svd.matrixV()) - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - svd.compute(a, ComputeThinV); - svd.matrixV(); - VERIFY_RAISES_ASSERT(svd.matrixU()) - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - } - else - { - VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinU)) - VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinV)) - } -} - -// work around stupid msvc error when constructing at compile time an expression that involves -// a division by zero, even if the numeric type has floating point -template<typename Scalar> -EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); } - -// workaround aggressive optimization in ICC -template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } - - -template<typename MatrixType, typename SVD> -void svd_inf_nan() -{ - // all this function does is verify we don't iterate infinitely on nan/inf values - - SVD svd; - typedef typename MatrixType::Scalar Scalar; - Scalar some_inf = Scalar(1) / zero<Scalar>(); - VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf)); - svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV); - - Scalar some_nan = zero<Scalar> () / zero<Scalar> (); - VERIFY(some_nan != some_nan); - svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV); - - MatrixType m = MatrixType::Zero(10,10); - m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf; - svd.compute(m, ComputeFullU | ComputeFullV); - - m = MatrixType::Zero(10,10); - m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_nan; - svd.compute(m, ComputeFullU | ComputeFullV); -} - - -template<typename SVD> -void svd_preallocate() -{ - Vector3f v(3.f, 2.f, 1.f); - MatrixXf m = v.asDiagonal(); - - internal::set_is_malloc_allowed(false); - VERIFY_RAISES_ASSERT(VectorXf v(10);) - SVD svd; - internal::set_is_malloc_allowed(true); - svd.compute(m); - VERIFY_IS_APPROX(svd.singularValues(), v); - - SVD svd2(3,3); - internal::set_is_malloc_allowed(false); - svd2.compute(m); - internal::set_is_malloc_allowed(true); - VERIFY_IS_APPROX(svd2.singularValues(), v); - VERIFY_RAISES_ASSERT(svd2.matrixU()); - VERIFY_RAISES_ASSERT(svd2.matrixV()); - svd2.compute(m, ComputeFullU | ComputeFullV); - VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); - VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); - internal::set_is_malloc_allowed(false); - svd2.compute(m); - internal::set_is_malloc_allowed(true); - - SVD svd3(3,3,ComputeFullU|ComputeFullV); - internal::set_is_malloc_allowed(false); - svd2.compute(m); - internal::set_is_malloc_allowed(true); - VERIFY_IS_APPROX(svd2.singularValues(), v); - VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); - VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); - internal::set_is_malloc_allowed(false); - svd2.compute(m, ComputeFullU|ComputeFullV); - internal::set_is_malloc_allowed(true); -} - - - - - |