mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-07 18:27:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

Extend unit tests of sefladjoint-eigensolver

Browse Source
This commit is contained in:
Gael Guennebaud 2015-05-07 15:54:07 +02:00
parent ebf8ca4fa8
commit c2107d30ce
3 changed files with 92 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
test/eigensolver_selfadjoint.cpp
View File

@ -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 <limits>
#include <Eigen/Eigenvalues>
template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
{
typedef typename MatrixType::Index Index;
@ -31,17 +34,8 @@ template<typename MatrixType> 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<Index>(-cols,cols);
for(Index k=0; k<count; ++k)
{
Index i = internal::random<Index>(0,cols-1);
symmA.row(i).setZero();
symmA.col(i).setZero();
}
}
svd_fill_random(symmA,Symmetric);
symmA.template triangularView<StrictlyUpper>().setZero();
symmC.template triangularView<StrictlyUpper>().setZero();

61
test/svd_common.h
View File

@ -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<typename SvdType, typename MatrixType>
@ -257,65 +259,6 @@ void svd_test_all_computation_options(const MatrixType& m, bool full_only)
}
}
template<typename MatrixType>
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<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)*std::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();
}
m = U * d.asDiagonal() * VT;
// (partly) cancel some coeffs
if(!(dup && unit_uv))
{
Matrix<Scalar,Dynamic,1> samples(7);
samples << 0, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -1./NumTraits<RealScalar>::highest(), 1./NumTraits<RealScalar>::highest();
Index n = internal::random<Index>(0,m.size()-1);
for(Index i=0; i<n; ++i)
m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = samples(internal::random<Index>(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

85
test/svd_fill.h Normal file
View File

@ -0,0 +1,85 @@
// 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 MatrixType>
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<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)*std::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();
}
if(Option==Symmetric)
{
m = U * d.asDiagonal() * U.transpose();
// randomly nullify some rows/columns
{
Index count = internal::random<Index>(-1,1);
for(Index k=0; k<count; ++k)
{
Index i = internal::random<Index>(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<Scalar,Dynamic,1> samples(7);
samples << 0, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -1./NumTraits<RealScalar>::highest(), 1./NumTraits<RealScalar>::highest();
Index n = internal::random<Index>(0,m.size()-1);
for(Index i=0; i<n; ++i)
m(internal::random<Index>(0,m.rows()-1), internal::random<Index>(0,m.cols()-1)) = samples(internal::random<Index>(0,6));
}
}
}