svd_fill.h
Go to the documentation of this file.
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));
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));
113  *(&numext::real_ref(m(i,j))+1) = samples.real()(internal::random<Index>(0,samples.size()-1));
114  }
115  }
116  }
117 }
118 
Matrix3f m
SCALAR Scalar
Definition: bench_gemm.cpp:46
void svd_fill_random(MatrixType &m, int Option=0)
Definition: svd_fill.h:21
EIGEN_DEVICE_FUNC internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) >::type real_ref(const Scalar &x)
#define min(a, b)
Definition: datatypes.h:19
int n
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_pow_op< typename Derived::Scalar, typename ExponentDerived::Scalar >, const Derived, const ExponentDerived > pow(const Eigen::ArrayBase< Derived > &x, const Eigen::ArrayBase< ExponentDerived > &exponents)
MatrixXf MatrixType
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Array< double, 1, 3 > e(1./3., 0.5, 2.)
RealScalar s
Array< T, 4, 1 > four_denorms()
Definition: svd_fill.h:18
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
void createRandomPIMatrixOfRank(Index desired_rank, Index rows, Index cols, MatrixType &m)
Definition: main.h:653
Jet< T, N > pow(const Jet< T, N > &f, double g)
Definition: jet.h:570
std::ptrdiff_t j
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy y set format x g set format y g set format x2 g set format y2 g set format z g set angles radians set nogrid set key title set key left top Right noreverse box linetype linewidth samplen spacing width set nolabel set noarrow set nologscale set logscale x set set pointsize set encoding default set nopolar set noparametric set set samples


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:36:29