src/Polynomials/Companion.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2010 Manuel Yguel <manuel.yguel@gmail.com>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_COMPANION_H
0011 #define EIGEN_COMPANION_H
0012
0013 // This file requires the user to include
0014 // * Eigen/Core
0015 // * Eigen/src/PolynomialSolver.h
0016
0017 namespace Eigen {
0018
0019 namespace internal {
0020
0021 #ifndef EIGEN_PARSED_BY_DOXYGEN
0022
0023 template<int Size>
0024 struct decrement_if_fixed_size
0025 {
0026   enum {
0027     ret = (Size == Dynamic) ? Dynamic : Size-1 };
0028 };
0029
0030 #endif
0031
0032 template< typename _Scalar, int _Deg >
0033 class companion
0034 {
0035   public:
0036     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
0037
0038     enum {
0039       Deg = _Deg,
0040       Deg_1=decrement_if_fixed_size<Deg>::ret
0041     };
0042
0043     typedef _Scalar                                Scalar;
0044     typedef typename NumTraits<Scalar>::Real       RealScalar;
0045     typedef Matrix<Scalar, Deg, 1>                 RightColumn;
0046     //typedef DiagonalMatrix< Scalar, Deg_1, Deg_1 > BottomLeftDiagonal;
0047     typedef Matrix<Scalar, Deg_1, 1>               BottomLeftDiagonal;
0048
0049     typedef Matrix<Scalar, Deg, Deg>               DenseCompanionMatrixType;
0050     typedef Matrix< Scalar, _Deg, Deg_1 >          LeftBlock;
0051     typedef Matrix< Scalar, Deg_1, Deg_1 >         BottomLeftBlock;
0052     typedef Matrix< Scalar, 1, Deg_1 >             LeftBlockFirstRow;
0053
0054     typedef DenseIndex Index;
0055
0056   public:
0057     EIGEN_STRONG_INLINE const _Scalar operator()(Index row, Index col ) const
0058     {
0059       if( m_bl_diag.rows() > col )
0060       {
0061         if( 0 < row ){ return m_bl_diag[col]; }
0062         else{ return 0; }
0063       }
0064       else{ return m_monic[row]; }
0065     }
0066
0067   public:
0068     template<typename VectorType>
0069     void setPolynomial( const VectorType& poly )
0070     {
0071       const Index deg = poly.size()-1;
0072       m_monic = -poly.head(deg)/poly[deg];
0073       m_bl_diag.setOnes(deg-1);
0074     }
0075
0076     template<typename VectorType>
0077     companion( const VectorType& poly ){
0078       setPolynomial( poly ); }
0079
0080   public:
0081     DenseCompanionMatrixType denseMatrix() const
0082     {
0083       const Index deg   = m_monic.size();
0084       const Index deg_1 = deg-1;
0085       DenseCompanionMatrixType companMat(deg,deg);
0086       companMat <<
0087         ( LeftBlock(deg,deg_1)
0088           << LeftBlockFirstRow::Zero(1,deg_1),
0089           BottomLeftBlock::Identity(deg-1,deg-1)*m_bl_diag.asDiagonal() ).finished()
0090         , m_monic;
0091       return companMat;
0092     }
0093
0094
0095
0096   protected:
0097     /** Helper function for the balancing algorithm.
0098      * \returns true if the row and the column, having colNorm and rowNorm
0099      * as norms, are balanced, false otherwise.
0100      * colB and rowB are respectively the multipliers for
0101      * the column and the row in order to balance them.
0102      * */
0103     bool balanced( RealScalar colNorm, RealScalar rowNorm,
0104         bool& isBalanced, RealScalar& colB, RealScalar& rowB );
0105
0106     /** Helper function for the balancing algorithm.
0107      * \returns true if the row and the column, having colNorm and rowNorm
0108      * as norms, are balanced, false otherwise.
0109      * colB and rowB are respectively the multipliers for
0110      * the column and the row in order to balance them.
0111      * */
0112     bool balancedR( RealScalar colNorm, RealScalar rowNorm,
0113         bool& isBalanced, RealScalar& colB, RealScalar& rowB );
0114
0115   public:
0116     /**
0117      * Balancing algorithm from B. N. PARLETT and C. REINSCH (1969)
0118      * "Balancing a matrix for calculation of eigenvalues and eigenvectors"
0119      * adapted to the case of companion matrices.
0120      * A matrix with non zero row and non zero column is balanced
0121      * for a certain norm if the i-th row and the i-th column
0122      * have same norm for all i.
0123      */
0124     void balance();
0125
0126   protected:
0127       RightColumn                m_monic;
0128       BottomLeftDiagonal         m_bl_diag;
0129 };
0130
0131
0132
0133 template< typename _Scalar, int _Deg >
0134 inline
0135 bool companion<_Scalar,_Deg>::balanced( RealScalar colNorm, RealScalar rowNorm,
0136     bool& isBalanced, RealScalar& colB, RealScalar& rowB )
0137 {
0138   if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm
0139       || !(numext::isfinite)(colNorm) || !(numext::isfinite)(rowNorm)){
0140     return true;
0141   }
0142   else
0143   {
0144     //To find the balancing coefficients, if the radix is 2,
0145     //one finds \f$ \sigma \f$ such that
0146     // \f$ 2^{2\sigma-1} < rowNorm / colNorm \le 2^{2\sigma+1} \f$
0147     // then the balancing coefficient for the row is \f$ 1/2^{\sigma} \f$
0148     // and the balancing coefficient for the column is \f$ 2^{\sigma} \f$
0149     const RealScalar radix = RealScalar(2);
0150     const RealScalar radix2 = RealScalar(4);
0151
0152     rowB = rowNorm / radix;
0153     colB = RealScalar(1);
0154     const RealScalar s = colNorm + rowNorm;
0155
0156     // Find sigma s.t. rowNorm / 2 <= 2^(2*sigma) * colNorm
0157     RealScalar scout = colNorm;
0158     while (scout < rowB)
0159     {
0160       colB *= radix;
0161       scout *= radix2;
0162     }
0163
0164     // We now have an upper-bound for sigma, try to lower it.
0165     // Find sigma s.t. 2^(2*sigma) * colNorm / 2 < rowNorm
0166     scout = colNorm * (colB / radix) * colB;  // Avoid overflow.
0167     while (scout >= rowNorm)
0168     {
0169       colB /= radix;
0170       scout /= radix2;
0171     }
0172
0173     // This line is used to avoid insubstantial balancing.
0174     if ((rowNorm + radix * scout) < RealScalar(0.95) * s * colB)
0175     {
0176       isBalanced = false;
0177       rowB = RealScalar(1) / colB;
0178       return false;
0179     }
0180     else
0181     {
0182       return true;
0183     }
0184   }
0185 }
0186
0187 template< typename _Scalar, int _Deg >
0188 inline
0189 bool companion<_Scalar,_Deg>::balancedR( RealScalar colNorm, RealScalar rowNorm,
0190     bool& isBalanced, RealScalar& colB, RealScalar& rowB )
0191 {
0192   if( RealScalar(0) == colNorm || RealScalar(0) == rowNorm ){ return true; }
0193   else
0194   {
0195     /**
0196      * Set the norm of the column and the row to the geometric mean
0197      * of the row and column norm
0198      */
0199     const RealScalar q = colNorm/rowNorm;
0200     if( !isApprox( q, _Scalar(1) ) )
0201     {
0202       rowB = sqrt( colNorm/rowNorm );
0203       colB = RealScalar(1)/rowB;
0204
0205       isBalanced = false;
0206       return false;
0207     }
0208     else{
0209       return true; }
0210   }
0211 }
0212
0213
0214 template< typename _Scalar, int _Deg >
0215 void companion<_Scalar,_Deg>::balance()
0216 {
0217   using std::abs;
0218   EIGEN_STATIC_ASSERT( Deg == Dynamic || 1 < Deg, YOU_MADE_A_PROGRAMMING_MISTAKE );
0219   const Index deg   = m_monic.size();
0220   const Index deg_1 = deg-1;
0221
0222   bool hasConverged=false;
0223   while( !hasConverged )
0224   {
0225     hasConverged = true;
0226     RealScalar colNorm,rowNorm;
0227     RealScalar colB,rowB;
0228
0229     //First row, first column excluding the diagonal
0230     //==============================================
0231     colNorm = abs(m_bl_diag[0]);
0232     rowNorm = abs(m_monic[0]);
0233
0234     //Compute balancing of the row and the column
0235     if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
0236     {
0237       m_bl_diag[0] *= colB;
0238       m_monic[0] *= rowB;
0239     }
0240
0241     //Middle rows and columns excluding the diagonal
0242     //==============================================
0243     for( Index i=1; i<deg_1; ++i )
0244     {
0245       // column norm, excluding the diagonal
0246       colNorm = abs(m_bl_diag[i]);
0247
0248       // row norm, excluding the diagonal
0249       rowNorm = abs(m_bl_diag[i-1]) + abs(m_monic[i]);
0250
0251       //Compute balancing of the row and the column
0252       if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
0253       {
0254         m_bl_diag[i]   *= colB;
0255         m_bl_diag[i-1] *= rowB;
0256         m_monic[i]     *= rowB;
0257       }
0258     }
0259
0260     //Last row, last column excluding the diagonal
0261     //============================================
0262     const Index ebl = m_bl_diag.size()-1;
0263     VectorBlock<RightColumn,Deg_1> headMonic( m_monic, 0, deg_1 );
0264     colNorm = headMonic.array().abs().sum();
0265     rowNorm = abs( m_bl_diag[ebl] );
0266
0267     //Compute balancing of the row and the column
0268     if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
0269     {
0270       headMonic      *= colB;
0271       m_bl_diag[ebl] *= rowB;
0272     }
0273   }
0274 }
0275
0276 } // end namespace internal
0277
0278 } // end namespace Eigen
0279
0280 #endif // EIGEN_COMPANION_H