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 (at) inria.fr> 6 // Copyright (C) 2010 Jitse Niesen <jitse (at) 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 typename MatrixType::Index Index; 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 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 {} 101 102 /** \brief Constructor; computes Schur decomposition of given matrix. 103 * 104 * \param[in] matrix Square matrix whose Schur decomposition is to be computed. 105 * \param[in] computeU If true, both T and U are computed; if false, only T is computed. 106 * 107 * This constructor calls compute() to compute the Schur decomposition. 108 * 109 * \sa matrixT() and matrixU() for examples. 110 */ 111 ComplexSchur(const MatrixType& matrix, bool computeU = true) 112 : m_matT(matrix.rows(),matrix.cols()), 113 m_matU(matrix.rows(),matrix.cols()), 114 m_hess(matrix.rows()), 115 m_isInitialized(false), 116 m_matUisUptodate(false) 117 { 118 compute(matrix, computeU); 119 } 120 121 /** \brief Returns the unitary matrix in the Schur decomposition. 122 * 123 * \returns A const reference to the matrix U. 124 * 125 * It is assumed that either the constructor 126 * ComplexSchur(const MatrixType& matrix, bool computeU) or the 127 * member function compute(const MatrixType& matrix, bool computeU) 128 * has been called before to compute the Schur decomposition of a 129 * matrix, and that \p computeU was set to true (the default 130 * value). 131 * 132 * Example: \include ComplexSchur_matrixU.cpp 133 * Output: \verbinclude ComplexSchur_matrixU.out 134 */ 135 const ComplexMatrixType& matrixU() const 136 { 137 eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 138 eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition."); 139 return m_matU; 140 } 141 142 /** \brief Returns the triangular matrix in the Schur decomposition. 143 * 144 * \returns A const reference to the matrix T. 145 * 146 * It is assumed that either the constructor 147 * ComplexSchur(const MatrixType& matrix, bool computeU) or the 148 * member function compute(const MatrixType& matrix, bool computeU) 149 * has been called before to compute the Schur decomposition of a 150 * matrix. 151 * 152 * Note that this function returns a plain square matrix. If you want to reference 153 * only the upper triangular part, use: 154 * \code schur.matrixT().triangularView<Upper>() \endcode 155 * 156 * Example: \include ComplexSchur_matrixT.cpp 157 * Output: \verbinclude ComplexSchur_matrixT.out 158 */ 159 const ComplexMatrixType& matrixT() const 160 { 161 eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 162 return m_matT; 163 } 164 165 /** \brief Computes Schur decomposition of given matrix. 166 * 167 * \param[in] matrix Square matrix whose Schur decomposition is to be computed. 168 * \param[in] computeU If true, both T and U are computed; if false, only T is computed. 169 * \returns Reference to \c *this 170 * 171 * The Schur decomposition is computed by first reducing the 172 * matrix to Hessenberg form using the class 173 * HessenbergDecomposition. The Hessenberg matrix is then reduced 174 * to triangular form by performing QR iterations with a single 175 * shift. The cost of computing the Schur decomposition depends 176 * on the number of iterations; as a rough guide, it may be taken 177 * on the number of iterations; as a rough guide, it may be taken 178 * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops 179 * if \a computeU is false. 180 * 181 * Example: \include ComplexSchur_compute.cpp 182 * Output: \verbinclude ComplexSchur_compute.out 183 */ 184 ComplexSchur& compute(const MatrixType& matrix, bool computeU = true); 185 186 /** \brief Reports whether previous computation was successful. 187 * 188 * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 189 */ 190 ComputationInfo info() const 191 { 192 eigen_assert(m_isInitialized && "RealSchur is not initialized."); 193 return m_info; 194 } 195 196 /** \brief Maximum number of iterations. 197 * 198 * Maximum number of iterations allowed for an eigenvalue to converge. 199 */ 200 static const int m_maxIterations = 30; 201 202 protected: 203 ComplexMatrixType m_matT, m_matU; 204 HessenbergDecomposition<MatrixType> m_hess; 205 ComputationInfo m_info; 206 bool m_isInitialized; 207 bool m_matUisUptodate; 208 209 private: 210 bool subdiagonalEntryIsNeglegible(Index i); 211 ComplexScalar computeShift(Index iu, Index iter); 212 void reduceToTriangularForm(bool computeU); 213 friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>; 214 }; 215 216 /** If m_matT(i+1,i) is neglegible in floating point arithmetic 217 * compared to m_matT(i,i) and m_matT(j,j), then set it to zero and 218 * return true, else return false. */ 219 template<typename MatrixType> 220 inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i) 221 { 222 RealScalar d = internal::norm1(m_matT.coeff(i,i)) + internal::norm1(m_matT.coeff(i+1,i+1)); 223 RealScalar sd = internal::norm1(m_matT.coeff(i+1,i)); 224 if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon())) 225 { 226 m_matT.coeffRef(i+1,i) = ComplexScalar(0); 227 return true; 228 } 229 return false; 230 } 231 232 233 /** Compute the shift in the current QR iteration. */ 234 template<typename MatrixType> 235 typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter) 236 { 237 if (iter == 10 || iter == 20) 238 { 239 // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f 240 return internal::abs(internal::real(m_matT.coeff(iu,iu-1))) + internal::abs(internal::real(m_matT.coeff(iu-1,iu-2))); 241 } 242 243 // compute the shift as one of the eigenvalues of t, the 2x2 244 // diagonal block on the bottom of the active submatrix 245 Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1); 246 RealScalar normt = t.cwiseAbs().sum(); 247 t /= normt; // the normalization by sf is to avoid under/overflow 248 249 ComplexScalar b = t.coeff(0,1) * t.coeff(1,0); 250 ComplexScalar c = t.coeff(0,0) - t.coeff(1,1); 251 ComplexScalar disc = sqrt(c*c + RealScalar(4)*b); 252 ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b; 253 ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1); 254 ComplexScalar eival1 = (trace + disc) / RealScalar(2); 255 ComplexScalar eival2 = (trace - disc) / RealScalar(2); 256 257 if(internal::norm1(eival1) > internal::norm1(eival2)) 258 eival2 = det / eival1; 259 else 260 eival1 = det / eival2; 261 262 // choose the eigenvalue closest to the bottom entry of the diagonal 263 if(internal::norm1(eival1-t.coeff(1,1)) < internal::norm1(eival2-t.coeff(1,1))) 264 return normt * eival1; 265 else 266 return normt * eival2; 267 } 268 269 270 template<typename MatrixType> 271 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU) 272 { 273 m_matUisUptodate = false; 274 eigen_assert(matrix.cols() == matrix.rows()); 275 276 if(matrix.cols() == 1) 277 { 278 m_matT = matrix.template cast<ComplexScalar>(); 279 if(computeU) m_matU = ComplexMatrixType::Identity(1,1); 280 m_info = Success; 281 m_isInitialized = true; 282 m_matUisUptodate = computeU; 283 return *this; 284 } 285 286 internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU); 287 reduceToTriangularForm(computeU); 288 return *this; 289 } 290 291 namespace internal { 292 293 /* Reduce given matrix to Hessenberg form */ 294 template<typename MatrixType, bool IsComplex> 295 struct complex_schur_reduce_to_hessenberg 296 { 297 // this is the implementation for the case IsComplex = true 298 static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU) 299 { 300 _this.m_hess.compute(matrix); 301 _this.m_matT = _this.m_hess.matrixH(); 302 if(computeU) _this.m_matU = _this.m_hess.matrixQ(); 303 } 304 }; 305 306 template<typename MatrixType> 307 struct complex_schur_reduce_to_hessenberg<MatrixType, false> 308 { 309 static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU) 310 { 311 typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar; 312 typedef typename ComplexSchur<MatrixType>::ComplexMatrixType ComplexMatrixType; 313 314 // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar 315 _this.m_hess.compute(matrix); 316 _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>(); 317 if(computeU) 318 { 319 // This may cause an allocation which seems to be avoidable 320 MatrixType Q = _this.m_hess.matrixQ(); 321 _this.m_matU = Q.template cast<ComplexScalar>(); 322 } 323 } 324 }; 325 326 } // end namespace internal 327 328 // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration. 329 template<typename MatrixType> 330 void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU) 331 { 332 // The matrix m_matT is divided in three parts. 333 // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 334 // Rows il,...,iu is the part we are working on (the active submatrix). 335 // Rows iu+1,...,end are already brought in triangular form. 336 Index iu = m_matT.cols() - 1; 337 Index il; 338 Index iter = 0; // number of iterations we are working on the (iu,iu) element 339 340 while(true) 341 { 342 // find iu, the bottom row of the active submatrix 343 while(iu > 0) 344 { 345 if(!subdiagonalEntryIsNeglegible(iu-1)) break; 346 iter = 0; 347 --iu; 348 } 349 350 // if iu is zero then we are done; the whole matrix is triangularized 351 if(iu==0) break; 352 353 // if we spent too many iterations on the current element, we give up 354 iter++; 355 if(iter > m_maxIterations * m_matT.cols()) break; 356 357 // find il, the top row of the active submatrix 358 il = iu-1; 359 while(il > 0 && !subdiagonalEntryIsNeglegible(il-1)) 360 { 361 --il; 362 } 363 364 /* perform the QR step using Givens rotations. The first rotation 365 creates a bulge; the (il+2,il) element becomes nonzero. This 366 bulge is chased down to the bottom of the active submatrix. */ 367 368 ComplexScalar shift = computeShift(iu, iter); 369 JacobiRotation<ComplexScalar> rot; 370 rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il)); 371 m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint()); 372 m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot); 373 if(computeU) m_matU.applyOnTheRight(il, il+1, rot); 374 375 for(Index i=il+1 ; i<iu ; i++) 376 { 377 rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1)); 378 m_matT.coeffRef(i+1,i-1) = ComplexScalar(0); 379 m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint()); 380 m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot); 381 if(computeU) m_matU.applyOnTheRight(i, i+1, rot); 382 } 383 } 384 385 if(iter <= m_maxIterations * m_matT.cols()) 386 m_info = Success; 387 else 388 m_info = NoConvergence; 389 390 m_isInitialized = true; 391 m_matUisUptodate = computeU; 392 } 393 394 } // end namespace Eigen 395 396 #endif // EIGEN_COMPLEX_SCHUR_H 397
