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<>
four_denorms()14 Array4f four_denorms() { return Array4f(5.60844e-39f, -5.60844e-39f, 4.94e-44f, -4.94e-44f); }
15 template<>
four_denorms()16 Array4d four_denorms() { return Array4d(5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324); }
17 template<typename T>
four_denorms()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 typedef typename MatrixType::Index Index;
27 Index diagSize = (std::min)(m.rows(), m.cols());
28 RealScalar s = std::numeric_limits<RealScalar>::max_exponent10/4;
29 s = internal::random<RealScalar>(1,s);
30 Matrix<RealScalar,Dynamic,1> d = Matrix<RealScalar,Dynamic,1>::Random(diagSize);
31 for(Index k=0; k<diagSize; ++k)
32 d(k) = d(k)*pow(RealScalar(10),internal::random<RealScalar>(-s,s));
33
34 bool dup = internal::random<int>(0,10) < 3;
35 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
36
37 // duplicate some singular values
38 if(dup)
39 {
40 Index n = internal::random<Index>(0,d.size()-1);
41 for(Index i=0; i<n; ++i)
42 d(internal::random<Index>(0,d.size()-1)) = d(internal::random<Index>(0,d.size()-1));
43 }
44
45 Matrix<Scalar,Dynamic,Dynamic> U(m.rows(),diagSize);
46 Matrix<Scalar,Dynamic,Dynamic> VT(diagSize,m.cols());
47 if(unit_uv)
48 {
49 // in very rare cases let's try with a pure diagonal matrix
50 if(internal::random<int>(0,10) < 1)
51 {
52 U.setIdentity();
53 VT.setIdentity();
54 }
55 else
56 {
57 createRandomPIMatrixOfRank(diagSize,U.rows(), U.cols(), U);
58 createRandomPIMatrixOfRank(diagSize,VT.rows(), VT.cols(), VT);
59 }
60 }
61 else
62 {
63 U.setRandom();
64 VT.setRandom();
65 }
66
67 Matrix<Scalar,Dynamic,1> samples(9);
68 samples << 0, four_denorms<RealScalar>(),
69 -RealScalar(1)/NumTraits<RealScalar>::highest(), RealScalar(1)/NumTraits<RealScalar>::highest(), (std::numeric_limits<RealScalar>::min)(), pow((std::numeric_limits<RealScalar>::min)(),0.8);
70
71 if(Option==Symmetric)
72 {
73 m = U * d.asDiagonal() * U.transpose();
74
75 // randomly nullify some rows/columns
76 {
77 Index count = internal::random<Index>(-diagSize,diagSize);
78 for(Index k=0; k<count; ++k)
79 {
80 Index i = internal::random<Index>(0,diagSize-1);
81 m.row(i).setZero();
82 m.col(i).setZero();
83 }
84 if(count<0)
85 // (partly) cancel some coeffs
86 if(!(dup && unit_uv))
87 {
88
89 Index n = internal::random<Index>(0,m.size()-1);
90 for(Index k=0; k<n; ++k)
91 {
92 Index i = internal::random<Index>(0,m.rows()-1);
93 Index j = internal::random<Index>(0,m.cols()-1);
94 m(j,i) = m(i,j) = samples(internal::random<Index>(0,samples.size()-1));
95 if(NumTraits<Scalar>::IsComplex)
96 *(&numext::real_ref(m(j,i))+1) = *(&numext::real_ref(m(i,j))+1) = samples.real()(internal::random<Index>(0,samples.size()-1));
97 }
98 }
99 }
100 }
101 else
102 {
103 m = U * d.asDiagonal() * VT;
104 // (partly) cancel some coeffs
105 if(!(dup && unit_uv))
106 {
107 Index n = internal::random<Index>(0,m.size()-1);
108 for(Index k=0; k<n; ++k)
109 {
110 Index i = internal::random<Index>(0,m.rows()-1);
111 Index j = internal::random<Index>(0,m.cols()-1);
112 m(i,j) = samples(internal::random<Index>(0,samples.size()-1));
113 if(NumTraits<Scalar>::IsComplex)
114 *(&numext::real_ref(m(i,j))+1) = samples.real()(internal::random<Index>(0,samples.size()-1));
115 }
116 }
117 }
118 }
119
120