$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_POLYNOMIAL_SOLVER_H 00011 #define EIGEN_POLYNOMIAL_SOLVER_H 00012 00013 namespace Eigen { 00014 00028 template< typename _Scalar, int _Deg > 00029 class PolynomialSolverBase 00030 { 00031 public: 00032 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg) 00033 00034 typedef _Scalar Scalar; 00035 typedef typename NumTraits<Scalar>::Real RealScalar; 00036 typedef std::complex<RealScalar> RootType; 00037 typedef Matrix<RootType,_Deg,1> RootsType; 00038 00039 typedef DenseIndex Index; 00040 00041 protected: 00042 template< typename OtherPolynomial > 00043 inline void setPolynomial( const OtherPolynomial& poly ){ 00044 m_roots.resize(poly.size()); } 00045 00046 public: 00047 template< typename OtherPolynomial > 00048 inline PolynomialSolverBase( const OtherPolynomial& poly ){ 00049 setPolynomial( poly() ); } 00050 00051 inline PolynomialSolverBase(){} 00052 00053 public: 00055 inline const RootsType& roots() const { return m_roots; } 00056 00057 public: 00068 template<typename Stl_back_insertion_sequence> 00069 inline void realRoots( Stl_back_insertion_sequence& bi_seq, 00070 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 00071 { 00072 using std::abs; 00073 bi_seq.clear(); 00074 for(Index i=0; i<m_roots.size(); ++i ) 00075 { 00076 if( abs( m_roots[i].imag() ) < absImaginaryThreshold ){ 00077 bi_seq.push_back( m_roots[i].real() ); } 00078 } 00079 } 00080 00081 protected: 00082 template<typename squaredNormBinaryPredicate> 00083 inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const 00084 { 00085 Index res=0; 00086 RealScalar norm2 = numext::abs2( m_roots[0] ); 00087 for( Index i=1; i<m_roots.size(); ++i ) 00088 { 00089 const RealScalar currNorm2 = numext::abs2( m_roots[i] ); 00090 if( pred( currNorm2, norm2 ) ){ 00091 res=i; norm2=currNorm2; } 00092 } 00093 return m_roots[res]; 00094 } 00095 00096 public: 00100 inline const RootType& greatestRoot() const 00101 { 00102 std::greater<Scalar> greater; 00103 return selectComplexRoot_withRespectToNorm( greater ); 00104 } 00105 00109 inline const RootType& smallestRoot() const 00110 { 00111 std::less<Scalar> less; 00112 return selectComplexRoot_withRespectToNorm( less ); 00113 } 00114 00115 protected: 00116 template<typename squaredRealPartBinaryPredicate> 00117 inline const RealScalar& selectRealRoot_withRespectToAbsRealPart( 00118 squaredRealPartBinaryPredicate& pred, 00119 bool& hasArealRoot, 00120 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 00121 { 00122 using std::abs; 00123 hasArealRoot = false; 00124 Index res=0; 00125 RealScalar abs2(0); 00126 00127 for( Index i=0; i<m_roots.size(); ++i ) 00128 { 00129 if( abs( m_roots[i].imag() ) < absImaginaryThreshold ) 00130 { 00131 if( !hasArealRoot ) 00132 { 00133 hasArealRoot = true; 00134 res = i; 00135 abs2 = m_roots[i].real() * m_roots[i].real(); 00136 } 00137 else 00138 { 00139 const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real(); 00140 if( pred( currAbs2, abs2 ) ) 00141 { 00142 abs2 = currAbs2; 00143 res = i; 00144 } 00145 } 00146 } 00147 else 00148 { 00149 if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){ 00150 res = i; } 00151 } 00152 } 00153 return numext::real_ref(m_roots[res]); 00154 } 00155 00156 00157 template<typename RealPartBinaryPredicate> 00158 inline const RealScalar& selectRealRoot_withRespectToRealPart( 00159 RealPartBinaryPredicate& pred, 00160 bool& hasArealRoot, 00161 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 00162 { 00163 using std::abs; 00164 hasArealRoot = false; 00165 Index res=0; 00166 RealScalar val(0); 00167 00168 for( Index i=0; i<m_roots.size(); ++i ) 00169 { 00170 if( abs( m_roots[i].imag() ) < absImaginaryThreshold ) 00171 { 00172 if( !hasArealRoot ) 00173 { 00174 hasArealRoot = true; 00175 res = i; 00176 val = m_roots[i].real(); 00177 } 00178 else 00179 { 00180 const RealScalar curr = m_roots[i].real(); 00181 if( pred( curr, val ) ) 00182 { 00183 val = curr; 00184 res = i; 00185 } 00186 } 00187 } 00188 else 00189 { 00190 if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){ 00191 res = i; } 00192 } 00193 } 00194 return numext::real_ref(m_roots[res]); 00195 } 00196 00197 public: 00212 inline const RealScalar& absGreatestRealRoot( 00213 bool& hasArealRoot, 00214 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 00215 { 00216 std::greater<Scalar> greater; 00217 return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold ); 00218 } 00219 00220 00235 inline const RealScalar& absSmallestRealRoot( 00236 bool& hasArealRoot, 00237 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 00238 { 00239 std::less<Scalar> less; 00240 return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold ); 00241 } 00242 00243 00258 inline const RealScalar& greatestRealRoot( 00259 bool& hasArealRoot, 00260 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 00261 { 00262 std::greater<Scalar> greater; 00263 return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold ); 00264 } 00265 00266 00281 inline const RealScalar& smallestRealRoot( 00282 bool& hasArealRoot, 00283 const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const 00284 { 00285 std::less<Scalar> less; 00286 return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold ); 00287 } 00288 00289 protected: 00290 RootsType m_roots; 00291 }; 00292 00293 #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE ) \ 00294 typedef typename BASE::Scalar Scalar; \ 00295 typedef typename BASE::RealScalar RealScalar; \ 00296 typedef typename BASE::RootType RootType; \ 00297 typedef typename BASE::RootsType RootsType; 00298 00299 00300 00330 template< typename _Scalar, int _Deg > 00331 class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> 00332 { 00333 public: 00334 EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg) 00335 00336 typedef PolynomialSolverBase<_Scalar,_Deg> PS_Base; 00337 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) 00338 00339 typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType; 00340 typedef EigenSolver<CompanionMatrixType> EigenSolverType; 00341 00342 public: 00344 template< typename OtherPolynomial > 00345 void compute( const OtherPolynomial& poly ) 00346 { 00347 eigen_assert( Scalar(0) != poly[poly.size()-1] ); 00348 internal::companion<Scalar,_Deg> companion( poly ); 00349 companion.balance(); 00350 m_eigenSolver.compute( companion.denseMatrix() ); 00351 m_roots = m_eigenSolver.eigenvalues(); 00352 } 00353 00354 public: 00355 template< typename OtherPolynomial > 00356 inline PolynomialSolver( const OtherPolynomial& poly ){ 00357 compute( poly ); } 00358 00359 inline PolynomialSolver(){} 00360 00361 protected: 00362 using PS_Base::m_roots; 00363 EigenSolverType m_eigenSolver; 00364 }; 00365 00366 00367 template< typename _Scalar > 00368 class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1> 00369 { 00370 public: 00371 typedef PolynomialSolverBase<_Scalar,1> PS_Base; 00372 EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) 00373 00374 public: 00376 template< typename OtherPolynomial > 00377 void compute( const OtherPolynomial& poly ) 00378 { 00379 eigen_assert( Scalar(0) != poly[poly.size()-1] ); 00380 m_roots[0] = -poly[0]/poly[poly.size()-1]; 00381 } 00382 00383 protected: 00384 using PS_Base::m_roots; 00385 }; 00386 00387 } // end namespace Eigen 00388 00389 #endif // EIGEN_POLYNOMIAL_SOLVER_H