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
262