26 #ifndef EIGEN_EIGENSOLVER_H
27 #define EIGEN_EIGENSOLVER_H
95 typedef typename MatrixType::Scalar
Scalar;
97 typedef typename MatrixType::Index
Index;
162 :
m_eivec(matrix.rows(), matrix.cols()),
167 m_matT(matrix.rows(), matrix.cols()),
170 compute(matrix, computeEigenvectors);
300 void doComputeEigenvectors();
314 template<
typename MatrixType>
317 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
318 Index n = m_eivalues.rows();
320 for (
Index i=0; i<n; ++i)
334 template<
typename MatrixType>
337 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
338 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
339 Index n = m_eivec.cols();
341 for (
Index j=0; j<n; ++j)
346 matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
347 matV.col(j).normalize();
352 for (
Index i=0; i<n; ++i)
357 matV.col(j).normalize();
358 matV.col(j+1).normalize();
365 template<
typename MatrixType>
368 assert(matrix.cols() == matrix.rows());
371 m_realSchur.
compute(matrix, computeEigenvectors);
372 if (m_realSchur.info() ==
Success)
374 m_matT = m_realSchur.matrixT();
375 if (computeEigenvectors)
376 m_eivec = m_realSchur.matrixU();
379 m_eivalues.resize(matrix.cols());
381 while (i < matrix.cols())
383 if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) ==
Scalar(0))
385 m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
390 Scalar p =
Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
392 m_eivalues.coeffRef(i) =
ComplexScalar(m_matT.coeff(i+1, i+1) +
p, z);
393 m_eivalues.coeffRef(i+1) =
ComplexScalar(m_matT.coeff(i+1, i+1) +
p, -z);
399 if (computeEigenvectors)
400 doComputeEigenvectors();
403 m_isInitialized =
true;
404 m_eigenvectorsOk = computeEigenvectors;
410 template<
typename Scalar>
411 std::complex<Scalar>
cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
418 return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
424 return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
429 template<
typename MatrixType>
430 void EigenSolver<MatrixType>::doComputeEigenvectors()
432 const Index size = m_eivec.cols();
433 const Scalar eps = NumTraits<Scalar>::epsilon();
437 for (Index j = 0; j < size; ++j)
439 norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
448 for (Index n = size-1; n >= 0; n--)
450 Scalar
p = m_eivalues.coeff(n).real();
451 Scalar
q = m_eivalues.coeff(n).imag();
456 Scalar lastr(0), lastw(0);
459 m_matT.coeffRef(n,n) = 1.0;
460 for (Index i = n-1; i >= 0; i--)
462 Scalar w = m_matT.coeff(i,i) -
p;
463 Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
465 if (m_eivalues.coeff(i).imag() < 0.0)
473 if (m_eivalues.coeff(i).imag() == 0.0)
476 m_matT.coeffRef(i,n) = -r / w;
478 m_matT.coeffRef(i,n) = -r / (eps * norm);
482 Scalar x = m_matT.coeff(i,i+1);
483 Scalar
y = m_matT.coeff(i+1,i);
484 Scalar denom = (m_eivalues.coeff(i).real() -
p) * (m_eivalues.coeff(i).real() -
p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
485 Scalar t = (x * lastr - lastw * r) / denom;
486 m_matT.coeffRef(i,n) = t;
488 m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
490 m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
495 if ((eps * t) * t > Scalar(1))
496 m_matT.col(n).tail(size-i) /= t;
500 else if (q < Scalar(0) && n > 0)
502 Scalar lastra(0), lastsa(0), lastw(0);
508 m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
509 m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) -
p) / m_matT.coeff(n,n-1);
513 std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-
p,
q);
517 m_matT.coeffRef(n,n-1) = 0.0;
518 m_matT.coeffRef(n,n) = 1.0;
519 for (Index i = n-2; i >= 0; i--)
521 Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
522 Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
523 Scalar w = m_matT.coeff(i,i) -
p;
525 if (m_eivalues.coeff(i).imag() < 0.0)
534 if (m_eivalues.coeff(i).imag() == RealScalar(0))
536 std::complex<Scalar> cc =
cdiv(-ra,-sa,w,q);
543 Scalar x = m_matT.coeff(i,i+1);
544 Scalar y = m_matT.coeff(i+1,i);
545 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;
546 Scalar vi = (m_eivalues.coeff(i).real() -
p) * Scalar(2) *
q;
547 if ((vr == 0.0) && (vi == 0.0))
550 std::complex<Scalar> cc =
cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
555 m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
556 m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
560 cc =
cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,
q);
569 if ((eps * t) * t > Scalar(1))
570 m_matT.block(i, n-1, size-i, 2) /= t;
585 for (Index j = size-1; j >= 0; j--)
587 m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
588 m_eivec.col(j) = m_tmp;
594 #endif // EIGEN_EIGENSOLVER_H