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 "./EigenvaluesCommon.h"
00031 #include "./HessenbergDecomposition.h"
00032
00033 namespace internal {
00034 template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg;
00035 }
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065 template<typename _MatrixType> class ComplexSchur
00066 {
00067 public:
00068 typedef _MatrixType MatrixType;
00069 enum {
00070 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00071 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00072 Options = MatrixType::Options,
00073 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00074 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00075 };
00076
00077
00078 typedef typename MatrixType::Scalar Scalar;
00079 typedef typename NumTraits<Scalar>::Real RealScalar;
00080 typedef typename MatrixType::Index Index;
00081
00082
00083
00084
00085
00086
00087
00088 typedef std::complex<RealScalar> ComplexScalar;
00089
00090
00091
00092
00093
00094
00095 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType;
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108 ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
00109 : m_matT(size,size),
00110 m_matU(size,size),
00111 m_hess(size),
00112 m_isInitialized(false),
00113 m_matUisUptodate(false)
00114 {}
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125 ComplexSchur(const MatrixType& matrix, bool computeU = true)
00126 : m_matT(matrix.rows(),matrix.cols()),
00127 m_matU(matrix.rows(),matrix.cols()),
00128 m_hess(matrix.rows()),
00129 m_isInitialized(false),
00130 m_matUisUptodate(false)
00131 {
00132 compute(matrix, computeU);
00133 }
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149 const ComplexMatrixType& matrixU() const
00150 {
00151 eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
00152 eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
00153 return m_matU;
00154 }
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173 const ComplexMatrixType& matrixT() const
00174 {
00175 eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
00176 return m_matT;
00177 }
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198 ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
00199
00200
00201
00202
00203
00204 ComputationInfo info() const
00205 {
00206 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
00207 return m_info;
00208 }
00209
00210
00211
00212
00213
00214 static const int m_maxIterations = 30;
00215
00216 protected:
00217 ComplexMatrixType m_matT, m_matU;
00218 HessenbergDecomposition<MatrixType> m_hess;
00219 ComputationInfo m_info;
00220 bool m_isInitialized;
00221 bool m_matUisUptodate;
00222
00223 private:
00224 bool subdiagonalEntryIsNeglegible(Index i);
00225 ComplexScalar computeShift(Index iu, Index iter);
00226 void reduceToTriangularForm(bool computeU);
00227 friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
00228 };
00229
00230 namespace internal {
00231
00232
00233 template<typename RealScalar>
00234 std::complex<RealScalar> sqrt(const std::complex<RealScalar> &z)
00235 {
00236 RealScalar t, tre, tim;
00237
00238 t = abs(z);
00239
00240 if (abs(real(z)) <= abs(imag(z)))
00241 {
00242
00243 tre = sqrt(RealScalar(0.5)*(t + real(z)));
00244 tim = sqrt(RealScalar(0.5)*(t - real(z)));
00245 }
00246 else
00247 {
00248
00249 if (z.real() > RealScalar(0))
00250 {
00251 tre = t + z.real();
00252 tim = abs(imag(z))*sqrt(RealScalar(0.5)/tre);
00253 tre = sqrt(RealScalar(0.5)*tre);
00254 }
00255 else
00256 {
00257 tim = t - z.real();
00258 tre = abs(imag(z))*sqrt(RealScalar(0.5)/tim);
00259 tim = sqrt(RealScalar(0.5)*tim);
00260 }
00261 }
00262 if(z.imag() < RealScalar(0))
00263 tim = -tim;
00264
00265 return (std::complex<RealScalar>(tre,tim));
00266 }
00267 }
00268
00269
00270
00271
00272
00273 template<typename MatrixType>
00274 inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
00275 {
00276 RealScalar d = internal::norm1(m_matT.coeff(i,i)) + internal::norm1(m_matT.coeff(i+1,i+1));
00277 RealScalar sd = internal::norm1(m_matT.coeff(i+1,i));
00278 if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
00279 {
00280 m_matT.coeffRef(i+1,i) = ComplexScalar(0);
00281 return true;
00282 }
00283 return false;
00284 }
00285
00286
00287
00288 template<typename MatrixType>
00289 typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
00290 {
00291 if (iter == 10 || iter == 20)
00292 {
00293
00294 return internal::abs(internal::real(m_matT.coeff(iu,iu-1))) + internal::abs(internal::real(m_matT.coeff(iu-1,iu-2)));
00295 }
00296
00297
00298
00299 Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
00300 RealScalar normt = t.cwiseAbs().sum();
00301 t /= normt;
00302
00303 ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
00304 ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
00305 ComplexScalar disc = internal::sqrt(c*c + RealScalar(4)*b);
00306 ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
00307 ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
00308 ComplexScalar eival1 = (trace + disc) / RealScalar(2);
00309 ComplexScalar eival2 = (trace - disc) / RealScalar(2);
00310
00311 if(internal::norm1(eival1) > internal::norm1(eival2))
00312 eival2 = det / eival1;
00313 else
00314 eival1 = det / eival2;
00315
00316
00317 if(internal::norm1(eival1-t.coeff(1,1)) < internal::norm1(eival2-t.coeff(1,1)))
00318 return normt * eival1;
00319 else
00320 return normt * eival2;
00321 }
00322
00323
00324 template<typename MatrixType>
00325 ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
00326 {
00327 m_matUisUptodate = false;
00328 eigen_assert(matrix.cols() == matrix.rows());
00329
00330 if(matrix.cols() == 1)
00331 {
00332 m_matT = matrix.template cast<ComplexScalar>();
00333 if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
00334 m_info = Success;
00335 m_isInitialized = true;
00336 m_matUisUptodate = computeU;
00337 return *this;
00338 }
00339
00340 internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU);
00341 reduceToTriangularForm(computeU);
00342 return *this;
00343 }
00344
00345 namespace internal {
00346
00347
00348 template<typename MatrixType, bool IsComplex>
00349 struct complex_schur_reduce_to_hessenberg
00350 {
00351
00352 static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
00353 {
00354 _this.m_hess.compute(matrix);
00355 _this.m_matT = _this.m_hess.matrixH();
00356 if(computeU) _this.m_matU = _this.m_hess.matrixQ();
00357 }
00358 };
00359
00360 template<typename MatrixType>
00361 struct complex_schur_reduce_to_hessenberg<MatrixType, false>
00362 {
00363 static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
00364 {
00365 typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
00366 typedef typename ComplexSchur<MatrixType>::ComplexMatrixType ComplexMatrixType;
00367
00368
00369 _this.m_hess.compute(matrix);
00370 _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
00371 if(computeU)
00372 {
00373
00374 MatrixType Q = _this.m_hess.matrixQ();
00375 _this.m_matU = Q.template cast<ComplexScalar>();
00376 }
00377 }
00378 };
00379
00380 }
00381
00382
00383 template<typename MatrixType>
00384 void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
00385 {
00386
00387
00388
00389
00390 Index iu = m_matT.cols() - 1;
00391 Index il;
00392 Index iter = 0;
00393
00394 while(true)
00395 {
00396
00397 while(iu > 0)
00398 {
00399 if(!subdiagonalEntryIsNeglegible(iu-1)) break;
00400 iter = 0;
00401 --iu;
00402 }
00403
00404
00405 if(iu==0) break;
00406
00407
00408 iter++;
00409 if(iter > m_maxIterations) break;
00410
00411
00412 il = iu-1;
00413 while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
00414 {
00415 --il;
00416 }
00417
00418
00419
00420
00421
00422 ComplexScalar shift = computeShift(iu, iter);
00423 JacobiRotation<ComplexScalar> rot;
00424 rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
00425 m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
00426 m_matT.topRows(std::min(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
00427 if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
00428
00429 for(Index i=il+1 ; i<iu ; i++)
00430 {
00431 rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
00432 m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
00433 m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
00434 m_matT.topRows(std::min(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
00435 if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
00436 }
00437 }
00438
00439 if(iter <= m_maxIterations)
00440 m_info = Success;
00441 else
00442 m_info = NoConvergence;
00443
00444 m_isInitialized = true;
00445 m_matUisUptodate = computeU;
00446 }
00447
00448 #endif // EIGEN_COMPLEX_SCHUR_H