00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027 #ifndef EIGEN_COMPLEX_SCHUR_H
00028 #define EIGEN_COMPLEX_SCHUR_H
00029
00030 #include "./HessenbergDecomposition.h"
00031
00032 namespace Eigen {
00033
00034 namespace internal {
00035 template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg;
00036 }
00037
00066 template<typename _MatrixType> class ComplexSchur
00067 {
00068 public:
00069 typedef _MatrixType MatrixType;
00070 enum {
00071 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00072 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00073 Options = MatrixType::Options,
00074 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00075 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00076 };
00077
00079 typedef typename MatrixType::Scalar Scalar;
00080 typedef typename NumTraits<Scalar>::Real RealScalar;
00081 typedef typename MatrixType::Index Index;
00082
00089 typedef std::complex<RealScalar> ComplexScalar;
00090
00096 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType;
00097
00109 ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
00110 : m_matT(size,size),
00111 m_matU(size,size),
00112 m_hess(size),
00113 m_isInitialized(false),
00114 m_matUisUptodate(false)
00115 {}
00116
00126 ComplexSchur(const MatrixType& matrix, bool computeU = true)
00127 : m_matT(matrix.rows(),matrix.cols()),
00128 m_matU(matrix.rows(),matrix.cols()),
00129 m_hess(matrix.rows()),
00130 m_isInitialized(false),
00131 m_matUisUptodate(false)
00132 {
00133 compute(matrix, computeU);
00134 }
00135
00150 const ComplexMatrixType& matrixU() const
00151 {
00152 eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
00153 eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
00154 return m_matU;
00155 }
00156
00174 const ComplexMatrixType& matrixT() const
00175 {
00176 eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
00177 return m_matT;
00178 }
00179
00199 ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
00200
00205 ComputationInfo info() const
00206 {
00207 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
00208 return m_info;
00209 }
00210
00215 static const int m_maxIterations = 30;
00216
00217 protected:
00218 ComplexMatrixType m_matT, m_matU;
00219 HessenbergDecomposition<MatrixType> m_hess;
00220 ComputationInfo m_info;
00221 bool m_isInitialized;
00222 bool m_matUisUptodate;
00223
00224 private:
00225 bool subdiagonalEntryIsNeglegible(Index i);
00226 ComplexScalar computeShift(Index iu, Index iter);
00227 void reduceToTriangularForm(bool computeU);
00228 friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
00229 };
00230
00234 template<typename MatrixType>
00235 inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
00236 {
00237 RealScalar d = internal::norm1(m_matT.coeff(i,i)) + internal::norm1(m_matT.coeff(i+1,i+1));
00238 RealScalar sd = internal::norm1(m_matT.coeff(i+1,i));
00239 if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
00240 {
00241 m_matT.coeffRef(i+1,i) = ComplexScalar(0);
00242 return true;
00243 }
00244 return false;
00245 }
00246
00247
00249 template<typename MatrixType>
00250 typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
00251 {
00252 if (iter == 10 || iter == 20)
00253 {
00254
00255 return internal::abs(internal::real(m_matT.coeff(iu,iu-1))) + internal::abs(internal::real(m_matT.coeff(iu-1,iu-2)));
00256 }
00257
00258
00259
00260 Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
00261 RealScalar normt = t.cwiseAbs().sum();
00262 t /= normt;
00263
00264 ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
00265 ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
00266 ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
00267 ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
00268 ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
00269 ComplexScalar eival1 = (trace + disc) / RealScalar(2);
00270 ComplexScalar eival2 = (trace - disc) / RealScalar(2);
00271
00272 if(internal::norm1(eival1) > internal::norm1(eival2))
00273 eival2 = det / eival1;
00274 else
00275 eival1 = det / eival2;
00276
00277
00278 if(internal::norm1(eival1-t.coeff(1,1)) < internal::norm1(eival2-t.coeff(1,1)))
00279 return normt * eival1;
00280 else
00281 return normt * eival2;
00282 }
00283
00284
00285 template<typename MatrixType>
00286 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
00287 {
00288 m_matUisUptodate = false;
00289 eigen_assert(matrix.cols() == matrix.rows());
00290
00291 if(matrix.cols() == 1)
00292 {
00293 m_matT = matrix.template cast<ComplexScalar>();
00294 if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
00295 m_info = Success;
00296 m_isInitialized = true;
00297 m_matUisUptodate = computeU;
00298 return *this;
00299 }
00300
00301 internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU);
00302 reduceToTriangularForm(computeU);
00303 return *this;
00304 }
00305
00306 namespace internal {
00307
00308
00309 template<typename MatrixType, bool IsComplex>
00310 struct complex_schur_reduce_to_hessenberg
00311 {
00312
00313 static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
00314 {
00315 _this.m_hess.compute(matrix);
00316 _this.m_matT = _this.m_hess.matrixH();
00317 if(computeU) _this.m_matU = _this.m_hess.matrixQ();
00318 }
00319 };
00320
00321 template<typename MatrixType>
00322 struct complex_schur_reduce_to_hessenberg<MatrixType, false>
00323 {
00324 static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
00325 {
00326 typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
00327 typedef typename ComplexSchur<MatrixType>::ComplexMatrixType ComplexMatrixType;
00328
00329
00330 _this.m_hess.compute(matrix);
00331 _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
00332 if(computeU)
00333 {
00334
00335 MatrixType Q = _this.m_hess.matrixQ();
00336 _this.m_matU = Q.template cast<ComplexScalar>();
00337 }
00338 }
00339 };
00340
00341 }
00342
00343
00344 template<typename MatrixType>
00345 void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
00346 {
00347
00348
00349
00350
00351 Index iu = m_matT.cols() - 1;
00352 Index il;
00353 Index iter = 0;
00354
00355 while(true)
00356 {
00357
00358 while(iu > 0)
00359 {
00360 if(!subdiagonalEntryIsNeglegible(iu-1)) break;
00361 iter = 0;
00362 --iu;
00363 }
00364
00365
00366 if(iu==0) break;
00367
00368
00369 iter++;
00370 if(iter > m_maxIterations) break;
00371
00372
00373 il = iu-1;
00374 while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
00375 {
00376 --il;
00377 }
00378
00379
00380
00381
00382
00383 ComplexScalar shift = computeShift(iu, iter);
00384 JacobiRotation<ComplexScalar> rot;
00385 rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
00386 m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
00387 m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
00388 if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
00389
00390 for(Index i=il+1 ; i<iu ; i++)
00391 {
00392 rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
00393 m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
00394 m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
00395 m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
00396 if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
00397 }
00398 }
00399
00400 if(iter <= m_maxIterations)
00401 m_info = Success;
00402 else
00403 m_info = NoConvergence;
00404
00405 m_isInitialized = true;
00406 m_matUisUptodate = computeU;
00407 }
00408
00409 }
00410
00411 #endif // EIGEN_COMPLEX_SCHUR_H