11 #ifndef EIGEN_EIGENSOLVER_H
12 #define EIGEN_EIGENSOLVER_H
14 #include "./RealSchur.h"
72 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
73 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
74 Options = MatrixType::Options,
75 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
76 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
80 typedef typename MatrixType::Scalar
Scalar;
82 typedef typename MatrixType::Index Index;
113 EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
122 : m_eivec(size, size),
124 m_isInitialized(false),
125 m_eigenvectorsOk(false),
147 : m_eivec(matrix.rows(), matrix.cols()),
148 m_eivalues(matrix.cols()),
149 m_isInitialized(false),
150 m_eigenvectorsOk(false),
151 m_realSchur(matrix.cols()),
152 m_matT(matrix.rows(), matrix.cols()),
155 compute(matrix, computeEigenvectors);
200 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
201 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
245 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
280 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
281 return m_realSchur.
info();
298 void doComputeEigenvectors();
303 bool m_isInitialized;
304 bool m_eigenvectorsOk;
309 ColumnVectorType m_tmp;
312 template<
typename MatrixType>
315 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
316 Index n = m_eivalues.rows();
318 for (Index i=0; i<n; ++i)
320 if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i))))
321 matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i));
324 matD.template block<2,2>(i,i) << numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)),
325 -numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i));
332 template<
typename MatrixType>
335 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
336 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
337 Index n = m_eivec.cols();
339 for (Index j=0; j<n; ++j)
341 if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j))) || j+1==n)
344 matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
345 matV.col(j).normalize();
350 for (Index i=0; i<n; ++i)
352 matV.coeffRef(i,j) =
ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1));
353 matV.coeffRef(i,j+1) =
ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
355 matV.col(j).normalize();
356 matV.col(j+1).normalize();
363 template<
typename MatrixType>
369 eigen_assert(matrix.cols() == matrix.rows());
372 m_realSchur.compute(matrix, computeEigenvectors);
374 if (m_realSchur.info() ==
Success)
376 m_matT = m_realSchur.matrixT();
377 if (computeEigenvectors)
378 m_eivec = m_realSchur.matrixU();
381 m_eivalues.resize(matrix.cols());
383 while (i < matrix.cols())
385 if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) ==
Scalar(0))
387 m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
392 Scalar p =
Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
393 Scalar z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
394 m_eivalues.coeffRef(i) =
ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
395 m_eivalues.coeffRef(i+1) =
ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
401 if (computeEigenvectors)
402 doComputeEigenvectors();
405 m_isInitialized =
true;
406 m_eigenvectorsOk = computeEigenvectors;
412 template<
typename Scalar>
413 std::complex<Scalar> cdiv(
const Scalar& xr,
const Scalar& xi,
const Scalar& yr,
const Scalar& yi)
417 if (abs(yr) > abs(yi))
421 return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
427 return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
432 template<
typename MatrixType>
433 void EigenSolver<MatrixType>::doComputeEigenvectors()
436 const Index size = m_eivec.cols();
437 const Scalar eps = NumTraits<Scalar>::epsilon();
441 for (Index j = 0; j < size; ++j)
443 norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
452 for (Index n = size-1; n >= 0; n--)
454 Scalar p = m_eivalues.coeff(n).real();
455 Scalar q = m_eivalues.coeff(n).imag();
460 Scalar lastr(0), lastw(0);
463 m_matT.coeffRef(n,n) = 1.0;
464 for (Index i = n-1; i >= 0; i--)
466 Scalar w = m_matT.coeff(i,i) - p;
467 Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
469 if (m_eivalues.coeff(i).imag() < 0.0)
477 if (m_eivalues.coeff(i).imag() == 0.0)
480 m_matT.coeffRef(i,n) = -r / w;
482 m_matT.coeffRef(i,n) = -r / (eps * norm);
486 Scalar x = m_matT.coeff(i,i+1);
487 Scalar y = m_matT.coeff(i+1,i);
488 Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
489 Scalar t = (x * lastr - lastw * r) / denom;
490 m_matT.coeffRef(i,n) = t;
491 if (abs(x) > abs(lastw))
492 m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
494 m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
498 Scalar t = abs(m_matT.coeff(i,n));
499 if ((eps * t) * t > Scalar(1))
500 m_matT.col(n).tail(size-i) /= t;
504 else if (q < Scalar(0) && n > 0)
506 Scalar lastra(0), lastsa(0), lastw(0);
510 if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n)))
512 m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
513 m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
517 std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
518 m_matT.coeffRef(n-1,n-1) = numext::real(cc);
519 m_matT.coeffRef(n-1,n) = numext::imag(cc);
521 m_matT.coeffRef(n,n-1) = 0.0;
522 m_matT.coeffRef(n,n) = 1.0;
523 for (Index i = n-2; i >= 0; i--)
525 Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
526 Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
527 Scalar w = m_matT.coeff(i,i) - p;
529 if (m_eivalues.coeff(i).imag() < 0.0)
538 if (m_eivalues.coeff(i).imag() == RealScalar(0))
540 std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
541 m_matT.coeffRef(i,n-1) = numext::real(cc);
542 m_matT.coeffRef(i,n) = numext::imag(cc);
547 Scalar x = m_matT.coeff(i,i+1);
548 Scalar y = m_matT.coeff(i+1,i);
549 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;
550 Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
551 if ((vr == 0.0) && (vi == 0.0))
552 vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw));
554 std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
555 m_matT.coeffRef(i,n-1) = numext::real(cc);
556 m_matT.coeffRef(i,n) = numext::imag(cc);
557 if (abs(x) > (abs(lastw) + abs(q)))
559 m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
560 m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
564 cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
565 m_matT.coeffRef(i+1,n-1) = numext::real(cc);
566 m_matT.coeffRef(i+1,n) = numext::imag(cc);
572 Scalar t = (max)(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
573 if ((eps * t) * t > Scalar(1))
574 m_matT.block(i, n-1, size-i, 2) /= t;
584 eigen_assert(0 &&
"Internal bug in EigenSolver");
589 for (Index j = size-1; j >= 0; j--)
591 m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
592 m_eivec.col(j) = m_tmp;
598 #endif // EIGEN_EIGENSOLVER_H