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