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
17template<typename MatrixType, int QRPreconditioner>
18void 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
46template<typename MatrixType, int QRPreconditioner>
47void 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
69template<typename MatrixType, int QRPreconditioner>
70void 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
92template<typename MatrixType, int QRPreconditioner>
93void 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
128template<typename MatrixType>
129void 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
139template<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
193template<typename MatrixType>
194void 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
208template<typename Scalar>
209EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); }
210
211// workaround aggressive optimization in ICC
212template<typename T> EIGEN_DONT_INLINE  T sub(T a, T b) { return a - b; }
213
214template<typename MatrixType>
215void 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.
240void 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
257void 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
295void 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