PolynomialSolver.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_POLYNOMIAL_SOLVER_H
11 #define EIGEN_POLYNOMIAL_SOLVER_H
12 
13 namespace Eigen {
14 
28 template< typename _Scalar, int _Deg >
30 {
31  public:
33 
34  typedef _Scalar Scalar;
35  typedef typename NumTraits<Scalar>::Real RealScalar;
36  typedef std::complex<RealScalar> RootType;
37  typedef Matrix<RootType,_Deg,1> RootsType;
38 
39  typedef DenseIndex Index;
40 
41  protected:
42  template< typename OtherPolynomial >
43  inline void setPolynomial( const OtherPolynomial& poly ){
44  m_roots.resize(poly.size()); }
45 
46  public:
47  template< typename OtherPolynomial >
48  inline PolynomialSolverBase( const OtherPolynomial& poly ){
49  setPolynomial( poly() ); }
50 
52 
53  public:
55  inline const RootsType& roots() const { return m_roots; }
56 
57  public:
68  template<typename Stl_back_insertion_sequence>
69  inline void realRoots( Stl_back_insertion_sequence& bi_seq,
70  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
71  {
72  using std::abs;
73  bi_seq.clear();
74  for(Index i=0; i<m_roots.size(); ++i )
75  {
76  if( abs( m_roots[i].imag() ) < absImaginaryThreshold ){
77  bi_seq.push_back( m_roots[i].real() ); }
78  }
79  }
80 
81  protected:
82  template<typename squaredNormBinaryPredicate>
83  inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred ) const
84  {
85  Index res=0;
86  RealScalar norm2 = numext::abs2( m_roots[0] );
87  for( Index i=1; i<m_roots.size(); ++i )
88  {
89  const RealScalar currNorm2 = numext::abs2( m_roots[i] );
90  if( pred( currNorm2, norm2 ) ){
91  res=i; norm2=currNorm2; }
92  }
93  return m_roots[res];
94  }
95 
96  public:
100  inline const RootType& greatestRoot() const
101  {
102  std::greater<Scalar> greater;
103  return selectComplexRoot_withRespectToNorm( greater );
104  }
105 
109  inline const RootType& smallestRoot() const
110  {
111  std::less<Scalar> less;
113  }
114 
115  protected:
116  template<typename squaredRealPartBinaryPredicate>
118  squaredRealPartBinaryPredicate& pred,
119  bool& hasArealRoot,
120  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
121  {
122  using std::abs;
123  hasArealRoot = false;
124  Index res=0;
125  RealScalar abs2(0);
126 
127  for( Index i=0; i<m_roots.size(); ++i )
128  {
129  if( abs( m_roots[i].imag() ) < absImaginaryThreshold )
130  {
131  if( !hasArealRoot )
132  {
133  hasArealRoot = true;
134  res = i;
135  abs2 = m_roots[i].real() * m_roots[i].real();
136  }
137  else
138  {
139  const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
140  if( pred( currAbs2, abs2 ) )
141  {
142  abs2 = currAbs2;
143  res = i;
144  }
145  }
146  }
147  else
148  {
149  if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){
150  res = i; }
151  }
152  }
153  return numext::real_ref(m_roots[res]);
154  }
155 
156 
157  template<typename RealPartBinaryPredicate>
159  RealPartBinaryPredicate& pred,
160  bool& hasArealRoot,
161  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
162  {
163  using std::abs;
164  hasArealRoot = false;
165  Index res=0;
166  RealScalar val(0);
167 
168  for( Index i=0; i<m_roots.size(); ++i )
169  {
170  if( abs( m_roots[i].imag() ) < absImaginaryThreshold )
171  {
172  if( !hasArealRoot )
173  {
174  hasArealRoot = true;
175  res = i;
176  val = m_roots[i].real();
177  }
178  else
179  {
180  const RealScalar curr = m_roots[i].real();
181  if( pred( curr, val ) )
182  {
183  val = curr;
184  res = i;
185  }
186  }
187  }
188  else
189  {
190  if( abs( m_roots[i].imag() ) < abs( m_roots[res].imag() ) ){
191  res = i; }
192  }
193  }
194  return numext::real_ref(m_roots[res]);
195  }
196 
197  public:
213  bool& hasArealRoot,
214  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
215  {
216  std::greater<Scalar> greater;
217  return selectRealRoot_withRespectToAbsRealPart( greater, hasArealRoot, absImaginaryThreshold );
218  }
219 
220 
236  bool& hasArealRoot,
237  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
238  {
239  std::less<Scalar> less;
240  return selectRealRoot_withRespectToAbsRealPart( less, hasArealRoot, absImaginaryThreshold );
241  }
242 
243 
259  bool& hasArealRoot,
260  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
261  {
262  std::greater<Scalar> greater;
263  return selectRealRoot_withRespectToRealPart( greater, hasArealRoot, absImaginaryThreshold );
264  }
265 
266 
282  bool& hasArealRoot,
283  const RealScalar& absImaginaryThreshold = NumTraits<Scalar>::dummy_precision() ) const
284  {
285  std::less<Scalar> less;
286  return selectRealRoot_withRespectToRealPart( less, hasArealRoot, absImaginaryThreshold );
287  }
288 
289  protected:
291 };
292 
293 #define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( BASE ) \
294  typedef typename BASE::Scalar Scalar; \
295  typedef typename BASE::RealScalar RealScalar; \
296  typedef typename BASE::RootType RootType; \
297  typedef typename BASE::RootsType RootsType;
298 
299 
300 
330 template< typename _Scalar, int _Deg >
331 class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
332 {
333  public:
335 
336  typedef PolynomialSolverBase<_Scalar,_Deg> PS_Base;
338 
339  typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType;
340  typedef EigenSolver<CompanionMatrixType> EigenSolverType;
341 
342  public:
344  template< typename OtherPolynomial >
345  void compute( const OtherPolynomial& poly )
346  {
347  eigen_assert( Scalar(0) != poly[poly.size()-1] );
348  internal::companion<Scalar,_Deg> companion( poly );
349  companion.balance();
350  m_eigenSolver.compute( companion.denseMatrix() );
351  m_roots = m_eigenSolver.eigenvalues();
352  }
353 
354  public:
355  template< typename OtherPolynomial >
356  inline PolynomialSolver( const OtherPolynomial& poly ){
357  compute( poly ); }
358 
359  inline PolynomialSolver(){}
360 
361  protected:
362  using PS_Base::m_roots;
364 };
365 
366 
367 template< typename _Scalar >
368 class PolynomialSolver<_Scalar,1> : public PolynomialSolverBase<_Scalar,1>
369 {
370  public:
373 
374  public:
376  template< typename OtherPolynomial >
377  void compute( const OtherPolynomial& poly )
378  {
379  eigen_assert( Scalar(0) != poly[poly.size()-1] );
380  m_roots[0] = -poly[0]/poly[poly.size()-1];
381  }
382 
383  protected:
384  using PS_Base::m_roots;
385 };
386 
387 } // end namespace Eigen
388 
389 #endif // EIGEN_POLYNOMIAL_SOLVER_H
#define inline
Definition: Types.h:109
const RealScalar & greatestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
void realRoots(Stl_back_insertion_sequence &bi_seq, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
PolynomialSolverBase< _Scalar, 1 > PS_Base
DerType::Scalar imag(const AutoDiffScalar< DerType > &)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
Defined to be inherited by polynomial solvers: it provides convenient methods such as...
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs2_op< Scalar >, const Derived > abs2() const
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
const RootType & smallestRoot() const
const RealScalar & smallestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
const RealScalar & absSmallestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
EIGEN_STRONG_INLINE void resize(Index nbRows, Index nbCols)
#define EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES(BASE)
PolynomialSolver(const OtherPolynomial &poly)
const RootType & greatestRoot() const
EIGEN_DEFAULT_DENSE_INDEX_TYPE DenseIndex
Definition: XprHelper.h:27
void setPolynomial(const OtherPolynomial &poly)
const RealScalar & absGreatestRealRoot(bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
PolynomialSolverBase(const OtherPolynomial &poly)
const RootsType & roots() const
internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) >::type real_ref(const Scalar &x)
EigenSolverType m_eigenSolver
Computes eigenvalues and eigenvectors of general matrices.
Definition: EigenSolver.h:64
#define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(Scalar, Size)
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:127
const RealScalar & selectRealRoot_withRespectToAbsRealPart(squaredRealPartBinaryPredicate &pred, bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
#define eigen_assert(x)
const AutoDiffScalar< DerType > & real(const AutoDiffScalar< DerType > &x)
A polynomial solver.
const RealScalar & selectRealRoot_withRespectToRealPart(RealPartBinaryPredicate &pred, bool &hasArealRoot, const RealScalar &absImaginaryThreshold=NumTraits< Scalar >::dummy_precision()) const
const RootType & selectComplexRoot_withRespectToNorm(squaredNormBinaryPredicate &pred) const


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:59