$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_EIGENSOLVER_H 00012 #define EIGEN_EIGENSOLVER_H 00013 00014 #include "./RealSchur.h" 00015 00016 namespace Eigen { 00017 00064 template<typename _MatrixType> class EigenSolver 00065 { 00066 public: 00067 00069 typedef _MatrixType MatrixType; 00070 00071 enum { 00072 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 00073 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00074 Options = MatrixType::Options, 00075 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00076 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00077 }; 00078 00080 typedef typename MatrixType::Scalar Scalar; 00081 typedef typename NumTraits<Scalar>::Real RealScalar; 00082 typedef typename MatrixType::Index Index; 00083 00090 typedef std::complex<RealScalar> ComplexScalar; 00091 00097 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType; 00098 00104 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType; 00105 00113 EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {} 00114 00121 EigenSolver(Index size) 00122 : m_eivec(size, size), 00123 m_eivalues(size), 00124 m_isInitialized(false), 00125 m_eigenvectorsOk(false), 00126 m_realSchur(size), 00127 m_matT(size, size), 00128 m_tmp(size) 00129 {} 00130 00146 EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) 00147 : m_eivec(matrix.rows(), matrix.cols()), 00148 m_eivalues(matrix.cols()), 00149 m_isInitialized(false), 00150 m_eigenvectorsOk(false), 00151 m_realSchur(matrix.cols()), 00152 m_matT(matrix.rows(), matrix.cols()), 00153 m_tmp(matrix.cols()) 00154 { 00155 compute(matrix, computeEigenvectors); 00156 } 00157 00178 EigenvectorsType eigenvectors() const; 00179 00198 const MatrixType& pseudoEigenvectors() const 00199 { 00200 eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 00201 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 00202 return m_eivec; 00203 } 00204 00223 MatrixType pseudoEigenvalueMatrix() const; 00224 00243 const EigenvalueType& eigenvalues() const 00244 { 00245 eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 00246 return m_eivalues; 00247 } 00248 00276 EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); 00277 00278 ComputationInfo info() const 00279 { 00280 eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 00281 return m_realSchur.info(); 00282 } 00283 00285 EigenSolver& setMaxIterations(Index maxIters) 00286 { 00287 m_realSchur.setMaxIterations(maxIters); 00288 return *this; 00289 } 00290 00292 Index getMaxIterations() 00293 { 00294 return m_realSchur.getMaxIterations(); 00295 } 00296 00297 private: 00298 void doComputeEigenvectors(); 00299 00300 protected: 00301 00302 static void check_template_parameters() 00303 { 00304 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00305 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL); 00306 } 00307 00308 MatrixType m_eivec; 00309 EigenvalueType m_eivalues; 00310 bool m_isInitialized; 00311 bool m_eigenvectorsOk; 00312 RealSchur<MatrixType> m_realSchur; 00313 MatrixType m_matT; 00314 00315 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; 00316 ColumnVectorType m_tmp; 00317 }; 00318 00319 template<typename MatrixType> 00320 MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const 00321 { 00322 eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 00323 Index n = m_eivalues.rows(); 00324 MatrixType matD = MatrixType::Zero(n,n); 00325 for (Index i=0; i<n; ++i) 00326 { 00327 if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)))) 00328 matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i)); 00329 else 00330 { 00331 matD.template block<2,2>(i,i) << numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)), 00332 -numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)); 00333 ++i; 00334 } 00335 } 00336 return matD; 00337 } 00338 00339 template<typename MatrixType> 00340 typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const 00341 { 00342 eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 00343 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 00344 Index n = m_eivec.cols(); 00345 EigenvectorsType matV(n,n); 00346 for (Index j=0; j<n; ++j) 00347 { 00348 if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j))) || j+1==n) 00349 { 00350 // we have a real eigen value 00351 matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>(); 00352 matV.col(j).normalize(); 00353 } 00354 else 00355 { 00356 // we have a pair of complex eigen values 00357 for (Index i=0; i<n; ++i) 00358 { 00359 matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1)); 00360 matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1)); 00361 } 00362 matV.col(j).normalize(); 00363 matV.col(j+1).normalize(); 00364 ++j; 00365 } 00366 } 00367 return matV; 00368 } 00369 00370 template<typename MatrixType> 00371 EigenSolver<MatrixType>& 00372 EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) 00373 { 00374 check_template_parameters(); 00375 00376 using std::sqrt; 00377 using std::abs; 00378 eigen_assert(matrix.cols() == matrix.rows()); 00379 00380 // Reduce to real Schur form. 00381 m_realSchur.compute(matrix, computeEigenvectors); 00382 00383 if (m_realSchur.info() == Success) 00384 { 00385 m_matT = m_realSchur.matrixT(); 00386 if (computeEigenvectors) 00387 m_eivec = m_realSchur.matrixU(); 00388 00389 // Compute eigenvalues from matT 00390 m_eivalues.resize(matrix.cols()); 00391 Index i = 0; 00392 while (i < matrix.cols()) 00393 { 00394 if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 00395 { 00396 m_eivalues.coeffRef(i) = m_matT.coeff(i, i); 00397 ++i; 00398 } 00399 else 00400 { 00401 Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1)); 00402 Scalar z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1))); 00403 m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z); 00404 m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z); 00405 i += 2; 00406 } 00407 } 00408 00409 // Compute eigenvectors. 00410 if (computeEigenvectors) 00411 doComputeEigenvectors(); 00412 } 00413 00414 m_isInitialized = true; 00415 m_eigenvectorsOk = computeEigenvectors; 00416 00417 return *this; 00418 } 00419 00420 // Complex scalar division. 00421 template<typename Scalar> 00422 std::complex<Scalar> cdiv(const Scalar& xr, const Scalar& xi, const Scalar& yr, const Scalar& yi) 00423 { 00424 using std::abs; 00425 Scalar r,d; 00426 if (abs(yr) > abs(yi)) 00427 { 00428 r = yi/yr; 00429 d = yr + r*yi; 00430 return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d); 00431 } 00432 else 00433 { 00434 r = yr/yi; 00435 d = yi + r*yr; 00436 return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d); 00437 } 00438 } 00439 00440 00441 template<typename MatrixType> 00442 void EigenSolver<MatrixType>::doComputeEigenvectors() 00443 { 00444 using std::abs; 00445 const Index size = m_eivec.cols(); 00446 const Scalar eps = NumTraits<Scalar>::epsilon(); 00447 00448 // inefficient! this is already computed in RealSchur 00449 Scalar norm(0); 00450 for (Index j = 0; j < size; ++j) 00451 { 00452 norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum(); 00453 } 00454 00455 // Backsubstitute to find vectors of upper triangular form 00456 if (norm == 0.0) 00457 { 00458 return; 00459 } 00460 00461 for (Index n = size-1; n >= 0; n--) 00462 { 00463 Scalar p = m_eivalues.coeff(n).real(); 00464 Scalar q = m_eivalues.coeff(n).imag(); 00465 00466 // Scalar vector 00467 if (q == Scalar(0)) 00468 { 00469 Scalar lastr(0), lastw(0); 00470 Index l = n; 00471 00472 m_matT.coeffRef(n,n) = 1.0; 00473 for (Index i = n-1; i >= 0; i--) 00474 { 00475 Scalar w = m_matT.coeff(i,i) - p; 00476 Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1)); 00477 00478 if (m_eivalues.coeff(i).imag() < 0.0) 00479 { 00480 lastw = w; 00481 lastr = r; 00482 } 00483 else 00484 { 00485 l = i; 00486 if (m_eivalues.coeff(i).imag() == 0.0) 00487 { 00488 if (w != 0.0) 00489 m_matT.coeffRef(i,n) = -r / w; 00490 else 00491 m_matT.coeffRef(i,n) = -r / (eps * norm); 00492 } 00493 else // Solve real equations 00494 { 00495 Scalar x = m_matT.coeff(i,i+1); 00496 Scalar y = m_matT.coeff(i+1,i); 00497 Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag(); 00498 Scalar t = (x * lastr - lastw * r) / denom; 00499 m_matT.coeffRef(i,n) = t; 00500 if (abs(x) > abs(lastw)) 00501 m_matT.coeffRef(i+1,n) = (-r - w * t) / x; 00502 else 00503 m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw; 00504 } 00505 00506 // Overflow control 00507 Scalar t = abs(m_matT.coeff(i,n)); 00508 if ((eps * t) * t > Scalar(1)) 00509 m_matT.col(n).tail(size-i) /= t; 00510 } 00511 } 00512 } 00513 else if (q < Scalar(0) && n > 0) // Complex vector 00514 { 00515 Scalar lastra(0), lastsa(0), lastw(0); 00516 Index l = n-1; 00517 00518 // Last vector component imaginary so matrix is triangular 00519 if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n))) 00520 { 00521 m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1); 00522 m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1); 00523 } 00524 else 00525 { 00526 std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q); 00527 m_matT.coeffRef(n-1,n-1) = numext::real(cc); 00528 m_matT.coeffRef(n-1,n) = numext::imag(cc); 00529 } 00530 m_matT.coeffRef(n,n-1) = 0.0; 00531 m_matT.coeffRef(n,n) = 1.0; 00532 for (Index i = n-2; i >= 0; i--) 00533 { 00534 Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1)); 00535 Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1)); 00536 Scalar w = m_matT.coeff(i,i) - p; 00537 00538 if (m_eivalues.coeff(i).imag() < 0.0) 00539 { 00540 lastw = w; 00541 lastra = ra; 00542 lastsa = sa; 00543 } 00544 else 00545 { 00546 l = i; 00547 if (m_eivalues.coeff(i).imag() == RealScalar(0)) 00548 { 00549 std::complex<Scalar> cc = cdiv(-ra,-sa,w,q); 00550 m_matT.coeffRef(i,n-1) = numext::real(cc); 00551 m_matT.coeffRef(i,n) = numext::imag(cc); 00552 } 00553 else 00554 { 00555 // Solve complex equations 00556 Scalar x = m_matT.coeff(i,i+1); 00557 Scalar y = m_matT.coeff(i+1,i); 00558 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; 00559 Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q; 00560 if ((vr == 0.0) && (vi == 0.0)) 00561 vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw)); 00562 00563 std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi); 00564 m_matT.coeffRef(i,n-1) = numext::real(cc); 00565 m_matT.coeffRef(i,n) = numext::imag(cc); 00566 if (abs(x) > (abs(lastw) + abs(q))) 00567 { 00568 m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x; 00569 m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x; 00570 } 00571 else 00572 { 00573 cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q); 00574 m_matT.coeffRef(i+1,n-1) = numext::real(cc); 00575 m_matT.coeffRef(i+1,n) = numext::imag(cc); 00576 } 00577 } 00578 00579 // Overflow control 00580 using std::max; 00581 Scalar t = (max)(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n))); 00582 if ((eps * t) * t > Scalar(1)) 00583 m_matT.block(i, n-1, size-i, 2) /= t; 00584 00585 } 00586 } 00587 00588 // We handled a pair of complex conjugate eigenvalues, so need to skip them both 00589 n--; 00590 } 00591 else 00592 { 00593 eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen 00594 } 00595 } 00596 00597 // Back transformation to get eigenvectors of original matrix 00598 for (Index j = size-1; j >= 0; j--) 00599 { 00600 m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1); 00601 m_eivec.col(j) = m_tmp; 00602 } 00603 } 00604 00605 } // end namespace Eigen 00606 00607 #endif // EIGEN_EIGENSOLVER_H