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_EIGENSOLVER_H
00027 #define EIGEN_EIGENSOLVER_H
00028
00029 #include "./EigenvaluesCommon.h"
00030 #include "./RealSchur.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
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078 template<typename _MatrixType> class EigenSolver
00079 {
00080 public:
00081
00082
00083 typedef _MatrixType MatrixType;
00084
00085 enum {
00086 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00087 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00088 Options = MatrixType::Options,
00089 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00090 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00091 };
00092
00093
00094 typedef typename MatrixType::Scalar Scalar;
00095 typedef typename NumTraits<Scalar>::Real RealScalar;
00096 typedef typename MatrixType::Index Index;
00097
00098
00099
00100
00101
00102
00103
00104 typedef std::complex<RealScalar> ComplexScalar;
00105
00106
00107
00108
00109
00110
00111 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
00112
00113
00114
00115
00116
00117
00118 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
00119
00120
00121
00122
00123
00124
00125
00126
00127 EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
00128
00129
00130
00131
00132
00133
00134
00135 EigenSolver(Index size)
00136 : m_eivec(size, size),
00137 m_eivalues(size),
00138 m_isInitialized(false),
00139 m_eigenvectorsOk(false),
00140 m_realSchur(size),
00141 m_matT(size, size),
00142 m_tmp(size)
00143 {}
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160 EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
00161 : m_eivec(matrix.rows(), matrix.cols()),
00162 m_eivalues(matrix.cols()),
00163 m_isInitialized(false),
00164 m_eigenvectorsOk(false),
00165 m_realSchur(matrix.cols()),
00166 m_matT(matrix.rows(), matrix.cols()),
00167 m_tmp(matrix.cols())
00168 {
00169 compute(matrix, computeEigenvectors);
00170 }
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192 EigenvectorsType eigenvectors() const;
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212 const MatrixType& pseudoEigenvectors() const
00213 {
00214 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00215 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00216 return m_eivec;
00217 }
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236 MatrixType pseudoEigenvalueMatrix() const;
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255 const EigenvalueType& eigenvalues() const
00256 {
00257 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00258 return m_eivalues;
00259 }
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288 EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
00289
00290 ComputationInfo info() const
00291 {
00292 eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00293 return m_realSchur.info();
00294 }
00295
00296 private:
00297 void doComputeEigenvectors();
00298
00299 protected:
00300 MatrixType m_eivec;
00301 EigenvalueType m_eivalues;
00302 bool m_isInitialized;
00303 bool m_eigenvectorsOk;
00304 RealSchur<MatrixType> m_realSchur;
00305 MatrixType m_matT;
00306
00307 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
00308 ColumnVectorType m_tmp;
00309 };
00310
00311 template<typename MatrixType>
00312 MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
00313 {
00314 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00315 Index n = m_eivalues.rows();
00316 MatrixType matD = MatrixType::Zero(n,n);
00317 for (Index i=0; i<n; ++i)
00318 {
00319 if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i))))
00320 matD.coeffRef(i,i) = internal::real(m_eivalues.coeff(i));
00321 else
00322 {
00323 matD.template block<2,2>(i,i) << internal::real(m_eivalues.coeff(i)), internal::imag(m_eivalues.coeff(i)),
00324 -internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i));
00325 ++i;
00326 }
00327 }
00328 return matD;
00329 }
00330
00331 template<typename MatrixType>
00332 typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const
00333 {
00334 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00335 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00336 Index n = m_eivec.cols();
00337 EigenvectorsType matV(n,n);
00338 for (Index j=0; j<n; ++j)
00339 {
00340 if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(j)), internal::real(m_eivalues.coeff(j))))
00341 {
00342
00343 matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
00344 }
00345 else
00346 {
00347
00348 for (Index i=0; i<n; ++i)
00349 {
00350 matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1));
00351 matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
00352 }
00353 matV.col(j).normalize();
00354 matV.col(j+1).normalize();
00355 ++j;
00356 }
00357 }
00358 return matV;
00359 }
00360
00361 template<typename MatrixType>
00362 EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
00363 {
00364 assert(matrix.cols() == matrix.rows());
00365
00366
00367 m_realSchur.compute(matrix, computeEigenvectors);
00368 if (m_realSchur.info() == Success)
00369 {
00370 m_matT = m_realSchur.matrixT();
00371 if (computeEigenvectors)
00372 m_eivec = m_realSchur.matrixU();
00373
00374
00375 m_eivalues.resize(matrix.cols());
00376 Index i = 0;
00377 while (i < matrix.cols())
00378 {
00379 if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
00380 {
00381 m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
00382 ++i;
00383 }
00384 else
00385 {
00386 Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
00387 Scalar z = internal::sqrt(internal::abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
00388 m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
00389 m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
00390 i += 2;
00391 }
00392 }
00393
00394
00395 if (computeEigenvectors)
00396 doComputeEigenvectors();
00397 }
00398
00399 m_isInitialized = true;
00400 m_eigenvectorsOk = computeEigenvectors;
00401
00402 return *this;
00403 }
00404
00405
00406 template<typename Scalar>
00407 std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
00408 {
00409 Scalar r,d;
00410 if (internal::abs(yr) > internal::abs(yi))
00411 {
00412 r = yi/yr;
00413 d = yr + r*yi;
00414 return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
00415 }
00416 else
00417 {
00418 r = yr/yi;
00419 d = yi + r*yr;
00420 return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
00421 }
00422 }
00423
00424
00425 template<typename MatrixType>
00426 void EigenSolver<MatrixType>::doComputeEigenvectors()
00427 {
00428 const Index size = m_eivec.cols();
00429 const Scalar eps = NumTraits<Scalar>::epsilon();
00430
00431
00432 Scalar norm = 0.0;
00433 for (Index j = 0; j < size; ++j)
00434 {
00435 norm += m_matT.row(j).segment(std::max(j-1,Index(0)), size-std::max(j-1,Index(0))).cwiseAbs().sum();
00436 }
00437
00438
00439 if (norm == 0.0)
00440 {
00441 return;
00442 }
00443
00444 for (Index n = size-1; n >= 0; n--)
00445 {
00446 Scalar p = m_eivalues.coeff(n).real();
00447 Scalar q = m_eivalues.coeff(n).imag();
00448
00449
00450 if (q == 0)
00451 {
00452 Scalar lastr=0, lastw=0;
00453 Index l = n;
00454
00455 m_matT.coeffRef(n,n) = 1.0;
00456 for (Index i = n-1; i >= 0; i--)
00457 {
00458 Scalar w = m_matT.coeff(i,i) - p;
00459 Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
00460
00461 if (m_eivalues.coeff(i).imag() < 0.0)
00462 {
00463 lastw = w;
00464 lastr = r;
00465 }
00466 else
00467 {
00468 l = i;
00469 if (m_eivalues.coeff(i).imag() == 0.0)
00470 {
00471 if (w != 0.0)
00472 m_matT.coeffRef(i,n) = -r / w;
00473 else
00474 m_matT.coeffRef(i,n) = -r / (eps * norm);
00475 }
00476 else
00477 {
00478 Scalar x = m_matT.coeff(i,i+1);
00479 Scalar y = m_matT.coeff(i+1,i);
00480 Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
00481 Scalar t = (x * lastr - lastw * r) / denom;
00482 m_matT.coeffRef(i,n) = t;
00483 if (internal::abs(x) > internal::abs(lastw))
00484 m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
00485 else
00486 m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
00487 }
00488
00489
00490 Scalar t = internal::abs(m_matT.coeff(i,n));
00491 if ((eps * t) * t > 1)
00492 m_matT.col(n).tail(size-i) /= t;
00493 }
00494 }
00495 }
00496 else if (q < 0)
00497 {
00498 Scalar lastra=0, lastsa=0, lastw=0;
00499 Index l = n-1;
00500
00501
00502 if (internal::abs(m_matT.coeff(n,n-1)) > internal::abs(m_matT.coeff(n-1,n)))
00503 {
00504 m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
00505 m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
00506 }
00507 else
00508 {
00509 std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
00510 m_matT.coeffRef(n-1,n-1) = internal::real(cc);
00511 m_matT.coeffRef(n-1,n) = internal::imag(cc);
00512 }
00513 m_matT.coeffRef(n,n-1) = 0.0;
00514 m_matT.coeffRef(n,n) = 1.0;
00515 for (Index i = n-2; i >= 0; i--)
00516 {
00517 Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
00518 Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
00519 Scalar w = m_matT.coeff(i,i) - p;
00520
00521 if (m_eivalues.coeff(i).imag() < 0.0)
00522 {
00523 lastw = w;
00524 lastra = ra;
00525 lastsa = sa;
00526 }
00527 else
00528 {
00529 l = i;
00530 if (m_eivalues.coeff(i).imag() == 0)
00531 {
00532 std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
00533 m_matT.coeffRef(i,n-1) = internal::real(cc);
00534 m_matT.coeffRef(i,n) = internal::imag(cc);
00535 }
00536 else
00537 {
00538
00539 Scalar x = m_matT.coeff(i,i+1);
00540 Scalar y = m_matT.coeff(i+1,i);
00541 Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
00542 Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
00543 if ((vr == 0.0) && (vi == 0.0))
00544 vr = eps * norm * (internal::abs(w) + internal::abs(q) + internal::abs(x) + internal::abs(y) + internal::abs(lastw));
00545
00546 std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
00547 m_matT.coeffRef(i,n-1) = internal::real(cc);
00548 m_matT.coeffRef(i,n) = internal::imag(cc);
00549 if (internal::abs(x) > (internal::abs(lastw) + internal::abs(q)))
00550 {
00551 m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
00552 m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
00553 }
00554 else
00555 {
00556 cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
00557 m_matT.coeffRef(i+1,n-1) = internal::real(cc);
00558 m_matT.coeffRef(i+1,n) = internal::imag(cc);
00559 }
00560 }
00561
00562
00563 Scalar t = std::max(internal::abs(m_matT.coeff(i,n-1)),internal::abs(m_matT.coeff(i,n)));
00564 if ((eps * t) * t > 1)
00565 m_matT.block(i, n-1, size-i, 2) /= t;
00566
00567 }
00568 }
00569 }
00570 }
00571
00572
00573 for (Index j = size-1; j >= 0; j--)
00574 {
00575 m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
00576 m_eivec.col(j) = m_tmp;
00577 }
00578 }
00579
00580 #endif // EIGEN_EIGENSOLVER_H