diff options
Diffstat (limited to 'Eigen/src/Eigenvalues')
| -rw-r--r-- | Eigen/src/Eigenvalues/CMakeLists.txt | 6 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/ComplexEigenSolver.h | 25 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/ComplexSchur.h | 19 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h (renamed from Eigen/src/Eigenvalues/ComplexSchur_MKL.h) | 41 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/EigenSolver.h | 113 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/GeneralizedEigenSolver.h | 197 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h | 3 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/HessenbergDecomposition.h | 15 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/RealQZ.h | 48 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/RealSchur.h | 41 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/RealSchur_LAPACKE.h (renamed from Eigen/src/Eigenvalues/RealSchur_MKL.h) | 42 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 280 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h (renamed from Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h) | 42 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/Tridiagonalization.h | 39 |
14 files changed, 554 insertions, 357 deletions
diff --git a/Eigen/src/Eigenvalues/CMakeLists.txt b/Eigen/src/Eigenvalues/CMakeLists.txt deleted file mode 100644 index 193e02685..000000000 --- a/Eigen/src/Eigenvalues/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_EIGENVALUES_SRCS "*.h") - -INSTALL(FILES - ${Eigen_EIGENVALUES_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigenvalues COMPONENT Devel - ) diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h index 417c72944..dc5fae06a 100644 --- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h +++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h @@ -60,7 +60,7 @@ template<typename _MatrixType> class ComplexEigenSolver /** \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. * @@ -104,7 +104,7 @@ template<typename _MatrixType> class ComplexEigenSolver * according to the specified problem \a size. * \sa ComplexEigenSolver() */ - ComplexEigenSolver(Index size) + explicit ComplexEigenSolver(Index size) : m_eivec(size, size), m_eivalues(size), m_schur(size), @@ -122,7 +122,8 @@ template<typename _MatrixType> class ComplexEigenSolver * * This constructor calls compute() to compute the eigendecomposition. */ - ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) + template<typename InputType> + explicit ComplexEigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true) : m_eivec(matrix.rows(),matrix.cols()), m_eivalues(matrix.cols()), m_schur(matrix.rows()), @@ -130,7 +131,7 @@ template<typename _MatrixType> class ComplexEigenSolver m_eigenvectorsOk(false), m_matX(matrix.rows(),matrix.cols()) { - compute(matrix, computeEigenvectors); + compute(matrix.derived(), computeEigenvectors); } /** \brief Returns the eigenvectors of given matrix. @@ -208,7 +209,8 @@ template<typename _MatrixType> class ComplexEigenSolver * Example: \include ComplexEigenSolver_compute.cpp * Output: \verbinclude ComplexEigenSolver_compute.out */ - ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); + template<typename InputType> + ComplexEigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true); /** \brief Reports whether previous computation was successful. * @@ -248,14 +250,15 @@ template<typename _MatrixType> class ComplexEigenSolver EigenvectorType m_matX; private: - void doComputeEigenvectors(const RealScalar& matrixnorm); + void doComputeEigenvectors(RealScalar matrixnorm); void sortEigenvalues(bool computeEigenvectors); }; template<typename MatrixType> +template<typename InputType> ComplexEigenSolver<MatrixType>& -ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) +ComplexEigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors) { check_template_parameters(); @@ -264,13 +267,13 @@ ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEi // Do a complex Schur decomposition, A = U T U^* // The eigenvalues are on the diagonal of T. - m_schur.compute(matrix, computeEigenvectors); + m_schur.compute(matrix.derived(), computeEigenvectors); if(m_schur.info() == Success) { m_eivalues = m_schur.matrixT().diagonal(); if(computeEigenvectors) - doComputeEigenvectors(matrix.norm()); + doComputeEigenvectors(m_schur.matrixT().norm()); sortEigenvalues(computeEigenvectors); } @@ -281,10 +284,12 @@ ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEi template<typename MatrixType> -void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(const RealScalar& matrixnorm) +void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm) { const Index n = m_eivalues.size(); + matrixnorm = numext::maxi(matrixnorm,(std::numeric_limits<RealScalar>::min)()); + // Compute X such that T = X D X^(-1), where D is the diagonal of T. // The matrix X is unit triangular. m_matX = EigenvectorType::Zero(n, n); diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h index 89e6cade3..7f38919f7 100644 --- a/Eigen/src/Eigenvalues/ComplexSchur.h +++ b/Eigen/src/Eigenvalues/ComplexSchur.h @@ -63,7 +63,7 @@ template<typename _MatrixType> class ComplexSchur /** \brief Scalar type for matrices of type \p _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 \p _MatrixType. * @@ -91,7 +91,7 @@ template<typename _MatrixType> class ComplexSchur * * \sa compute() for an example. */ - ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) + explicit ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) : m_matT(size,size), m_matU(size,size), m_hess(size), @@ -109,7 +109,8 @@ template<typename _MatrixType> class ComplexSchur * * \sa matrixT() and matrixU() for examples. */ - ComplexSchur(const MatrixType& matrix, bool computeU = true) + template<typename InputType> + explicit ComplexSchur(const EigenBase<InputType>& matrix, bool computeU = true) : m_matT(matrix.rows(),matrix.cols()), m_matU(matrix.rows(),matrix.cols()), m_hess(matrix.rows()), @@ -117,7 +118,7 @@ template<typename _MatrixType> class ComplexSchur m_matUisUptodate(false), m_maxIters(-1) { - compute(matrix, computeU); + compute(matrix.derived(), computeU); } /** \brief Returns the unitary matrix in the Schur decomposition. @@ -186,7 +187,8 @@ template<typename _MatrixType> class ComplexSchur * * \sa compute(const MatrixType&, bool, Index) */ - ComplexSchur& compute(const MatrixType& matrix, bool computeU = true); + template<typename InputType> + ComplexSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true); /** \brief Compute Schur decomposition from a given Hessenberg matrix * \param[in] matrixH Matrix in Hessenberg form H @@ -313,14 +315,15 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu template<typename MatrixType> -ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU) +template<typename InputType> +ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU) { m_matUisUptodate = false; eigen_assert(matrix.cols() == matrix.rows()); if(matrix.cols() == 1) { - m_matT = matrix.template cast<ComplexScalar>(); + m_matT = matrix.derived().template cast<ComplexScalar>(); if(computeU) m_matU = ComplexMatrixType::Identity(1,1); m_info = Success; m_isInitialized = true; @@ -328,7 +331,7 @@ ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& ma return *this; } - internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU); + internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix.derived(), computeU); computeFromHessenberg(m_matT, m_matU, computeU); return *this; } diff --git a/Eigen/src/Eigenvalues/ComplexSchur_MKL.h b/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h index 91496ae5b..4980a3ede 100644 --- a/Eigen/src/Eigenvalues/ComplexSchur_MKL.h +++ b/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h @@ -25,27 +25,24 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ******************************************************************************** - * Content : Eigen bindings to Intel(R) MKL + * Content : Eigen bindings to LAPACKe * Complex Schur needed to complex unsymmetrical eigenvalues/eigenvectors. ******************************************************************************** */ -#ifndef EIGEN_COMPLEX_SCHUR_MKL_H -#define EIGEN_COMPLEX_SCHUR_MKL_H - -#include "Eigen/src/Core/util/MKL_support.h" +#ifndef EIGEN_COMPLEX_SCHUR_LAPACKE_H +#define EIGEN_COMPLEX_SCHUR_LAPACKE_H namespace Eigen { -/** \internal Specialization for the data types supported by MKL */ +/** \internal Specialization for the data types supported by LAPACKe */ -#define EIGEN_MKL_SCHUR_COMPLEX(EIGTYPE, MKLTYPE, MKLPREFIX, MKLPREFIX_U, EIGCOLROW, MKLCOLROW) \ -template<> inline \ +#define EIGEN_LAPACKE_SCHUR_COMPLEX(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX, LAPACKE_PREFIX_U, EIGCOLROW, LAPACKE_COLROW) \ +template<> template<typename InputType> inline \ ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >& \ -ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW>& matrix, bool computeU) \ +ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const EigenBase<InputType>& matrix, bool computeU) \ { \ typedef Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> MatrixType; \ - typedef MatrixType::Scalar Scalar; \ typedef MatrixType::RealScalar RealScalar; \ typedef std::complex<RealScalar> ComplexScalar; \ \ @@ -54,25 +51,25 @@ ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Matri m_matUisUptodate = false; \ if(matrix.cols() == 1) \ { \ - m_matT = matrix.cast<ComplexScalar>(); \ + m_matT = matrix.derived().template cast<ComplexScalar>(); \ if(computeU) m_matU = ComplexMatrixType::Identity(1,1); \ m_info = Success; \ m_isInitialized = true; \ m_matUisUptodate = computeU; \ return *this; \ } \ - lapack_int n = matrix.cols(), sdim, info; \ - lapack_int lda = matrix.outerStride(); \ - lapack_int matrix_order = MKLCOLROW; \ + lapack_int n = internal::convert_index<lapack_int>(matrix.cols()), sdim, info; \ + lapack_int matrix_order = LAPACKE_COLROW; \ char jobvs, sort='N'; \ - LAPACK_##MKLPREFIX_U##_SELECT1 select = 0; \ + LAPACK_##LAPACKE_PREFIX_U##_SELECT1 select = 0; \ jobvs = (computeU) ? 'V' : 'N'; \ m_matU.resize(n, n); \ - lapack_int ldvs = m_matU.outerStride(); \ + lapack_int ldvs = internal::convert_index<lapack_int>(m_matU.outerStride()); \ m_matT = matrix; \ + lapack_int lda = internal::convert_index<lapack_int>(m_matT.outerStride()); \ Matrix<EIGTYPE, Dynamic, Dynamic> w; \ w.resize(n, 1);\ - info = LAPACKE_##MKLPREFIX##gees( matrix_order, jobvs, sort, select, n, (MKLTYPE*)m_matT.data(), lda, &sdim, (MKLTYPE*)w.data(), (MKLTYPE*)m_matU.data(), ldvs ); \ + info = LAPACKE_##LAPACKE_PREFIX##gees( matrix_order, jobvs, sort, select, n, (LAPACKE_TYPE*)m_matT.data(), lda, &sdim, (LAPACKE_TYPE*)w.data(), (LAPACKE_TYPE*)m_matU.data(), ldvs ); \ if(info == 0) \ m_info = Success; \ else \ @@ -84,11 +81,11 @@ ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Matri \ } -EIGEN_MKL_SCHUR_COMPLEX(dcomplex, MKL_Complex16, z, Z, ColMajor, LAPACK_COL_MAJOR) -EIGEN_MKL_SCHUR_COMPLEX(scomplex, MKL_Complex8, c, C, ColMajor, LAPACK_COL_MAJOR) -EIGEN_MKL_SCHUR_COMPLEX(dcomplex, MKL_Complex16, z, Z, RowMajor, LAPACK_ROW_MAJOR) -EIGEN_MKL_SCHUR_COMPLEX(scomplex, MKL_Complex8, c, C, RowMajor, LAPACK_ROW_MAJOR) +EIGEN_LAPACKE_SCHUR_COMPLEX(dcomplex, lapack_complex_double, z, Z, ColMajor, LAPACK_COL_MAJOR) +EIGEN_LAPACKE_SCHUR_COMPLEX(scomplex, lapack_complex_float, c, C, ColMajor, LAPACK_COL_MAJOR) +EIGEN_LAPACKE_SCHUR_COMPLEX(dcomplex, lapack_complex_double, z, Z, RowMajor, LAPACK_ROW_MAJOR) +EIGEN_LAPACKE_SCHUR_COMPLEX(scomplex, lapack_complex_float, c, C, RowMajor, LAPACK_ROW_MAJOR) } // end namespace Eigen -#endif // EIGEN_COMPLEX_SCHUR_MKL_H +#endif // EIGEN_COMPLEX_SCHUR_LAPACKE_H 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 } } diff --git a/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h b/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h index 956e80d9e..36a91dffc 100644 --- a/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h +++ b/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h @@ -1,8 +1,9 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> +// Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud@inria.fr> // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> +// Copyright (C) 2016 Tobias Wood <tobias@spinicist.org.uk> // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -72,7 +73,7 @@ template<typename _MatrixType> class GeneralizedEigenSolver /** \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. * @@ -89,7 +90,7 @@ template<typename _MatrixType> class GeneralizedEigenSolver */ typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType; - /** \brief Type for vector of complex scalar values eigenvalues as returned by betas(). + /** \brief Type for vector of complex scalar values eigenvalues as returned by alphas(). * * This is a column vector with entries of type #ComplexScalar. * The length of the vector is the size of #MatrixType. @@ -114,7 +115,14 @@ template<typename _MatrixType> class GeneralizedEigenSolver * * \sa compute() for an example. */ - GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {} + GeneralizedEigenSolver() + : m_eivec(), + m_alphas(), + m_betas(), + m_valuesOkay(false), + m_vectorsOkay(false), + m_realQZ() + {} /** \brief Default constructor with memory preallocation * @@ -122,14 +130,13 @@ template<typename _MatrixType> class GeneralizedEigenSolver * according to the specified problem \a size. * \sa GeneralizedEigenSolver() */ - GeneralizedEigenSolver(Index size) + explicit GeneralizedEigenSolver(Index size) : m_eivec(size, size), m_alphas(size), m_betas(size), - m_isInitialized(false), - m_eigenvectorsOk(false), + m_valuesOkay(false), + m_vectorsOkay(false), m_realQZ(size), - m_matS(size, size), m_tmp(size) {} @@ -149,10 +156,9 @@ template<typename _MatrixType> class GeneralizedEigenSolver : m_eivec(A.rows(), A.cols()), m_alphas(A.cols()), m_betas(A.cols()), - m_isInitialized(false), - m_eigenvectorsOk(false), + m_valuesOkay(false), + m_vectorsOkay(false), m_realQZ(A.cols()), - m_matS(A.rows(), A.cols()), m_tmp(A.cols()) { compute(A, B, computeEigenvectors); @@ -160,22 +166,20 @@ template<typename _MatrixType> class GeneralizedEigenSolver /* \brief Returns the computed generalized eigenvectors. * - * \returns %Matrix whose columns are the (possibly complex) eigenvectors. + * \returns %Matrix whose columns are the (possibly complex) right eigenvectors. + * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues. * * \pre Either the constructor * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function * compute(const MatrixType&, const MatrixType& bool) has been called before, and * \p computeEigenvectors was set to true (the default). * - * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding - * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The - * eigenvectors are normalized to have (Euclidean) norm equal to one. The - * matrix returned by this function is the matrix \f$ V \f$ in the - * generalized eigendecomposition \f$ A = B V D V^{-1} \f$, if it exists. - * * \sa eigenvalues() */ -// EigenvectorsType eigenvectors() const; + EigenvectorsType eigenvectors() const { + eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated."); + return m_eivec; + } /** \brief Returns an expression of the computed generalized eigenvalues. * @@ -197,7 +201,7 @@ template<typename _MatrixType> class GeneralizedEigenSolver */ EigenvalueType eigenvalues() const { - eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); + eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); return EigenvalueType(m_alphas,m_betas); } @@ -208,7 +212,7 @@ template<typename _MatrixType> class GeneralizedEigenSolver * \sa betas(), eigenvalues() */ ComplexVectorType alphas() const { - eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); + eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); return m_alphas; } @@ -219,7 +223,7 @@ template<typename _MatrixType> class GeneralizedEigenSolver * \sa alphas(), eigenvalues() */ VectorType betas() const { - eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized."); + eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); return m_betas; } @@ -250,7 +254,7 @@ template<typename _MatrixType> class GeneralizedEigenSolver ComputationInfo info() const { - eigen_assert(m_isInitialized && "EigenSolver is not initialized."); + eigen_assert(m_valuesOkay && "EigenSolver is not initialized."); return m_realQZ.info(); } @@ -270,29 +274,14 @@ template<typename _MatrixType> class GeneralizedEigenSolver EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL); } - MatrixType m_eivec; + EigenvectorsType m_eivec; ComplexVectorType m_alphas; VectorType m_betas; - bool m_isInitialized; - bool m_eigenvectorsOk; + bool m_valuesOkay, m_vectorsOkay; RealQZ<MatrixType> m_realQZ; - MatrixType m_matS; - - typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; - ColumnVectorType m_tmp; + ComplexVectorType m_tmp; }; -//template<typename MatrixType> -//typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const -//{ -// eigen_assert(m_isInitialized && "EigenSolver is not initialized."); -// eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); -// Index n = m_eivec.cols(); -// EigenvectorsType matV(n,n); -// // TODO -// return matV; -//} - template<typename MatrixType> GeneralizedEigenSolver<MatrixType>& GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) @@ -302,46 +291,126 @@ GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixTyp using std::sqrt; using std::abs; eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows()); - + Index size = A.cols(); + m_valuesOkay = false; + m_vectorsOkay = false; // Reduce to generalized real Schur form: // A = Q S Z and B = Q T Z m_realQZ.compute(A, B, computeEigenvectors); - if (m_realQZ.info() == Success) { - m_matS = m_realQZ.matrixS(); + // Resize storage + m_alphas.resize(size); + m_betas.resize(size); if (computeEigenvectors) - m_eivec = m_realQZ.matrixZ().transpose(); - - // Compute eigenvalues from matS - m_alphas.resize(A.cols()); - m_betas.resize(A.cols()); + { + m_eivec.resize(size,size); + m_tmp.resize(size); + } + + // Aliases: + Map<VectorType> v(reinterpret_cast<Scalar*>(m_tmp.data()), size); + ComplexVectorType &cv = m_tmp; + const MatrixType &mZ = m_realQZ.matrixZ(); + const MatrixType &mS = m_realQZ.matrixS(); + const MatrixType &mT = m_realQZ.matrixT(); + Index i = 0; - while (i < A.cols()) + while (i < size) { - if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0)) + if (i == size - 1 || mS.coeff(i+1, i) == Scalar(0)) { - m_alphas.coeffRef(i) = m_matS.coeff(i, i); - m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i); + // Real eigenvalue + m_alphas.coeffRef(i) = mS.diagonal().coeff(i); + m_betas.coeffRef(i) = mT.diagonal().coeff(i); + if (computeEigenvectors) + { + v.setConstant(Scalar(0.0)); + v.coeffRef(i) = Scalar(1.0); + // For singular eigenvalues do nothing more + if(abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)()) + { + // Non-singular eigenvalue + const Scalar alpha = real(m_alphas.coeffRef(i)); + const Scalar beta = m_betas.coeffRef(i); + for (Index j = i-1; j >= 0; j--) + { + const Index st = j+1; + const Index sz = i-j; + if (j > 0 && mS.coeff(j, j-1) != Scalar(0)) + { + // 2x2 block + Matrix<Scalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( v.segment(st,sz) ); + Matrix<Scalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1); + v.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs); + j--; + } + else + { + v.coeffRef(j) = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() / (beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j)); + } + } + } + m_eivec.col(i).real().noalias() = mZ.transpose() * v; + m_eivec.col(i).real().normalize(); + m_eivec.col(i).imag().setConstant(0); + } ++i; } else { - Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1)); - Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1))); - m_alphas.coeffRef(i) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z); - m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z); - - m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i); - m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i); + // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T + // Then taking beta=T_00*T_11, we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00): + + // T = [a 0] + // [0 b] + RealScalar a = mT.diagonal().coeff(i), + b = mT.diagonal().coeff(i+1); + const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b; + + // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug. + Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal(); + + Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1)); + Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1))); + const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z); + m_alphas.coeffRef(i) = conj(alpha); + m_alphas.coeffRef(i+1) = alpha; + + if (computeEigenvectors) { + // Compute eigenvector in position (i+1) and then position (i) is just the conjugate + cv.setZero(); + cv.coeffRef(i+1) = Scalar(1.0); + // here, the "static_cast" workaound expression template issues. + cv.coeffRef(i) = -(static_cast<Scalar>(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1)) + / (static_cast<Scalar>(beta*mS.coeffRef(i,i)) - alpha*mT.coeffRef(i,i)); + for (Index j = i-1; j >= 0; j--) + { + const Index st = j+1; + const Index sz = i+1-j; + if (j > 0 && mS.coeff(j, j-1) != Scalar(0)) + { + // 2x2 block + Matrix<ComplexScalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( cv.segment(st,sz) ); + Matrix<ComplexScalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1); + cv.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs); + j--; + } else { + cv.coeffRef(j) = cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() + / (alpha*mT.coeffRef(j,j) - static_cast<Scalar>(beta*mS.coeffRef(j,j))); + } + } + m_eivec.col(i+1).noalias() = (mZ.transpose() * cv); + m_eivec.col(i+1).normalize(); + m_eivec.col(i) = m_eivec.col(i+1).conjugate(); + } i += 2; } } - } - - m_isInitialized = true; - m_eigenvectorsOk = false;//computeEigenvectors; + m_valuesOkay = true; + m_vectorsOkay = computeEigenvectors; + } return *this; } diff --git a/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h index 07bf1ea09..5f6bb8289 100644 --- a/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h @@ -50,7 +50,6 @@ class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT typedef SelfAdjointEigenSolver<_MatrixType> Base; public: - typedef typename Base::Index Index; typedef _MatrixType MatrixType; /** \brief Default constructor for fixed-size matrices. @@ -74,7 +73,7 @@ class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT * * \sa compute() for an example */ - GeneralizedSelfAdjointEigenSolver(Index size) + explicit GeneralizedSelfAdjointEigenSolver(Index size) : Base(size) {} diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h index 3db0c0106..f647f69b0 100644 --- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h +++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h @@ -71,7 +71,7 @@ template<typename _MatrixType> class HessenbergDecomposition /** \brief Scalar type for matrices of type #MatrixType. */ typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::Index Index; + typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 /** \brief Type for vector of Householder coefficients. * @@ -97,7 +97,7 @@ template<typename _MatrixType> class HessenbergDecomposition * * \sa compute() for an example. */ - HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size) + explicit HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size) : m_matrix(size,size), m_temp(size), m_isInitialized(false) @@ -115,8 +115,9 @@ template<typename _MatrixType> class HessenbergDecomposition * * \sa matrixH() for an example. */ - HessenbergDecomposition(const MatrixType& matrix) - : m_matrix(matrix), + template<typename InputType> + explicit HessenbergDecomposition(const EigenBase<InputType>& matrix) + : m_matrix(matrix.derived()), m_temp(matrix.rows()), m_isInitialized(false) { @@ -147,9 +148,10 @@ template<typename _MatrixType> class HessenbergDecomposition * Example: \include HessenbergDecomposition_compute.cpp * Output: \verbinclude HessenbergDecomposition_compute.out */ - HessenbergDecomposition& compute(const MatrixType& matrix) + template<typename InputType> + HessenbergDecomposition& compute(const EigenBase<InputType>& matrix) { - m_matrix = matrix; + m_matrix = matrix.derived(); if(matrix.rows()<2) { m_isInitialized = true; @@ -337,7 +339,6 @@ namespace internal { template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> > { - typedef typename MatrixType::Index Index; public: /** \brief Constructor. * diff --git a/Eigen/src/Eigenvalues/RealQZ.h b/Eigen/src/Eigenvalues/RealQZ.h index fba6f1d77..b3a910dd9 100644 --- a/Eigen/src/Eigenvalues/RealQZ.h +++ b/Eigen/src/Eigenvalues/RealQZ.h @@ -67,7 +67,7 @@ namespace Eigen { }; typedef typename MatrixType::Scalar Scalar; typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; - typedef typename MatrixType::Index Index; + typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType; typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; @@ -83,7 +83,7 @@ namespace Eigen { * * \sa compute() for an example. */ - RealQZ(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) : + explicit RealQZ(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) : m_S(size, size), m_T(size, size), m_Q(size, size), @@ -244,7 +244,7 @@ namespace Eigen { if (m_computeQZ) m_Q.applyOnTheRight(i-1,i,G); } - // update Q + // kill T(i,i-1) if(m_T.coeff(i,i-1)!=Scalar(0)) { G.makeGivens(m_T.coeff(i,i), m_T.coeff(i,i-1), &m_T.coeffRef(i,i)); @@ -276,7 +276,7 @@ namespace Eigen { /** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0) */ template<typename MatrixType> - inline typename MatrixType::Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu) + inline Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu) { using std::abs; Index res = iu; @@ -294,7 +294,7 @@ namespace Eigen { /** \internal Look for single small diagonal element T(res, res) for res between f and l, and return res (or f-1) */ template<typename MatrixType> - inline typename MatrixType::Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l) + inline Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l) { using std::abs; Index res = l; @@ -315,8 +315,8 @@ namespace Eigen { const Index dim=m_S.cols(); if (abs(m_S.coeff(i+1,i))==Scalar(0)) return; - Index z = findSmallDiagEntry(i,i+1); - if (z==i-1) + Index j = findSmallDiagEntry(i,i+1); + if (j==i-1) { // block of (S T^{-1}) Matrix2s STi = m_T.template block<2,2>(i,i).template triangularView<Upper>(). @@ -352,7 +352,7 @@ namespace Eigen { } else { - pushDownZero(z,i,i+1); + pushDownZero(j,i,i+1); } } @@ -552,7 +552,6 @@ namespace Eigen { m_T.coeffRef(l,l-1) = Scalar(0.0); } - template<typename MatrixType> RealQZ<MatrixType>& RealQZ<MatrixType>::compute(const MatrixType& A_in, const MatrixType& B_in, bool computeQZ) { @@ -616,6 +615,37 @@ namespace Eigen { } // check if we converged before reaching iterations limit m_info = (local_iter<m_maxIters) ? Success : NoConvergence; + + // For each non triangular 2x2 diagonal block of S, + // reduce the respective 2x2 diagonal block of T to positive diagonal form using 2x2 SVD. + // This step is not mandatory for QZ, but it does help further extraction of eigenvalues/eigenvectors, + // and is in par with Lapack/Matlab QZ. + if(m_info==Success) + { + for(Index i=0; i<dim-1; ++i) + { + if(m_S.coeff(i+1, i) != Scalar(0)) + { + JacobiRotation<Scalar> j_left, j_right; + internal::real_2x2_jacobi_svd(m_T, i, i+1, &j_left, &j_right); + + // Apply resulting Jacobi rotations + m_S.applyOnTheLeft(i,i+1,j_left); + m_S.applyOnTheRight(i,i+1,j_right); + m_T.applyOnTheLeft(i,i+1,j_left); + m_T.applyOnTheRight(i,i+1,j_right); + m_T(i+1,i) = m_T(i,i+1) = Scalar(0); + + if(m_computeQZ) { + m_Q.applyOnTheRight(i,i+1,j_left.transpose()); + m_Z.applyOnTheLeft(i,i+1,j_right.transpose()); + } + + i++; + } + } + } + return *this; } // end compute diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h index 16d387537..f5c86041d 100644 --- a/Eigen/src/Eigenvalues/RealSchur.h +++ b/Eigen/src/Eigenvalues/RealSchur.h @@ -64,7 +64,7 @@ template<typename _MatrixType> class RealSchur }; typedef typename MatrixType::Scalar Scalar; typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; - typedef typename MatrixType::Index Index; + typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType; typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; @@ -80,7 +80,7 @@ template<typename _MatrixType> class RealSchur * * \sa compute() for an example. */ - RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) + explicit RealSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) : m_matT(size, size), m_matU(size, size), m_workspaceVector(size), @@ -100,7 +100,8 @@ template<typename _MatrixType> class RealSchur * Example: \include RealSchur_RealSchur_MatrixType.cpp * Output: \verbinclude RealSchur_RealSchur_MatrixType.out */ - RealSchur(const MatrixType& matrix, bool computeU = true) + template<typename InputType> + explicit RealSchur(const EigenBase<InputType>& matrix, bool computeU = true) : m_matT(matrix.rows(),matrix.cols()), m_matU(matrix.rows(),matrix.cols()), m_workspaceVector(matrix.rows()), @@ -109,7 +110,7 @@ template<typename _MatrixType> class RealSchur m_matUisUptodate(false), m_maxIters(-1) { - compute(matrix, computeU); + compute(matrix.derived(), computeU); } /** \brief Returns the orthogonal matrix in the Schur decomposition. @@ -165,7 +166,8 @@ template<typename _MatrixType> class RealSchur * * \sa compute(const MatrixType&, bool, Index) */ - RealSchur& compute(const MatrixType& matrix, bool computeU = true); + template<typename InputType> + RealSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true); /** \brief Computes Schur decomposition of a Hessenberg matrix H = Z T Z^T * \param[in] matrixH Matrix in Hessenberg form H @@ -243,26 +245,45 @@ template<typename _MatrixType> class RealSchur template<typename MatrixType> -RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU) +template<typename InputType> +RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU) { + const Scalar considerAsZero = (std::numeric_limits<Scalar>::min)(); + eigen_assert(matrix.cols() == matrix.rows()); Index maxIters = m_maxIters; if (maxIters == -1) maxIters = m_maxIterationsPerRow * matrix.rows(); + Scalar scale = matrix.derived().cwiseAbs().maxCoeff(); + if(scale<considerAsZero) + { + m_matT.setZero(matrix.rows(),matrix.cols()); + if(computeU) + m_matU.setIdentity(matrix.rows(),matrix.cols()); + m_info = Success; + m_isInitialized = true; + m_matUisUptodate = computeU; + return *this; + } + // Step 1. Reduce to Hessenberg form - m_hess.compute(matrix); + m_hess.compute(matrix.derived()/scale); // Step 2. Reduce to real Schur form computeFromHessenberg(m_hess.matrixH(), m_hess.matrixQ(), computeU); + + m_matT *= scale; return *this; } template<typename MatrixType> template<typename HessMatrixType, typename OrthMatrixType> RealSchur<MatrixType>& RealSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU) -{ - m_matT = matrixH; +{ + using std::abs; + + m_matT = matrixH; if(computeU) m_matU = matrixQ; @@ -343,7 +364,7 @@ inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT() /** \internal Look for single small sub-diagonal element and returns its index */ template<typename MatrixType> -inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu) +inline Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu) { using std::abs; Index res = iu; diff --git a/Eigen/src/Eigenvalues/RealSchur_MKL.h b/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h index ad9736460..2c2251715 100644 --- a/Eigen/src/Eigenvalues/RealSchur_MKL.h +++ b/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h @@ -25,43 +25,37 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ******************************************************************************** - * Content : Eigen bindings to Intel(R) MKL + * Content : Eigen bindings to LAPACKe * Real Schur needed to real unsymmetrical eigenvalues/eigenvectors. ******************************************************************************** */ -#ifndef EIGEN_REAL_SCHUR_MKL_H -#define EIGEN_REAL_SCHUR_MKL_H - -#include "Eigen/src/Core/util/MKL_support.h" +#ifndef EIGEN_REAL_SCHUR_LAPACKE_H +#define EIGEN_REAL_SCHUR_LAPACKE_H namespace Eigen { -/** \internal Specialization for the data types supported by MKL */ +/** \internal Specialization for the data types supported by LAPACKe */ -#define EIGEN_MKL_SCHUR_REAL(EIGTYPE, MKLTYPE, MKLPREFIX, MKLPREFIX_U, EIGCOLROW, MKLCOLROW) \ -template<> inline \ +#define EIGEN_LAPACKE_SCHUR_REAL(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX, LAPACKE_PREFIX_U, EIGCOLROW, LAPACKE_COLROW) \ +template<> template<typename InputType> inline \ RealSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >& \ -RealSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW>& matrix, bool computeU) \ +RealSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const EigenBase<InputType>& matrix, bool computeU) \ { \ - typedef Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> MatrixType; \ - typedef MatrixType::Scalar Scalar; \ - typedef MatrixType::RealScalar RealScalar; \ -\ eigen_assert(matrix.cols() == matrix.rows()); \ \ - lapack_int n = matrix.cols(), sdim, info; \ - lapack_int lda = matrix.outerStride(); \ - lapack_int matrix_order = MKLCOLROW; \ + lapack_int n = internal::convert_index<lapack_int>(matrix.cols()), sdim, info; \ + lapack_int matrix_order = LAPACKE_COLROW; \ char jobvs, sort='N'; \ - LAPACK_##MKLPREFIX_U##_SELECT2 select = 0; \ + LAPACK_##LAPACKE_PREFIX_U##_SELECT2 select = 0; \ jobvs = (computeU) ? 'V' : 'N'; \ m_matU.resize(n, n); \ - lapack_int ldvs = m_matU.outerStride(); \ + lapack_int ldvs = internal::convert_index<lapack_int>(m_matU.outerStride()); \ m_matT = matrix; \ + lapack_int lda = internal::convert_index<lapack_int>(m_matT.outerStride()); \ Matrix<EIGTYPE, Dynamic, Dynamic> wr, wi; \ wr.resize(n, 1); wi.resize(n, 1); \ - info = LAPACKE_##MKLPREFIX##gees( matrix_order, jobvs, sort, select, n, (MKLTYPE*)m_matT.data(), lda, &sdim, (MKLTYPE*)wr.data(), (MKLTYPE*)wi.data(), (MKLTYPE*)m_matU.data(), ldvs ); \ + info = LAPACKE_##LAPACKE_PREFIX##gees( matrix_order, jobvs, sort, select, n, (LAPACKE_TYPE*)m_matT.data(), lda, &sdim, (LAPACKE_TYPE*)wr.data(), (LAPACKE_TYPE*)wi.data(), (LAPACKE_TYPE*)m_matU.data(), ldvs ); \ if(info == 0) \ m_info = Success; \ else \ @@ -73,11 +67,11 @@ RealSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Matrix<E \ } -EIGEN_MKL_SCHUR_REAL(double, double, d, D, ColMajor, LAPACK_COL_MAJOR) -EIGEN_MKL_SCHUR_REAL(float, float, s, S, ColMajor, LAPACK_COL_MAJOR) -EIGEN_MKL_SCHUR_REAL(double, double, d, D, RowMajor, LAPACK_ROW_MAJOR) -EIGEN_MKL_SCHUR_REAL(float, float, s, S, RowMajor, LAPACK_ROW_MAJOR) +EIGEN_LAPACKE_SCHUR_REAL(double, double, d, D, ColMajor, LAPACK_COL_MAJOR) +EIGEN_LAPACKE_SCHUR_REAL(float, float, s, S, ColMajor, LAPACK_COL_MAJOR) +EIGEN_LAPACKE_SCHUR_REAL(double, double, d, D, RowMajor, LAPACK_ROW_MAJOR) +EIGEN_LAPACKE_SCHUR_REAL(float, float, s, S, RowMajor, LAPACK_ROW_MAJOR) } // end namespace Eigen -#endif // EIGEN_REAL_SCHUR_MKL_H +#endif // EIGEN_REAL_SCHUR_LAPACKE_H diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index c2e76c884..9ddd553f2 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -20,6 +20,8 @@ class GeneralizedSelfAdjointEigenSolver; namespace internal { template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues; +template<typename MatrixType, typename DiagType, typename SubDiagType> +ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec); } /** \eigenvalues_module \ingroup Eigenvalues_Module @@ -79,7 +81,9 @@ template<typename _MatrixType> class SelfAdjointEigenSolver /** \brief Scalar type for matrices of type \p _MatrixType. */ typedef typename MatrixType::Scalar Scalar; - typedef typename MatrixType::Index Index; + typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 + + typedef Matrix<Scalar,Size,Size,ColMajor,MaxColsAtCompileTime,MaxColsAtCompileTime> EigenvectorsType; /** \brief Real scalar type for \p _MatrixType. * @@ -98,6 +102,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver */ typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType; typedef Tridiagonalization<MatrixType> TridiagonalizationType; + typedef typename TridiagonalizationType::SubDiagonalType SubDiagonalType; /** \brief Default constructor for fixed-size matrices. * @@ -109,6 +114,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out */ + EIGEN_DEVICE_FUNC SelfAdjointEigenSolver() : m_eivec(), m_eivalues(), @@ -128,7 +134,8 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa compute() for an example */ - SelfAdjointEigenSolver(Index size) + EIGEN_DEVICE_FUNC + explicit SelfAdjointEigenSolver(Index size) : m_eivec(size, size), m_eivalues(size), m_subdiag(size > 1 ? size - 1 : 1), @@ -150,13 +157,15 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa compute(const MatrixType&, int) */ - SelfAdjointEigenSolver(const MatrixType& matrix, int options = ComputeEigenvectors) + template<typename InputType> + EIGEN_DEVICE_FUNC + explicit SelfAdjointEigenSolver(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors) : m_eivec(matrix.rows(), matrix.cols()), m_eivalues(matrix.cols()), m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1), m_isInitialized(false) { - compute(matrix, options); + compute(matrix.derived(), options); } /** \brief Computes eigendecomposition of given matrix. @@ -189,24 +198,45 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa SelfAdjointEigenSolver(const MatrixType&, int) */ - SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors); + template<typename InputType> + EIGEN_DEVICE_FUNC + SelfAdjointEigenSolver& compute(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors); - /** \brief Computes eigendecomposition of given matrix using a direct algorithm + /** \brief Computes eigendecomposition of given matrix using a closed-form algorithm * * This is a variant of compute(const MatrixType&, int options) which * directly solves the underlying polynomial equation. * - * Currently only 3x3 matrices for which the sizes are known at compile time are supported (e.g., Matrix3d). + * Currently only 2x2 and 3x3 matrices for which the sizes are known at compile time are supported (e.g., Matrix3d). * - * This method is usually significantly faster than the QR algorithm + * This method is usually significantly faster than the QR iterative algorithm * but it might also be less accurate. It is also worth noting that * for 3x3 matrices it involves trigonometric operations which are * not necessarily available for all scalar types. + * + * For the 3x3 case, we observed the following worst case relative error regarding the eigenvalues: + * - double: 1e-8 + * - float: 1e-3 * * \sa compute(const MatrixType&, int options) */ + EIGEN_DEVICE_FUNC SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors); + /** + *\brief Computes the eigen decomposition from a tridiagonal symmetric matrix + * + * \param[in] diag The vector containing the diagonal of the matrix. + * \param[in] subdiag The subdiagonal of the matrix. + * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. + * \returns Reference to \c *this + * + * This function assumes that the matrix has been reduced to tridiagonal form. + * + * \sa compute(const MatrixType&, int) for more information + */ + SelfAdjointEigenSolver& computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options=ComputeEigenvectors); + /** \brief Returns the eigenvectors of given matrix. * * \returns A const reference to the matrix whose columns are the eigenvectors. @@ -225,7 +255,8 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa eigenvalues() */ - const MatrixType& eigenvectors() const + EIGEN_DEVICE_FUNC + const EigenvectorsType& eigenvectors() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); @@ -247,6 +278,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa eigenvectors(), MatrixBase::eigenvalues() */ + EIGEN_DEVICE_FUNC const RealVectorType& eigenvalues() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -268,9 +300,9 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * Example: \include SelfAdjointEigenSolver_operatorSqrt.cpp * Output: \verbinclude SelfAdjointEigenSolver_operatorSqrt.out * - * \sa operatorInverseSqrt(), - * \ref MatrixFunctions_Module "MatrixFunctions Module" + * \sa operatorInverseSqrt(), <a href="unsupported/group__MatrixFunctions__Module.html">MatrixFunctions Module</a> */ + EIGEN_DEVICE_FUNC MatrixType operatorSqrt() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -293,9 +325,9 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * Example: \include SelfAdjointEigenSolver_operatorInverseSqrt.cpp * Output: \verbinclude SelfAdjointEigenSolver_operatorInverseSqrt.out * - * \sa operatorSqrt(), MatrixBase::inverse(), - * \ref MatrixFunctions_Module "MatrixFunctions Module" + * \sa operatorSqrt(), MatrixBase::inverse(), <a href="unsupported/group__MatrixFunctions__Module.html">MatrixFunctions Module</a> */ + EIGEN_DEVICE_FUNC MatrixType operatorInverseSqrt() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -307,6 +339,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \returns \c Success if computation was succesful, \c NoConvergence otherwise. */ + EIGEN_DEVICE_FUNC ComputationInfo info() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -320,43 +353,13 @@ template<typename _MatrixType> class SelfAdjointEigenSolver */ static const int m_maxIterations = 30; - #ifdef EIGEN2_SUPPORT - SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors) - : m_eivec(matrix.rows(), matrix.cols()), - m_eivalues(matrix.cols()), - m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1), - m_isInitialized(false) - { - compute(matrix, computeEigenvectors); - } - - SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) - : m_eivec(matA.cols(), matA.cols()), - m_eivalues(matA.cols()), - m_subdiag(matA.cols() > 1 ? matA.cols() - 1 : 1), - m_isInitialized(false) - { - static_cast<GeneralizedSelfAdjointEigenSolver<MatrixType>*>(this)->compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly); - } - - void compute(const MatrixType& matrix, bool computeEigenvectors) - { - compute(matrix, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly); - } - - void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) - { - compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly); - } - #endif // EIGEN2_SUPPORT - protected: static void check_template_parameters() { EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); } - MatrixType m_eivec; + EigenvectorsType m_eivec; RealVectorType m_eivalues; typename TridiagonalizationType::SubDiagonalType m_subdiag; ComputationInfo m_info; @@ -364,6 +367,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver bool m_eigenvectorsOk; }; +namespace internal { /** \internal * * \eigenvalues_module \ingroup Eigenvalues_Module @@ -371,8 +375,12 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * Performs a QR step on a tridiagonal symmetric matrix represented as a * pair of two vectors \a diag and \a subdiag. * - * \param matA the input selfadjoint matrix - * \param hCoeffs returned Householder coefficients + * \param diag the diagonal part of the input selfadjoint tridiagonal matrix + * \param subdiag the sub-diagonal part of the input selfadjoint tridiagonal matrix + * \param start starting index of the submatrix to work on + * \param end last+1 index of the submatrix to work on + * \param matrixQ pointer to the column-major matrix holding the eigenvectors, can be 0 + * \param n size of the input matrix * * For compilation efficiency reasons, this procedure does not use eigen expression * for its arguments. @@ -380,17 +388,21 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * Implemented from Golub's "Matrix Computations", algorithm 8.3.2: * "implicit symmetric QR step with Wilkinson shift" */ -namespace internal { template<int StorageOrder,typename RealScalar, typename Scalar, typename Index> +EIGEN_DEVICE_FUNC static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n); } template<typename MatrixType> +template<typename InputType> +EIGEN_DEVICE_FUNC SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> -::compute(const MatrixType& matrix, int options) +::compute(const EigenBase<InputType>& a_matrix, int options) { check_template_parameters(); + const InputType &matrix(a_matrix.derived()); + using std::abs; eigen_assert(matrix.cols() == matrix.rows()); eigen_assert((options&~(EigVecMask|GenEigMask))==0 @@ -402,7 +414,8 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> if(n==1) { - m_eivalues.coeffRef(0,0) = numext::real(matrix.coeff(0,0)); + m_eivec = matrix; + m_eivalues.coeffRef(0,0) = numext::real(m_eivec.coeff(0,0)); if(computeEigenvectors) m_eivec.setOnes(n,n); m_info = Success; @@ -413,7 +426,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> // declare some aliases RealVectorType& diag = m_eivalues; - MatrixType& mat = m_eivec; + EigenvectorsType& mat = m_eivec; // map the matrix coefficients to [-1:1] to avoid over- and underflow. mat = matrix.template triangularView<Lower>(); @@ -422,19 +435,74 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> mat.template triangularView<Lower>() /= scale; m_subdiag.resize(n-1); internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors); + + m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec); + // scale back the eigen values + m_eivalues *= scale; + + m_isInitialized = true; + m_eigenvectorsOk = computeEigenvectors; + return *this; +} + +template<typename MatrixType> +SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> +::computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options) +{ + //TODO : Add an option to scale the values beforehand + bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; + + m_eivalues = diag; + m_subdiag = subdiag; + if (computeEigenvectors) + { + m_eivec.setIdentity(diag.size(), diag.size()); + } + m_info = internal::computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec); + + m_isInitialized = true; + m_eigenvectorsOk = computeEigenvectors; + return *this; +} + +namespace internal { +/** + * \internal + * \brief Compute the eigendecomposition from a tridiagonal matrix + * + * \param[in,out] diag : On input, the diagonal of the matrix, on output the eigenvalues + * \param[in,out] subdiag : The subdiagonal part of the matrix (entries are modified during the decomposition) + * \param[in] maxIterations : the maximum number of iterations + * \param[in] computeEigenvectors : whether the eigenvectors have to be computed or not + * \param[out] eivec : The matrix to store the eigenvectors if computeEigenvectors==true. Must be allocated on input. + * \returns \c Success or \c NoConvergence + */ +template<typename MatrixType, typename DiagType, typename SubDiagType> +ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec) +{ + using std::abs; + + ComputationInfo info; + typedef typename MatrixType::Scalar Scalar; + + Index n = diag.size(); Index end = n-1; Index start = 0; Index iter = 0; // total number of iterations - + + typedef typename DiagType::RealScalar RealScalar; + const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)(); + const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon(); + while (end>0) { for (Index i = start; i<end; ++i) - if (internal::isMuchSmallerThan(abs(m_subdiag[i]),(abs(diag[i])+abs(diag[i+1])))) - m_subdiag[i] = 0; + if (internal::isMuchSmallerThan(abs(subdiag[i]),(abs(diag[i])+abs(diag[i+1])),precision) || abs(subdiag[i]) <= considerAsZero) + subdiag[i] = 0; // find the largest unreduced block - while (end>0 && m_subdiag[end-1]==0) + while (end>0 && subdiag[end-1]==RealScalar(0)) { end--; } @@ -443,51 +511,42 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> // if we spent too many iterations, we give up iter++; - if(iter > m_maxIterations * n) break; + if(iter > maxIterations * n) break; start = end - 1; - while (start>0 && m_subdiag[start-1]!=0) + while (start>0 && subdiag[start-1]!=0) start--; - internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n); + internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n); } - - if (iter <= m_maxIterations * n) - m_info = Success; + if (iter <= maxIterations * n) + info = Success; else - m_info = NoConvergence; + info = NoConvergence; // Sort eigenvalues and corresponding vectors. // TODO make the sort optional ? // TODO use a better sort algorithm !! - if (m_info == Success) + if (info == Success) { for (Index i = 0; i < n-1; ++i) { Index k; - m_eivalues.segment(i,n-i).minCoeff(&k); + diag.segment(i,n-i).minCoeff(&k); if (k > 0) { - std::swap(m_eivalues[i], m_eivalues[k+i]); + std::swap(diag[i], diag[k+i]); if(computeEigenvectors) - m_eivec.col(i).swap(m_eivec.col(k+i)); + eivec.col(i).swap(eivec.col(k+i)); } } } - - // scale back the eigen values - m_eivalues *= scale; - - m_isInitialized = true; - m_eigenvectorsOk = computeEigenvectors; - return *this; + return info; } - - -namespace internal { template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues { + EIGEN_DEVICE_FUNC static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options) { eig.compute(A,options); } }; @@ -497,20 +556,22 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 typedef typename SolverType::MatrixType MatrixType; typedef typename SolverType::RealVectorType VectorType; typedef typename SolverType::Scalar Scalar; - typedef typename MatrixType::Index Index; + typedef typename SolverType::EigenvectorsType EigenvectorsType; + /** \internal * Computes the roots of the characteristic polynomial of \a m. * For numerical stability m.trace() should be near zero and to avoid over- or underflow m should be normalized. */ + EIGEN_DEVICE_FUNC static inline void computeRoots(const MatrixType& m, VectorType& roots) { - using std::sqrt; - using std::atan2; - using std::cos; - using std::sin; - const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0); - const Scalar s_sqrt3 = sqrt(Scalar(3.0)); + EIGEN_USING_STD_MATH(sqrt) + EIGEN_USING_STD_MATH(atan2) + EIGEN_USING_STD_MATH(cos) + EIGEN_USING_STD_MATH(sin) + const Scalar s_inv3 = Scalar(1)/Scalar(3); + const Scalar s_sqrt3 = sqrt(Scalar(3)); // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0. The // eigenvalues are the roots to this equation, all guaranteed to be @@ -523,14 +584,12 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 // and in solving the equation for the roots in closed form. Scalar c2_over_3 = c2*s_inv3; Scalar a_over_3 = (c2*c2_over_3 - c1)*s_inv3; - if(a_over_3<Scalar(0)) - a_over_3 = Scalar(0); + a_over_3 = numext::maxi(a_over_3, Scalar(0)); Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1)); Scalar q = a_over_3*a_over_3*a_over_3 - half_b*half_b; - if(q<Scalar(0)) - q = Scalar(0); + q = numext::maxi(q, Scalar(0)); // Compute the eigenvalues by solving for the roots of the polynomial. Scalar rho = sqrt(a_over_3); @@ -543,6 +602,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 roots(2) = c2_over_3 + Scalar(2)*rho*cos_theta; } + EIGEN_DEVICE_FUNC static inline bool extract_kernel(MatrixType& mat, Ref<VectorType> res, Ref<VectorType> representative) { using std::abs; @@ -562,6 +622,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 return true; } + EIGEN_DEVICE_FUNC static inline void run(SolverType& solver, const MatrixType& mat, int options) { eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows()); @@ -570,7 +631,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 && "invalid option parameter"); bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; - MatrixType& eivecs = solver.m_eivec; + EigenvectorsType& eivecs = solver.m_eivec; VectorType& eivals = solver.m_eivalues; // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow. @@ -603,7 +664,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 Index k(0), l(2); if(d0 > d1) { - std::swap(k,l); + numext::swap(k,l); d0 = d1; } @@ -647,12 +708,15 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 }; // 2x2 direct eigenvalues decomposition, code from Hauke Heibel -template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2,false> +template<typename SolverType> +struct direct_selfadjoint_eigenvalues<SolverType,2,false> { typedef typename SolverType::MatrixType MatrixType; typedef typename SolverType::RealVectorType VectorType; typedef typename SolverType::Scalar Scalar; + typedef typename SolverType::EigenvectorsType EigenvectorsType; + EIGEN_DEVICE_FUNC static inline void computeRoots(const MatrixType& m, VectorType& roots) { using std::sqrt; @@ -662,28 +726,33 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2 roots(1) = t1 + t0; } + EIGEN_DEVICE_FUNC static inline void run(SolverType& solver, const MatrixType& mat, int options) { - using std::sqrt; - using std::abs; - + EIGEN_USING_STD_MATH(sqrt); + EIGEN_USING_STD_MATH(abs); + eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows()); eigen_assert((options&~(EigVecMask|GenEigMask))==0 && (options&EigVecMask)!=EigVecMask && "invalid option parameter"); bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; - MatrixType& eivecs = solver.m_eivec; + EigenvectorsType& eivecs = solver.m_eivec; VectorType& eivals = solver.m_eivalues; - // map the matrix coefficients to [-1:1] to avoid over- and underflow. - Scalar scale = mat.cwiseAbs().maxCoeff(); - scale = (std::max)(scale,Scalar(1)); - MatrixType scaledMat = mat / scale; - + // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow. + Scalar shift = mat.trace() / Scalar(2); + MatrixType scaledMat = mat; + scaledMat.coeffRef(0,1) = mat.coeff(1,0); + scaledMat.diagonal().array() -= shift; + Scalar scale = scaledMat.cwiseAbs().maxCoeff(); + if(scale > Scalar(0)) + scaledMat /= scale; + // Compute the eigenvalues computeRoots(scaledMat,eivals); - + // compute the eigen vectors if(computeEigenvectors) { @@ -711,10 +780,11 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2 eivecs.col(0) << eivecs.col(1).unitOrthogonal(); } } - + // Rescale back to the original size. eivals *= scale; - + eivals.array() += shift; + solver.m_info = Success; solver.m_isInitialized = true; solver.m_eigenvectorsOk = computeEigenvectors; @@ -724,6 +794,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2 } template<typename MatrixType> +EIGEN_DEVICE_FUNC SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> ::computeDirect(const MatrixType& matrix, int options) { @@ -733,6 +804,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> namespace internal { template<int StorageOrder,typename RealScalar, typename Scalar, typename Index> +EIGEN_DEVICE_FUNC static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n) { using std::abs; @@ -744,14 +816,14 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta // RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2)); // This explain the following, somewhat more complicated, version: RealScalar mu = diag[end]; - if(td==0) + if(td==RealScalar(0)) mu -= abs(e); else { RealScalar e2 = numext::abs2(subdiag[end-1]); RealScalar h = numext::hypot(td,e); - if(e2==0) mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h); - else mu -= e2 / (td + (td>0 ? h : -h)); + if(e2==RealScalar(0)) mu -= (e / (td + (td>RealScalar(0) ? RealScalar(1) : RealScalar(-1)))) * (e / h); + else mu -= e2 / (td + (td>RealScalar(0) ? h : -h)); } RealScalar x = diag[start] - mu; diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h index 17c0dadd2..3891cf883 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h @@ -25,38 +25,36 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ******************************************************************************** - * Content : Eigen bindings to Intel(R) MKL + * Content : Eigen bindings to LAPACKe * Self-adjoint eigenvalues/eigenvectors. ******************************************************************************** */ -#ifndef EIGEN_SAEIGENSOLVER_MKL_H -#define EIGEN_SAEIGENSOLVER_MKL_H - -#include "Eigen/src/Core/util/MKL_support.h" +#ifndef EIGEN_SAEIGENSOLVER_LAPACKE_H +#define EIGEN_SAEIGENSOLVER_LAPACKE_H namespace Eigen { -/** \internal Specialization for the data types supported by MKL */ +/** \internal Specialization for the data types supported by LAPACKe */ -#define EIGEN_MKL_EIG_SELFADJ(EIGTYPE, MKLTYPE, MKLRTYPE, MKLNAME, EIGCOLROW, MKLCOLROW ) \ -template<> inline \ +#define EIGEN_LAPACKE_EIG_SELFADJ(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_NAME, EIGCOLROW, LAPACKE_COLROW ) \ +template<> template<typename InputType> inline \ SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >& \ -SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW>& matrix, int options) \ +SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const EigenBase<InputType>& matrix, int options) \ { \ eigen_assert(matrix.cols() == matrix.rows()); \ eigen_assert((options&~(EigVecMask|GenEigMask))==0 \ && (options&EigVecMask)!=EigVecMask \ && "invalid option parameter"); \ bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; \ - lapack_int n = matrix.cols(), lda, matrix_order, info; \ + lapack_int n = internal::convert_index<lapack_int>(matrix.cols()), lda, matrix_order, info; \ m_eivalues.resize(n,1); \ m_subdiag.resize(n-1); \ m_eivec = matrix; \ \ if(n==1) \ { \ - m_eivalues.coeffRef(0,0) = numext::real(matrix.coeff(0,0)); \ + m_eivalues.coeffRef(0,0) = numext::real(m_eivec.coeff(0,0)); \ if(computeEigenvectors) m_eivec.setOnes(n,n); \ m_info = Success; \ m_isInitialized = true; \ @@ -64,12 +62,12 @@ SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(c return *this; \ } \ \ - lda = matrix.outerStride(); \ - matrix_order=MKLCOLROW; \ + lda = internal::convert_index<lapack_int>(m_eivec.outerStride()); \ + matrix_order=LAPACKE_COLROW; \ char jobz, uplo='L'/*, range='A'*/; \ jobz = computeEigenvectors ? 'V' : 'N'; \ \ - info = LAPACKE_##MKLNAME( matrix_order, jobz, uplo, n, (MKLTYPE*)m_eivec.data(), lda, (MKLRTYPE*)m_eivalues.data() ); \ + info = LAPACKE_##LAPACKE_NAME( matrix_order, jobz, uplo, n, (LAPACKE_TYPE*)m_eivec.data(), lda, (LAPACKE_RTYPE*)m_eivalues.data() ); \ m_info = (info==0) ? Success : NoConvergence; \ m_isInitialized = true; \ m_eigenvectorsOk = computeEigenvectors; \ @@ -77,15 +75,15 @@ SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(c } -EIGEN_MKL_EIG_SELFADJ(double, double, double, dsyev, ColMajor, LAPACK_COL_MAJOR) -EIGEN_MKL_EIG_SELFADJ(float, float, float, ssyev, ColMajor, LAPACK_COL_MAJOR) -EIGEN_MKL_EIG_SELFADJ(dcomplex, MKL_Complex16, double, zheev, ColMajor, LAPACK_COL_MAJOR) -EIGEN_MKL_EIG_SELFADJ(scomplex, MKL_Complex8, float, cheev, ColMajor, LAPACK_COL_MAJOR) +EIGEN_LAPACKE_EIG_SELFADJ(double, double, double, dsyev, ColMajor, LAPACK_COL_MAJOR) +EIGEN_LAPACKE_EIG_SELFADJ(float, float, float, ssyev, ColMajor, LAPACK_COL_MAJOR) +EIGEN_LAPACKE_EIG_SELFADJ(dcomplex, lapack_complex_double, double, zheev, ColMajor, LAPACK_COL_MAJOR) +EIGEN_LAPACKE_EIG_SELFADJ(scomplex, lapack_complex_float, float, cheev, ColMajor, LAPACK_COL_MAJOR) -EIGEN_MKL_EIG_SELFADJ(double, double, double, dsyev, RowMajor, LAPACK_ROW_MAJOR) -EIGEN_MKL_EIG_SELFADJ(float, float, float, ssyev, RowMajor, LAPACK_ROW_MAJOR) -EIGEN_MKL_EIG_SELFADJ(dcomplex, MKL_Complex16, double, zheev, RowMajor, LAPACK_ROW_MAJOR) -EIGEN_MKL_EIG_SELFADJ(scomplex, MKL_Complex8, float, cheev, RowMajor, LAPACK_ROW_MAJOR) +EIGEN_LAPACKE_EIG_SELFADJ(double, double, double, dsyev, RowMajor, LAPACK_ROW_MAJOR) +EIGEN_LAPACKE_EIG_SELFADJ(float, float, float, ssyev, RowMajor, LAPACK_ROW_MAJOR) +EIGEN_LAPACKE_EIG_SELFADJ(dcomplex, lapack_complex_double, double, zheev, RowMajor, LAPACK_ROW_MAJOR) +EIGEN_LAPACKE_EIG_SELFADJ(scomplex, lapack_complex_float, float, cheev, RowMajor, LAPACK_ROW_MAJOR) } // end namespace Eigen diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h index 192278d68..1d102c17b 100644 --- a/Eigen/src/Eigenvalues/Tridiagonalization.h +++ b/Eigen/src/Eigenvalues/Tridiagonalization.h @@ -18,8 +18,10 @@ namespace internal { template<typename MatrixType> struct TridiagonalizationMatrixTReturnType; template<typename MatrixType> struct traits<TridiagonalizationMatrixTReturnType<MatrixType> > + : public traits<typename MatrixType::PlainObject> { - typedef typename MatrixType::PlainObject ReturnType; + typedef typename MatrixType::PlainObject ReturnType; // FIXME shall it be a BandMatrix? + enum { Flags = 0 }; }; template<typename MatrixType, typename CoeffVectorType> @@ -67,7 +69,7 @@ template<typename _MatrixType> class Tridiagonalization 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 enum { Size = MatrixType::RowsAtCompileTime, @@ -89,10 +91,8 @@ template<typename _MatrixType> class Tridiagonalization >::type DiagonalReturnType; typedef typename internal::conditional<NumTraits<Scalar>::IsComplex, - typename internal::add_const_on_value_type<typename Diagonal< - Block<const MatrixType,SizeMinusOne,SizeMinusOne> >::RealReturnType>::type, - const Diagonal< - Block<const MatrixType,SizeMinusOne,SizeMinusOne> > + typename internal::add_const_on_value_type<typename Diagonal<const MatrixType, -1>::RealReturnType>::type, + const Diagonal<const MatrixType, -1> >::type SubDiagonalReturnType; /** \brief Return type of matrixQ() */ @@ -110,7 +110,7 @@ template<typename _MatrixType> class Tridiagonalization * * \sa compute() for an example. */ - Tridiagonalization(Index size = Size==Dynamic ? 2 : Size) + explicit Tridiagonalization(Index size = Size==Dynamic ? 2 : Size) : m_matrix(size,size), m_hCoeffs(size > 1 ? size-1 : 1), m_isInitialized(false) @@ -126,8 +126,9 @@ template<typename _MatrixType> class Tridiagonalization * Example: \include Tridiagonalization_Tridiagonalization_MatrixType.cpp * Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out */ - Tridiagonalization(const MatrixType& matrix) - : m_matrix(matrix), + template<typename InputType> + explicit Tridiagonalization(const EigenBase<InputType>& matrix) + : m_matrix(matrix.derived()), m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1), m_isInitialized(false) { @@ -152,9 +153,10 @@ template<typename _MatrixType> class Tridiagonalization * Example: \include Tridiagonalization_compute.cpp * Output: \verbinclude Tridiagonalization_compute.out */ - Tridiagonalization& compute(const MatrixType& matrix) + template<typename InputType> + Tridiagonalization& compute(const EigenBase<InputType>& matrix) { - m_matrix = matrix; + m_matrix = matrix.derived(); m_hCoeffs.resize(matrix.rows()-1, 1); internal::tridiagonalization_inplace(m_matrix, m_hCoeffs); m_isInitialized = true; @@ -305,7 +307,7 @@ typename Tridiagonalization<MatrixType>::DiagonalReturnType Tridiagonalization<MatrixType>::diagonal() const { eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); - return m_matrix.diagonal(); + return m_matrix.diagonal().real(); } template<typename MatrixType> @@ -313,8 +315,7 @@ typename Tridiagonalization<MatrixType>::SubDiagonalReturnType Tridiagonalization<MatrixType>::subDiagonal() const { eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); - Index n = m_matrix.rows(); - return Block<const MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal(); + return m_matrix.template diagonal<-1>().real(); } namespace internal { @@ -346,7 +347,6 @@ template<typename MatrixType, typename CoeffVectorType> void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs) { using numext::conj; - typedef typename MatrixType::Index Index; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; Index n = matA.rows(); @@ -367,10 +367,10 @@ void tridiagonalization_inplace(MatrixType& matA, CoeffVectorType& hCoeffs) hCoeffs.tail(n-i-1).noalias() = (matA.bottomRightCorner(remainingSize,remainingSize).template selfadjointView<Lower>() * (conj(h) * matA.col(i).tail(remainingSize))); - hCoeffs.tail(n-i-1) += (conj(h)*Scalar(-0.5)*(hCoeffs.tail(remainingSize).dot(matA.col(i).tail(remainingSize)))) * matA.col(i).tail(n-i-1); + hCoeffs.tail(n-i-1) += (conj(h)*RealScalar(-0.5)*(hCoeffs.tail(remainingSize).dot(matA.col(i).tail(remainingSize)))) * matA.col(i).tail(n-i-1); matA.bottomRightCorner(remainingSize, remainingSize).template selfadjointView<Lower>() - .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), -1); + .rankUpdate(matA.col(i).tail(remainingSize), hCoeffs.tail(remainingSize), Scalar(-1)); matA.col(i).coeffRef(i+1) = beta; hCoeffs.coeffRef(i) = h; @@ -438,7 +438,6 @@ struct tridiagonalization_inplace_selector { typedef typename Tridiagonalization<MatrixType>::CoeffVectorType CoeffVectorType; typedef typename Tridiagonalization<MatrixType>::HouseholderSequenceType HouseholderSequenceType; - typedef typename MatrixType::Index Index; template<typename DiagonalType, typename SubDiagonalType> static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) { @@ -467,9 +466,10 @@ struct tridiagonalization_inplace_selector<MatrixType,3,false> static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) { using std::sqrt; + const RealScalar tol = (std::numeric_limits<RealScalar>::min)(); diag[0] = mat(0,0); RealScalar v1norm2 = numext::abs2(mat(2,0)); - if(v1norm2 == RealScalar(0)) + if(v1norm2 <= tol) { diag[1] = mat(1,1); diag[2] = mat(2,2); @@ -526,7 +526,6 @@ struct tridiagonalization_inplace_selector<MatrixType,1,IsComplex> template<typename MatrixType> struct TridiagonalizationMatrixTReturnType : public ReturnByValue<TridiagonalizationMatrixTReturnType<MatrixType> > { - typedef typename MatrixType::Index Index; public: /** \brief Constructor. * |