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>
bdcsvd(const MatrixType & a=MatrixType (),bool pickrandom=true)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>
bdcsvd_method()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>
compare_bdc_jacobi(const MatrixType & a=MatrixType (),unsigned int computationOptions=0)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
test_bdcsvd()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