GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-14 04:35:57 +08:00
Code Issues Packages Projects Releases Wiki Activity

Extend polynomial solver unit tests to complexes

Browse Source 
(grafted from f12b368417992f0974678646f2fb7fa2db44b633
)
This commit is contained in:
Gael Guennebaud 2016-11-23 16:05:45 +01:00
parent 222ce4b49d
commit e777674a87
1 changed files with 16 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
unsupported/test/polynomialsolver.cpp
View File

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