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 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 // discard stack allocation as that too bypasses malloc
12 #define EIGEN_STACK_ALLOCATION_LIMIT 0
13 #define EIGEN_RUNTIME_NO_MALLOC
14 #include "main.h"
15 #include <Eigen/SVD>
16
17 template<typename MatrixType, int QRPreconditioner>
jacobisvd_check_full(const MatrixType & m,const JacobiSVD<MatrixType,QRPreconditioner> & svd)18 void jacobisvd_check_full(const MatrixType& m, const JacobiSVD<MatrixType, QRPreconditioner>& svd)
19 {
20 typedef typename MatrixType::Index Index;
21 Index rows = m.rows();
22 Index cols = m.cols();
23
24 enum {
25 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
26 ColsAtCompileTime = MatrixType::ColsAtCompileTime
27 };
28
29 typedef typename MatrixType::Scalar Scalar;
30 typedef typename NumTraits<Scalar>::Real RealScalar;
31 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType;
32 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType;
33 typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType;
34 typedef Matrix<Scalar, ColsAtCompileTime, 1> InputVectorType;
35
36 MatrixType sigma = MatrixType::Zero(rows,cols);
37 sigma.diagonal() = svd.singularValues().template cast<Scalar>();
38 MatrixUType u = svd.matrixU();
39 MatrixVType v = svd.matrixV();
40
41 VERIFY_IS_APPROX(m, u * sigma * v.adjoint());
42 VERIFY_IS_UNITARY(u);
43 VERIFY_IS_UNITARY(v);
44 }
45
46 template<typename MatrixType, int QRPreconditioner>
jacobisvd_compare_to_full(const MatrixType & m,unsigned int computationOptions,const JacobiSVD<MatrixType,QRPreconditioner> & referenceSvd)47 void jacobisvd_compare_to_full(const MatrixType& m,
48 unsigned int computationOptions,
49 const JacobiSVD<MatrixType, QRPreconditioner>& referenceSvd)
50 {
51 typedef typename MatrixType::Index Index;
52 Index rows = m.rows();
53 Index cols = m.cols();
54 Index diagSize = (std::min)(rows, cols);
55
56 JacobiSVD<MatrixType, QRPreconditioner> svd(m, computationOptions);
57
58 VERIFY_IS_APPROX(svd.singularValues(), referenceSvd.singularValues());
59 if(computationOptions & ComputeFullU)
60 VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU());
61 if(computationOptions & ComputeThinU)
62 VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize));
63 if(computationOptions & ComputeFullV)
64 VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV());
65 if(computationOptions & ComputeThinV)
66 VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize));
67 }
68
69 template<typename MatrixType, int QRPreconditioner>
jacobisvd_solve(const MatrixType & m,unsigned int computationOptions)70 void jacobisvd_solve(const MatrixType& m, unsigned int computationOptions)
71 {
72 typedef typename MatrixType::Scalar Scalar;
73 typedef typename MatrixType::Index Index;
74 Index rows = m.rows();
75 Index cols = m.cols();
76
77 enum {
78 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
79 ColsAtCompileTime = MatrixType::ColsAtCompileTime
80 };
81
82 typedef Matrix<Scalar, RowsAtCompileTime, Dynamic> RhsType;
83 typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType;
84
85 RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols));
86 JacobiSVD<MatrixType, QRPreconditioner> svd(m, computationOptions);
87 SolutionType x = svd.solve(rhs);
88 // evaluate normal equation which works also for least-squares solutions
89 VERIFY_IS_APPROX(m.adjoint()*m*x,m.adjoint()*rhs);
90 }
91
92 template<typename MatrixType, int QRPreconditioner>
jacobisvd_test_all_computation_options(const MatrixType & m)93 void jacobisvd_test_all_computation_options(const MatrixType& m)
94 {
95 if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols())
96 return;
97 JacobiSVD<MatrixType, QRPreconditioner> fullSvd(m, ComputeFullU|ComputeFullV);
98
99 jacobisvd_check_full(m, fullSvd);
100 jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeFullV);
101
102 if(QRPreconditioner == FullPivHouseholderQRPreconditioner)
103 return;
104
105 jacobisvd_compare_to_full(m, ComputeFullU, fullSvd);
106 jacobisvd_compare_to_full(m, ComputeFullV, fullSvd);
107 jacobisvd_compare_to_full(m, 0, fullSvd);
108
109 if (MatrixType::ColsAtCompileTime == Dynamic) {
110 // thin U/V are only available with dynamic number of columns
111 jacobisvd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd);
112 jacobisvd_compare_to_full(m, ComputeThinV, fullSvd);
113 jacobisvd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd);
114 jacobisvd_compare_to_full(m, ComputeThinU , fullSvd);
115 jacobisvd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd);
116 jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeFullU | ComputeThinV);
117 jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeFullV);
118 jacobisvd_solve<MatrixType, QRPreconditioner>(m, ComputeThinU | ComputeThinV);
119
120 // test reconstruction
121 typedef typename MatrixType::Index Index;
122 Index diagSize = (std::min)(m.rows(), m.cols());
123 JacobiSVD<MatrixType, QRPreconditioner> svd(m, ComputeThinU | ComputeThinV);
124 VERIFY_IS_APPROX(m, svd.matrixU().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint());
125 }
126 }
127
128 template<typename MatrixType>
jacobisvd(const MatrixType & a=MatrixType (),bool pickrandom=true)129 void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true)
130 {
131 MatrixType m = pickrandom ? MatrixType::Random(a.rows(), a.cols()) : a;
132
133 jacobisvd_test_all_computation_options<MatrixType, FullPivHouseholderQRPreconditioner>(m);
134 jacobisvd_test_all_computation_options<MatrixType, ColPivHouseholderQRPreconditioner>(m);
135 jacobisvd_test_all_computation_options<MatrixType, HouseholderQRPreconditioner>(m);
136 jacobisvd_test_all_computation_options<MatrixType, NoQRPreconditioner>(m);
137 }
138
jacobisvd_verify_assert(const MatrixType & m)139 template<typename MatrixType> void jacobisvd_verify_assert(const MatrixType& m)
140 {
141 typedef typename MatrixType::Scalar Scalar;
142 typedef typename MatrixType::Index Index;
143 Index rows = m.rows();
144 Index cols = m.cols();
145
146 enum {
147 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
148 ColsAtCompileTime = MatrixType::ColsAtCompileTime
149 };
150
151 typedef Matrix<Scalar, RowsAtCompileTime, 1> RhsType;
152
153 RhsType rhs(rows);
154
155 JacobiSVD<MatrixType> svd;
156 VERIFY_RAISES_ASSERT(svd.matrixU())
157 VERIFY_RAISES_ASSERT(svd.singularValues())
158 VERIFY_RAISES_ASSERT(svd.matrixV())
159 VERIFY_RAISES_ASSERT(svd.solve(rhs))
160
161 MatrixType a = MatrixType::Zero(rows, cols);
162 a.setZero();
163 svd.compute(a, 0);
164 VERIFY_RAISES_ASSERT(svd.matrixU())
165 VERIFY_RAISES_ASSERT(svd.matrixV())
166 svd.singularValues();
167 VERIFY_RAISES_ASSERT(svd.solve(rhs))
168
169 if (ColsAtCompileTime == Dynamic)
170 {
171 svd.compute(a, ComputeThinU);
172 svd.matrixU();
173 VERIFY_RAISES_ASSERT(svd.matrixV())
174 VERIFY_RAISES_ASSERT(svd.solve(rhs))
175
176 svd.compute(a, ComputeThinV);
177 svd.matrixV();
178 VERIFY_RAISES_ASSERT(svd.matrixU())
179 VERIFY_RAISES_ASSERT(svd.solve(rhs))
180
181 JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner> svd_fullqr;
182 VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeFullU|ComputeThinV))
183 VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeThinV))
184 VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeFullV))
185 }
186 else
187 {
188 VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinU))
189 VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinV))
190 }
191 }
192
193 template<typename MatrixType>
jacobisvd_method()194 void jacobisvd_method()
195 {
196 enum { Size = MatrixType::RowsAtCompileTime };
197 typedef typename MatrixType::RealScalar RealScalar;
198 typedef Matrix<RealScalar, Size, 1> RealVecType;
199 MatrixType m = MatrixType::Identity();
200 VERIFY_IS_APPROX(m.jacobiSvd().singularValues(), RealVecType::Ones());
201 VERIFY_RAISES_ASSERT(m.jacobiSvd().matrixU());
202 VERIFY_RAISES_ASSERT(m.jacobiSvd().matrixV());
203 VERIFY_IS_APPROX(m.jacobiSvd(ComputeFullU|ComputeFullV).solve(m), m);
204 }
205
206 // work around stupid msvc error when constructing at compile time an expression that involves
207 // a division by zero, even if the numeric type has floating point
208 template<typename Scalar>
zero()209 EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); }
210
211 // workaround aggressive optimization in ICC
sub(T a,T b)212 template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; }
213
214 template<typename MatrixType>
jacobisvd_inf_nan()215 void jacobisvd_inf_nan()
216 {
217 // all this function does is verify we don't iterate infinitely on nan/inf values
218
219 JacobiSVD<MatrixType> svd;
220 typedef typename MatrixType::Scalar Scalar;
221 Scalar some_inf = Scalar(1) / zero<Scalar>();
222 VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf));
223 svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV);
224
225 Scalar some_nan = zero<Scalar>() / zero<Scalar>();
226 VERIFY(some_nan != some_nan);
227 svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV);
228
229 MatrixType m = MatrixType::Zero(10,10);
230 m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf;
231 svd.compute(m, ComputeFullU | ComputeFullV);
232
233 m = MatrixType::Zero(10,10);
234 m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_nan;
235 svd.compute(m, ComputeFullU | ComputeFullV);
236 }
237
238 // Regression test for bug 286: JacobiSVD loops indefinitely with some
239 // matrices containing denormal numbers.
jacobisvd_bug286()240 void jacobisvd_bug286()
241 {
242 #if defined __INTEL_COMPILER
243 // shut up warning #239: floating point underflow
244 #pragma warning push
245 #pragma warning disable 239
246 #endif
247 Matrix2d M;
248 M << -7.90884e-313, -4.94e-324,
249 0, 5.60844e-313;
250 #if defined __INTEL_COMPILER
251 #pragma warning pop
252 #endif
253 JacobiSVD<Matrix2d> svd;
254 svd.compute(M); // just check we don't loop indefinitely
255 }
256
jacobisvd_preallocate()257 void jacobisvd_preallocate()
258 {
259 Vector3f v(3.f, 2.f, 1.f);
260 MatrixXf m = v.asDiagonal();
261
262 internal::set_is_malloc_allowed(false);
263 VERIFY_RAISES_ASSERT(VectorXf v(10);)
264 JacobiSVD<MatrixXf> svd;
265 internal::set_is_malloc_allowed(true);
266 svd.compute(m);
267 VERIFY_IS_APPROX(svd.singularValues(), v);
268
269 JacobiSVD<MatrixXf> svd2(3,3);
270 internal::set_is_malloc_allowed(false);
271 svd2.compute(m);
272 internal::set_is_malloc_allowed(true);
273 VERIFY_IS_APPROX(svd2.singularValues(), v);
274 VERIFY_RAISES_ASSERT(svd2.matrixU());
275 VERIFY_RAISES_ASSERT(svd2.matrixV());
276 svd2.compute(m, ComputeFullU | ComputeFullV);
277 VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
278 VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
279 internal::set_is_malloc_allowed(false);
280 svd2.compute(m);
281 internal::set_is_malloc_allowed(true);
282
283 JacobiSVD<MatrixXf> svd3(3,3,ComputeFullU|ComputeFullV);
284 internal::set_is_malloc_allowed(false);
285 svd2.compute(m);
286 internal::set_is_malloc_allowed(true);
287 VERIFY_IS_APPROX(svd2.singularValues(), v);
288 VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
289 VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
290 internal::set_is_malloc_allowed(false);
291 svd2.compute(m, ComputeFullU|ComputeFullV);
292 internal::set_is_malloc_allowed(true);
293 }
294
test_jacobisvd()295 void test_jacobisvd()
296 {
297 CALL_SUBTEST_3(( jacobisvd_verify_assert(Matrix3f()) ));
298 CALL_SUBTEST_4(( jacobisvd_verify_assert(Matrix4d()) ));
299 CALL_SUBTEST_7(( jacobisvd_verify_assert(MatrixXf(10,12)) ));
300 CALL_SUBTEST_8(( jacobisvd_verify_assert(MatrixXcd(7,5)) ));
301
302 for(int i = 0; i < g_repeat; i++) {
303 Matrix2cd m;
304 m << 0, 1,
305 0, 1;
306 CALL_SUBTEST_1(( jacobisvd(m, false) ));
307 m << 1, 0,
308 1, 0;
309 CALL_SUBTEST_1(( jacobisvd(m, false) ));
310
311 Matrix2d n;
312 n << 0, 0,
313 0, 0;
314 CALL_SUBTEST_2(( jacobisvd(n, false) ));
315 n << 0, 0,
316 0, 1;
317 CALL_SUBTEST_2(( jacobisvd(n, false) ));
318
319 CALL_SUBTEST_3(( jacobisvd<Matrix3f>() ));
320 CALL_SUBTEST_4(( jacobisvd<Matrix4d>() ));
321 CALL_SUBTEST_5(( jacobisvd<Matrix<float,3,5> >() ));
322 CALL_SUBTEST_6(( jacobisvd<Matrix<double,Dynamic,2> >(Matrix<double,Dynamic,2>(10,2)) ));
323
324 int r = internal::random<int>(1, 30),
325 c = internal::random<int>(1, 30);
326 CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(r,c)) ));
327 CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(r,c)) ));
328 (void) r;
329 (void) c;
330
331 // Test on inf/nan matrix
332 CALL_SUBTEST_7( jacobisvd_inf_nan<MatrixXf>() );
333 }
334
335 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))) ));
336 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))) ));
337
338 // test matrixbase method
339 CALL_SUBTEST_1(( jacobisvd_method<Matrix2cd>() ));
340 CALL_SUBTEST_3(( jacobisvd_method<Matrix3f>() ));
341
342 // Test problem size constructors
343 CALL_SUBTEST_7( JacobiSVD<MatrixXf>(10,10) );
344
345 // Check that preallocation avoids subsequent mallocs
346 CALL_SUBTEST_9( jacobisvd_preallocate() );
347
348 // Regression check for bug 286
349 CALL_SUBTEST_2( jacobisvd_bug286() );
350 }
351