1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
6 //
7 // Copyright (C) 2013 Gauthier Brun <brun.gauthier@gmail.com>
8 // Copyright (C) 2013 Nicolas Carre <nicolas.carre@ensimag.fr>
9 // Copyright (C) 2013 Jean Ceccato <jean.ceccato@ensimag.fr>
10 // Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli@ensimag.fr>
11 //
12 // This Source Code Form is subject to the terms of the Mozilla
13 // Public License v. 2.0. If a copy of the MPL was not distributed
14 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
15
16 // discard stack allocation as that too bypasses malloc
17 #define EIGEN_STACK_ALLOCATION_LIMIT 0
18 #define EIGEN_RUNTIME_NO_MALLOC
19
20 #include "main.h"
21 #include <unsupported/Eigen/SVD>
22 #include <Eigen/LU>
23
24
25 // check if "svd" is the good image of "m"
26 template<typename MatrixType, typename SVD>
svd_check_full(const MatrixType & m,const SVD & svd)27 void svd_check_full(const MatrixType& m, const SVD& svd)
28 {
29 typedef typename MatrixType::Index Index;
30 Index rows = m.rows();
31 Index cols = m.cols();
32 enum {
33 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
34 ColsAtCompileTime = MatrixType::ColsAtCompileTime
35 };
36
37 typedef typename MatrixType::Scalar Scalar;
38 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType;
39 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType;
40
41
42 MatrixType sigma = MatrixType::Zero(rows, cols);
43 sigma.diagonal() = svd.singularValues().template cast<Scalar>();
44 MatrixUType u = svd.matrixU();
45 MatrixVType v = svd.matrixV();
46 VERIFY_IS_APPROX(m, u * sigma * v.adjoint());
47 VERIFY_IS_UNITARY(u);
48 VERIFY_IS_UNITARY(v);
49 } // end svd_check_full
50
51
52
53 // Compare to a reference value
54 template<typename MatrixType, typename SVD>
svd_compare_to_full(const MatrixType & m,unsigned int computationOptions,const SVD & referenceSvd)55 void svd_compare_to_full(const MatrixType& m,
56 unsigned int computationOptions,
57 const SVD& referenceSvd)
58 {
59 typedef typename MatrixType::Index Index;
60 Index rows = m.rows();
61 Index cols = m.cols();
62 Index diagSize = (std::min)(rows, cols);
63
64 SVD svd(m, computationOptions);
65
66 VERIFY_IS_APPROX(svd.singularValues(), referenceSvd.singularValues());
67 if(computationOptions & ComputeFullU)
68 VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU());
69 if(computationOptions & ComputeThinU)
70 VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize));
71 if(computationOptions & ComputeFullV)
72 VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV());
73 if(computationOptions & ComputeThinV)
74 VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize));
75 } // end svd_compare_to_full
76
77
78
79 template<typename MatrixType, typename SVD>
svd_solve(const MatrixType & m,unsigned int computationOptions)80 void svd_solve(const MatrixType& m, unsigned int computationOptions)
81 {
82 typedef typename MatrixType::Scalar Scalar;
83 typedef typename MatrixType::Index Index;
84 Index rows = m.rows();
85 Index cols = m.cols();
86
87 enum {
88 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
89 ColsAtCompileTime = MatrixType::ColsAtCompileTime
90 };
91
92 typedef Matrix<Scalar, RowsAtCompileTime, Dynamic> RhsType;
93 typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType;
94
95 RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols));
96 SVD svd(m, computationOptions);
97 SolutionType x = svd.solve(rhs);
98 // evaluate normal equation which works also for least-squares solutions
99 VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs);
100 } // end svd_solve
101
102
103 // test computations options
104 // 2 functions because Jacobisvd can return before the second function
105 template<typename MatrixType, typename SVD>
svd_test_computation_options_1(const MatrixType & m,const SVD & fullSvd)106 void svd_test_computation_options_1(const MatrixType& m, const SVD& fullSvd)
107 {
108 svd_check_full< MatrixType, SVD >(m, fullSvd);
109 svd_solve< MatrixType, SVD >(m, ComputeFullU | ComputeFullV);
110 }
111
112
113 template<typename MatrixType, typename SVD>
svd_test_computation_options_2(const MatrixType & m,const SVD & fullSvd)114 void svd_test_computation_options_2(const MatrixType& m, const SVD& fullSvd)
115 {
116 svd_compare_to_full< MatrixType, SVD >(m, ComputeFullU, fullSvd);
117 svd_compare_to_full< MatrixType, SVD >(m, ComputeFullV, fullSvd);
118 svd_compare_to_full< MatrixType, SVD >(m, 0, fullSvd);
119
120 if (MatrixType::ColsAtCompileTime == Dynamic) {
121 // thin U/V are only available with dynamic number of columns
122
123 svd_compare_to_full< MatrixType, SVD >(m, ComputeFullU|ComputeThinV, fullSvd);
124 svd_compare_to_full< MatrixType, SVD >(m, ComputeThinV, fullSvd);
125 svd_compare_to_full< MatrixType, SVD >(m, ComputeThinU|ComputeFullV, fullSvd);
126 svd_compare_to_full< MatrixType, SVD >(m, ComputeThinU , fullSvd);
127 svd_compare_to_full< MatrixType, SVD >(m, ComputeThinU|ComputeThinV, fullSvd);
128 svd_solve<MatrixType, SVD>(m, ComputeFullU | ComputeThinV);
129 svd_solve<MatrixType, SVD>(m, ComputeThinU | ComputeFullV);
130 svd_solve<MatrixType, SVD>(m, ComputeThinU | ComputeThinV);
131
132 typedef typename MatrixType::Index Index;
133 Index diagSize = (std::min)(m.rows(), m.cols());
134 SVD svd(m, ComputeThinU | ComputeThinV);
135 VERIFY_IS_APPROX(m, svd.matrixU().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint());
136 }
137 }
138
139 template<typename MatrixType, typename SVD>
svd_verify_assert(const MatrixType & m)140 void svd_verify_assert(const MatrixType& m)
141 {
142 typedef typename MatrixType::Scalar Scalar;
143 typedef typename MatrixType::Index Index;
144 Index rows = m.rows();
145 Index cols = m.cols();
146
147 enum {
148 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
149 ColsAtCompileTime = MatrixType::ColsAtCompileTime
150 };
151
152 typedef Matrix<Scalar, RowsAtCompileTime, 1> RhsType;
153 RhsType rhs(rows);
154 SVD svd;
155 VERIFY_RAISES_ASSERT(svd.matrixU())
156 VERIFY_RAISES_ASSERT(svd.singularValues())
157 VERIFY_RAISES_ASSERT(svd.matrixV())
158 VERIFY_RAISES_ASSERT(svd.solve(rhs))
159 MatrixType a = MatrixType::Zero(rows, cols);
160 a.setZero();
161 svd.compute(a, 0);
162 VERIFY_RAISES_ASSERT(svd.matrixU())
163 VERIFY_RAISES_ASSERT(svd.matrixV())
164 svd.singularValues();
165 VERIFY_RAISES_ASSERT(svd.solve(rhs))
166
167 if (ColsAtCompileTime == Dynamic)
168 {
169 svd.compute(a, ComputeThinU);
170 svd.matrixU();
171 VERIFY_RAISES_ASSERT(svd.matrixV())
172 VERIFY_RAISES_ASSERT(svd.solve(rhs))
173 svd.compute(a, ComputeThinV);
174 svd.matrixV();
175 VERIFY_RAISES_ASSERT(svd.matrixU())
176 VERIFY_RAISES_ASSERT(svd.solve(rhs))
177 }
178 else
179 {
180 VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinU))
181 VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinV))
182 }
183 }
184
185 // work around stupid msvc error when constructing at compile time an expression that involves
186 // a division by zero, even if the numeric type has floating point
187 template<typename Scalar>
zero()188 EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); }
189
190 // workaround aggressive optimization in ICC
sub(T a,T b)191 template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; }
192
193
194 template<typename MatrixType, typename SVD>
svd_inf_nan()195 void svd_inf_nan()
196 {
197 // all this function does is verify we don't iterate infinitely on nan/inf values
198
199 SVD svd;
200 typedef typename MatrixType::Scalar Scalar;
201 Scalar some_inf = Scalar(1) / zero<Scalar>();
202 VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf));
203 svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV);
204
205 Scalar some_nan = zero<Scalar> () / zero<Scalar> ();
206 VERIFY(some_nan != some_nan);
207 svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV);
208
209 MatrixType m = MatrixType::Zero(10,10);
210 m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf;
211 svd.compute(m, ComputeFullU | ComputeFullV);
212
213 m = MatrixType::Zero(10,10);
214 m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_nan;
215 svd.compute(m, ComputeFullU | ComputeFullV);
216 }
217
218
219 template<typename SVD>
svd_preallocate()220 void svd_preallocate()
221 {
222 Vector3f v(3.f, 2.f, 1.f);
223 MatrixXf m = v.asDiagonal();
224
225 internal::set_is_malloc_allowed(false);
226 VERIFY_RAISES_ASSERT(VectorXf v(10);)
227 SVD svd;
228 internal::set_is_malloc_allowed(true);
229 svd.compute(m);
230 VERIFY_IS_APPROX(svd.singularValues(), v);
231
232 SVD svd2(3,3);
233 internal::set_is_malloc_allowed(false);
234 svd2.compute(m);
235 internal::set_is_malloc_allowed(true);
236 VERIFY_IS_APPROX(svd2.singularValues(), v);
237 VERIFY_RAISES_ASSERT(svd2.matrixU());
238 VERIFY_RAISES_ASSERT(svd2.matrixV());
239 svd2.compute(m, ComputeFullU | ComputeFullV);
240 VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
241 VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
242 internal::set_is_malloc_allowed(false);
243 svd2.compute(m);
244 internal::set_is_malloc_allowed(true);
245
246 SVD svd3(3,3,ComputeFullU|ComputeFullV);
247 internal::set_is_malloc_allowed(false);
248 svd2.compute(m);
249 internal::set_is_malloc_allowed(true);
250 VERIFY_IS_APPROX(svd2.singularValues(), v);
251 VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
252 VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
253 internal::set_is_malloc_allowed(false);
254 svd2.compute(m, ComputeFullU|ComputeFullV);
255 internal::set_is_malloc_allowed(true);
256 }
257
258
259
260
261
262