bdcsvd.cpp (4443B)
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@gmail.com> 5 // Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr> 6 // Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr> 7 // Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@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; 32 if(pickrandom) { 33 m.resizeLike(a); 34 svd_fill_random(m); 35 } 36 else 37 m = a; 38 39 CALL_SUBTEST(( svd_test_all_computation_options<BDCSVD<MatrixType> >(m, false) )); 40 } 41 42 template<typename MatrixType> 43 void bdcsvd_method() 44 { 45 enum { Size = MatrixType::RowsAtCompileTime }; 46 typedef typename MatrixType::RealScalar RealScalar; 47 typedef Matrix<RealScalar, Size, 1> RealVecType; 48 MatrixType m = MatrixType::Identity(); 49 VERIFY_IS_APPROX(m.bdcSvd().singularValues(), RealVecType::Ones()); 50 VERIFY_RAISES_ASSERT(m.bdcSvd().matrixU()); 51 VERIFY_RAISES_ASSERT(m.bdcSvd().matrixV()); 52 VERIFY_IS_APPROX(m.bdcSvd(ComputeFullU|ComputeFullV).solve(m), m); 53 VERIFY_IS_APPROX(m.bdcSvd(ComputeFullU|ComputeFullV).transpose().solve(m), m); 54 VERIFY_IS_APPROX(m.bdcSvd(ComputeFullU|ComputeFullV).adjoint().solve(m), m); 55 } 56 57 // compare the Singular values returned with Jacobi and Bdc 58 template<typename MatrixType> 59 void compare_bdc_jacobi(const MatrixType& a = MatrixType(), unsigned int computationOptions = 0) 60 { 61 MatrixType m = MatrixType::Random(a.rows(), a.cols()); 62 BDCSVD<MatrixType> bdc_svd(m); 63 JacobiSVD<MatrixType> jacobi_svd(m); 64 VERIFY_IS_APPROX(bdc_svd.singularValues(), jacobi_svd.singularValues()); 65 if(computationOptions & ComputeFullU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU()); 66 if(computationOptions & ComputeThinU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU()); 67 if(computationOptions & ComputeFullV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV()); 68 if(computationOptions & ComputeThinV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV()); 69 } 70 71 EIGEN_DECLARE_TEST(bdcsvd) 72 { 73 CALL_SUBTEST_3(( svd_verify_assert<BDCSVD<Matrix3f> >(Matrix3f()) )); 74 CALL_SUBTEST_4(( svd_verify_assert<BDCSVD<Matrix4d> >(Matrix4d()) )); 75 CALL_SUBTEST_7(( svd_verify_assert<BDCSVD<MatrixXf> >(MatrixXf(10,12)) )); 76 CALL_SUBTEST_8(( svd_verify_assert<BDCSVD<MatrixXcd> >(MatrixXcd(7,5)) )); 77 78 CALL_SUBTEST_101(( svd_all_trivial_2x2(bdcsvd<Matrix2cd>) )); 79 CALL_SUBTEST_102(( svd_all_trivial_2x2(bdcsvd<Matrix2d>) )); 80 81 for(int i = 0; i < g_repeat; i++) { 82 CALL_SUBTEST_3(( bdcsvd<Matrix3f>() )); 83 CALL_SUBTEST_4(( bdcsvd<Matrix4d>() )); 84 CALL_SUBTEST_5(( bdcsvd<Matrix<float,3,5> >() )); 85 86 int r = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2), 87 c = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2); 88 89 TEST_SET_BUT_UNUSED_VARIABLE(r) 90 TEST_SET_BUT_UNUSED_VARIABLE(c) 91 92 CALL_SUBTEST_6(( bdcsvd(Matrix<double,Dynamic,2>(r,2)) )); 93 CALL_SUBTEST_7(( bdcsvd(MatrixXf(r,c)) )); 94 CALL_SUBTEST_7(( compare_bdc_jacobi(MatrixXf(r,c)) )); 95 CALL_SUBTEST_10(( bdcsvd(MatrixXd(r,c)) )); 96 CALL_SUBTEST_10(( compare_bdc_jacobi(MatrixXd(r,c)) )); 97 CALL_SUBTEST_8(( bdcsvd(MatrixXcd(r,c)) )); 98 CALL_SUBTEST_8(( compare_bdc_jacobi(MatrixXcd(r,c)) )); 99 100 // Test on inf/nan matrix 101 CALL_SUBTEST_7( (svd_inf_nan<BDCSVD<MatrixXf>, MatrixXf>()) ); 102 CALL_SUBTEST_10( (svd_inf_nan<BDCSVD<MatrixXd>, MatrixXd>()) ); 103 } 104 105 // test matrixbase method 106 CALL_SUBTEST_1(( bdcsvd_method<Matrix2cd>() )); 107 CALL_SUBTEST_3(( bdcsvd_method<Matrix3f>() )); 108 109 // Test problem size constructors 110 CALL_SUBTEST_7( BDCSVD<MatrixXf>(10,10) ); 111 112 // Check that preallocation avoids subsequent mallocs 113 // Disabled because not supported by BDCSVD 114 // CALL_SUBTEST_9( svd_preallocate<void>() ); 115 116 CALL_SUBTEST_2( svd_underoverflow<void>() ); 117 } 118