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