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 // Copyright (C) 2013 Gauthier Brun <brun.gauthier (at) gmail.com> 8 // Copyright (C) 2013 Nicolas Carre <nicolas.carre (at) ensimag.fr> 9 // Copyright (C) 2013 Jean Ceccato <jean.ceccato (at) ensimag.fr> 10 // Copyright (C) 2013 Pierre Zoppitelli <pierre.zoppitelli (at) 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" 26 template<typename MatrixType, typename SVD> 27 void 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 54 template<typename MatrixType, typename SVD> 55 void 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 79 template<typename MatrixType, typename SVD> 80 void 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 105 template<typename MatrixType, typename SVD> 106 void 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 113 template<typename MatrixType, typename SVD> 114 void 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 139 template<typename MatrixType, typename SVD> 140 void 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 187 template<typename Scalar> 188 EIGEN_DONT_INLINE Scalar zero() { return Scalar(0); } 189 190 // workaround aggressive optimization in ICC 191 template<typename T> EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } 192 193 194 template<typename MatrixType, typename SVD> 195 void 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 219 template<typename SVD> 220 void 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
