source: pacpussensors/trunk/Vislab/lib3dv/eigen/unsupported/test/svd_common.h@ 136

Last change on this file since 136 was 136, checked in by ldecherf, 7 years ago

Doc
File size: 8.7 KB
Line 
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"
26template<typename MatrixType, typename SVD>
27void 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
54template<typename MatrixType, typename SVD>
55void 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
79template<typename MatrixType, typename SVD>
80void 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
105template<typename MatrixType, typename SVD>
106void 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
113template<typename MatrixType, typename SVD>
114void 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
139template<typename MatrixType, typename SVD>
140void 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
187template<typename Scalar>
188EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); }
189
190// workaround aggressive optimization in ICC
191template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; }
192
193
194template<typename MatrixType, typename SVD>
195void 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
219template<typename SVD>
220void 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
Note: See TracBrowser for help on using the repository browser.