ComplexSchur.h (17274B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2009 Claire Maurice 5 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 6 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> 7 // 8 // This Source Code Form is subject to the terms of the Mozilla 9 // Public License v. 2.0. If a copy of the MPL was not distributed 10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 11 12 #ifndef EIGEN_COMPLEX_SCHUR_H 13 #define EIGEN_COMPLEX_SCHUR_H 14 15 #include "./HessenbergDecomposition.h" 16 17 namespace Eigen { 18 19 namespace internal { 20 template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg; 21 } 22 23 /** \eigenvalues_module \ingroup Eigenvalues_Module 24 * 25 * 26 * \class ComplexSchur 27 * 28 * \brief Performs a complex Schur decomposition of a real or complex square matrix 29 * 30 * \tparam _MatrixType the type of the matrix of which we are 31 * computing the Schur decomposition; this is expected to be an 32 * instantiation of the Matrix class template. 33 * 34 * Given a real or complex square matrix A, this class computes the 35 * Schur decomposition: \f$ A = U T U^*\f$ where U is a unitary 36 * complex matrix, and T is a complex upper triangular matrix. The 37 * diagonal of the matrix T corresponds to the eigenvalues of the 38 * matrix A. 39 * 40 * Call the function compute() to compute the Schur decomposition of 41 * a given matrix. Alternatively, you can use the 42 * ComplexSchur(const MatrixType&, bool) constructor which computes 43 * the Schur decomposition at construction time. Once the 44 * decomposition is computed, you can use the matrixU() and matrixT() 45 * functions to retrieve the matrices U and V in the decomposition. 46 * 47 * \note This code is inspired from Jampack 48 * 49 * \sa class RealSchur, class EigenSolver, class ComplexEigenSolver 50 */ 51 template<typename _MatrixType> class ComplexSchur 52 { 53 public: 54 typedef _MatrixType MatrixType; 55 enum { 56 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 57 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 58 Options = MatrixType::Options, 59 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 60 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 61 }; 62 63 /** \brief Scalar type for matrices of type \p _MatrixType. */ 64 typedef typename MatrixType::Scalar Scalar; 65 typedef typename NumTraits<Scalar>::Real RealScalar; 66 typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 67 68 /** \brief Complex scalar type for \p _MatrixType. 69 * 70 * This is \c std::complex<Scalar> if #Scalar is real (e.g., 71 * \c float or \c double) and just \c Scalar if #Scalar is 72 * complex. 73 */ 74 typedef std::complex<RealScalar> ComplexScalar; 75 76 /** \brief Type for the matrices in the Schur decomposition. 77 * 78 * This is a square matrix with entries of type #ComplexScalar. 79 * The size is the same as the size of \p _MatrixType. 80 */ 81 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType; 82 83 /** \brief Default constructor. 84 * 85 * \param [in] size Positive integer, size of the matrix whose Schur decomposition will be computed. 86 * 87 * The default constructor is useful in cases in which the user 88 * intends to perform decompositions via compute(). The \p size 89 * parameter is only used as a hint. It is not an error to give a 90 * wrong \p size, but it may impair performance. 91 * 92 * \sa compute() for an example. 93 */ 94 explicit ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) 95 : m_matT(size,size), 96 m_matU(size,size), 97 m_hess(size), 98 m_isInitialized(false), 99 m_matUisUptodate(false), 100 m_maxIters(-1) 101 {} 102 103 /** \brief Constructor; computes Schur decomposition of given matrix. 104 * 105 * \param[in] matrix Square matrix whose Schur decomposition is to be computed. 106 * \param[in] computeU If true, both T and U are computed; if false, only T is computed. 107 * 108 * This constructor calls compute() to compute the Schur decomposition. 109 * 110 * \sa matrixT() and matrixU() for examples. 111 */ 112 template<typename InputType> 113 explicit ComplexSchur(const EigenBase<InputType>& matrix, bool computeU = true) 114 : m_matT(matrix.rows(),matrix.cols()), 115 m_matU(matrix.rows(),matrix.cols()), 116 m_hess(matrix.rows()), 117 m_isInitialized(false), 118 m_matUisUptodate(false), 119 m_maxIters(-1) 120 { 121 compute(matrix.derived(), computeU); 122 } 123 124 /** \brief Returns the unitary matrix in the Schur decomposition. 125 * 126 * \returns A const reference to the matrix U. 127 * 128 * It is assumed that either the constructor 129 * ComplexSchur(const MatrixType& matrix, bool computeU) or the 130 * member function compute(const MatrixType& matrix, bool computeU) 131 * has been called before to compute the Schur decomposition of a 132 * matrix, and that \p computeU was set to true (the default 133 * value). 134 * 135 * Example: \include ComplexSchur_matrixU.cpp 136 * Output: \verbinclude ComplexSchur_matrixU.out 137 */ 138 const ComplexMatrixType& matrixU() const 139 { 140 eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 141 eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition."); 142 return m_matU; 143 } 144 145 /** \brief Returns the triangular matrix in the Schur decomposition. 146 * 147 * \returns A const reference to the matrix T. 148 * 149 * It is assumed that either the constructor 150 * ComplexSchur(const MatrixType& matrix, bool computeU) or the 151 * member function compute(const MatrixType& matrix, bool computeU) 152 * has been called before to compute the Schur decomposition of a 153 * matrix. 154 * 155 * Note that this function returns a plain square matrix. If you want to reference 156 * only the upper triangular part, use: 157 * \code schur.matrixT().triangularView<Upper>() \endcode 158 * 159 * Example: \include ComplexSchur_matrixT.cpp 160 * Output: \verbinclude ComplexSchur_matrixT.out 161 */ 162 const ComplexMatrixType& matrixT() const 163 { 164 eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 165 return m_matT; 166 } 167 168 /** \brief Computes Schur decomposition of given matrix. 169 * 170 * \param[in] matrix Square matrix whose Schur decomposition is to be computed. 171 * \param[in] computeU If true, both T and U are computed; if false, only T is computed. 172 173 * \returns Reference to \c *this 174 * 175 * The Schur decomposition is computed by first reducing the 176 * matrix to Hessenberg form using the class 177 * HessenbergDecomposition. The Hessenberg matrix is then reduced 178 * to triangular form by performing QR iterations with a single 179 * shift. The cost of computing the Schur decomposition depends 180 * on the number of iterations; as a rough guide, it may be taken 181 * on the number of iterations; as a rough guide, it may be taken 182 * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops 183 * if \a computeU is false. 184 * 185 * Example: \include ComplexSchur_compute.cpp 186 * Output: \verbinclude ComplexSchur_compute.out 187 * 188 * \sa compute(const MatrixType&, bool, Index) 189 */ 190 template<typename InputType> 191 ComplexSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true); 192 193 /** \brief Compute Schur decomposition from a given Hessenberg matrix 194 * \param[in] matrixH Matrix in Hessenberg form H 195 * \param[in] matrixQ orthogonal matrix Q that transform a matrix A to H : A = Q H Q^T 196 * \param computeU Computes the matriX U of the Schur vectors 197 * \return Reference to \c *this 198 * 199 * This routine assumes that the matrix is already reduced in Hessenberg form matrixH 200 * using either the class HessenbergDecomposition or another mean. 201 * It computes the upper quasi-triangular matrix T of the Schur decomposition of H 202 * When computeU is true, this routine computes the matrix U such that 203 * A = U T U^T = (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix 204 * 205 * NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix 206 * is not available, the user should give an identity matrix (Q.setIdentity()) 207 * 208 * \sa compute(const MatrixType&, bool) 209 */ 210 template<typename HessMatrixType, typename OrthMatrixType> 211 ComplexSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU=true); 212 213 /** \brief Reports whether previous computation was successful. 214 * 215 * \returns \c Success if computation was successful, \c NoConvergence otherwise. 216 */ 217 ComputationInfo info() const 218 { 219 eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 220 return m_info; 221 } 222 223 /** \brief Sets the maximum number of iterations allowed. 224 * 225 * If not specified by the user, the maximum number of iterations is m_maxIterationsPerRow times the size 226 * of the matrix. 227 */ 228 ComplexSchur& setMaxIterations(Index maxIters) 229 { 230 m_maxIters = maxIters; 231 return *this; 232 } 233 234 /** \brief Returns the maximum number of iterations. */ 235 Index getMaxIterations() 236 { 237 return m_maxIters; 238 } 239 240 /** \brief Maximum number of iterations per row. 241 * 242 * If not otherwise specified, the maximum number of iterations is this number times the size of the 243 * matrix. It is currently set to 30. 244 */ 245 static const int m_maxIterationsPerRow = 30; 246 247 protected: 248 ComplexMatrixType m_matT, m_matU; 249 HessenbergDecomposition<MatrixType> m_hess; 250 ComputationInfo m_info; 251 bool m_isInitialized; 252 bool m_matUisUptodate; 253 Index m_maxIters; 254 255 private: 256 bool subdiagonalEntryIsNeglegible(Index i); 257 ComplexScalar computeShift(Index iu, Index iter); 258 void reduceToTriangularForm(bool computeU); 259 friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>; 260 }; 261 262 /** If m_matT(i+1,i) is neglegible in floating point arithmetic 263 * compared to m_matT(i,i) and m_matT(j,j), then set it to zero and 264 * return true, else return false. */ 265 template<typename MatrixType> 266 inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i) 267 { 268 RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1)); 269 RealScalar sd = numext::norm1(m_matT.coeff(i+1,i)); 270 if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon())) 271 { 272 m_matT.coeffRef(i+1,i) = ComplexScalar(0); 273 return true; 274 } 275 return false; 276 } 277 278 279 /** Compute the shift in the current QR iteration. */ 280 template<typename MatrixType> 281 typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter) 282 { 283 using std::abs; 284 if (iter == 10 || iter == 20) 285 { 286 // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f 287 return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2))); 288 } 289 290 // compute the shift as one of the eigenvalues of t, the 2x2 291 // diagonal block on the bottom of the active submatrix 292 Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1); 293 RealScalar normt = t.cwiseAbs().sum(); 294 t /= normt; // the normalization by sf is to avoid under/overflow 295 296 ComplexScalar b = t.coeff(0,1) * t.coeff(1,0); 297 ComplexScalar c = t.coeff(0,0) - t.coeff(1,1); 298 ComplexScalar disc = sqrt(c*c + RealScalar(4)*b); 299 ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b; 300 ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1); 301 ComplexScalar eival1 = (trace + disc) / RealScalar(2); 302 ComplexScalar eival2 = (trace - disc) / RealScalar(2); 303 RealScalar eival1_norm = numext::norm1(eival1); 304 RealScalar eival2_norm = numext::norm1(eival2); 305 // A division by zero can only occur if eival1==eival2==0. 306 // In this case, det==0, and all we have to do is checking that eival2_norm!=0 307 if(eival1_norm > eival2_norm) 308 eival2 = det / eival1; 309 else if(eival2_norm!=RealScalar(0)) 310 eival1 = det / eival2; 311 312 // choose the eigenvalue closest to the bottom entry of the diagonal 313 if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1))) 314 return normt * eival1; 315 else 316 return normt * eival2; 317 } 318 319 320 template<typename MatrixType> 321 template<typename InputType> 322 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU) 323 { 324 m_matUisUptodate = false; 325 eigen_assert(matrix.cols() == matrix.rows()); 326 327 if(matrix.cols() == 1) 328 { 329 m_matT = matrix.derived().template cast<ComplexScalar>(); 330 if(computeU) m_matU = ComplexMatrixType::Identity(1,1); 331 m_info = Success; 332 m_isInitialized = true; 333 m_matUisUptodate = computeU; 334 return *this; 335 } 336 337 internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix.derived(), computeU); 338 computeFromHessenberg(m_matT, m_matU, computeU); 339 return *this; 340 } 341 342 template<typename MatrixType> 343 template<typename HessMatrixType, typename OrthMatrixType> 344 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU) 345 { 346 m_matT = matrixH; 347 if(computeU) 348 m_matU = matrixQ; 349 reduceToTriangularForm(computeU); 350 return *this; 351 } 352 namespace internal { 353 354 /* Reduce given matrix to Hessenberg form */ 355 template<typename MatrixType, bool IsComplex> 356 struct complex_schur_reduce_to_hessenberg 357 { 358 // this is the implementation for the case IsComplex = true 359 static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU) 360 { 361 _this.m_hess.compute(matrix); 362 _this.m_matT = _this.m_hess.matrixH(); 363 if(computeU) _this.m_matU = _this.m_hess.matrixQ(); 364 } 365 }; 366 367 template<typename MatrixType> 368 struct complex_schur_reduce_to_hessenberg<MatrixType, false> 369 { 370 static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU) 371 { 372 typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar; 373 374 // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar 375 _this.m_hess.compute(matrix); 376 _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>(); 377 if(computeU) 378 { 379 // This may cause an allocation which seems to be avoidable 380 MatrixType Q = _this.m_hess.matrixQ(); 381 _this.m_matU = Q.template cast<ComplexScalar>(); 382 } 383 } 384 }; 385 386 } // end namespace internal 387 388 // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration. 389 template<typename MatrixType> 390 void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU) 391 { 392 Index maxIters = m_maxIters; 393 if (maxIters == -1) 394 maxIters = m_maxIterationsPerRow * m_matT.rows(); 395 396 // The matrix m_matT is divided in three parts. 397 // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 398 // Rows il,...,iu is the part we are working on (the active submatrix). 399 // Rows iu+1,...,end are already brought in triangular form. 400 Index iu = m_matT.cols() - 1; 401 Index il; 402 Index iter = 0; // number of iterations we are working on the (iu,iu) element 403 Index totalIter = 0; // number of iterations for whole matrix 404 405 while(true) 406 { 407 // find iu, the bottom row of the active submatrix 408 while(iu > 0) 409 { 410 if(!subdiagonalEntryIsNeglegible(iu-1)) break; 411 iter = 0; 412 --iu; 413 } 414 415 // if iu is zero then we are done; the whole matrix is triangularized 416 if(iu==0) break; 417 418 // if we spent too many iterations, we give up 419 iter++; 420 totalIter++; 421 if(totalIter > maxIters) break; 422 423 // find il, the top row of the active submatrix 424 il = iu-1; 425 while(il > 0 && !subdiagonalEntryIsNeglegible(il-1)) 426 { 427 --il; 428 } 429 430 /* perform the QR step using Givens rotations. The first rotation 431 creates a bulge; the (il+2,il) element becomes nonzero. This 432 bulge is chased down to the bottom of the active submatrix. */ 433 434 ComplexScalar shift = computeShift(iu, iter); 435 JacobiRotation<ComplexScalar> rot; 436 rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il)); 437 m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint()); 438 m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot); 439 if(computeU) m_matU.applyOnTheRight(il, il+1, rot); 440 441 for(Index i=il+1 ; i<iu ; i++) 442 { 443 rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1)); 444 m_matT.coeffRef(i+1,i-1) = ComplexScalar(0); 445 m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint()); 446 m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot); 447 if(computeU) m_matU.applyOnTheRight(i, i+1, rot); 448 } 449 } 450 451 if(totalIter <= maxIters) 452 m_info = Success; 453 else 454 m_info = NoConvergence; 455 456 m_isInitialized = true; 457 m_matUisUptodate = computeU; 458 } 459 460 } // end namespace Eigen 461 462 #endif // EIGEN_COMPLEX_SCHUR_H