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>
|
---|
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>
|
---|
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>
|
---|
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>
|
---|
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>
|
---|
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>
|
---|
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>
|
---|
188 | EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); }
|
---|
189 |
|
---|
190 | // workaround aggressive optimization in ICC
|
---|
191 | template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; }
|
---|
192 |
|
---|
193 |
|
---|
194 | template<typename MatrixType, typename SVD>
|
---|
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>
|
---|
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 |
|
---|