diff options
Diffstat (limited to 'eigen/unsupported/Eigen/src/Polynomials')
4 files changed, 61 insertions, 48 deletions
diff --git a/eigen/unsupported/Eigen/src/Polynomials/CMakeLists.txt b/eigen/unsupported/Eigen/src/Polynomials/CMakeLists.txt deleted file mode 100644 index 51f13f3..0000000 --- a/eigen/unsupported/Eigen/src/Polynomials/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Polynomials_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Polynomials_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/Polynomials COMPONENT Devel - ) diff --git a/eigen/unsupported/Eigen/src/Polynomials/Companion.h b/eigen/unsupported/Eigen/src/Polynomials/Companion.h index b515c29..e0af6eb 100644 --- a/eigen/unsupported/Eigen/src/Polynomials/Companion.h +++ b/eigen/unsupported/Eigen/src/Polynomials/Companion.h @@ -75,7 +75,7 @@ class companion void setPolynomial( const VectorType& poly ) { const Index deg = poly.size()-1; - m_monic = -1/poly[deg] * poly.head(deg); + m_monic = Scalar(-1)/poly[deg] * poly.head(deg); //m_bl_diag.setIdentity( deg-1 ); m_bl_diag.setOnes(deg-1); } @@ -107,8 +107,8 @@ class companion * colB and rowB are repectively the multipliers for * the column and the row in order to balance them. * */ - bool balanced( Scalar colNorm, Scalar rowNorm, - bool& isBalanced, Scalar& colB, Scalar& rowB ); + bool balanced( RealScalar colNorm, RealScalar rowNorm, + bool& isBalanced, RealScalar& colB, RealScalar& rowB ); /** Helper function for the balancing algorithm. * \returns true if the row and the column, having colNorm and rowNorm @@ -116,8 +116,8 @@ class companion * colB and rowB are repectively the multipliers for * the column and the row in order to balance them. * */ - bool balancedR( Scalar colNorm, Scalar rowNorm, - bool& isBalanced, Scalar& colB, Scalar& rowB ); + bool balancedR( RealScalar colNorm, RealScalar rowNorm, + bool& isBalanced, RealScalar& colB, RealScalar& rowB ); public: /** @@ -139,10 +139,10 @@ class companion template< typename _Scalar, int _Deg > inline -bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm, - bool& isBalanced, Scalar& colB, Scalar& rowB ) +bool companion<_Scalar,_Deg>::balanced( RealScalar colNorm, RealScalar rowNorm, + bool& isBalanced, RealScalar& colB, RealScalar& rowB ) { - if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; } + if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm ){ return true; } else { //To find the balancing coefficients, if the radix is 2, @@ -150,29 +150,29 @@ bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm, // \f$ 2^{2\sigma-1} < rowNorm / colNorm \le 2^{2\sigma+1} \f$ // then the balancing coefficient for the row is \f$ 1/2^{\sigma} \f$ // and the balancing coefficient for the column is \f$ 2^{\sigma} \f$ - rowB = rowNorm / radix<Scalar>(); - colB = Scalar(1); - const Scalar s = colNorm + rowNorm; + rowB = rowNorm / radix<RealScalar>(); + colB = RealScalar(1); + const RealScalar s = colNorm + rowNorm; while (colNorm < rowB) { - colB *= radix<Scalar>(); - colNorm *= radix2<Scalar>(); + colB *= radix<RealScalar>(); + colNorm *= radix2<RealScalar>(); } - rowB = rowNorm * radix<Scalar>(); + rowB = rowNorm * radix<RealScalar>(); while (colNorm >= rowB) { - colB /= radix<Scalar>(); - colNorm /= radix2<Scalar>(); + colB /= radix<RealScalar>(); + colNorm /= radix2<RealScalar>(); } //This line is used to avoid insubstantial balancing - if ((rowNorm + colNorm) < Scalar(0.95) * s * colB) + if ((rowNorm + colNorm) < RealScalar(0.95) * s * colB) { isBalanced = false; - rowB = Scalar(1) / colB; + rowB = RealScalar(1) / colB; return false; } else{ @@ -182,21 +182,21 @@ bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm, template< typename _Scalar, int _Deg > inline -bool companion<_Scalar,_Deg>::balancedR( Scalar colNorm, Scalar rowNorm, - bool& isBalanced, Scalar& colB, Scalar& rowB ) +bool companion<_Scalar,_Deg>::balancedR( RealScalar colNorm, RealScalar rowNorm, + bool& isBalanced, RealScalar& colB, RealScalar& rowB ) { - if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; } + if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm ){ return true; } else { /** * Set the norm of the column and the row to the geometric mean * of the row and column norm */ - const _Scalar q = colNorm/rowNorm; + const RealScalar q = colNorm/rowNorm; if( !isApprox( q, _Scalar(1) ) ) { rowB = sqrt( colNorm/rowNorm ); - colB = Scalar(1)/rowB; + colB = RealScalar(1)/rowB; isBalanced = false; return false; @@ -219,8 +219,8 @@ void companion<_Scalar,_Deg>::balance() while( !hasConverged ) { hasConverged = true; - Scalar colNorm,rowNorm; - Scalar colB,rowB; + RealScalar colNorm,rowNorm; + RealScalar colB,rowB; //First row, first column excluding the diagonal //============================================== diff --git a/eigen/unsupported/Eigen/src/Polynomials/PolynomialSolver.h b/eigen/unsupported/Eigen/src/Polynomials/PolynomialSolver.h index cd5c04b..7885942 100644 --- a/eigen/unsupported/Eigen/src/Polynomials/PolynomialSolver.h +++ b/eigen/unsupported/Eigen/src/Polynomials/PolynomialSolver.h @@ -41,7 +41,7 @@ class PolynomialSolverBase protected: template< typename OtherPolynomial > inline void setPolynomial( const OtherPolynomial& poly ){ - m_roots.resize(poly.size()); } + m_roots.resize(poly.size()-1); } public: template< typename OtherPolynomial > @@ -99,7 +99,7 @@ class PolynomialSolverBase */ inline const RootType& greatestRoot() const { - std::greater<Scalar> greater; + std::greater<RealScalar> greater; return selectComplexRoot_withRespectToNorm( greater ); } @@ -108,7 +108,7 @@ class PolynomialSolverBase */ inline const RootType& smallestRoot() const { - std::less<Scalar> less; + std::less<RealScalar> less; return selectComplexRoot_withRespectToNorm( less ); } @@ -213,7 +213,7 @@ class PolynomialSolverBase bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const { - std::greater<Scalar> greater; + std::greater<RealScalar> greater; return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold ); } @@ -236,7 +236,7 @@ class PolynomialSolverBase bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const { - std::less<Scalar> less; + std::less<RealScalar> less; return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold ); } @@ -259,7 +259,7 @@ class PolynomialSolverBase bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const { - std::greater<Scalar> greater; + std::greater<RealScalar> greater; return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold ); } @@ -282,7 +282,7 @@ class PolynomialSolverBase bool& hasArealRoot, const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const { - std::less<Scalar> less; + std::less<RealScalar> less; return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold ); } @@ -316,7 +316,7 @@ class PolynomialSolverBase * - real roots with greatest, smallest absolute real value. * - greatest, smallest real roots. * - * WARNING: this polynomial solver is experimental, part of the unsuported Eigen modules. + * WARNING: this polynomial solver is experimental, part of the unsupported Eigen modules. * * * Currently a QR algorithm is used to compute the eigenvalues of the companion matrix of @@ -327,7 +327,7 @@ class PolynomialSolverBase * However, almost always, correct accuracy is reached even in these cases for 64bit * (double) floating types and small polynomial degree (<20). */ -template< typename _Scalar, int _Deg > +template<typename _Scalar, int _Deg> class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> { public: @@ -337,7 +337,9 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType; - typedef EigenSolver<CompanionMatrixType> EigenSolverType; + typedef typename internal::conditional<NumTraits<Scalar>::IsComplex, + ComplexEigenSolver<CompanionMatrixType>, + EigenSolver<CompanionMatrixType> >::type EigenSolverType; public: /** Computes the complex roots of a new polynomial. */ @@ -345,10 +347,19 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> void compute( const OtherPolynomial& poly ) { eigen_assert( Scalar(0) != poly[poly.size()-1] ); - internal::companion<Scalar,_Deg> companion( poly ); - companion.balance(); - m_eigenSolver.compute( companion.denseMatrix() ); - m_roots = m_eigenSolver.eigenvalues(); + eigen_assert( poly.size() > 1 ); + if(poly.size() > 2 ) + { + internal::companion<Scalar,_Deg> companion( poly ); + companion.balance(); + m_eigenSolver.compute( companion.denseMatrix() ); + m_roots = m_eigenSolver.eigenvalues(); + } + else if(poly.size () == 2) + { + m_roots.resize(1); + m_roots[0] = -poly[0]/poly[1]; + } } public: @@ -376,10 +387,18 @@ class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1> template< typename OtherPolynomial > void compute( const OtherPolynomial& poly ) { - eigen_assert( Scalar(0) != poly[poly.size()-1] ); - m_roots[0] = -poly[0]/poly[poly.size()-1]; + eigen_assert( poly.size() == 2 ); + eigen_assert( Scalar(0) != poly[1] ); + m_roots[0] = -poly[0]/poly[1]; } + public: + template< typename OtherPolynomial > + inline PolynomialSolver( const OtherPolynomial& poly ){ + compute( poly ); } + + inline PolynomialSolver(){} + protected: using PS_Base::m_roots; }; diff --git a/eigen/unsupported/Eigen/src/Polynomials/PolynomialUtils.h b/eigen/unsupported/Eigen/src/Polynomials/PolynomialUtils.h index 2bb8bc8..40ba65b 100644 --- a/eigen/unsupported/Eigen/src/Polynomials/PolynomialUtils.h +++ b/eigen/unsupported/Eigen/src/Polynomials/PolynomialUtils.h @@ -56,7 +56,7 @@ T poly_eval( const Polynomials& poly, const T& x ) for( DenseIndex i=1; i<poly.size(); ++i ){ val = val*inv_x + poly[i]; } - return std::pow(x,(T)(poly.size()-1)) * val; + return numext::pow(x,(T)(poly.size()-1)) * val; } } |