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
17namespace Eigen {
18
19namespace internal {
20
21#ifndef EIGEN_PARSED_BY_DOXYGEN
22
23template <typename T>
24T radix(){ return 2; }
25
26template <typename T>
27T radix2(){ return radix<T>()*radix<T>(); }
28
29template<int Size>
30struct decrement_if_fixed_size
31{
32  enum {
33    ret = (Size == Dynamic) ? Dynamic : Size-1 };
34};
35
36#endif
37
38template< typename _Scalar, int _Deg >
39class 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:
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>
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>
84    companion( const VectorType& poly ){
85      setPolynomial( poly ); }
86
87  public:
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
140template< typename _Scalar, int _Deg >
141inline
142bool 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
183template< typename _Scalar, int _Deg >
184inline
185bool 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
210template< typename _Scalar, int _Deg >
211void 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