1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2013 Gauthier Brun <brun.gauthier (at) gmail.com> 5 // Copyright (C) 2013 Nicolas Carre <nicolas.carre (at) ensimag.fr> 6 // Copyright (C) 2013 Jean Ceccato <jean.ceccato (at) ensimag.fr> 7 // Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli (at) ensimag.fr> 8 // 9 // This Source Code Form is subject to the terms of the Mozilla 10 // Public License v. 2.0. If a copy of the MPL was not distributed 11 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/ 12 13 // discard stack allocation as that too bypasses malloc 14 #define EIGEN_STACK_ALLOCATION_LIMIT 0 15 #define EIGEN_RUNTIME_NO_MALLOC 16 17 #include "main.h" 18 #include <Eigen/SVD> 19 #include <iostream> 20 #include <Eigen/LU> 21 22 23 #define SVD_DEFAULT(M) BDCSVD<M> 24 #define SVD_FOR_MIN_NORM(M) BDCSVD<M> 25 #include "svd_common.h" 26 27 // Check all variants of JacobiSVD 28 template<typename MatrixType> 29 void bdcsvd(const MatrixType& a = MatrixType(), bool pickrandom = true) 30 { 31 MatrixType m = a; 32 if(pickrandom) 33 svd_fill_random(m); 34 35 CALL_SUBTEST(( svd_test_all_computation_options<BDCSVD<MatrixType> >(m, false) )); 36 } 37 38 template<typename MatrixType> 39 void bdcsvd_method() 40 { 41 enum { Size = MatrixType::RowsAtCompileTime }; 42 typedef typename MatrixType::RealScalar RealScalar; 43 typedef Matrix<RealScalar, Size, 1> RealVecType; 44 MatrixType m = MatrixType::Identity(); 45 VERIFY_IS_APPROX(m.bdcSvd().singularValues(), RealVecType::Ones()); 46 VERIFY_RAISES_ASSERT(m.bdcSvd().matrixU()); 47 VERIFY_RAISES_ASSERT(m.bdcSvd().matrixV()); 48 VERIFY_IS_APPROX(m.bdcSvd(ComputeFullU|ComputeFullV).solve(m), m); 49 } 50 51 // compare the Singular values returned with Jacobi and Bdc 52 template<typename MatrixType> 53 void compare_bdc_jacobi(const MatrixType& a = MatrixType(), unsigned int computationOptions = 0) 54 { 55 MatrixType m = MatrixType::Random(a.rows(), a.cols()); 56 BDCSVD<MatrixType> bdc_svd(m); 57 JacobiSVD<MatrixType> jacobi_svd(m); 58 VERIFY_IS_APPROX(bdc_svd.singularValues(), jacobi_svd.singularValues()); 59 if(computationOptions & ComputeFullU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU()); 60 if(computationOptions & ComputeThinU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU()); 61 if(computationOptions & ComputeFullV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV()); 62 if(computationOptions & ComputeThinV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV()); 63 } 64 65 void test_bdcsvd() 66 { 67 CALL_SUBTEST_3(( svd_verify_assert<BDCSVD<Matrix3f> >(Matrix3f()) )); 68 CALL_SUBTEST_4(( svd_verify_assert<BDCSVD<Matrix4d> >(Matrix4d()) )); 69 CALL_SUBTEST_7(( svd_verify_assert<BDCSVD<MatrixXf> >(MatrixXf(10,12)) )); 70 CALL_SUBTEST_8(( svd_verify_assert<BDCSVD<MatrixXcd> >(MatrixXcd(7,5)) )); 71 72 CALL_SUBTEST_101(( svd_all_trivial_2x2(bdcsvd<Matrix2cd>) )); 73 CALL_SUBTEST_102(( svd_all_trivial_2x2(bdcsvd<Matrix2d>) )); 74 75 for(int i = 0; i < g_repeat; i++) { 76 CALL_SUBTEST_3(( bdcsvd<Matrix3f>() )); 77 CALL_SUBTEST_4(( bdcsvd<Matrix4d>() )); 78 CALL_SUBTEST_5(( bdcsvd<Matrix<float,3,5> >() )); 79 80 int r = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2), 81 c = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2); 82 83 TEST_SET_BUT_UNUSED_VARIABLE(r) 84 TEST_SET_BUT_UNUSED_VARIABLE(c) 85 86 CALL_SUBTEST_6(( bdcsvd(Matrix<double,Dynamic,2>(r,2)) )); 87 CALL_SUBTEST_7(( bdcsvd(MatrixXf(r,c)) )); 88 CALL_SUBTEST_7(( compare_bdc_jacobi(MatrixXf(r,c)) )); 89 CALL_SUBTEST_10(( bdcsvd(MatrixXd(r,c)) )); 90 CALL_SUBTEST_10(( compare_bdc_jacobi(MatrixXd(r,c)) )); 91 CALL_SUBTEST_8(( bdcsvd(MatrixXcd(r,c)) )); 92 CALL_SUBTEST_8(( compare_bdc_jacobi(MatrixXcd(r,c)) )); 93 94 // Test on inf/nan matrix 95 CALL_SUBTEST_7( (svd_inf_nan<BDCSVD<MatrixXf>, MatrixXf>()) ); 96 CALL_SUBTEST_10( (svd_inf_nan<BDCSVD<MatrixXd>, MatrixXd>()) ); 97 } 98 99 // test matrixbase method 100 CALL_SUBTEST_1(( bdcsvd_method<Matrix2cd>() )); 101 CALL_SUBTEST_3(( bdcsvd_method<Matrix3f>() )); 102 103 // Test problem size constructors 104 CALL_SUBTEST_7( BDCSVD<MatrixXf>(10,10) ); 105 106 // Check that preallocation avoids subsequent mallocs 107 CALL_SUBTEST_9( svd_preallocate<void>() ); 108 109 CALL_SUBTEST_2( svd_underoverflow<void>() ); 110 } 111 112
