Re: [eigen] Unsupported Polynomial solver |
[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]
Sorry for late reply, but thanks for the patches.Regarding the eigen-solver choice issue, it's possible to automatically use the right one based on NumTraits<Scalar>::IsComplex and conditional<>.We would also need to extended the unit test to cover complexes.gaelOn Tue, Nov 8, 2016 at 3:41 PM, Oleg Shirokobrod <oleg.shirokobrod@xxxxxxxxx> wrote:Oleg ShirokobrodBest regards,Matlab for the same coefficients gives following roots (remember that matlab array of coefficients and Eigen one are reversed to each other)Output:TestWith this code I have run test against Matlab and it gave similar results.I have to replace in number of places in the file Companion.h (mainly in functions balanced() and balance()) Scalar with RealScalar, where variables are really real.So I have to add additional template parameter EigenSolverType with default value for real coefficients:Dear list,Polynomial solver currently finds roots only for real polynomial coefficients. Having Eigen::ComplexEigenSolver lets Polynomial module to compute roots for complex polynomial coefficients. I have made necessary modifications. I have attached corresponding patch files. Unfortunately module cannot deduce type of the solver from partial specialization of corresponding matrix type, such that
template<typename RealScalar> class EigenSolver {};
template<typenameRealScalar> class EigenSolver<Eigen::Matrix<RealScalar,_Deg,_Deg> > {// EigenSolver };
template<typename RealScalar> class EigenSolver<Eigen::Matrix<std::complex<RealScalar>,_Deg,_Deg > > { //ComplexEigenSolver};
template< typename _Scalar, int _Deg , typename EigenSolverType = EigenSolver<Matrix<_Scalar,_Deg,_Deg> > >
class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>{};
VectorXcd roots = VectorXcd::Random(4);
VectorXcd polynomialCoefficients;
roots_to_monicPolynomial(roots, polynomialCoefficients);
cout << "roots : " << endl;
cout << setprecision(14) << roots << endl;
cout << "polynomialCoefficients : " << endl;
cout << setprecision(14) << polynomialCoefficients << endl;
PolynomialSolver<std::complex<double>, Dynamic, ComplexEigenSolver<MatrixXcd> > psolve(polynomialCoefficients) ;
VectorXcd computedRoots = psolve.roots();
cout << "computedRoots : " << endl;
cout << setprecision(14) << computedRoots << endl;
for(int index = 0; index < computedRoots.size(); ++index)
{
cout << setprecision(14) << "polynom at computedRoots * 1.0 e-11: " << poly_eval(polynomialCoefficients, computedRoots(index))*1..0e11 << endl;
}
roots :
(0.12717062898648,-0.99749748222297)
(0.61748100222785,-0.61339152195807)
(-0.04025391399884,0.17001861629078)
(0.79192480239265,-0.29941709646901)
polynomialCoefficients :
(0.091649023983597,-0.091441416775918)
(0.24020598856634,0.37401934653925)
(-0.16301627948124,-1.8544616197629)
(-1.4963225196081,1.7402874843593)
(1,0)
computedRoots :
(-0.04025391399884,0.17001861629078)
(0.79192480239265,-0.29941709646901)
(0.61748100222785,-0.61339152195807)
(0.12717062898648,-0.99749748222297)
polynom at computedRoots * 1.0 e-11: (8.3266726846887e-006,-1.2490009027033e-005)
polynom at computedRoots * 1.0 e-11: (4.7184478546569e-005,-1.6653345369377e-005)
polynom at computedRoots * 1.0 e-11: (2.2204460492503e-005,-1.3877787807814e-005)
polynom at computedRoots * 1.0 e-11: (1.5163809063462e-005,-2.7286889655471e-005)
computedRoots :
0.127170628986480 - 0.997497482222969i
0..617481002227849 - 0.613391521958071i
0.791924802392650 - 0.299417096469009i
-0.040253913998840 + 0.170018616290780i
Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |