diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 7b0077a6d..4748bdd0b 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -9,9 +9,12 @@ // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #include "main.h" +#include "svd_fill.h" #include #include + + template void selfadjointeigensolver(const MatrixType& m) { typedef typename MatrixType::Index Index; @@ -31,17 +34,8 @@ template void selfadjointeigensolver(const MatrixType& m) MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1; MatrixType symmC = symmA; - // randomly nullify some rows/columns - { - Index count = 1;//internal::random(-cols,cols); - for(Index k=0; k(0,cols-1); - symmA.row(i).setZero(); - symmA.col(i).setZero(); - } - } - + svd_fill_random(symmA,Symmetric); + symmA.template triangularView().setZero(); symmC.template triangularView().setZero(); diff --git a/test/svd_common.h b/test/svd_common.h index b44b79124..b06a8a0f2 100644 --- a/test/svd_common.h +++ b/test/svd_common.h @@ -16,6 +16,8 @@ #error a macro SVD_FOR_MIN_NORM(MatrixType) must be defined prior to including svd_common.h #endif +#include "svd_fill.h" + // Check that the matrix m is properly reconstructed and that the U and V factors are unitary // The SVD must have already been computed. template @@ -257,65 +259,6 @@ void svd_test_all_computation_options(const MatrixType& m, bool full_only) } } -template -void svd_fill_random(MatrixType &m) -{ - typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - typedef typename MatrixType::Index Index; - Index diagSize = (std::min)(m.rows(), m.cols()); - RealScalar s = std::numeric_limits::max_exponent10/4; - s = internal::random(1,s); - Matrix d = Matrix::Random(diagSize); - for(Index k=0; k(-s,s)); - - bool dup = internal::random(0,10) < 3; - bool unit_uv = internal::random(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(0,d.size()-1); - for(Index i=0; i(0,d.size()-1)) = d(internal::random(0,d.size()-1)); - } - - Matrix U(m.rows(),diagSize); - Matrix VT(diagSize,m.cols()); - if(unit_uv) - { - // in very rare cases let's try with a pure diagonal matrix - if(internal::random(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(); - } - - m = U * d.asDiagonal() * VT; - - // (partly) cancel some coeffs - if(!(dup && unit_uv)) - { - Matrix samples(7); - samples << 0, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -1./NumTraits::highest(), 1./NumTraits::highest(); - Index n = internal::random(0,m.size()-1); - for(Index i=0; i(0,m.rows()-1), internal::random(0,m.cols()-1)) = samples(internal::random(0,6)); - } -} - // work around stupid msvc error when constructing at compile time an expression that involves // a division by zero, even if the numeric type has floating point diff --git a/test/svd_fill.h b/test/svd_fill.h new file mode 100644 index 000000000..aa75326d9 --- /dev/null +++ b/test/svd_fill.h @@ -0,0 +1,85 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2014-2015 Gael Guennebaud +// +// 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 +void svd_fill_random(MatrixType &m, int Option = 0) +{ + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + Index diagSize = (std::min)(m.rows(), m.cols()); + RealScalar s = std::numeric_limits::max_exponent10/4; + s = internal::random(1,s); + Matrix d = Matrix::Random(diagSize); + for(Index k=0; k(-s,s)); + + bool dup = internal::random(0,10) < 3; + bool unit_uv = internal::random(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(0,d.size()-1); + for(Index i=0; i(0,d.size()-1)) = d(internal::random(0,d.size()-1)); + } + + Matrix U(m.rows(),diagSize); + Matrix VT(diagSize,m.cols()); + if(unit_uv) + { + // in very rare cases let's try with a pure diagonal matrix + if(internal::random(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(); + } + + if(Option==Symmetric) + { + m = U * d.asDiagonal() * U.transpose(); + + // randomly nullify some rows/columns + { + Index count = internal::random(-1,1); + for(Index k=0; k(0,diagSize-1); + m.row(i).setZero(); + m.col(i).setZero(); + } + } + } + else + { + m = U * d.asDiagonal() * VT; + // (partly) cancel some coeffs + if(!(dup && unit_uv)) + { + Matrix samples(7); + samples << 0, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -1./NumTraits::highest(), 1./NumTraits::highest(); + Index n = internal::random(0,m.size()-1); + for(Index i=0; i(0,m.rows()-1), internal::random(0,m.cols()-1)) = samples(internal::random(0,6)); + } + } +} +