PolynomialSolver1.cpp
Go to the documentation of this file.
00001 #include <unsupported/Eigen/Polynomials>
00002 #include <vector>
00003 #include <iostream>
00004 
00005 using namespace Eigen;
00006 using namespace std;
00007 
00008 int main()
00009 {
00010   typedef Matrix<double,5,1> Vector5d;
00011 
00012   Vector5d roots = Vector5d::Random();
00013   cout << "Roots: " << roots.transpose() << endl;
00014   Eigen::Matrix<double,6,1> polynomial;
00015   roots_to_monicPolynomial( roots, polynomial );
00016 
00017   PolynomialSolver<double,5> psolve( polynomial );
00018   cout << "Complex roots: " << psolve.roots().transpose() << endl;
00019 
00020   std::vector<double> realRoots;
00021   psolve.realRoots( realRoots );
00022   Map<Vector5d> mapRR( &realRoots[0] );
00023   cout << "Real roots: " << mapRR.transpose() << endl;
00024 
00025   cout << endl;
00026   cout << "Illustration of the convergence problem with the QR algorithm: " << endl;
00027   cout << "---------------------------------------------------------------" << endl;
00028   Eigen::Matrix<float,7,1> hardCase_polynomial;
00029   hardCase_polynomial <<
00030   -0.957, 0.9219, 0.3516, 0.9453, -0.4023, -0.5508, -0.03125;
00031   cout << "Hard case polynomial defined by floats: " << hardCase_polynomial.transpose() << endl;
00032   PolynomialSolver<float,6> psolvef( hardCase_polynomial );
00033   cout << "Complex roots: " << psolvef.roots().transpose() << endl;
00034   Eigen::Matrix<float,6,1> evals;
00035   for( int i=0; i<6; ++i ){ evals[i] = std::abs( poly_eval( hardCase_polynomial, psolvef.roots()[i] ) ); }
00036   cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl;
00037 
00038   cout << "Using double's almost always solves the problem for small degrees: " << endl;
00039   cout << "-------------------------------------------------------------------" << endl;
00040   PolynomialSolver<double,6> psolve6d( hardCase_polynomial.cast<double>() );
00041   cout << "Complex roots: " << psolve6d.roots().transpose() << endl;
00042   for( int i=0; i<6; ++i )
00043   {
00044     std::complex<float> castedRoot( psolve6d.roots()[i].real(), psolve6d.roots()[i].imag() );
00045     evals[i] = std::abs( poly_eval( hardCase_polynomial, castedRoot ) );
00046   }
00047   cout << "Norms of the evaluations of the polynomial at the roots: " << evals.transpose() << endl << endl;
00048 
00049   cout.precision(10);
00050   cout << "The last root in float then in double: " << psolvef.roots()[5] << "\t" << psolve6d.roots()[5] << endl;
00051   std::complex<float> castedRoot( psolve6d.roots()[5].real(), psolve6d.roots()[5].imag() );
00052   cout << "Norm of the difference: " << internal::abs( psolvef.roots()[5] - castedRoot ) << endl;
00053 }


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:33:11