1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2014 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 #ifndef SVD_DEFAULT
12 #error a macro SVD_DEFAULT(MatrixType) must be defined prior to including svd_common.h
13 #endif
14
15 #ifndef SVD_FOR_MIN_NORM
16 #error a macro SVD_FOR_MIN_NORM(MatrixType) must be defined prior to including svd_common.h
17 #endif
18
19 #include "svd_fill.h"
20 #include "solverbase.h"
21
22 // Check that the matrix m is properly reconstructed and that the U and V factors are unitary
23 // The SVD must have already been computed.
24 template<typename SvdType, typename MatrixType>
svd_check_full(const MatrixType & m,const SvdType & svd)25 void svd_check_full(const MatrixType& m, const SvdType& svd)
26 {
27 Index rows = m.rows();
28 Index cols = m.cols();
29
30 enum {
31 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
32 ColsAtCompileTime = MatrixType::ColsAtCompileTime
33 };
34
35 typedef typename MatrixType::Scalar Scalar;
36 typedef typename MatrixType::RealScalar RealScalar;
37 typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType;
38 typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType;
39
40 MatrixType sigma = MatrixType::Zero(rows,cols);
41 sigma.diagonal() = svd.singularValues().template cast<Scalar>();
42 MatrixUType u = svd.matrixU();
43 MatrixVType v = svd.matrixV();
44 RealScalar scaling = m.cwiseAbs().maxCoeff();
45 if(scaling<(std::numeric_limits<RealScalar>::min)())
46 {
47 VERIFY(sigma.cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)());
48 }
49 else
50 {
51 VERIFY_IS_APPROX(m/scaling, u * (sigma/scaling) * v.adjoint());
52 }
53 VERIFY_IS_UNITARY(u);
54 VERIFY_IS_UNITARY(v);
55 }
56
57 // Compare partial SVD defined by computationOptions to a full SVD referenceSvd
58 template<typename SvdType, typename MatrixType>
svd_compare_to_full(const MatrixType & m,unsigned int computationOptions,const SvdType & referenceSvd)59 void svd_compare_to_full(const MatrixType& m,
60 unsigned int computationOptions,
61 const SvdType& referenceSvd)
62 {
63 typedef typename MatrixType::RealScalar RealScalar;
64 Index rows = m.rows();
65 Index cols = m.cols();
66 Index diagSize = (std::min)(rows, cols);
67 RealScalar prec = test_precision<RealScalar>();
68
69 SvdType svd(m, computationOptions);
70
71 VERIFY_IS_APPROX(svd.singularValues(), referenceSvd.singularValues());
72
73 if(computationOptions & (ComputeFullV|ComputeThinV))
74 {
75 VERIFY( (svd.matrixV().adjoint()*svd.matrixV()).isIdentity(prec) );
76 VERIFY_IS_APPROX( svd.matrixV().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint(),
77 referenceSvd.matrixV().leftCols(diagSize) * referenceSvd.singularValues().asDiagonal() * referenceSvd.matrixV().leftCols(diagSize).adjoint());
78 }
79
80 if(computationOptions & (ComputeFullU|ComputeThinU))
81 {
82 VERIFY( (svd.matrixU().adjoint()*svd.matrixU()).isIdentity(prec) );
83 VERIFY_IS_APPROX( svd.matrixU().leftCols(diagSize) * svd.singularValues().cwiseAbs2().asDiagonal() * svd.matrixU().leftCols(diagSize).adjoint(),
84 referenceSvd.matrixU().leftCols(diagSize) * referenceSvd.singularValues().cwiseAbs2().asDiagonal() * referenceSvd.matrixU().leftCols(diagSize).adjoint());
85 }
86
87 // The following checks are not critical.
88 // For instance, with Dived&Conquer SVD, if only the factor 'V' is computedt then different matrix-matrix product implementation will be used
89 // and the resulting 'V' factor might be significantly different when the SVD decomposition is not unique, especially with single precision float.
90 ++g_test_level;
91 if(computationOptions & ComputeFullU) VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU());
92 if(computationOptions & ComputeThinU) VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize));
93 if(computationOptions & ComputeFullV) VERIFY_IS_APPROX(svd.matrixV().cwiseAbs(), referenceSvd.matrixV().cwiseAbs());
94 if(computationOptions & ComputeThinV) VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize));
95 --g_test_level;
96 }
97
98 //
99 template<typename SvdType, typename MatrixType>
svd_least_square(const MatrixType & m,unsigned int computationOptions)100 void svd_least_square(const MatrixType& m, unsigned int computationOptions)
101 {
102 typedef typename MatrixType::Scalar Scalar;
103 typedef typename MatrixType::RealScalar RealScalar;
104 Index rows = m.rows();
105 Index cols = m.cols();
106
107 enum {
108 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
109 ColsAtCompileTime = MatrixType::ColsAtCompileTime
110 };
111
112 typedef Matrix<Scalar, RowsAtCompileTime, Dynamic> RhsType;
113 typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType;
114
115 RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols));
116 SvdType svd(m, computationOptions);
117
118 if(internal::is_same<RealScalar,double>::value) svd.setThreshold(1e-8);
119 else if(internal::is_same<RealScalar,float>::value) svd.setThreshold(2e-4);
120
121 SolutionType x = svd.solve(rhs);
122
123 RealScalar residual = (m*x-rhs).norm();
124 RealScalar rhs_norm = rhs.norm();
125 if(!test_isMuchSmallerThan(residual,rhs.norm()))
126 {
127 // ^^^ If the residual is very small, then we have an exact solution, so we are already good.
128
129 // evaluate normal equation which works also for least-squares solutions
130 if(internal::is_same<RealScalar,double>::value || svd.rank()==m.diagonal().size())
131 {
132 using std::sqrt;
133 // This test is not stable with single precision.
134 // This is probably because squaring m signicantly affects the precision.
135 if(internal::is_same<RealScalar,float>::value) ++g_test_level;
136
137 VERIFY_IS_APPROX(m.adjoint()*(m*x),m.adjoint()*rhs);
138
139 if(internal::is_same<RealScalar,float>::value) --g_test_level;
140 }
141
142 // Check that there is no significantly better solution in the neighborhood of x
143 for(Index k=0;k<x.rows();++k)
144 {
145 using std::abs;
146
147 SolutionType y(x);
148 y.row(k) = (RealScalar(1)+2*NumTraits<RealScalar>::epsilon())*x.row(k);
149 RealScalar residual_y = (m*y-rhs).norm();
150 VERIFY( test_isMuchSmallerThan(abs(residual_y-residual), rhs_norm) || residual < residual_y );
151 if(internal::is_same<RealScalar,float>::value) ++g_test_level;
152 VERIFY( test_isApprox(residual_y,residual) || residual < residual_y );
153 if(internal::is_same<RealScalar,float>::value) --g_test_level;
154
155 y.row(k) = (RealScalar(1)-2*NumTraits<RealScalar>::epsilon())*x.row(k);
156 residual_y = (m*y-rhs).norm();
157 VERIFY( test_isMuchSmallerThan(abs(residual_y-residual), rhs_norm) || residual < residual_y );
158 if(internal::is_same<RealScalar,float>::value) ++g_test_level;
159 VERIFY( test_isApprox(residual_y,residual) || residual < residual_y );
160 if(internal::is_same<RealScalar,float>::value) --g_test_level;
161 }
162 }
163 }
164
165 // check minimal norm solutions, the inoput matrix m is only used to recover problem size
166 template<typename MatrixType>
svd_min_norm(const MatrixType & m,unsigned int computationOptions)167 void svd_min_norm(const MatrixType& m, unsigned int computationOptions)
168 {
169 typedef typename MatrixType::Scalar Scalar;
170 Index cols = m.cols();
171
172 enum {
173 ColsAtCompileTime = MatrixType::ColsAtCompileTime
174 };
175
176 typedef Matrix<Scalar, ColsAtCompileTime, Dynamic> SolutionType;
177
178 // generate a full-rank m x n problem with m<n
179 enum {
180 RankAtCompileTime2 = ColsAtCompileTime==Dynamic ? Dynamic : (ColsAtCompileTime)/2+1,
181 RowsAtCompileTime3 = ColsAtCompileTime==Dynamic ? Dynamic : ColsAtCompileTime+1
182 };
183 typedef Matrix<Scalar, RankAtCompileTime2, ColsAtCompileTime> MatrixType2;
184 typedef Matrix<Scalar, RankAtCompileTime2, 1> RhsType2;
185 typedef Matrix<Scalar, ColsAtCompileTime, RankAtCompileTime2> MatrixType2T;
186 Index rank = RankAtCompileTime2==Dynamic ? internal::random<Index>(1,cols) : Index(RankAtCompileTime2);
187 MatrixType2 m2(rank,cols);
188 int guard = 0;
189 do {
190 m2.setRandom();
191 } while(SVD_FOR_MIN_NORM(MatrixType2)(m2).setThreshold(test_precision<Scalar>()).rank()!=rank && (++guard)<10);
192 VERIFY(guard<10);
193
194 RhsType2 rhs2 = RhsType2::Random(rank);
195 // use QR to find a reference minimal norm solution
196 HouseholderQR<MatrixType2T> qr(m2.adjoint());
197 Matrix<Scalar,Dynamic,1> tmp = qr.matrixQR().topLeftCorner(rank,rank).template triangularView<Upper>().adjoint().solve(rhs2);
198 tmp.conservativeResize(cols);
199 tmp.tail(cols-rank).setZero();
200 SolutionType x21 = qr.householderQ() * tmp;
201 // now check with SVD
202 SVD_FOR_MIN_NORM(MatrixType2) svd2(m2, computationOptions);
203 SolutionType x22 = svd2.solve(rhs2);
204 VERIFY_IS_APPROX(m2*x21, rhs2);
205 VERIFY_IS_APPROX(m2*x22, rhs2);
206 VERIFY_IS_APPROX(x21, x22);
207
208 // Now check with a rank deficient matrix
209 typedef Matrix<Scalar, RowsAtCompileTime3, ColsAtCompileTime> MatrixType3;
210 typedef Matrix<Scalar, RowsAtCompileTime3, 1> RhsType3;
211 Index rows3 = RowsAtCompileTime3==Dynamic ? internal::random<Index>(rank+1,2*cols) : Index(RowsAtCompileTime3);
212 Matrix<Scalar,RowsAtCompileTime3,Dynamic> C = Matrix<Scalar,RowsAtCompileTime3,Dynamic>::Random(rows3,rank);
213 MatrixType3 m3 = C * m2;
214 RhsType3 rhs3 = C * rhs2;
215 SVD_FOR_MIN_NORM(MatrixType3) svd3(m3, computationOptions);
216 SolutionType x3 = svd3.solve(rhs3);
217 VERIFY_IS_APPROX(m3*x3, rhs3);
218 VERIFY_IS_APPROX(m3*x21, rhs3);
219 VERIFY_IS_APPROX(m2*x3, rhs2);
220 VERIFY_IS_APPROX(x21, x3);
221 }
222
223 template<typename MatrixType, typename SolverType>
svd_test_solvers(const MatrixType & m,const SolverType & solver)224 void svd_test_solvers(const MatrixType& m, const SolverType& solver) {
225 Index rows, cols, cols2;
226
227 rows = m.rows();
228 cols = m.cols();
229
230 if(MatrixType::ColsAtCompileTime==Dynamic)
231 {
232 cols2 = internal::random<int>(2,EIGEN_TEST_MAX_SIZE);
233 }
234 else
235 {
236 cols2 = cols;
237 }
238 typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> CMatrixType;
239 check_solverbase<CMatrixType, MatrixType>(m, solver, rows, cols, cols2);
240 }
241
242 // Check full, compare_to_full, least_square, and min_norm for all possible compute-options
243 template<typename SvdType, typename MatrixType>
svd_test_all_computation_options(const MatrixType & m,bool full_only)244 void svd_test_all_computation_options(const MatrixType& m, bool full_only)
245 {
246 // if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols())
247 // return;
248 STATIC_CHECK(( internal::is_same<typename SvdType::StorageIndex,int>::value ));
249
250 SvdType fullSvd(m, ComputeFullU|ComputeFullV);
251 CALL_SUBTEST(( svd_check_full(m, fullSvd) ));
252 CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeFullU | ComputeFullV) ));
253 CALL_SUBTEST(( svd_min_norm(m, ComputeFullU | ComputeFullV) ));
254
255 #if defined __INTEL_COMPILER
256 // remark #111: statement is unreachable
257 #pragma warning disable 111
258 #endif
259
260 svd_test_solvers(m, fullSvd);
261
262 if(full_only)
263 return;
264
265 CALL_SUBTEST(( svd_compare_to_full(m, ComputeFullU, fullSvd) ));
266 CALL_SUBTEST(( svd_compare_to_full(m, ComputeFullV, fullSvd) ));
267 CALL_SUBTEST(( svd_compare_to_full(m, 0, fullSvd) ));
268
269 if (MatrixType::ColsAtCompileTime == Dynamic) {
270 // thin U/V are only available with dynamic number of columns
271 CALL_SUBTEST(( svd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd) ));
272 CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinV, fullSvd) ));
273 CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd) ));
274 CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinU , fullSvd) ));
275 CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd) ));
276
277 CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeFullU | ComputeThinV) ));
278 CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeThinU | ComputeFullV) ));
279 CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeThinU | ComputeThinV) ));
280
281 CALL_SUBTEST(( svd_min_norm(m, ComputeFullU | ComputeThinV) ));
282 CALL_SUBTEST(( svd_min_norm(m, ComputeThinU | ComputeFullV) ));
283 CALL_SUBTEST(( svd_min_norm(m, ComputeThinU | ComputeThinV) ));
284
285 // test reconstruction
286 Index diagSize = (std::min)(m.rows(), m.cols());
287 SvdType svd(m, ComputeThinU | ComputeThinV);
288 VERIFY_IS_APPROX(m, svd.matrixU().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint());
289 }
290 }
291
292
293 // work around stupid msvc error when constructing at compile time an expression that involves
294 // a division by zero, even if the numeric type has floating point
295 template<typename Scalar>
zero()296 EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); }
297
298 // workaround aggressive optimization in ICC
sub(T a,T b)299 template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; }
300
301 // This function verifies we don't iterate infinitely on nan/inf values,
302 // and that info() returns InvalidInput.
303 template<typename SvdType, typename MatrixType>
svd_inf_nan()304 void svd_inf_nan()
305 {
306 SvdType svd;
307 typedef typename MatrixType::Scalar Scalar;
308 Scalar some_inf = Scalar(1) / zero<Scalar>();
309 VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf));
310 svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV);
311 VERIFY(svd.info() == InvalidInput);
312
313 Scalar nan = std::numeric_limits<Scalar>::quiet_NaN();
314 VERIFY(nan != nan);
315 svd.compute(MatrixType::Constant(10,10,nan), ComputeFullU | ComputeFullV);
316 VERIFY(svd.info() == InvalidInput);
317
318 MatrixType m = MatrixType::Zero(10,10);
319 m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf;
320 svd.compute(m, ComputeFullU | ComputeFullV);
321 VERIFY(svd.info() == InvalidInput);
322
323 m = MatrixType::Zero(10,10);
324 m(internal::random<int>(0,9), internal::random<int>(0,9)) = nan;
325 svd.compute(m, ComputeFullU | ComputeFullV);
326 VERIFY(svd.info() == InvalidInput);
327
328 // regression test for bug 791
329 m.resize(3,3);
330 m << 0, 2*NumTraits<Scalar>::epsilon(), 0.5,
331 0, -0.5, 0,
332 nan, 0, 0;
333 svd.compute(m, ComputeFullU | ComputeFullV);
334 VERIFY(svd.info() == InvalidInput);
335
336 m.resize(4,4);
337 m << 1, 0, 0, 0,
338 0, 3, 1, 2e-308,
339 1, 0, 1, nan,
340 0, nan, nan, 0;
341 svd.compute(m, ComputeFullU | ComputeFullV);
342 VERIFY(svd.info() == InvalidInput);
343 }
344
345 // Regression test for bug 286: JacobiSVD loops indefinitely with some
346 // matrices containing denormal numbers.
347 template<typename>
svd_underoverflow()348 void svd_underoverflow()
349 {
350 #if defined __INTEL_COMPILER
351 // shut up warning #239: floating point underflow
352 #pragma warning push
353 #pragma warning disable 239
354 #endif
355 Matrix2d M;
356 M << -7.90884e-313, -4.94e-324,
357 0, 5.60844e-313;
358 SVD_DEFAULT(Matrix2d) svd;
359 svd.compute(M,ComputeFullU|ComputeFullV);
360 CALL_SUBTEST( svd_check_full(M,svd) );
361
362 // Check all 2x2 matrices made with the following coefficients:
363 VectorXd value_set(9);
364 value_set << 0, 1, -1, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -4.94e-223, 4.94e-223;
365 Array4i id(0,0,0,0);
366 int k = 0;
367 do
368 {
369 M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3));
370 svd.compute(M,ComputeFullU|ComputeFullV);
371 CALL_SUBTEST( svd_check_full(M,svd) );
372
373 id(k)++;
374 if(id(k)>=value_set.size())
375 {
376 while(k<3 && id(k)>=value_set.size()) id(++k)++;
377 id.head(k).setZero();
378 k=0;
379 }
380
381 } while((id<int(value_set.size())).all());
382
383 #if defined __INTEL_COMPILER
384 #pragma warning pop
385 #endif
386
387 // Check for overflow:
388 Matrix3d M3;
389 M3 << 4.4331978442502944e+307, -5.8585363752028680e+307, 6.4527017443412964e+307,
390 3.7841695601406358e+307, 2.4331702789740617e+306, -3.5235707140272905e+307,
391 -8.7190887618028355e+307, -7.3453213709232193e+307, -2.4367363684472105e+307;
392
393 SVD_DEFAULT(Matrix3d) svd3;
394 svd3.compute(M3,ComputeFullU|ComputeFullV); // just check we don't loop indefinitely
395 CALL_SUBTEST( svd_check_full(M3,svd3) );
396 }
397
398 // void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true)
399
400 template<typename MatrixType>
svd_all_trivial_2x2(void (* cb)(const MatrixType &,bool))401 void svd_all_trivial_2x2( void (*cb)(const MatrixType&,bool) )
402 {
403 MatrixType M;
404 VectorXd value_set(3);
405 value_set << 0, 1, -1;
406 Array4i id(0,0,0,0);
407 int k = 0;
408 do
409 {
410 M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3));
411
412 cb(M,false);
413
414 id(k)++;
415 if(id(k)>=value_set.size())
416 {
417 while(k<3 && id(k)>=value_set.size()) id(++k)++;
418 id.head(k).setZero();
419 k=0;
420 }
421
422 } while((id<int(value_set.size())).all());
423 }
424
425 template<typename>
svd_preallocate()426 void svd_preallocate()
427 {
428 Vector3f v(3.f, 2.f, 1.f);
429 MatrixXf m = v.asDiagonal();
430
431 internal::set_is_malloc_allowed(false);
432 VERIFY_RAISES_ASSERT(VectorXf tmp(10);)
433 SVD_DEFAULT(MatrixXf) svd;
434 internal::set_is_malloc_allowed(true);
435 svd.compute(m);
436 VERIFY_IS_APPROX(svd.singularValues(), v);
437
438 SVD_DEFAULT(MatrixXf) svd2(3,3);
439 internal::set_is_malloc_allowed(false);
440 svd2.compute(m);
441 internal::set_is_malloc_allowed(true);
442 VERIFY_IS_APPROX(svd2.singularValues(), v);
443 VERIFY_RAISES_ASSERT(svd2.matrixU());
444 VERIFY_RAISES_ASSERT(svd2.matrixV());
445 svd2.compute(m, ComputeFullU | ComputeFullV);
446 VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
447 VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
448 internal::set_is_malloc_allowed(false);
449 svd2.compute(m);
450 internal::set_is_malloc_allowed(true);
451
452 SVD_DEFAULT(MatrixXf) svd3(3,3,ComputeFullU|ComputeFullV);
453 internal::set_is_malloc_allowed(false);
454 svd2.compute(m);
455 internal::set_is_malloc_allowed(true);
456 VERIFY_IS_APPROX(svd2.singularValues(), v);
457 VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
458 VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
459 internal::set_is_malloc_allowed(false);
460 svd2.compute(m, ComputeFullU|ComputeFullV);
461 internal::set_is_malloc_allowed(true);
462 }
463
464 template<typename SvdType,typename MatrixType>
465 void svd_verify_assert(const MatrixType& m, bool fullOnly = false)
466 {
467 typedef typename MatrixType::Scalar Scalar;
468 Index rows = m.rows();
469 Index cols = m.cols();
470
471 enum {
472 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
473 ColsAtCompileTime = MatrixType::ColsAtCompileTime
474 };
475
476 typedef Matrix<Scalar, RowsAtCompileTime, 1> RhsType;
477 RhsType rhs(rows);
478 SvdType svd;
479 VERIFY_RAISES_ASSERT(svd.matrixU())
480 VERIFY_RAISES_ASSERT(svd.singularValues())
481 VERIFY_RAISES_ASSERT(svd.matrixV())
482 VERIFY_RAISES_ASSERT(svd.solve(rhs))
483 VERIFY_RAISES_ASSERT(svd.transpose().solve(rhs))
484 VERIFY_RAISES_ASSERT(svd.adjoint().solve(rhs))
485 MatrixType a = MatrixType::Zero(rows, cols);
486 a.setZero();
487 svd.compute(a, 0);
488 VERIFY_RAISES_ASSERT(svd.matrixU())
489 VERIFY_RAISES_ASSERT(svd.matrixV())
490 svd.singularValues();
491 VERIFY_RAISES_ASSERT(svd.solve(rhs))
492
493 svd.compute(a, ComputeFullU);
494 svd.matrixU();
495 VERIFY_RAISES_ASSERT(svd.matrixV())
496 VERIFY_RAISES_ASSERT(svd.solve(rhs))
497 svd.compute(a, ComputeFullV);
498 svd.matrixV();
499 VERIFY_RAISES_ASSERT(svd.matrixU())
500 VERIFY_RAISES_ASSERT(svd.solve(rhs))
501
502 if (!fullOnly && ColsAtCompileTime == Dynamic)
503 {
504 svd.compute(a, ComputeThinU);
505 svd.matrixU();
506 VERIFY_RAISES_ASSERT(svd.matrixV())
507 VERIFY_RAISES_ASSERT(svd.solve(rhs))
508 svd.compute(a, ComputeThinV);
509 svd.matrixV();
510 VERIFY_RAISES_ASSERT(svd.matrixU())
511 VERIFY_RAISES_ASSERT(svd.solve(rhs))
512 }
513 else
514 {
515 VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinU))
516 VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinV))
517 }
518 }
519
520 #undef SVD_DEFAULT
521 #undef SVD_FOR_MIN_NORM
522