svd_common.h (19317B)
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> 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> 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> 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> 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> 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> 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> 296 EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); } 297 298 // workaround aggressive optimization in ICC 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> 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> 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> 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> 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