Home | History | Annotate | Download | only in test 
      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 (at) 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 
    120