• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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_COMPANION_H
11 #define EIGEN_COMPANION_H
12 
13 // This file requires the user to include
14 // * Eigen/Core
15 // * Eigen/src/PolynomialSolver.h
16 
17 namespace Eigen {
18 
19 namespace internal {
20 
21 #ifndef EIGEN_PARSED_BY_DOXYGEN
22 
23 template <typename T>
radix()24 T radix(){ return 2; }
25 
26 template <typename T>
radix2()27 T radix2(){ return radix<T>()*radix<T>(); }
28 
29 template<int Size>
30 struct decrement_if_fixed_size
31 {
32   enum {
33     ret = (Size == Dynamic) ? Dynamic : Size-1 };
34 };
35 
36 #endif
37 
38 template< typename _Scalar, int _Deg >
39 class companion
40 {
41   public:
42     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)
43 
44     enum {
45       Deg = _Deg,
46       Deg_1=decrement_if_fixed_size<Deg>::ret
47     };
48 
49     typedef _Scalar                                Scalar;
50     typedef typename NumTraits<Scalar>::Real       RealScalar;
51     typedef Matrix<Scalar, Deg, 1>                 RightColumn;
52     //typedef DiagonalMatrix< Scalar, Deg_1, Deg_1 > BottomLeftDiagonal;
53     typedef Matrix<Scalar, Deg_1, 1>               BottomLeftDiagonal;
54 
55     typedef Matrix<Scalar, Deg, Deg>               DenseCompanionMatrixType;
56     typedef Matrix< Scalar, _Deg, Deg_1 >          LeftBlock;
57     typedef Matrix< Scalar, Deg_1, Deg_1 >         BottomLeftBlock;
58     typedef Matrix< Scalar, 1, Deg_1 >             LeftBlockFirstRow;
59 
60     typedef DenseIndex Index;
61 
62   public:
operator()63     EIGEN_STRONG_INLINE const _Scalar operator()(Index row, Index col ) const
64     {
65       if( m_bl_diag.rows() > col )
66       {
67         if( 0 < row ){ return m_bl_diag[col]; }
68         else{ return 0; }
69       }
70       else{ return m_monic[row]; }
71     }
72 
73   public:
74     template<typename VectorType>
setPolynomial(const VectorType & poly)75     void setPolynomial( const VectorType& poly )
76     {
77       const Index deg = poly.size()-1;
78       m_monic = -1/poly[deg] * poly.head(deg);
79       //m_bl_diag.setIdentity( deg-1 );
80       m_bl_diag.setOnes(deg-1);
81     }
82 
83     template<typename VectorType>
companion(const VectorType & poly)84     companion( const VectorType& poly ){
85       setPolynomial( poly ); }
86 
87   public:
denseMatrix()88     DenseCompanionMatrixType denseMatrix() const
89     {
90       const Index deg   = m_monic.size();
91       const Index deg_1 = deg-1;
92       DenseCompanionMatrixType companion(deg,deg);
93       companion <<
94         ( LeftBlock(deg,deg_1)
95           << LeftBlockFirstRow::Zero(1,deg_1),
96           BottomLeftBlock::Identity(deg-1,deg-1)*m_bl_diag.asDiagonal() ).finished()
97         , m_monic;
98       return companion;
99     }
100 
101 
102 
103   protected:
104     /** Helper function for the balancing algorithm.
105      * \returns true if the row and the column, having colNorm and rowNorm
106      * as norms, are balanced, false otherwise.
107      * colB and rowB are repectively the multipliers for
108      * the column and the row in order to balance them.
109      * */
110     bool balanced( Scalar colNorm, Scalar rowNorm,
111         bool& isBalanced, Scalar& colB, Scalar& rowB );
112 
113     /** Helper function for the balancing algorithm.
114      * \returns true if the row and the column, having colNorm and rowNorm
115      * as norms, are balanced, false otherwise.
116      * colB and rowB are repectively the multipliers for
117      * the column and the row in order to balance them.
118      * */
119     bool balancedR( Scalar colNorm, Scalar rowNorm,
120         bool& isBalanced, Scalar& colB, Scalar& rowB );
121 
122   public:
123     /**
124      * Balancing algorithm from B. N. PARLETT and C. REINSCH (1969)
125      * "Balancing a matrix for calculation of eigenvalues and eigenvectors"
126      * adapted to the case of companion matrices.
127      * A matrix with non zero row and non zero column is balanced
128      * for a certain norm if the i-th row and the i-th column
129      * have same norm for all i.
130      */
131     void balance();
132 
133   protected:
134       RightColumn                m_monic;
135       BottomLeftDiagonal         m_bl_diag;
136 };
137 
138 
139 
140 template< typename _Scalar, int _Deg >
141 inline
balanced(Scalar colNorm,Scalar rowNorm,bool & isBalanced,Scalar & colB,Scalar & rowB)142 bool companion<_Scalar,_Deg>::balanced( Scalar colNorm, Scalar rowNorm,
143     bool& isBalanced, Scalar& colB, Scalar& rowB )
144 {
145   if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; }
146   else
147   {
148     //To find the balancing coefficients, if the radix is 2,
149     //one finds \f$ \sigma \f$ such that
150     // \f$ 2^{2\sigma-1} < rowNorm / colNorm \le 2^{2\sigma+1} \f$
151     // then the balancing coefficient for the row is \f$ 1/2^{\sigma} \f$
152     // and the balancing coefficient for the column is \f$ 2^{\sigma} \f$
153     rowB = rowNorm / radix<Scalar>();
154     colB = Scalar(1);
155     const Scalar s = colNorm + rowNorm;
156 
157     while (colNorm < rowB)
158     {
159       colB *= radix<Scalar>();
160       colNorm *= radix2<Scalar>();
161     }
162 
163     rowB = rowNorm * radix<Scalar>();
164 
165     while (colNorm >= rowB)
166     {
167       colB /= radix<Scalar>();
168       colNorm /= radix2<Scalar>();
169     }
170 
171     //This line is used to avoid insubstantial balancing
172     if ((rowNorm + colNorm) < Scalar(0.95) * s * colB)
173     {
174       isBalanced = false;
175       rowB = Scalar(1) / colB;
176       return false;
177     }
178     else{
179       return true; }
180   }
181 }
182 
183 template< typename _Scalar, int _Deg >
184 inline
balancedR(Scalar colNorm,Scalar rowNorm,bool & isBalanced,Scalar & colB,Scalar & rowB)185 bool companion<_Scalar,_Deg>::balancedR( Scalar colNorm, Scalar rowNorm,
186     bool& isBalanced, Scalar& colB, Scalar& rowB )
187 {
188   if( Scalar(0) == colNorm || Scalar(0) == rowNorm ){ return true; }
189   else
190   {
191     /**
192      * Set the norm of the column and the row to the geometric mean
193      * of the row and column norm
194      */
195     const _Scalar q = colNorm/rowNorm;
196     if( !isApprox( q, _Scalar(1) ) )
197     {
198       rowB = sqrt( colNorm/rowNorm );
199       colB = Scalar(1)/rowB;
200 
201       isBalanced = false;
202       return false;
203     }
204     else{
205       return true; }
206   }
207 }
208 
209 
210 template< typename _Scalar, int _Deg >
balance()211 void companion<_Scalar,_Deg>::balance()
212 {
213   using std::abs;
214   EIGEN_STATIC_ASSERT( Deg == Dynamic || 1 < Deg, YOU_MADE_A_PROGRAMMING_MISTAKE );
215   const Index deg   = m_monic.size();
216   const Index deg_1 = deg-1;
217 
218   bool hasConverged=false;
219   while( !hasConverged )
220   {
221     hasConverged = true;
222     Scalar colNorm,rowNorm;
223     Scalar colB,rowB;
224 
225     //First row, first column excluding the diagonal
226     //==============================================
227     colNorm = abs(m_bl_diag[0]);
228     rowNorm = abs(m_monic[0]);
229 
230     //Compute balancing of the row and the column
231     if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
232     {
233       m_bl_diag[0] *= colB;
234       m_monic[0] *= rowB;
235     }
236 
237     //Middle rows and columns excluding the diagonal
238     //==============================================
239     for( Index i=1; i<deg_1; ++i )
240     {
241       // column norm, excluding the diagonal
242       colNorm = abs(m_bl_diag[i]);
243 
244       // row norm, excluding the diagonal
245       rowNorm = abs(m_bl_diag[i-1]) + abs(m_monic[i]);
246 
247       //Compute balancing of the row and the column
248       if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
249       {
250         m_bl_diag[i]   *= colB;
251         m_bl_diag[i-1] *= rowB;
252         m_monic[i]     *= rowB;
253       }
254     }
255 
256     //Last row, last column excluding the diagonal
257     //============================================
258     const Index ebl = m_bl_diag.size()-1;
259     VectorBlock<RightColumn,Deg_1> headMonic( m_monic, 0, deg_1 );
260     colNorm = headMonic.array().abs().sum();
261     rowNorm = abs( m_bl_diag[ebl] );
262 
263     //Compute balancing of the row and the column
264     if( !balanced( colNorm, rowNorm, hasConverged, colB, rowB ) )
265     {
266       headMonic      *= colB;
267       m_bl_diag[ebl] *= rowB;
268     }
269   }
270 }
271 
272 } // end namespace internal
273 
274 } // end namespace Eigen
275 
276 #endif // EIGEN_COMPANION_H
277