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 "./RealSchur.h"
00030
00031 namespace Eigen {
00032
00079 template<typename _MatrixType> class EigenSolver
00080 {
00081 public:
00082
00084 typedef _MatrixType MatrixType;
00085
00086 enum {
00087 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00088 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00089 Options = MatrixType::Options,
00090 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00091 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00092 };
00093
00095 typedef typename MatrixType::Scalar Scalar;
00096 typedef typename NumTraits<Scalar>::Real RealScalar;
00097 typedef typename MatrixType::Index Index;
00098
00105 typedef std::complex<RealScalar> ComplexScalar;
00106
00112 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
00113
00119 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
00120
00128 EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
00129
00136 EigenSolver(Index size)
00137 : m_eivec(size, size),
00138 m_eivalues(size),
00139 m_isInitialized(false),
00140 m_eigenvectorsOk(false),
00141 m_realSchur(size),
00142 m_matT(size, size),
00143 m_tmp(size)
00144 {}
00145
00161 EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
00162 : m_eivec(matrix.rows(), matrix.cols()),
00163 m_eivalues(matrix.cols()),
00164 m_isInitialized(false),
00165 m_eigenvectorsOk(false),
00166 m_realSchur(matrix.cols()),
00167 m_matT(matrix.rows(), matrix.cols()),
00168 m_tmp(matrix.cols())
00169 {
00170 compute(matrix, computeEigenvectors);
00171 }
00172
00193 EigenvectorsType eigenvectors() const;
00194
00213 const MatrixType& pseudoEigenvectors() const
00214 {
00215 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00216 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00217 return m_eivec;
00218 }
00219
00238 MatrixType pseudoEigenvalueMatrix() const;
00239
00258 const EigenvalueType& eigenvalues() const
00259 {
00260 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00261 return m_eivalues;
00262 }
00263
00291 EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
00292
00293 ComputationInfo info() const
00294 {
00295 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00296 return m_realSchur.info();
00297 }
00298
00299 private:
00300 void doComputeEigenvectors();
00301
00302 protected:
00303 MatrixType m_eivec;
00304 EigenvalueType m_eivalues;
00305 bool m_isInitialized;
00306 bool m_eigenvectorsOk;
00307 RealSchur<MatrixType> m_realSchur;
00308 MatrixType m_matT;
00309
00310 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
00311 ColumnVectorType m_tmp;
00312 };
00313
00314 template<typename MatrixType>
00315 MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
00316 {
00317 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00318 Index n = m_eivalues.rows();
00319 MatrixType matD = MatrixType::Zero(n,n);
00320 for (Index i=0; i<n; ++i)
00321 {
00322 if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i))))
00323 matD.coeffRef(i,i) = internal::real(m_eivalues.coeff(i));
00324 else
00325 {
00326 matD.template block<2,2>(i,i) << internal::real(m_eivalues.coeff(i)), internal::imag(m_eivalues.coeff(i)),
00327 -internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i));
00328 ++i;
00329 }
00330 }
00331 return matD;
00332 }
00333
00334 template<typename MatrixType>
00335 typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const
00336 {
00337 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00338 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00339 Index n = m_eivec.cols();
00340 EigenvectorsType matV(n,n);
00341 for (Index j=0; j<n; ++j)
00342 {
00343 if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(j)), internal::real(m_eivalues.coeff(j))) || j+1==n)
00344 {
00345
00346 matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
00347 matV.col(j).normalize();
00348 }
00349 else
00350 {
00351
00352 for (Index i=0; i<n; ++i)
00353 {
00354 matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1));
00355 matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
00356 }
00357 matV.col(j).normalize();
00358 matV.col(j+1).normalize();
00359 ++j;
00360 }
00361 }
00362 return matV;
00363 }
00364
00365 template<typename MatrixType>
00366 EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
00367 {
00368 assert(matrix.cols() == matrix.rows());
00369
00370
00371 m_realSchur.compute(matrix, computeEigenvectors);
00372 if (m_realSchur.info() == Success)
00373 {
00374 m_matT = m_realSchur.matrixT();
00375 if (computeEigenvectors)
00376 m_eivec = m_realSchur.matrixU();
00377
00378
00379 m_eivalues.resize(matrix.cols());
00380 Index i = 0;
00381 while (i < matrix.cols())
00382 {
00383 if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
00384 {
00385 m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
00386 ++i;
00387 }
00388 else
00389 {
00390 Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
00391 Scalar z = internal::sqrt(internal::abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
00392 m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
00393 m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
00394 i += 2;
00395 }
00396 }
00397
00398
00399 if (computeEigenvectors)
00400 doComputeEigenvectors();
00401 }
00402
00403 m_isInitialized = true;
00404 m_eigenvectorsOk = computeEigenvectors;
00405
00406 return *this;
00407 }
00408
00409
00410 template<typename Scalar>
00411 std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
00412 {
00413 Scalar r,d;
00414 if (internal::abs(yr) > internal::abs(yi))
00415 {
00416 r = yi/yr;
00417 d = yr + r*yi;
00418 return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
00419 }
00420 else
00421 {
00422 r = yr/yi;
00423 d = yi + r*yr;
00424 return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
00425 }
00426 }
00427
00428
00429 template<typename MatrixType>
00430 void EigenSolver<MatrixType>::doComputeEigenvectors()
00431 {
00432 const Index size = m_eivec.cols();
00433 const Scalar eps = NumTraits<Scalar>::epsilon();
00434
00435
00436 Scalar norm(0);
00437 for (Index j = 0; j < size; ++j)
00438 {
00439 norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
00440 }
00441
00442
00443 if (norm == 0.0)
00444 {
00445 return;
00446 }
00447
00448 for (Index n = size-1; n >= 0; n--)
00449 {
00450 Scalar p = m_eivalues.coeff(n).real();
00451 Scalar q = m_eivalues.coeff(n).imag();
00452
00453
00454 if (q == Scalar(0))
00455 {
00456 Scalar lastr(0), lastw(0);
00457 Index l = n;
00458
00459 m_matT.coeffRef(n,n) = 1.0;
00460 for (Index i = n-1; i >= 0; i--)
00461 {
00462 Scalar w = m_matT.coeff(i,i) - p;
00463 Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
00464
00465 if (m_eivalues.coeff(i).imag() < 0.0)
00466 {
00467 lastw = w;
00468 lastr = r;
00469 }
00470 else
00471 {
00472 l = i;
00473 if (m_eivalues.coeff(i).imag() == 0.0)
00474 {
00475 if (w != 0.0)
00476 m_matT.coeffRef(i,n) = -r / w;
00477 else
00478 m_matT.coeffRef(i,n) = -r / (eps * norm);
00479 }
00480 else
00481 {
00482 Scalar x = m_matT.coeff(i,i+1);
00483 Scalar y = m_matT.coeff(i+1,i);
00484 Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
00485 Scalar t = (x * lastr - lastw * r) / denom;
00486 m_matT.coeffRef(i,n) = t;
00487 if (internal::abs(x) > internal::abs(lastw))
00488 m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
00489 else
00490 m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
00491 }
00492
00493
00494 Scalar t = internal::abs(m_matT.coeff(i,n));
00495 if ((eps * t) * t > Scalar(1))
00496 m_matT.col(n).tail(size-i) /= t;
00497 }
00498 }
00499 }
00500 else if (q < Scalar(0) && n > 0)
00501 {
00502 Scalar lastra(0), lastsa(0), lastw(0);
00503 Index l = n-1;
00504
00505
00506 if (internal::abs(m_matT.coeff(n,n-1)) > internal::abs(m_matT.coeff(n-1,n)))
00507 {
00508 m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
00509 m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
00510 }
00511 else
00512 {
00513 std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
00514 m_matT.coeffRef(n-1,n-1) = internal::real(cc);
00515 m_matT.coeffRef(n-1,n) = internal::imag(cc);
00516 }
00517 m_matT.coeffRef(n,n-1) = 0.0;
00518 m_matT.coeffRef(n,n) = 1.0;
00519 for (Index i = n-2; i >= 0; i--)
00520 {
00521 Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
00522 Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
00523 Scalar w = m_matT.coeff(i,i) - p;
00524
00525 if (m_eivalues.coeff(i).imag() < 0.0)
00526 {
00527 lastw = w;
00528 lastra = ra;
00529 lastsa = sa;
00530 }
00531 else
00532 {
00533 l = i;
00534 if (m_eivalues.coeff(i).imag() == RealScalar(0))
00535 {
00536 std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
00537 m_matT.coeffRef(i,n-1) = internal::real(cc);
00538 m_matT.coeffRef(i,n) = internal::imag(cc);
00539 }
00540 else
00541 {
00542
00543 Scalar x = m_matT.coeff(i,i+1);
00544 Scalar y = m_matT.coeff(i+1,i);
00545 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;
00546 Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
00547 if ((vr == 0.0) && (vi == 0.0))
00548 vr = eps * norm * (internal::abs(w) + internal::abs(q) + internal::abs(x) + internal::abs(y) + internal::abs(lastw));
00549
00550 std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
00551 m_matT.coeffRef(i,n-1) = internal::real(cc);
00552 m_matT.coeffRef(i,n) = internal::imag(cc);
00553 if (internal::abs(x) > (internal::abs(lastw) + internal::abs(q)))
00554 {
00555 m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
00556 m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
00557 }
00558 else
00559 {
00560 cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
00561 m_matT.coeffRef(i+1,n-1) = internal::real(cc);
00562 m_matT.coeffRef(i+1,n) = internal::imag(cc);
00563 }
00564 }
00565
00566
00567 using std::max;
00568 Scalar t = (max)(internal::abs(m_matT.coeff(i,n-1)),internal::abs(m_matT.coeff(i,n)));
00569 if ((eps * t) * t > Scalar(1))
00570 m_matT.block(i, n-1, size-i, 2) /= t;
00571
00572 }
00573 }
00574
00575
00576 n--;
00577 }
00578 else
00579 {
00580 eigen_assert(0 && "Internal bug in EigenSolver");
00581 }
00582 }
00583
00584
00585 for (Index j = size-1; j >= 0; j--)
00586 {
00587 m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
00588 m_eivec.col(j) = m_tmp;
00589 }
00590 }
00591
00592 }
00593
00594 #endif // EIGEN_EIGENSOLVER_H