Eigne Matrix Class Library
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
ComplexSchur.h
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2009 Claire Maurice 00005 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 00006 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> 00007 // 00008 // This Source Code Form is subject to the terms of the Mozilla 00009 // Public License v. 2.0. If a copy of the MPL was not distributed 00010 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00011 00012 #ifndef EIGEN_COMPLEX_SCHUR_H 00013 #define EIGEN_COMPLEX_SCHUR_H 00014 00015 #include "./HessenbergDecomposition.h" 00016 00017 namespace Eigen { 00018 00019 namespace internal { 00020 template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg; 00021 } 00022 00023 /** \eigenvalues_module \ingroup Eigenvalues_Module 00024 * 00025 * 00026 * \class ComplexSchur 00027 * 00028 * \brief Performs a complex Schur decomposition of a real or complex square matrix 00029 * 00030 * \tparam _MatrixType the type of the matrix of which we are 00031 * computing the Schur decomposition; this is expected to be an 00032 * instantiation of the Matrix class template. 00033 * 00034 * Given a real or complex square matrix A, this class computes the 00035 * Schur decomposition: \f$ A = U T U^*\f$ where U is a unitary 00036 * complex matrix, and T is a complex upper triangular matrix. The 00037 * diagonal of the matrix T corresponds to the eigenvalues of the 00038 * matrix A. 00039 * 00040 * Call the function compute() to compute the Schur decomposition of 00041 * a given matrix. Alternatively, you can use the 00042 * ComplexSchur(const MatrixType&, bool) constructor which computes 00043 * the Schur decomposition at construction time. Once the 00044 * decomposition is computed, you can use the matrixU() and matrixT() 00045 * functions to retrieve the matrices U and V in the decomposition. 00046 * 00047 * \note This code is inspired from Jampack 00048 * 00049 * \sa class RealSchur, class EigenSolver, class ComplexEigenSolver 00050 */ 00051 template<typename _MatrixType> class ComplexSchur 00052 { 00053 public: 00054 typedef _MatrixType MatrixType; 00055 enum { 00056 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00057 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00058 Options = MatrixType::Options, 00059 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00060 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00061 }; 00062 00063 /** \brief Scalar type for matrices of type \p _MatrixType. */ 00064 typedef typename MatrixType::Scalar Scalar; 00065 typedef typename NumTraits<Scalar>::Real RealScalar; 00066 typedef typename MatrixType::Index Index; 00067 00068 /** \brief Complex scalar type for \p _MatrixType. 00069 * 00070 * This is \c std::complex<Scalar> if #Scalar is real (e.g., 00071 * \c float or \c double) and just \c Scalar if #Scalar is 00072 * complex. 00073 */ 00074 typedef std::complex<RealScalar> ComplexScalar; 00075 00076 /** \brief Type for the matrices in the Schur decomposition. 00077 * 00078 * This is a square matrix with entries of type #ComplexScalar. 00079 * The size is the same as the size of \p _MatrixType. 00080 */ 00081 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType; 00082 00083 /** \brief Default constructor. 00084 * 00085 * \param [in] size Positive integer, size of the matrix whose Schur decomposition will be computed. 00086 * 00087 * The default constructor is useful in cases in which the user 00088 * intends to perform decompositions via compute(). The \p size 00089 * parameter is only used as a hint. It is not an error to give a 00090 * wrong \p size, but it may impair performance. 00091 * 00092 * \sa compute() for an example. 00093 */ 00094 ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) 00095 : m_matT(size,size), 00096 m_matU(size,size), 00097 m_hess(size), 00098 m_isInitialized(false), 00099 m_matUisUptodate(false), 00100 m_maxIters(-1) 00101 {} 00102 00103 /** \brief Constructor; computes Schur decomposition of given matrix. 00104 * 00105 * \param[in] matrix Square matrix whose Schur decomposition is to be computed. 00106 * \param[in] computeU If true, both T and U are computed; if false, only T is computed. 00107 * 00108 * This constructor calls compute() to compute the Schur decomposition. 00109 * 00110 * \sa matrixT() and matrixU() for examples. 00111 */ 00112 ComplexSchur(const MatrixType& matrix, bool computeU = true) 00113 : m_matT(matrix.rows(),matrix.cols()), 00114 m_matU(matrix.rows(),matrix.cols()), 00115 m_hess(matrix.rows()), 00116 m_isInitialized(false), 00117 m_matUisUptodate(false), 00118 m_maxIters(-1) 00119 { 00120 compute(matrix, computeU); 00121 } 00122 00123 /** \brief Returns the unitary matrix in the Schur decomposition. 00124 * 00125 * \returns A const reference to the matrix U. 00126 * 00127 * It is assumed that either the constructor 00128 * ComplexSchur(const MatrixType& matrix, bool computeU) or the 00129 * member function compute(const MatrixType& matrix, bool computeU) 00130 * has been called before to compute the Schur decomposition of a 00131 * matrix, and that \p computeU was set to true (the default 00132 * value). 00133 * 00134 * Example: \include ComplexSchur_matrixU.cpp 00135 * Output: \verbinclude ComplexSchur_matrixU.out 00136 */ 00137 const ComplexMatrixType & matrixU() const 00138 { 00139 eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 00140 eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition."); 00141 return m_matU; 00142 } 00143 00144 /** \brief Returns the triangular matrix in the Schur decomposition. 00145 * 00146 * \returns A const reference to the matrix T. 00147 * 00148 * It is assumed that either the constructor 00149 * ComplexSchur(const MatrixType& matrix, bool computeU) or the 00150 * member function compute(const MatrixType& matrix, bool computeU) 00151 * has been called before to compute the Schur decomposition of a 00152 * matrix. 00153 * 00154 * Note that this function returns a plain square matrix. If you want to reference 00155 * only the upper triangular part, use: 00156 * \code schur.matrixT().triangularView<Upper>() \endcode 00157 * 00158 * Example: \include ComplexSchur_matrixT.cpp 00159 * Output: \verbinclude ComplexSchur_matrixT.out 00160 */ 00161 const ComplexMatrixType & matrixT() const 00162 { 00163 eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 00164 return m_matT; 00165 } 00166 00167 /** \brief Computes Schur decomposition of given matrix. 00168 * 00169 * \param[in] matrix Square matrix whose Schur decomposition is to be computed. 00170 * \param[in] computeU If true, both T and U are computed; if false, only T is computed. 00171 00172 * \returns Reference to \c *this 00173 * 00174 * The Schur decomposition is computed by first reducing the 00175 * matrix to Hessenberg form using the class 00176 * HessenbergDecomposition. The Hessenberg matrix is then reduced 00177 * to triangular form by performing QR iterations with a single 00178 * shift. The cost of computing the Schur decomposition depends 00179 * on the number of iterations; as a rough guide, it may be taken 00180 * on the number of iterations; as a rough guide, it may be taken 00181 * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops 00182 * if \a computeU is false. 00183 * 00184 * Example: \include ComplexSchur_compute.cpp 00185 * Output: \verbinclude ComplexSchur_compute.out 00186 * 00187 * \sa compute(const MatrixType&, bool, Index) 00188 */ 00189 ComplexSchur & compute(const MatrixType& matrix, bool computeU = true); 00190 00191 /** \brief Compute Schur decomposition from a given Hessenberg matrix 00192 * \param[in] matrixH Matrix in Hessenberg form H 00193 * \param[in] matrixQ orthogonal matrix Q that transform a matrix A to H : A = Q H Q^T 00194 * \param computeU Computes the matriX U of the Schur vectors 00195 * \return Reference to \c *this 00196 * 00197 * This routine assumes that the matrix is already reduced in Hessenberg form matrixH 00198 * using either the class HessenbergDecomposition or another mean. 00199 * It computes the upper quasi-triangular matrix T of the Schur decomposition of H 00200 * When computeU is true, this routine computes the matrix U such that 00201 * A = U T U^T = (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix 00202 * 00203 * NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix 00204 * is not available, the user should give an identity matrix (Q.setIdentity()) 00205 * 00206 * \sa compute(const MatrixType&, bool) 00207 */ 00208 template<typename HessMatrixType, typename OrthMatrixType> 00209 ComplexSchur & computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU=true); 00210 00211 /** \brief Reports whether previous computation was successful. 00212 * 00213 * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 00214 */ 00215 ComputationInfo info() const 00216 { 00217 eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 00218 return m_info; 00219 } 00220 00221 /** \brief Sets the maximum number of iterations allowed. 00222 * 00223 * If not specified by the user, the maximum number of iterations is m_maxIterationsPerRow times the size 00224 * of the matrix. 00225 */ 00226 ComplexSchur & setMaxIterations(Index maxIters) 00227 { 00228 m_maxIters = maxIters; 00229 return *this; 00230 } 00231 00232 /** \brief Returns the maximum number of iterations. */ 00233 Index getMaxIterations() 00234 { 00235 return m_maxIters; 00236 } 00237 00238 /** \brief Maximum number of iterations per row. 00239 * 00240 * If not otherwise specified, the maximum number of iterations is this number times the size of the 00241 * matrix. It is currently set to 30. 00242 */ 00243 static const int m_maxIterationsPerRow = 30; 00244 00245 protected: 00246 ComplexMatrixType m_matT, m_matU; 00247 HessenbergDecomposition<MatrixType> m_hess; 00248 ComputationInfo m_info; 00249 bool m_isInitialized; 00250 bool m_matUisUptodate; 00251 Index m_maxIters; 00252 00253 private: 00254 bool subdiagonalEntryIsNeglegible(Index i); 00255 ComplexScalar computeShift(Index iu, Index iter); 00256 void reduceToTriangularForm(bool computeU); 00257 friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>; 00258 }; 00259 00260 /** If m_matT(i+1,i) is neglegible in floating point arithmetic 00261 * compared to m_matT(i,i) and m_matT(j,j), then set it to zero and 00262 * return true, else return false. */ 00263 template<typename MatrixType> 00264 inline bool ComplexSchur <MatrixType>::subdiagonalEntryIsNeglegible(Index i) 00265 { 00266 RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1)); 00267 RealScalar sd = numext::norm1(m_matT.coeff(i+1,i)); 00268 if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon())) 00269 { 00270 m_matT.coeffRef(i+1,i) = ComplexScalar(0); 00271 return true; 00272 } 00273 return false; 00274 } 00275 00276 00277 /** Compute the shift in the current QR iteration. */ 00278 template<typename MatrixType> 00279 typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter) 00280 { 00281 using std::abs; 00282 if (iter == 10 || iter == 20) 00283 { 00284 // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f 00285 return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2))); 00286 } 00287 00288 // compute the shift as one of the eigenvalues of t, the 2x2 00289 // diagonal block on the bottom of the active submatrix 00290 Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1); 00291 RealScalar normt = t.cwiseAbs().sum(); 00292 t /= normt; // the normalization by sf is to avoid under/overflow 00293 00294 ComplexScalar b = t.coeff(0,1) * t.coeff(1,0); 00295 ComplexScalar c = t.coeff(0,0) - t.coeff(1,1); 00296 ComplexScalar disc = sqrt(c*c + RealScalar(4)*b); 00297 ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b; 00298 ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1); 00299 ComplexScalar eival1 = (trace + disc) / RealScalar(2); 00300 ComplexScalar eival2 = (trace - disc) / RealScalar(2); 00301 00302 if(numext::norm1(eival1) > numext::norm1(eival2)) 00303 eival2 = det / eival1; 00304 else 00305 eival1 = det / eival2; 00306 00307 // choose the eigenvalue closest to the bottom entry of the diagonal 00308 if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1))) 00309 return normt * eival1; 00310 else 00311 return normt * eival2; 00312 } 00313 00314 00315 template<typename MatrixType> 00316 ComplexSchur<MatrixType> & ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU) 00317 { 00318 m_matUisUptodate = false; 00319 eigen_assert(matrix.cols() == matrix.rows()); 00320 00321 if(matrix.cols() == 1) 00322 { 00323 m_matT = matrix.template cast<ComplexScalar>(); 00324 if(computeU) m_matU = ComplexMatrixType::Identity(1,1); 00325 m_info = Success; 00326 m_isInitialized = true; 00327 m_matUisUptodate = computeU; 00328 return *this; 00329 } 00330 00331 internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU); 00332 computeFromHessenberg(m_matT, m_matU, computeU); 00333 return *this; 00334 } 00335 00336 template<typename MatrixType> 00337 template<typename HessMatrixType, typename OrthMatrixType> 00338 ComplexSchur<MatrixType> & ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU) 00339 { 00340 m_matT = matrixH; 00341 if(computeU) 00342 m_matU = matrixQ; 00343 reduceToTriangularForm(computeU); 00344 return *this; 00345 } 00346 namespace internal { 00347 00348 /* Reduce given matrix to Hessenberg form */ 00349 template<typename MatrixType, bool IsComplex> 00350 struct complex_schur_reduce_to_hessenberg 00351 { 00352 // this is the implementation for the case IsComplex = true 00353 static void run(ComplexSchur<MatrixType> & _this, const MatrixType& matrix, bool computeU) 00354 { 00355 _this.m_hess.compute(matrix); 00356 _this.m_matT = _this.m_hess.matrixH(); 00357 if(computeU) _this.m_matU = _this.m_hess.matrixQ(); 00358 } 00359 }; 00360 00361 template<typename MatrixType> 00362 struct complex_schur_reduce_to_hessenberg<MatrixType, false> 00363 { 00364 static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU) 00365 { 00366 typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar; 00367 00368 // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar 00369 _this.m_hess.compute(matrix); 00370 _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>(); 00371 if(computeU) 00372 { 00373 // This may cause an allocation which seems to be avoidable 00374 MatrixType Q = _this.m_hess.matrixQ(); 00375 _this.m_matU = Q.template cast<ComplexScalar>(); 00376 } 00377 } 00378 }; 00379 00380 } // end namespace internal 00381 00382 // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration. 00383 template<typename MatrixType> 00384 void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU) 00385 { 00386 Index maxIters = m_maxIters; 00387 if (maxIters == -1) 00388 maxIters = m_maxIterationsPerRow * m_matT.rows(); 00389 00390 // The matrix m_matT is divided in three parts. 00391 // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 00392 // Rows il,...,iu is the part we are working on (the active submatrix). 00393 // Rows iu+1,...,end are already brought in triangular form. 00394 Index iu = m_matT.cols() - 1; 00395 Index il; 00396 Index iter = 0; // number of iterations we are working on the (iu,iu) element 00397 Index totalIter = 0; // number of iterations for whole matrix 00398 00399 while(true) 00400 { 00401 // find iu, the bottom row of the active submatrix 00402 while(iu > 0) 00403 { 00404 if(!subdiagonalEntryIsNeglegible(iu-1)) break; 00405 iter = 0; 00406 --iu; 00407 } 00408 00409 // if iu is zero then we are done; the whole matrix is triangularized 00410 if(iu==0) break; 00411 00412 // if we spent too many iterations, we give up 00413 iter++; 00414 totalIter++; 00415 if(totalIter > maxIters) break; 00416 00417 // find il, the top row of the active submatrix 00418 il = iu-1; 00419 while(il > 0 && !subdiagonalEntryIsNeglegible(il-1)) 00420 { 00421 --il; 00422 } 00423 00424 /* perform the QR step using Givens rotations. The first rotation 00425 creates a bulge; the (il+2,il) element becomes nonzero. This 00426 bulge is chased down to the bottom of the active submatrix. */ 00427 00428 ComplexScalar shift = computeShift(iu, iter); 00429 JacobiRotation<ComplexScalar> rot; 00430 rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il)); 00431 m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint()); 00432 m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot); 00433 if(computeU) m_matU.applyOnTheRight(il, il+1, rot); 00434 00435 for(Index i=il+1 ; i<iu ; i++) 00436 { 00437 rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1)); 00438 m_matT.coeffRef(i+1,i-1) = ComplexScalar(0); 00439 m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint()); 00440 m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot); 00441 if(computeU) m_matU.applyOnTheRight(i, i+1, rot); 00442 } 00443 } 00444 00445 if(totalIter <= maxIters) 00446 m_info = Success; 00447 else 00448 m_info = NoConvergence; 00449 00450 m_isInitialized = true; 00451 m_matUisUptodate = computeU; 00452 } 00453 00454 } // end namespace Eigen 00455 00456 #endif // EIGEN_COMPLEX_SCHUR_H
Generated on Tue Jul 12 2022 17:46:51 by 1.7.2