ground/gcs/src/libs/eigen/test/svd_fill.h

   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   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