svd_fill.h (4050B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2014-2015 Gael Guennebaud <gael.guennebaud@inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 template<typename T> 11 Array<T,4,1> four_denorms(); 12 13 template<> 14 Array4f four_denorms() { return Array4f(5.60844e-39f, -5.60844e-39f, 4.94e-44f, -4.94e-44f); } 15 template<> 16 Array4d four_denorms() { return Array4d(5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324); } 17 template<typename T> 18 Array<T,4,1> four_denorms() { return four_denorms<double>().cast<T>(); } 19 20 template<typename MatrixType> 21 void svd_fill_random(MatrixType &m, int Option = 0) 22 { 23 using std::pow; 24 typedef typename MatrixType::Scalar Scalar; 25 typedef typename MatrixType::RealScalar RealScalar; 26 Index diagSize = (std::min)(m.rows(), m.cols()); 27 RealScalar s = std::numeric_limits<RealScalar>::max_exponent10/4; 28 s = internal::random<RealScalar>(1,s); 29 Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(diagSize); 30 for(Index k=0; k<diagSize; ++k) 31 d(k) = d(k)*pow(RealScalar(10),internal::random<RealScalar>(-s,s)); 32 33 bool dup = internal::random<int>(0,10) < 3; 34 bool unit_uv = internal::random<int>(0,10) < (dup?7:3); // if we duplicate some diagonal entries, then increase the chance to preserve them using unitary U and V factors 35 36 // duplicate some singular values 37 if(dup) 38 { 39 Index n = internal::random<Index>(0,d.size()-1); 40 for(Index i=0; i<n; ++i) 41 d(internal::random<Index>(0,d.size()-1)) = d(internal::random<Index>(0,d.size()-1)); 42 } 43 44 Matrix<Scalar,Dynamic,Dynamic> U(m.rows(),diagSize); 45 Matrix<Scalar,Dynamic,Dynamic> VT(diagSize,m.cols()); 46 if(unit_uv) 47 { 48 // in very rare cases let's try with a pure diagonal matrix 49 if(internal::random<int>(0,10) < 1) 50 { 51 U.setIdentity(); 52 VT.setIdentity(); 53 } 54 else 55 { 56 createRandomPIMatrixOfRank(diagSize,U.rows(), U.cols(), U); 57 createRandomPIMatrixOfRank(diagSize,VT.rows(), VT.cols(), VT); 58 } 59 } 60 else 61 { 62 U.setRandom(); 63 VT.setRandom(); 64 } 65 66 Matrix<Scalar,Dynamic,1> samples(9); 67 samples << 0, four_denorms<RealScalar>(), 68 -RealScalar(1)/NumTraits<RealScalar>::highest(), RealScalar(1)/NumTraits<RealScalar>::highest(), (std::numeric_limits<RealScalar>::min)(), pow((std::numeric_limits<RealScalar>::min)(),0.8); 69 70 if(Option==Symmetric) 71 { 72 m = U * d.asDiagonal() * U.transpose(); 73 74 // randomly nullify some rows/columns 75 { 76 Index count = internal::random<Index>(-diagSize,diagSize); 77 for(Index k=0; k<count; ++k) 78 { 79 Index i = internal::random<Index>(0,diagSize-1); 80 m.row(i).setZero(); 81 m.col(i).setZero(); 82 } 83 if(count<0) 84 // (partly) cancel some coeffs 85 if(!(dup && unit_uv)) 86 { 87 88 Index n = internal::random<Index>(0,m.size()-1); 89 for(Index k=0; k<n; ++k) 90 { 91 Index i = internal::random<Index>(0,m.rows()-1); 92 Index j = internal::random<Index>(0,m.cols()-1); 93 m(j,i) = m(i,j) = samples(internal::random<Index>(0,samples.size()-1)); 94 if(NumTraits<Scalar>::IsComplex) 95 *(&numext::real_ref(m(j,i))+1) = *(&numext::real_ref(m(i,j))+1) = samples.real()(internal::random<Index>(0,samples.size()-1)); 96 } 97 } 98 } 99 } 100 else 101 { 102 m = U * d.asDiagonal() * VT; 103 // (partly) cancel some coeffs 104 if(!(dup && unit_uv)) 105 { 106 Index n = internal::random<Index>(0,m.size()-1); 107 for(Index k=0; k<n; ++k) 108 { 109 Index i = internal::random<Index>(0,m.rows()-1); 110 Index j = internal::random<Index>(0,m.cols()-1); 111 m(i,j) = samples(internal::random<Index>(0,samples.size()-1)); 112 if(NumTraits<Scalar>::IsComplex) 113 *(&numext::real_ref(m(i,j))+1) = samples.real()(internal::random<Index>(0,samples.size()-1)); 114 } 115 } 116 } 117 } 118