GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-05-18 23:57:39 +08:00
Code Issues Packages Projects Releases Wiki Activity

Patch from Oleg Shirokobrod to extend polynomial solver to complexes

Browse Source
This commit is contained in:
Gael Guennebaud 2016-11-23 15:42:26 +01:00
parent e340866c81
commit 9246587122
2 changed files with 27 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
unsupported/Eigen/src/Polynomials/Companion.h
View File

@ -75,7 +75,7 @@ class companion
void setPolynomial( const VectorType& poly ) void setPolynomial( const VectorType& poly )
{ {
const Index deg = poly.size()-1; 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.setIdentity( deg-1 );
m_bl_diag.setOnes(deg-1); m_bl_diag.setOnes(deg-1);
} }
@ -107,8 +107,8 @@ class companion
* colB and rowB are repectively the multipliers for * colB and rowB are repectively the multipliers for
* the column and the row in order to balance them. * the column and the row in order to balance them.
* */ * */
bool balanced( Scalar colNorm, Scalar rowNorm, bool balanced( RealScalar colNorm, RealScalar rowNorm,
bool& isBalanced, Scalar& colB, Scalar& rowB ); bool& isBalanced, RealScalar& colB, RealScalar& rowB );
/** Helper function for the balancing algorithm. /** Helper function for the balancing algorithm.
* \returns true if the row and the column, having colNorm and rowNorm * \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 * colB and rowB are repectively the multipliers for
* the column and the row in order to balance them. * the column and the row in order to balance them.
* */ * */
bool balancedR( Scalar colNorm, Scalar rowNorm, bool balancedR( RealScalar colNorm, RealScalar rowNorm,
bool& isBalanced, Scalar& colB, Scalar& rowB ); bool& isBalanced, RealScalar& colB, RealScalar& rowB );
public: public:
/** /**
@ -139,10 +139,10 @@ class companion
template< typename _Scalar, int _Deg > template< typename _Scalar, int _Deg >
inline inline
bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm, bool companion<_Scalar,_Deg>::balanced( RealScalar colNorm, RealScalar rowNorm,
bool& isBalanced, Scalar& colB, Scalar& rowB ) 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 else
{ {
//To find the balancing coefficients, if the radix is 2, //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$ // \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$ // 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$ // and the balancing coefficient for the column is \f$ 2^{\sigma} \f$
rowB = rowNorm / radix<Scalar>(); rowB = rowNorm / radix<RealScalar>();
colB = Scalar(1); colB = RealScalar(1);
const Scalar s = colNorm + rowNorm; const RealScalar s = colNorm + rowNorm;
while (colNorm < rowB) while (colNorm < rowB)
{ {
colB *= radix<Scalar>(); colB *= radix<RealScalar>();
colNorm *= radix2<Scalar>(); colNorm *= radix2<RealScalar>();
} }
rowB = rowNorm * radix<Scalar>(); rowB = rowNorm * radix<RealScalar>();
while (colNorm >= rowB) while (colNorm >= rowB)
{ {
colB /= radix<Scalar>(); colB /= radix<RealScalar>();
colNorm /= radix2<Scalar>(); colNorm /= radix2<RealScalar>();
} }
//This line is used to avoid insubstantial balancing //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; isBalanced = false;
rowB = Scalar(1) / colB; rowB = RealScalar(1) / colB;
return false; return false;
} }
else{ else{
@ -182,21 +182,21 @@ bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm,
template< typename _Scalar, int _Deg > template< typename _Scalar, int _Deg >
inline inline
bool companion<_Scalar,_Deg>::balancedR( Scalar colNorm, Scalar rowNorm, bool companion<_Scalar,_Deg>::balancedR( RealScalar colNorm, RealScalar rowNorm,
bool& isBalanced, Scalar& colB, Scalar& rowB ) 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 else
{ {
/** /**
* Set the norm of the column and the row to the geometric mean * Set the norm of the column and the row to the geometric mean
* of the row and column norm * of the row and column norm
*/ */
const _Scalar q = colNorm/rowNorm; const RealScalar q = colNorm/rowNorm;
if( !isApprox( q, _Scalar(1) ) ) if( !isApprox( q, _Scalar(1) ) )
{ {
rowB = sqrt( colNorm/rowNorm ); rowB = sqrt( colNorm/rowNorm );
colB = Scalar(1)/rowB; colB = RealScalar(1)/rowB;
isBalanced = false; isBalanced = false;
return false; return false;
@ -219,8 +219,8 @@ void companion<_Scalar,_Deg>::balance()
while( !hasConverged ) while( !hasConverged )
{ {
hasConverged = true; hasConverged = true;
Scalar colNorm,rowNorm; RealScalar colNorm,rowNorm;
Scalar colB,rowB; RealScalar colB,rowB;
//First row, first column excluding the diagonal //First row, first column excluding the diagonal
//============================================== //==============================================

4
unsupported/Eigen/src/Polynomials/PolynomialSolver.h
View File

@ -327,7 +327,7 @@ class PolynomialSolverBase
* However, almost always, correct accuracy is reached even in these cases for 64bit * However, almost always, correct accuracy is reached even in these cases for 64bit
* (double) floating types and small polynomial degree (<20). * (double) floating types and small polynomial degree (<20).
*/ */
template< typename _Scalar, int _Deg > template< typename _Scalar, int _Deg , typename EigenSolverType = EigenSolver<Matrix<_Scalar,_Deg,_Deg> > >
class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg> class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
{ {
public: public:
@ -337,7 +337,7 @@ class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>
EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base ) EIGEN_POLYNOMIAL_SOLVER_BASE_INHERITED_TYPES( PS_Base )
typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType; typedef Matrix<Scalar,_Deg,_Deg> CompanionMatrixType;
typedef EigenSolver<CompanionMatrixType> EigenSolverType; //typedef EigenSolver<CompanionMatrixType> EigenSolverType;
public: public:
/** Computes the complex roots of a new polynomial. */ /** Computes the complex roots of a new polynomial. */