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