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