diff options
Diffstat (limited to 'Eigen/src/Eigenvalues/EigenSolver.h')
| -rw-r--r-- | Eigen/src/Eigenvalues/EigenSolver.h | 113 |
1 files changed, 64 insertions, 49 deletions
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index 20c59a7a2..f205b185d 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -79,7 +79,7 @@ template<typename _MatrixType> class EigenSolver /** \brief Scalar type for matrices of type #MatrixType. */ typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; - typedef typename MatrixType::Index Index; + typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 /** \brief Complex scalar type for #MatrixType. * @@ -110,7 +110,7 @@ template<typename _MatrixType> class EigenSolver * * \sa compute() for an example. */ - EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {} + EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {} /** \brief Default constructor with memory preallocation * @@ -118,7 +118,7 @@ template<typename _MatrixType> class EigenSolver * according to the specified problem \a size. * \sa EigenSolver() */ - EigenSolver(Index size) + explicit EigenSolver(Index size) : m_eivec(size, size), m_eivalues(size), m_isInitialized(false), @@ -143,7 +143,8 @@ template<typename _MatrixType> class EigenSolver * * \sa compute() */ - EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) + template<typename InputType> + explicit EigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true) : m_eivec(matrix.rows(), matrix.cols()), m_eivalues(matrix.cols()), m_isInitialized(false), @@ -152,7 +153,7 @@ template<typename _MatrixType> class EigenSolver m_matT(matrix.rows(), matrix.cols()), m_tmp(matrix.cols()) { - compute(matrix, computeEigenvectors); + compute(matrix.derived(), computeEigenvectors); } /** \brief Returns the eigenvectors of given matrix. @@ -273,12 +274,14 @@ template<typename _MatrixType> class EigenSolver * Example: \include EigenSolver_compute.cpp * Output: \verbinclude EigenSolver_compute.out */ - EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); + template<typename InputType> + EigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true); + /** \returns NumericalIssue if the input contains INF or NaN values or overflow occured. Returns Success otherwise. */ ComputationInfo info() const { eigen_assert(m_isInitialized && "EigenSolver is not initialized."); - return m_realSchur.info(); + return m_info; } /** \brief Sets the maximum number of iterations allowed. */ @@ -309,6 +312,7 @@ template<typename _MatrixType> class EigenSolver EigenvalueType m_eivalues; bool m_isInitialized; bool m_eigenvectorsOk; + ComputationInfo m_info; RealSchur<MatrixType> m_realSchur; MatrixType m_matT; @@ -320,11 +324,12 @@ template<typename MatrixType> MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const { eigen_assert(m_isInitialized && "EigenSolver is not initialized."); + const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon(); Index n = m_eivalues.rows(); MatrixType matD = MatrixType::Zero(n,n); for (Index i=0; i<n; ++i) { - if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)))) + if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)), precision)) matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i)); else { @@ -341,11 +346,12 @@ typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eige { eigen_assert(m_isInitialized && "EigenSolver is not initialized."); eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); + const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon(); Index n = m_eivec.cols(); EigenvectorsType matV(n,n); for (Index j=0; j<n; ++j) { - if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j))) || j+1==n) + if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j)), precision) || j+1==n) { // we have a real eigen value matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>(); @@ -368,19 +374,23 @@ typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eige } template<typename MatrixType> +template<typename InputType> EigenSolver<MatrixType>& -EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) +EigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors) { check_template_parameters(); using std::sqrt; using std::abs; + using numext::isfinite; eigen_assert(matrix.cols() == matrix.rows()); // Reduce to real Schur form. - m_realSchur.compute(matrix, computeEigenvectors); + m_realSchur.compute(matrix.derived(), computeEigenvectors); + + m_info = m_realSchur.info(); - if (m_realSchur.info() == Success) + if (m_info == Success) { m_matT = m_realSchur.matrixT(); if (computeEigenvectors) @@ -394,14 +404,40 @@ EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvect if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) { m_eivalues.coeffRef(i) = m_matT.coeff(i, i); + if(!(isfinite)(m_eivalues.coeffRef(i))) + { + m_isInitialized = true; + m_eigenvectorsOk = false; + m_info = NumericalIssue; + return *this; + } ++i; } else { Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1)); - Scalar z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1))); + Scalar z; + // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1))); + // without overflow + { + Scalar t0 = m_matT.coeff(i+1, i); + Scalar t1 = m_matT.coeff(i, i+1); + Scalar maxval = numext::maxi<Scalar>(abs(p),numext::maxi<Scalar>(abs(t0),abs(t1))); + t0 /= maxval; + t1 /= maxval; + Scalar p0 = p/maxval; + z = maxval * sqrt(abs(p0 * p0 + t0 * t1)); + } + m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z); m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z); + if(!((isfinite)(m_eivalues.coeffRef(i)) && (isfinite)(m_eivalues.coeffRef(i+1)))) + { + m_isInitialized = true; + m_eigenvectorsOk = false; + m_info = NumericalIssue; + return *this; + } i += 2; } } @@ -417,26 +453,6 @@ EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvect return *this; } -// Complex scalar division. -template<typename Scalar> -std::complex<Scalar> cdiv(const Scalar& xr, const Scalar& xi, const Scalar& yr, const Scalar& yi) -{ - using std::abs; - Scalar r,d; - if (abs(yr) > abs(yi)) - { - r = yi/yr; - d = yr + r*yi; - return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d); - } - else - { - r = yr/yi; - d = yi + r*yr; - return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d); - } -} - template<typename MatrixType> void EigenSolver<MatrixType>::doComputeEigenvectors() @@ -453,7 +469,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors() } // Backsubstitute to find vectors of upper triangular form - if (norm == 0.0) + if (norm == Scalar(0)) { return; } @@ -469,13 +485,13 @@ void EigenSolver<MatrixType>::doComputeEigenvectors() Scalar lastr(0), lastw(0); Index l = n; - m_matT.coeffRef(n,n) = 1.0; + m_matT.coeffRef(n,n) = Scalar(1); for (Index i = n-1; i >= 0; i--) { Scalar w = m_matT.coeff(i,i) - p; Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1)); - if (m_eivalues.coeff(i).imag() < 0.0) + if (m_eivalues.coeff(i).imag() < Scalar(0)) { lastw = w; lastr = r; @@ -483,9 +499,9 @@ void EigenSolver<MatrixType>::doComputeEigenvectors() else { l = i; - if (m_eivalues.coeff(i).imag() == 0.0) + if (m_eivalues.coeff(i).imag() == Scalar(0)) { - if (w != 0.0) + if (w != Scalar(0)) m_matT.coeffRef(i,n) = -r / w; else m_matT.coeffRef(i,n) = -r / (eps * norm); @@ -523,19 +539,19 @@ void EigenSolver<MatrixType>::doComputeEigenvectors() } else { - std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q); + ComplexScalar cc = ComplexScalar(Scalar(0),-m_matT.coeff(n-1,n)) / ComplexScalar(m_matT.coeff(n-1,n-1)-p,q); m_matT.coeffRef(n-1,n-1) = numext::real(cc); m_matT.coeffRef(n-1,n) = numext::imag(cc); } - m_matT.coeffRef(n,n-1) = 0.0; - m_matT.coeffRef(n,n) = 1.0; + m_matT.coeffRef(n,n-1) = Scalar(0); + m_matT.coeffRef(n,n) = Scalar(1); for (Index i = n-2; i >= 0; i--) { Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1)); Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1)); Scalar w = m_matT.coeff(i,i) - p; - if (m_eivalues.coeff(i).imag() < 0.0) + if (m_eivalues.coeff(i).imag() < Scalar(0)) { lastw = w; lastra = ra; @@ -546,7 +562,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors() l = i; if (m_eivalues.coeff(i).imag() == RealScalar(0)) { - std::complex<Scalar> cc = cdiv(-ra,-sa,w,q); + ComplexScalar cc = ComplexScalar(-ra,-sa) / ComplexScalar(w,q); m_matT.coeffRef(i,n-1) = numext::real(cc); m_matT.coeffRef(i,n) = numext::imag(cc); } @@ -557,10 +573,10 @@ void EigenSolver<MatrixType>::doComputeEigenvectors() Scalar y = m_matT.coeff(i+1,i); 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; Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q; - if ((vr == 0.0) && (vi == 0.0)) + if ((vr == Scalar(0)) && (vi == Scalar(0))) vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw)); - std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi); + ComplexScalar cc = ComplexScalar(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra) / ComplexScalar(vr,vi); m_matT.coeffRef(i,n-1) = numext::real(cc); m_matT.coeffRef(i,n) = numext::imag(cc); if (abs(x) > (abs(lastw) + abs(q))) @@ -570,15 +586,14 @@ void EigenSolver<MatrixType>::doComputeEigenvectors() } else { - cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q); + cc = ComplexScalar(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n)) / ComplexScalar(lastw,q); m_matT.coeffRef(i+1,n-1) = numext::real(cc); m_matT.coeffRef(i+1,n) = numext::imag(cc); } } // Overflow control - using std::max; - Scalar t = (max)(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n))); + Scalar t = numext::maxi<Scalar>(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n))); if ((eps * t) * t > Scalar(1)) m_matT.block(i, n-1, size-i, 2) /= t; @@ -590,7 +605,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors() } else { - eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen + eigen_assert(0 && "Internal bug in EigenSolver (INF or NaN has not been detected)"); // this should not happen } } |