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 #ifndef EIGEN_REAL_SCHUR_H
00027 #define EIGEN_REAL_SCHUR_H
00028
00029 #include "./EigenvaluesCommon.h"
00030 #include "./HessenbergDecomposition.h"
00031
00032
00033
00034
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
00066
00067
00068 template<typename _MatrixType> class RealSchur
00069 {
00070 public:
00071 typedef _MatrixType MatrixType;
00072 enum {
00073 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00074 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00075 Options = MatrixType::Options,
00076 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00077 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00078 };
00079 typedef typename MatrixType::Scalar Scalar;
00080 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
00081 typedef typename MatrixType::Index Index;
00082
00083 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
00084 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097 RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
00098 : m_matT(size, size),
00099 m_matU(size, size),
00100 m_workspaceVector(size),
00101 m_hess(size),
00102 m_isInitialized(false),
00103 m_matUisUptodate(false)
00104 { }
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115
00116 RealSchur(const MatrixType& matrix, bool computeU = true)
00117 : m_matT(matrix.rows(),matrix.cols()),
00118 m_matU(matrix.rows(),matrix.cols()),
00119 m_workspaceVector(matrix.rows()),
00120 m_hess(matrix.rows()),
00121 m_isInitialized(false),
00122 m_matUisUptodate(false)
00123 {
00124 compute(matrix, computeU);
00125 }
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137
00138 const MatrixType& matrixU() const
00139 {
00140 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
00141 eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the RealSchur decomposition.");
00142 return m_matU;
00143 }
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 const MatrixType& matrixT() const
00156 {
00157 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
00158 return m_matT;
00159 }
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178 RealSchur& compute(const MatrixType& matrix, bool computeU = true);
00179
00180
00181
00182
00183
00184 ComputationInfo info() const
00185 {
00186 eigen_assert(m_isInitialized && "RealSchur is not initialized.");
00187 return m_info;
00188 }
00189
00190
00191
00192
00193
00194 static const int m_maxIterations = 40;
00195
00196 private:
00197
00198 MatrixType m_matT;
00199 MatrixType m_matU;
00200 ColumnVectorType m_workspaceVector;
00201 HessenbergDecomposition<MatrixType> m_hess;
00202 ComputationInfo m_info;
00203 bool m_isInitialized;
00204 bool m_matUisUptodate;
00205
00206 typedef Matrix<Scalar,3,1> Vector3s;
00207
00208 Scalar computeNormOfT();
00209 Index findSmallSubdiagEntry(Index iu, Scalar norm);
00210 void splitOffTwoRows(Index iu, bool computeU, Scalar exshift);
00211 void computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo);
00212 void initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector);
00213 void performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace);
00214 };
00215
00216
00217 template<typename MatrixType>
00218 RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU)
00219 {
00220 assert(matrix.cols() == matrix.rows());
00221
00222
00223 m_hess.compute(matrix);
00224 m_matT = m_hess.matrixH();
00225 if (computeU)
00226 m_matU = m_hess.matrixQ();
00227
00228
00229 m_workspaceVector.resize(m_matT.cols());
00230 Scalar* workspace = &m_workspaceVector.coeffRef(0);
00231
00232
00233
00234
00235
00236 Index iu = m_matT.cols() - 1;
00237 Index iter = 0;
00238 Scalar exshift = 0.0;
00239 Scalar norm = computeNormOfT();
00240
00241 while (iu >= 0)
00242 {
00243 Index il = findSmallSubdiagEntry(iu, norm);
00244
00245
00246 if (il == iu)
00247 {
00248 m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
00249 if (iu > 0)
00250 m_matT.coeffRef(iu, iu-1) = Scalar(0);
00251 iu--;
00252 iter = 0;
00253 }
00254 else if (il == iu-1)
00255 {
00256 splitOffTwoRows(iu, computeU, exshift);
00257 iu -= 2;
00258 iter = 0;
00259 }
00260 else
00261 {
00262
00263 Vector3s firstHouseholderVector(0,0,0), shiftInfo;
00264 computeShift(iu, iter, exshift, shiftInfo);
00265 iter = iter + 1;
00266 if (iter > m_maxIterations) break;
00267 Index im;
00268 initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
00269 performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
00270 }
00271 }
00272
00273 if(iter <= m_maxIterations)
00274 m_info = Success;
00275 else
00276 m_info = NoConvergence;
00277
00278 m_isInitialized = true;
00279 m_matUisUptodate = computeU;
00280 return *this;
00281 }
00282
00283
00284 template<typename MatrixType>
00285 inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
00286 {
00287 const Index size = m_matT.cols();
00288
00289
00290
00291 Scalar norm = 0.0;
00292 for (Index j = 0; j < size; ++j)
00293 norm += m_matT.row(j).segment(std::max(j-1,Index(0)), size-std::max(j-1,Index(0))).cwiseAbs().sum();
00294 return norm;
00295 }
00296
00297
00298 template<typename MatrixType>
00299 inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, Scalar norm)
00300 {
00301 Index res = iu;
00302 while (res > 0)
00303 {
00304 Scalar s = internal::abs(m_matT.coeff(res-1,res-1)) + internal::abs(m_matT.coeff(res,res));
00305 if (s == 0.0)
00306 s = norm;
00307 if (internal::abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
00308 break;
00309 res--;
00310 }
00311 return res;
00312 }
00313
00314
00315 template<typename MatrixType>
00316 inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scalar exshift)
00317 {
00318 const Index size = m_matT.cols();
00319
00320
00321
00322 Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
00323 Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
00324 m_matT.coeffRef(iu,iu) += exshift;
00325 m_matT.coeffRef(iu-1,iu-1) += exshift;
00326
00327 if (q >= 0)
00328 {
00329 Scalar z = internal::sqrt(internal::abs(q));
00330 JacobiRotation<Scalar> rot;
00331 if (p >= 0)
00332 rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
00333 else
00334 rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
00335
00336 m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
00337 m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
00338 m_matT.coeffRef(iu, iu-1) = Scalar(0);
00339 if (computeU)
00340 m_matU.applyOnTheRight(iu-1, iu, rot);
00341 }
00342
00343 if (iu > 1)
00344 m_matT.coeffRef(iu-1, iu-2) = Scalar(0);
00345 }
00346
00347
00348 template<typename MatrixType>
00349 inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
00350 {
00351 shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
00352 shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
00353 shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
00354
00355
00356 if (iter == 10)
00357 {
00358 exshift += shiftInfo.coeff(0);
00359 for (Index i = 0; i <= iu; ++i)
00360 m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
00361 Scalar s = internal::abs(m_matT.coeff(iu,iu-1)) + internal::abs(m_matT.coeff(iu-1,iu-2));
00362 shiftInfo.coeffRef(0) = Scalar(0.75) * s;
00363 shiftInfo.coeffRef(1) = Scalar(0.75) * s;
00364 shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
00365 }
00366
00367
00368 if (iter == 30)
00369 {
00370 Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
00371 s = s * s + shiftInfo.coeff(2);
00372 if (s > 0)
00373 {
00374 s = internal::sqrt(s);
00375 if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
00376 s = -s;
00377 s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
00378 s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
00379 exshift += s;
00380 for (Index i = 0; i <= iu; ++i)
00381 m_matT.coeffRef(i,i) -= s;
00382 shiftInfo.setConstant(Scalar(0.964));
00383 }
00384 }
00385 }
00386
00387
00388 template<typename MatrixType>
00389 inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu, const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
00390 {
00391 Vector3s& v = firstHouseholderVector;
00392
00393 for (im = iu-2; im >= il; --im)
00394 {
00395 const Scalar Tmm = m_matT.coeff(im,im);
00396 const Scalar r = shiftInfo.coeff(0) - Tmm;
00397 const Scalar s = shiftInfo.coeff(1) - Tmm;
00398 v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
00399 v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s;
00400 v.coeffRef(2) = m_matT.coeff(im+2,im+1);
00401 if (im == il) {
00402 break;
00403 }
00404 const Scalar lhs = m_matT.coeff(im,im-1) * (internal::abs(v.coeff(1)) + internal::abs(v.coeff(2)));
00405 const Scalar rhs = v.coeff(0) * (internal::abs(m_matT.coeff(im-1,im-1)) + internal::abs(Tmm) + internal::abs(m_matT.coeff(im+1,im+1)));
00406 if (internal::abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
00407 {
00408 break;
00409 }
00410 }
00411 }
00412
00413
00414 template<typename MatrixType>
00415 inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu, bool computeU, const Vector3s& firstHouseholderVector, Scalar* workspace)
00416 {
00417 assert(im >= il);
00418 assert(im <= iu-2);
00419
00420 const Index size = m_matT.cols();
00421
00422 for (Index k = im; k <= iu-2; ++k)
00423 {
00424 bool firstIteration = (k == im);
00425
00426 Vector3s v;
00427 if (firstIteration)
00428 v = firstHouseholderVector;
00429 else
00430 v = m_matT.template block<3,1>(k,k-1);
00431
00432 Scalar tau, beta;
00433 Matrix<Scalar, 2, 1> ess;
00434 v.makeHouseholder(ess, tau, beta);
00435
00436 if (beta != Scalar(0))
00437 {
00438 if (firstIteration && k > il)
00439 m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
00440 else if (!firstIteration)
00441 m_matT.coeffRef(k,k-1) = beta;
00442
00443
00444 m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
00445 m_matT.block(0, k, std::min(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
00446 if (computeU)
00447 m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
00448 }
00449 }
00450
00451 Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2);
00452 Scalar tau, beta;
00453 Matrix<Scalar, 1, 1> ess;
00454 v.makeHouseholder(ess, tau, beta);
00455
00456 if (beta != Scalar(0))
00457 {
00458 m_matT.coeffRef(iu-1, iu-2) = beta;
00459 m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
00460 m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
00461 if (computeU)
00462 m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
00463 }
00464
00465
00466 for (Index i = im+2; i <= iu; ++i)
00467 {
00468 m_matT.coeffRef(i,i-2) = Scalar(0);
00469 if (i > im+2)
00470 m_matT.coeffRef(i,i-3) = Scalar(0);
00471 }
00472 }
00473
00474 #endif // EIGEN_REAL_SCHUR_H