mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-30 17:40:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

Extend polynomial solver unit tests to complexes

Browse Source
This commit is contained in:
Gael Guennebaud 2016-11-23 16:05:45 +01:00
parent 56e5ec07c6
commit f12b368417
1 changed files with 16 additions and 18 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
unsupported/test/polynomialsolver.cpp
View File

@ -32,9 +32,10 @@ bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve )
{
typedef typename POLYNOMIAL::Index Index;
typedef typename POLYNOMIAL::Scalar Scalar;
typedef typename POLYNOMIAL::RealScalar RealScalar;
typedef typename SOLVER::RootsType RootsType;
typedef Matrix<Scalar,Deg,1> EvalRootsType;
typedef Matrix<RealScalar,Deg,1> EvalRootsType;
const Index deg = pols.size()-1;
@ -57,7 +58,7 @@ bool aux_evalSolver( const POLYNOMIAL& pols, SOLVER& psolve )
cerr << endl;
}
std::vector<Scalar> rootModuli( roots.size() );
std::vector<RealScalar> rootModuli( roots.size() );
Map< EvalRootsType > aux( &rootModuli[0], roots.size() );
aux = roots.array().abs();
std::sort( rootModuli.begin(), rootModuli.end() );
@ -83,7 +84,7 @@ void evalSolver( const POLYNOMIAL& pols )
{
typedef typename POLYNOMIAL::Scalar Scalar;
typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType;
typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType;
PolynomialSolverType psolve;
aux_evalSolver<Deg, POLYNOMIAL, PolynomialSolverType>( pols, psolve );
@ -97,6 +98,7 @@ void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const
{
using std::sqrt;
typedef typename POLYNOMIAL::Scalar Scalar;
typedef typename POLYNOMIAL::RealScalar RealScalar;
typedef PolynomialSolver<Scalar, Deg > PolynomialSolverType;
@ -107,15 +109,12 @@ void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const
// 1) the roots found are correct
// 2) the roots have distinct moduli
typedef typename POLYNOMIAL::Scalar Scalar;
typedef typename REAL_ROOTS::Scalar Real;
//Test realRoots
std::vector< Real > calc_realRoots;
psolve.realRoots( calc_realRoots );
VERIFY( calc_realRoots.size() == (size_t)real_roots.size() );
std::vector< RealScalar > calc_realRoots;
psolve.realRoots( calc_realRoots, test_precision<RealScalar>());
VERIFY_IS_EQUAL( calc_realRoots.size() , (size_t)real_roots.size() );
const Scalar psPrec = sqrt( test_precision<Scalar>() );
const RealScalar psPrec = sqrt( test_precision<RealScalar>() );
for( size_t i=0; i<calc_realRoots.size(); ++i )
{
@ -138,7 +137,7 @@ void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const
bool hasRealRoot;
//Test absGreatestRealRoot
Real r = psolve.absGreatestRealRoot( hasRealRoot );
RealScalar r = psolve.absGreatestRealRoot( hasRealRoot );
VERIFY( hasRealRoot == (real_roots.size() > 0 ) );
if( hasRealRoot ){
VERIFY( internal::isApprox( real_roots.array().abs().maxCoeff(), abs(r), psPrec ) ); }
@ -167,9 +166,11 @@ void evalSolverSugarFunction( const POLYNOMIAL& pols, const ROOTS& roots, const
template<typename _Scalar, int _Deg>
void polynomialsolver(int deg)
{
typedef internal::increment_if_fixed_size<_Deg> Dim;
typedef typename NumTraits<_Scalar>::Real RealScalar;
typedef internal::increment_if_fixed_size<_Deg> Dim;
typedef Matrix<_Scalar,Dim::ret,1> PolynomialType;
typedef Matrix<_Scalar,_Deg,1> EvalRootsType;
typedef Matrix<RealScalar,_Deg,1> RealRootsType;
cout << "Standard cases" << endl;
PolynomialType pols = PolynomialType::Random(deg+1);
@ -182,15 +183,11 @@ void polynomialsolver(int deg)
evalSolver<_Deg,PolynomialType>( pols );
cout << "Test sugar" << endl;
EvalRootsType realRoots = EvalRootsType::Random(deg);
RealRootsType realRoots = RealRootsType::Random(deg);
roots_to_monicPolynomial( realRoots, pols );
evalSolverSugarFunction<_Deg>(
pols,
realRoots.template cast <
std::complex<
typename NumTraits<_Scalar>::Real
>
>(),
realRoots.template cast <std::complex<RealScalar> >().eval(),
realRoots );
}
@ -214,5 +211,6 @@ void test_polynomialsolver()
internal::random<int>(9,13)
)) );
CALL_SUBTEST_11((polynomialsolver<float,Dynamic>(1)) );
CALL_SUBTEST_12((polynomialsolver<std::complex<double>,Dynamic>(internal::random<int>(2,13))) );
}
}