Companion.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_COMPANION_H
00026 #define EIGEN_COMPANION_H
00027 
00028 // This file requires the user to include
00029 // * Eigen/Core
00030 // * Eigen/src/PolynomialSolver.h
00031 
00032 #ifndef EIGEN_PARSED_BY_DOXYGEN
00033 
00034 namespace internal {
00035 
00036 template <typename T>
00037 T radix(){ return 2; }
00038 
00039 template <typename T>
00040 T radix2(){ return radix<T>()*radix<T>(); }
00041 
00042 template<int Size>
00043 struct decrement_if_fixed_size
00044 {
00045   enum {
00046     ret = (Size == Dynamic) ? Dynamic : Size-1 };
00047 };
00048 
00049 #endif
00050 
00051 template< typename _Scalar, int _Deg >
00052 class companion
00053 {
00054   public:
00055     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
00056 
00057     enum {
00058       Deg = _Deg,
00059       Deg_1=decrement_if_fixed_size<Deg>::ret
00060     };
00061 
00062     typedef _Scalar                                Scalar;
00063     typedef typename NumTraits<Scalar>::Real       RealScalar;
00064     typedef Matrix<Scalar, Deg, 1>                 RightColumn;
00065     //typedef DiagonalMatrix< Scalar, Deg_1, Deg_1 > BottomLeftDiagonal;
00066     typedef Matrix<Scalar, Deg_1, 1>               BottomLeftDiagonal;
00067 
00068     typedef Matrix<Scalar, Deg, Deg>               DenseCompanionMatrixType;
00069     typedef Matrix< Scalar, _Deg, Deg_1 >          LeftBlock;
00070     typedef Matrix< Scalar, Deg_1, Deg_1 >         BottomLeftBlock;
00071     typedef Matrix< Scalar, 1, Deg_1 >             LeftBlockFirstRow;
00072 
00073     typedef DenseIndex Index;
00074 
00075   public:
00076     EIGEN_STRONG_INLINE const _Scalar operator()(Index row, Index col ) const
00077     {
00078       if( m_bl_diag.rows() > col )
00079       {
00080         if( 0 < row ){ return m_bl_diag[col]; }
00081         else{ return 0; }
00082       }
00083       else{ return m_monic[row]; }
00084     }
00085 
00086   public:
00087     template<typename VectorType>
00088     void setPolynomial( const VectorType& poly )
00089     {
00090       const Index deg = poly.size()-1;
00091       m_monic = -1/poly[deg] * poly.head(deg);
00092       //m_bl_diag.setIdentity( deg-1 );
00093       m_bl_diag.setOnes(deg-1);
00094     }
00095 
00096     template<typename VectorType>
00097     companion( const VectorType& poly ){
00098       setPolynomial( poly ); }
00099 
00100   public:
00101     DenseCompanionMatrixType denseMatrix() const
00102     {
00103       const Index deg   = m_monic.size();
00104       const Index deg_1 = deg-1;
00105       DenseCompanionMatrixType companion(deg,deg);
00106       companion <<
00107         ( LeftBlock(deg,deg_1)
00108           << LeftBlockFirstRow::Zero(1,deg_1),
00109           BottomLeftBlock::Identity(deg-1,deg-1)*m_bl_diag.asDiagonal() ).finished()
00110         , m_monic;
00111       return companion;
00112     }
00113 
00114 
00115 
00116   protected:
00123     bool balanced( Scalar colNorm, Scalar rowNorm,
00124         bool& isBalanced, Scalar& colB, Scalar& rowB );
00125 
00132     bool balancedR( Scalar colNorm, Scalar rowNorm,
00133         bool& isBalanced, Scalar& colB, Scalar& rowB );
00134 
00135   public:
00144     void balance();
00145 
00146   protected:
00147       RightColumn                m_monic;
00148       BottomLeftDiagonal         m_bl_diag;
00149 };
00150 
00151 
00152 
00153 template< typename _Scalar, int _Deg >
00154 inline
00155 bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm,
00156     bool& isBalanced, Scalar& colB, Scalar& rowB )
00157 {
00158   if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; }
00159   else
00160   {
00161     //To find the balancing coefficients, if the radix is 2,
00162     //one finds \f$ \sigma \f$ such that
00163     // \f$ 2^{2\sigma-1} < rowNorm / colNorm \le 2^{2\sigma+1} \f$
00164     // then the balancing coefficient for the row is \f$ 1/2^{\sigma} \f$
00165     // and the balancing coefficient for the column is \f$ 2^{\sigma} \f$
00166     rowB = rowNorm / radix<Scalar>();
00167     colB = Scalar(1);
00168     const Scalar s = colNorm + rowNorm;
00169 
00170     while (colNorm < rowB)
00171     {
00172       colB *= radix<Scalar>();
00173       colNorm *= radix2<Scalar>();
00174     }
00175 
00176     rowB = rowNorm * radix<Scalar>();
00177 
00178     while (colNorm >= rowB)
00179     {
00180       colB /= radix<Scalar>();
00181       colNorm /= radix2<Scalar>();
00182     }
00183 
00184     //This line is used to avoid insubstantial balancing
00185     if ((rowNorm + colNorm) < Scalar(0.95) * s * colB)
00186     {
00187       isBalanced = false;
00188       rowB = Scalar(1) / colB;
00189       return false;
00190     }
00191     else{
00192       return true; }
00193   }
00194 }
00195 
00196 template< typename _Scalar, int _Deg >
00197 inline
00198 bool companion<_Scalar,_Deg>::balancedR( Scalar colNorm, Scalar rowNorm,
00199     bool& isBalanced, Scalar& colB, Scalar& rowB )
00200 {
00201   if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; }
00202   else
00203   {
00208     const _Scalar q = colNorm/rowNorm;
00209     if( !isApprox( q, _Scalar(1) ) )
00210     {
00211       rowB = sqrt( colNorm/rowNorm );
00212       colB = Scalar(1)/rowB;
00213 
00214       isBalanced = false;
00215       return false;
00216     }
00217     else{
00218       return true; }
00219   }
00220 }
00221 
00222 
00223 template< typename _Scalar, int _Deg >
00224 void companion<_Scalar,_Deg>::balance()
00225 {
00226   EIGEN_STATIC_ASSERT( Deg == Dynamic || 1 < Deg, YOU_MADE_A_PROGRAMMING_MISTAKE );
00227   const Index deg   = m_monic.size();
00228   const Index deg_1 = deg-1;
00229 
00230   bool hasConverged=false;
00231   while( !hasConverged )
00232   {
00233     hasConverged = true;
00234     Scalar colNorm,rowNorm;
00235     Scalar colB,rowB;
00236 
00237     //First row, first column excluding the diagonal
00238     //==============================================
00239     colNorm = abs(m_bl_diag[0]);
00240     rowNorm = abs(m_monic[0]);
00241 
00242     //Compute balancing of the row and the column
00243     if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
00244     {
00245       m_bl_diag[0] *= colB;
00246       m_monic[0] *= rowB;
00247     }
00248 
00249     //Middle rows and columns excluding the diagonal
00250     //==============================================
00251     for( Index i=1; i<deg_1; ++i )
00252     {
00253       // column norm, excluding the diagonal
00254       colNorm = abs(m_bl_diag[i]);
00255 
00256       // row norm, excluding the diagonal
00257       rowNorm = abs(m_bl_diag[i-1]) + abs(m_monic[i]);
00258 
00259       //Compute balancing of the row and the column
00260       if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
00261       {
00262         m_bl_diag[i]   *= colB;
00263         m_bl_diag[i-1] *= rowB;
00264         m_monic[i]     *= rowB;
00265       }
00266     }
00267 
00268     //Last row, last column excluding the diagonal
00269     //============================================
00270     const Index ebl = m_bl_diag.size()-1;
00271     VectorBlock<RightColumn,Deg_1> headMonic( m_monic, 0, deg_1 );
00272     colNorm = headMonic.array().abs().sum();
00273     rowNorm = abs( m_bl_diag[ebl] );
00274 
00275     //Compute balancing of the row and the column
00276     if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
00277     {
00278       headMonic      *= colB;
00279       m_bl_diag[ebl] *= rowB;
00280     }
00281   }
00282 }
00283 
00284 } // end namespace internal
00285 
00286 #endif // EIGEN_COMPANION_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:35