Squashed 'third_party/eigen/' changes from 61d72f6..cf794d3
Change-Id: I9b814151b01f49af6337a8605d0c42a3a1ed4c72
git-subtree-dir: third_party/eigen
git-subtree-split: cf794d3b741a6278df169e58461f8529f43bce5d
diff --git a/test/svd_fill.h b/test/svd_fill.h
new file mode 100644
index 0000000..d68647e
--- /dev/null
+++ b/test/svd_fill.h
@@ -0,0 +1,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));
+ }
+ }
+ }
+}
+