From 35f7829af10c61e33dd2e2a7a015058e11a11ea0 Mon Sep 17 00:00:00 2001 From: Stanislaw Halik Date: Sat, 25 Mar 2017 14:17:07 +0100 Subject: update --- eigen/Eigen/src/Eigenvalues/CMakeLists.txt | 6 - eigen/Eigen/src/Eigenvalues/ComplexEigenSolver.h | 25 +- eigen/Eigen/src/Eigenvalues/ComplexSchur.h | 19 +- eigen/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h | 91 +++++++ eigen/Eigen/src/Eigenvalues/ComplexSchur_MKL.h | 93 ------- eigen/Eigen/src/Eigenvalues/EigenSolver.h | 113 +++++---- .../Eigen/src/Eigenvalues/GeneralizedEigenSolver.h | 217 +++++++++------- .../GeneralizedSelfAdjointEigenSolver.h | 3 +- .../src/Eigenvalues/HessenbergDecomposition.h | 15 +- .../Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h | 4 +- eigen/Eigen/src/Eigenvalues/RealQZ.h | 46 +++- eigen/Eigen/src/Eigenvalues/RealSchur.h | 41 +++- eigen/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h | 77 ++++++ eigen/Eigen/src/Eigenvalues/RealSchur_MKL.h | 79 ------ .../Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 273 +++++++++++++-------- .../Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h | 90 +++++++ .../src/Eigenvalues/SelfAdjointEigenSolver_MKL.h | 92 ------- eigen/Eigen/src/Eigenvalues/Tridiagonalization.h | 35 ++- 18 files changed, 749 insertions(+), 570 deletions(-) delete mode 100644 eigen/Eigen/src/Eigenvalues/CMakeLists.txt create mode 100644 eigen/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h delete mode 100644 eigen/Eigen/src/Eigenvalues/ComplexSchur_MKL.h create mode 100644 eigen/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h delete mode 100644 eigen/Eigen/src/Eigenvalues/RealSchur_MKL.h create mode 100644 eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h delete mode 100644 eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h (limited to 'eigen/Eigen/src/Eigenvalues') diff --git a/eigen/Eigen/src/Eigenvalues/CMakeLists.txt b/eigen/Eigen/src/Eigenvalues/CMakeLists.txt deleted file mode 100644 index 193e026..0000000 --- a/eigen/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/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/eigen/Eigen/src/Eigenvalues/ComplexEigenSolver.h index 417c729..dc5fae0 100644 --- a/eigen/Eigen/src/Eigenvalues/ComplexEigenSolver.h +++ b/eigen/Eigen/src/Eigenvalues/ComplexEigenSolver.h @@ -60,7 +60,7 @@ template class ComplexEigenSolver /** \brief Scalar type for matrices of type #MatrixType. */ typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::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 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 class ComplexEigenSolver * * This constructor calls compute() to compute the eigendecomposition. */ - ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) + template + explicit ComplexEigenSolver(const EigenBase& matrix, bool computeEigenvectors = true) : m_eivec(matrix.rows(),matrix.cols()), m_eivalues(matrix.cols()), m_schur(matrix.rows()), @@ -130,7 +131,7 @@ template 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 class ComplexEigenSolver * Example: \include ComplexEigenSolver_compute.cpp * Output: \verbinclude ComplexEigenSolver_compute.out */ - ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); + template + ComplexEigenSolver& compute(const EigenBase& matrix, bool computeEigenvectors = true); /** \brief Reports whether previous computation was successful. * @@ -248,14 +250,15 @@ template class ComplexEigenSolver EigenvectorType m_matX; private: - void doComputeEigenvectors(const RealScalar& matrixnorm); + void doComputeEigenvectors(RealScalar matrixnorm); void sortEigenvalues(bool computeEigenvectors); }; template +template ComplexEigenSolver& -ComplexEigenSolver::compute(const MatrixType& matrix, bool computeEigenvectors) +ComplexEigenSolver::compute(const EigenBase& matrix, bool computeEigenvectors) { check_template_parameters(); @@ -264,13 +267,13 @@ ComplexEigenSolver::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::compute(const MatrixType& matrix, bool computeEi template -void ComplexEigenSolver::doComputeEigenvectors(const RealScalar& matrixnorm) +void ComplexEigenSolver::doComputeEigenvectors(RealScalar matrixnorm) { const Index n = m_eivalues.size(); + matrixnorm = numext::maxi(matrixnorm,(std::numeric_limits::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/Eigen/src/Eigenvalues/ComplexSchur.h b/eigen/Eigen/src/Eigenvalues/ComplexSchur.h index 89e6cad..7f38919 100644 --- a/eigen/Eigen/src/Eigenvalues/ComplexSchur.h +++ b/eigen/Eigen/src/Eigenvalues/ComplexSchur.h @@ -63,7 +63,7 @@ template class ComplexSchur /** \brief Scalar type for matrices of type \p _MatrixType. */ typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::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 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 class ComplexSchur * * \sa matrixT() and matrixU() for examples. */ - ComplexSchur(const MatrixType& matrix, bool computeU = true) + template + explicit ComplexSchur(const EigenBase& 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 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 class ComplexSchur * * \sa compute(const MatrixType&, bool, Index) */ - ComplexSchur& compute(const MatrixType& matrix, bool computeU = true); + template + ComplexSchur& compute(const EigenBase& 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::ComplexScalar ComplexSchur::compu template -ComplexSchur& ComplexSchur::compute(const MatrixType& matrix, bool computeU) +template +ComplexSchur& ComplexSchur::compute(const EigenBase& matrix, bool computeU) { m_matUisUptodate = false; eigen_assert(matrix.cols() == matrix.rows()); if(matrix.cols() == 1) { - m_matT = matrix.template cast(); + m_matT = matrix.derived().template cast(); if(computeU) m_matU = ComplexMatrixType::Identity(1,1); m_info = Success; m_isInitialized = true; @@ -328,7 +331,7 @@ ComplexSchur& ComplexSchur::compute(const MatrixType& ma return *this; } - internal::complex_schur_reduce_to_hessenberg::IsComplex>::run(*this, matrix, computeU); + internal::complex_schur_reduce_to_hessenberg::IsComplex>::run(*this, matrix.derived(), computeU); computeFromHessenberg(m_matT, m_matU, computeU); return *this; } diff --git a/eigen/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h b/eigen/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h new file mode 100644 index 0000000..4980a3e --- /dev/null +++ b/eigen/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h @@ -0,0 +1,91 @@ +/* + Copyright (c) 2011, Intel Corporation. All rights reserved. + + Redistribution and use in source and binary forms, with or without modification, + are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of Intel Corporation nor the names of its contributors may + be used to endorse or promote products derived from this software without + specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + ******************************************************************************** + * Content : Eigen bindings to LAPACKe + * Complex Schur needed to complex unsymmetrical eigenvalues/eigenvectors. + ******************************************************************************** +*/ + +#ifndef EIGEN_COMPLEX_SCHUR_LAPACKE_H +#define EIGEN_COMPLEX_SCHUR_LAPACKE_H + +namespace Eigen { + +/** \internal Specialization for the data types supported by LAPACKe */ + +#define EIGEN_LAPACKE_SCHUR_COMPLEX(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX, LAPACKE_PREFIX_U, EIGCOLROW, LAPACKE_COLROW) \ +template<> template inline \ +ComplexSchur >& \ +ComplexSchur >::compute(const EigenBase& matrix, bool computeU) \ +{ \ + typedef Matrix MatrixType; \ + typedef MatrixType::RealScalar RealScalar; \ + typedef std::complex ComplexScalar; \ +\ + eigen_assert(matrix.cols() == matrix.rows()); \ +\ + m_matUisUptodate = false; \ + if(matrix.cols() == 1) \ + { \ + m_matT = matrix.derived().template cast(); \ + if(computeU) m_matU = ComplexMatrixType::Identity(1,1); \ + m_info = Success; \ + m_isInitialized = true; \ + m_matUisUptodate = computeU; \ + return *this; \ + } \ + lapack_int n = internal::convert_index(matrix.cols()), sdim, info; \ + lapack_int matrix_order = LAPACKE_COLROW; \ + char jobvs, sort='N'; \ + LAPACK_##LAPACKE_PREFIX_U##_SELECT1 select = 0; \ + jobvs = (computeU) ? 'V' : 'N'; \ + m_matU.resize(n, n); \ + lapack_int ldvs = internal::convert_index(m_matU.outerStride()); \ + m_matT = matrix; \ + lapack_int lda = internal::convert_index(m_matT.outerStride()); \ + Matrix w; \ + w.resize(n, 1);\ + 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 \ + m_info = NoConvergence; \ +\ + m_isInitialized = true; \ + m_matUisUptodate = computeU; \ + return *this; \ +\ +} + +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_LAPACKE_H diff --git a/eigen/Eigen/src/Eigenvalues/ComplexSchur_MKL.h b/eigen/Eigen/src/Eigenvalues/ComplexSchur_MKL.h deleted file mode 100644 index 27aed92..0000000 --- a/eigen/Eigen/src/Eigenvalues/ComplexSchur_MKL.h +++ /dev/null @@ -1,93 +0,0 @@ -/* - Copyright (c) 2011, Intel Corporation. All rights reserved. - - Redistribution and use in source and binary forms, with or without modification, - are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright notice, this - list of conditions and the following disclaimer. - * Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - * Neither the name of Intel Corporation nor the names of its contributors may - be used to endorse or promote products derived from this software without - specific prior written permission. - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND - ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED - WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE - DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR - ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES - (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; - LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON - ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS - SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - - ******************************************************************************** - * Content : Eigen bindings to Intel(R) MKL - * 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" - -namespace Eigen { - -/** \internal Specialization for the data types supported by MKL */ - -#define EIGEN_MKL_SCHUR_COMPLEX(EIGTYPE, MKLTYPE, MKLPREFIX, MKLPREFIX_U, EIGCOLROW, MKLCOLROW) \ -template<> inline \ -ComplexSchur >& \ -ComplexSchur >::compute(const Matrix& matrix, bool computeU) \ -{ \ - typedef Matrix MatrixType; \ - typedef MatrixType::RealScalar RealScalar; \ - typedef std::complex ComplexScalar; \ -\ - eigen_assert(matrix.cols() == matrix.rows()); \ -\ - m_matUisUptodate = false; \ - if(matrix.cols() == 1) \ - { \ - m_matT = matrix.cast(); \ - 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; \ - char jobvs, sort='N'; \ - LAPACK_##MKLPREFIX_U##_SELECT1 select = 0; \ - jobvs = (computeU) ? 'V' : 'N'; \ - m_matU.resize(n, n); \ - lapack_int ldvs = m_matU.outerStride(); \ - m_matT = matrix; \ - Matrix 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 ); \ - if(info == 0) \ - m_info = Success; \ - else \ - m_info = NoConvergence; \ -\ - m_isInitialized = true; \ - m_matUisUptodate = computeU; \ - return *this; \ -\ -} - -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) - -} // end namespace Eigen - -#endif // EIGEN_COMPLEX_SCHUR_MKL_H diff --git a/eigen/Eigen/src/Eigenvalues/EigenSolver.h b/eigen/Eigen/src/Eigenvalues/EigenSolver.h index 20c59a7..f205b18 100644 --- a/eigen/Eigen/src/Eigenvalues/EigenSolver.h +++ b/eigen/Eigen/src/Eigenvalues/EigenSolver.h @@ -79,7 +79,7 @@ template class EigenSolver /** \brief Scalar type for matrices of type #MatrixType. */ typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::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 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 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 class EigenSolver * * \sa compute() */ - EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) + template + explicit EigenSolver(const EigenBase& matrix, bool computeEigenvectors = true) : m_eivec(matrix.rows(), matrix.cols()), m_eivalues(matrix.cols()), m_isInitialized(false), @@ -152,7 +153,7 @@ template 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 class EigenSolver * Example: \include EigenSolver_compute.cpp * Output: \verbinclude EigenSolver_compute.out */ - EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); + template + EigenSolver& compute(const EigenBase& 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 class EigenSolver EigenvalueType m_eivalues; bool m_isInitialized; bool m_eigenvectorsOk; + ComputationInfo m_info; RealSchur m_realSchur; MatrixType m_matT; @@ -320,11 +324,12 @@ template MatrixType EigenSolver::pseudoEigenvalueMatrix() const { eigen_assert(m_isInitialized && "EigenSolver is not initialized."); + const RealScalar precision = RealScalar(2)*NumTraits::epsilon(); Index n = m_eivalues.rows(); MatrixType matD = MatrixType::Zero(n,n); for (Index i=0; i::EigenvectorsType EigenSolver::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::epsilon(); Index n = m_eivec.cols(); EigenvectorsType matV(n,n); for (Index j=0; j(); @@ -368,19 +374,23 @@ typename EigenSolver::EigenvectorsType EigenSolver::eige } template +template EigenSolver& -EigenSolver::compute(const MatrixType& matrix, bool computeEigenvectors) +EigenSolver::compute(const EigenBase& 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::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(abs(p),numext::maxi(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::compute(const MatrixType& matrix, bool computeEigenvect return *this; } -// Complex scalar division. -template -std::complex 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((xr + r*xi)/d, (xi - r*xr)/d); - } - else - { - r = yr/yi; - d = yi + r*yr; - return std::complex((r*xr + xi)/d, (r*xi - xr)/d); - } -} - template void EigenSolver::doComputeEigenvectors() @@ -453,7 +469,7 @@ void EigenSolver::doComputeEigenvectors() } // Backsubstitute to find vectors of upper triangular form - if (norm == 0.0) + if (norm == Scalar(0)) { return; } @@ -469,13 +485,13 @@ void EigenSolver::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::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::doComputeEigenvectors() } else { - std::complex cc = cdiv(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::doComputeEigenvectors() l = i; if (m_eivalues.coeff(i).imag() == RealScalar(0)) { - std::complex 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::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 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::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(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::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/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h b/eigen/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h index e5131d2..36a91df 100644 --- a/eigen/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h +++ b/eigen/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 +// Copyright (C) 2012-2016 Gael Guennebaud // Copyright (C) 2010,2012 Jitse Niesen +// Copyright (C) 2016 Tobias Wood // // 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 class GeneralizedEigenSolver /** \brief Scalar type for matrices of type #MatrixType. */ typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::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 class GeneralizedEigenSolver */ typedef Matrix 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 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 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 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 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 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 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 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 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 class GeneralizedEigenSolver EIGEN_STATIC_ASSERT(!NumTraits::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 m_realQZ; - MatrixType m_matS; - - typedef Matrix ColumnVectorType; - ColumnVectorType m_tmp; + ComplexVectorType m_tmp; }; -//template -//typename GeneralizedEigenSolver::EigenvectorsType GeneralizedEigenSolver::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 GeneralizedEigenSolver& GeneralizedEigenSolver::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) @@ -302,66 +291,126 @@ GeneralizedEigenSolver::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 v(reinterpret_cast(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::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 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 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 { - // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a triangular 2x2 block T - // From the eigen decomposition of T = U * E * U^-1, - // we can extract the eigenvalues of (U^-1 * S * U) / E - // Here, we can take advantage that E = diag(T), and U = [ 1 T_01 ; 0 T_11-T_00], and U^-1 = [1 -T_11/(T_11-T_00) ; 0 1/(T_11-T_00)]. - // Then taking beta=T_00*T_11*(T_11-T_00), we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00) * (T_11-T_00): - - // T = [a b ; 0 c] - // S = [e f ; g h] - RealScalar a = m_realQZ.matrixT().coeff(i, i), b = m_realQZ.matrixT().coeff(i, i+1), c = m_realQZ.matrixT().coeff(i+1, i+1); - RealScalar e = m_matS.coeff(i, i), f = m_matS.coeff(i, i+1), g = m_matS.coeff(i+1, i), h = m_matS.coeff(i+1, i+1); - RealScalar d = c-a; - RealScalar gb = g*b; - Matrix A; - A << (e*d-gb)*c, ((e*b+f*d-h*b)*d-gb*b)*a, - g*c , (gb+h*d)*a; - - // NOTE, we could also compute the SVD of T's block during the QZ factorization so that the respective T block is guaranteed to be diagonal, - // and then we could directly apply the formula below (while taking care of scaling S columns by T11,T00): - - Scalar p = Scalar(0.5) * (A.coeff(i, i) - A.coeff(i+1, i+1)); - Scalar z = sqrt(abs(p * p + A.coeff(i+1, i) * A.coeff(i, i+1))); - m_alphas.coeffRef(i) = ComplexScalar(A.coeff(i+1, i+1) + p, z); - m_alphas.coeffRef(i+1) = ComplexScalar(A.coeff(i+1, i+1) + p, -z); - - m_betas.coeffRef(i) = - m_betas.coeffRef(i+1) = a*c*d; - + // 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 S2 = mS.template block<2,2>(i,i) * Matrix(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(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1)) + / (static_cast(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 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 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(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/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h b/eigen/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h index 07bf1ea..5f6bb82 100644 --- a/eigen/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h +++ b/eigen/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/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/eigen/Eigen/src/Eigenvalues/HessenbergDecomposition.h index 3db0c01..f647f69 100644 --- a/eigen/Eigen/src/Eigenvalues/HessenbergDecomposition.h +++ b/eigen/Eigen/src/Eigenvalues/HessenbergDecomposition.h @@ -71,7 +71,7 @@ template 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 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 class HessenbergDecomposition * * \sa matrixH() for an example. */ - HessenbergDecomposition(const MatrixType& matrix) - : m_matrix(matrix), + template + explicit HessenbergDecomposition(const EigenBase& matrix) + : m_matrix(matrix.derived()), m_temp(matrix.rows()), m_isInitialized(false) { @@ -147,9 +148,10 @@ template class HessenbergDecomposition * Example: \include HessenbergDecomposition_compute.cpp * Output: \verbinclude HessenbergDecomposition_compute.out */ - HessenbergDecomposition& compute(const MatrixType& matrix) + template + HessenbergDecomposition& compute(const EigenBase& matrix) { - m_matrix = matrix; + m_matrix = matrix.derived(); if(matrix.rows()<2) { m_isInitialized = true; @@ -337,7 +339,6 @@ namespace internal { template struct HessenbergDecompositionMatrixHReturnType : public ReturnByValue > { - typedef typename MatrixType::Index Index; public: /** \brief Constructor. * diff --git a/eigen/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h b/eigen/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h index 4fec8af..dbbd480 100644 --- a/eigen/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h +++ b/eigen/Eigen/src/Eigenvalues/MatrixBaseEigenvalues.h @@ -85,7 +85,7 @@ MatrixBase::eigenvalues() const * \sa SelfAdjointEigenSolver::eigenvalues(), MatrixBase::eigenvalues() */ template -inline typename SelfAdjointView::EigenvaluesReturnType +EIGEN_DEVICE_FUNC inline typename SelfAdjointView::EigenvaluesReturnType SelfAdjointView::eigenvalues() const { typedef typename SelfAdjointView::PlainObject PlainObject; @@ -149,7 +149,7 @@ MatrixBase::operatorNorm() const * \sa eigenvalues(), MatrixBase::operatorNorm() */ template -inline typename SelfAdjointView::RealScalar +EIGEN_DEVICE_FUNC inline typename SelfAdjointView::RealScalar SelfAdjointView::operatorNorm() const { return eigenvalues().cwiseAbs().maxCoeff(); diff --git a/eigen/Eigen/src/Eigenvalues/RealQZ.h b/eigen/Eigen/src/Eigenvalues/RealQZ.h index aa3833e..b3a910d 100644 --- a/eigen/Eigen/src/Eigenvalues/RealQZ.h +++ b/eigen/Eigen/src/Eigenvalues/RealQZ.h @@ -67,7 +67,7 @@ namespace Eigen { }; typedef typename MatrixType::Scalar Scalar; typedef std::complex::Real> ComplexScalar; - typedef typename MatrixType::Index Index; + typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 typedef Matrix EigenvalueType; typedef Matrix 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), @@ -276,7 +276,7 @@ namespace Eigen { /** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0) */ template - inline typename MatrixType::Index RealQZ::findSmallSubdiagEntry(Index iu) + inline Index RealQZ::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 - inline typename MatrixType::Index RealQZ::findSmallDiagEntry(Index f, Index l) + inline Index RealQZ::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(). @@ -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 RealQZ& RealQZ::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 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/Eigen/src/Eigenvalues/RealSchur.h b/eigen/Eigen/src/Eigenvalues/RealSchur.h index 16d3875..f5c8604 100644 --- a/eigen/Eigen/src/Eigenvalues/RealSchur.h +++ b/eigen/Eigen/src/Eigenvalues/RealSchur.h @@ -64,7 +64,7 @@ template class RealSchur }; typedef typename MatrixType::Scalar Scalar; typedef std::complex::Real> ComplexScalar; - typedef typename MatrixType::Index Index; + typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 typedef Matrix EigenvalueType; typedef Matrix ColumnVectorType; @@ -80,7 +80,7 @@ template 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 class RealSchur * Example: \include RealSchur_RealSchur_MatrixType.cpp * Output: \verbinclude RealSchur_RealSchur_MatrixType.out */ - RealSchur(const MatrixType& matrix, bool computeU = true) + template + explicit RealSchur(const EigenBase& 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 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 class RealSchur * * \sa compute(const MatrixType&, bool, Index) */ - RealSchur& compute(const MatrixType& matrix, bool computeU = true); + template + RealSchur& compute(const EigenBase& 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 class RealSchur template -RealSchur& RealSchur::compute(const MatrixType& matrix, bool computeU) +template +RealSchur& RealSchur::compute(const EigenBase& matrix, bool computeU) { + const Scalar considerAsZero = (std::numeric_limits::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 template RealSchur& RealSchur::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::computeNormOfT() /** \internal Look for single small sub-diagonal element and returns its index */ template -inline typename MatrixType::Index RealSchur::findSmallSubdiagEntry(Index iu) +inline Index RealSchur::findSmallSubdiagEntry(Index iu) { using std::abs; Index res = iu; diff --git a/eigen/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h b/eigen/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h new file mode 100644 index 0000000..2c22517 --- /dev/null +++ b/eigen/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h @@ -0,0 +1,77 @@ +/* + Copyright (c) 2011, Intel Corporation. All rights reserved. + + Redistribution and use in source and binary forms, with or without modification, + are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of Intel Corporation nor the names of its contributors may + be used to endorse or promote products derived from this software without + specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + ******************************************************************************** + * Content : Eigen bindings to LAPACKe + * Real Schur needed to real unsymmetrical eigenvalues/eigenvectors. + ******************************************************************************** +*/ + +#ifndef EIGEN_REAL_SCHUR_LAPACKE_H +#define EIGEN_REAL_SCHUR_LAPACKE_H + +namespace Eigen { + +/** \internal Specialization for the data types supported by LAPACKe */ + +#define EIGEN_LAPACKE_SCHUR_REAL(EIGTYPE, LAPACKE_TYPE, LAPACKE_PREFIX, LAPACKE_PREFIX_U, EIGCOLROW, LAPACKE_COLROW) \ +template<> template inline \ +RealSchur >& \ +RealSchur >::compute(const EigenBase& matrix, bool computeU) \ +{ \ + eigen_assert(matrix.cols() == matrix.rows()); \ +\ + lapack_int n = internal::convert_index(matrix.cols()), sdim, info; \ + lapack_int matrix_order = LAPACKE_COLROW; \ + char jobvs, sort='N'; \ + LAPACK_##LAPACKE_PREFIX_U##_SELECT2 select = 0; \ + jobvs = (computeU) ? 'V' : 'N'; \ + m_matU.resize(n, n); \ + lapack_int ldvs = internal::convert_index(m_matU.outerStride()); \ + m_matT = matrix; \ + lapack_int lda = internal::convert_index(m_matT.outerStride()); \ + Matrix wr, wi; \ + wr.resize(n, 1); wi.resize(n, 1); \ + 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 \ + m_info = NoConvergence; \ +\ + m_isInitialized = true; \ + m_matUisUptodate = computeU; \ + return *this; \ +\ +} + +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_LAPACKE_H diff --git a/eigen/Eigen/src/Eigenvalues/RealSchur_MKL.h b/eigen/Eigen/src/Eigenvalues/RealSchur_MKL.h deleted file mode 100644 index c3089b4..0000000 --- a/eigen/Eigen/src/Eigenvalues/RealSchur_MKL.h +++ /dev/null @@ -1,79 +0,0 @@ -/* - Copyright (c) 2011, Intel Corporation. All rights reserved. - - Redistribution and use in source and binary forms, with or without modification, - are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright notice, this - list of conditions and the following disclaimer. - * Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - * Neither the name of Intel Corporation nor the names of its contributors may - be used to endorse or promote products derived from this software without - specific prior written permission. - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND - ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED - WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE - DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR - ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES - (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; - LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON - ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS - SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - - ******************************************************************************** - * Content : Eigen bindings to Intel(R) MKL - * 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" - -namespace Eigen { - -/** \internal Specialization for the data types supported by MKL */ - -#define EIGEN_MKL_SCHUR_REAL(EIGTYPE, MKLTYPE, MKLPREFIX, MKLPREFIX_U, EIGCOLROW, MKLCOLROW) \ -template<> inline \ -RealSchur >& \ -RealSchur >::compute(const Matrix& matrix, bool computeU) \ -{ \ - eigen_assert(matrix.cols() == matrix.rows()); \ -\ - lapack_int n = matrix.cols(), sdim, info; \ - lapack_int lda = matrix.outerStride(); \ - lapack_int matrix_order = MKLCOLROW; \ - char jobvs, sort='N'; \ - LAPACK_##MKLPREFIX_U##_SELECT2 select = 0; \ - jobvs = (computeU) ? 'V' : 'N'; \ - m_matU.resize(n, n); \ - lapack_int ldvs = m_matU.outerStride(); \ - m_matT = matrix; \ - Matrix 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 ); \ - if(info == 0) \ - m_info = Success; \ - else \ - m_info = NoConvergence; \ -\ - m_isInitialized = true; \ - m_matUisUptodate = computeU; \ - return *this; \ -\ -} - -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) - -} // end namespace Eigen - -#endif // EIGEN_REAL_SCHUR_MKL_H diff --git a/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 1131c8a..9ddd553 100644 --- a/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -20,6 +20,8 @@ class GeneralizedSelfAdjointEigenSolver; namespace internal { template struct direct_selfadjoint_eigenvalues; +template +ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec); } /** \eigenvalues_module \ingroup Eigenvalues_Module @@ -79,7 +81,7 @@ template 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 EigenvectorsType; @@ -100,6 +102,7 @@ template class SelfAdjointEigenSolver */ typedef typename internal::plain_col_type::type RealVectorType; typedef Tridiagonalization TridiagonalizationType; + typedef typename TridiagonalizationType::SubDiagonalType SubDiagonalType; /** \brief Default constructor for fixed-size matrices. * @@ -111,6 +114,7 @@ template class SelfAdjointEigenSolver * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out */ + EIGEN_DEVICE_FUNC SelfAdjointEigenSolver() : m_eivec(), m_eivalues(), @@ -130,7 +134,8 @@ template 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), @@ -152,13 +157,15 @@ template class SelfAdjointEigenSolver * * \sa compute(const MatrixType&, int) */ - SelfAdjointEigenSolver(const MatrixType& matrix, int options = ComputeEigenvectors) + template + EIGEN_DEVICE_FUNC + explicit SelfAdjointEigenSolver(const EigenBase& 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. @@ -191,24 +198,45 @@ template class SelfAdjointEigenSolver * * \sa SelfAdjointEigenSolver(const MatrixType&, int) */ - SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors); + template + EIGEN_DEVICE_FUNC + SelfAdjointEigenSolver& compute(const EigenBase& 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. @@ -227,6 +255,7 @@ template class SelfAdjointEigenSolver * * \sa eigenvalues() */ + EIGEN_DEVICE_FUNC const EigenvectorsType& eigenvectors() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -249,6 +278,7 @@ template class SelfAdjointEigenSolver * * \sa eigenvectors(), MatrixBase::eigenvalues() */ + EIGEN_DEVICE_FUNC const RealVectorType& eigenvalues() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -270,9 +300,9 @@ template class SelfAdjointEigenSolver * Example: \include SelfAdjointEigenSolver_operatorSqrt.cpp * Output: \verbinclude SelfAdjointEigenSolver_operatorSqrt.out * - * \sa operatorInverseSqrt(), - * \ref MatrixFunctions_Module "MatrixFunctions Module" + * \sa operatorInverseSqrt(), MatrixFunctions Module */ + EIGEN_DEVICE_FUNC MatrixType operatorSqrt() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -295,9 +325,9 @@ template 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(), MatrixFunctions Module */ + EIGEN_DEVICE_FUNC MatrixType operatorInverseSqrt() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -309,6 +339,7 @@ template 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."); @@ -322,36 +353,6 @@ template 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*>(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() { @@ -366,6 +367,7 @@ template class SelfAdjointEigenSolver bool m_eigenvectorsOk; }; +namespace internal { /** \internal * * \eigenvalues_module \ingroup Eigenvalues_Module @@ -373,8 +375,12 @@ template 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. @@ -382,17 +388,21 @@ template class SelfAdjointEigenSolver * Implemented from Golub's "Matrix Computations", algorithm 8.3.2: * "implicit symmetric QR step with Wilkinson shift" */ -namespace internal { -template +template +EIGEN_DEVICE_FUNC static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n); } template +template +EIGEN_DEVICE_FUNC SelfAdjointEigenSolver& SelfAdjointEigenSolver -::compute(const MatrixType& matrix, int options) +::compute(const EigenBase& 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 @@ -404,7 +414,8 @@ SelfAdjointEigenSolver& SelfAdjointEigenSolver 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; @@ -424,19 +435,74 @@ SelfAdjointEigenSolver& SelfAdjointEigenSolver mat.template triangularView() /= 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 +SelfAdjointEigenSolver& SelfAdjointEigenSolver +::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 +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::min)(); + const RealScalar precision = RealScalar(2)*NumTraits::epsilon(); + while (end>0) { for (Index i = start; i0 && m_subdiag[end-1]==0) + while (end>0 && subdiag[end-1]==RealScalar(0)) { end--; } @@ -445,51 +511,42 @@ SelfAdjointEigenSolver& SelfAdjointEigenSolver // 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(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n); + internal::tridiagonal_qr_step(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 struct direct_selfadjoint_eigenvalues { + EIGEN_DEVICE_FUNC static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options) { eig.compute(A,options); } }; @@ -499,21 +556,22 @@ template struct direct_selfadjoint_eigenvalues struct direct_selfadjoint_eigenvalues struct direct_selfadjoint_eigenvalues res, Ref representative) { using std::abs; @@ -565,6 +622,7 @@ template struct direct_selfadjoint_eigenvalues struct direct_selfadjoint_eigenvalues d1) { - std::swap(k,l); + numext::swap(k,l); d0 = d1; } @@ -650,13 +708,15 @@ template struct direct_selfadjoint_eigenvalues struct direct_selfadjoint_eigenvalues +template +struct direct_selfadjoint_eigenvalues { 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; @@ -666,11 +726,12 @@ template struct direct_selfadjoint_eigenvalues struct direct_selfadjoint_eigenvalues Scalar(0)) + scaledMat /= scale; + // Compute the eigenvalues computeRoots(scaledMat,eivals); - + // compute the eigen vectors if(computeEigenvectors) { @@ -715,10 +780,11 @@ template struct direct_selfadjoint_eigenvalues struct direct_selfadjoint_eigenvalues +EIGEN_DEVICE_FUNC SelfAdjointEigenSolver& SelfAdjointEigenSolver ::computeDirect(const MatrixType& matrix, int options) { @@ -736,7 +803,8 @@ SelfAdjointEigenSolver& SelfAdjointEigenSolver } namespace internal { -template +template +EIGEN_DEVICE_FUNC static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n) { using std::abs; @@ -748,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; @@ -788,7 +856,8 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta // apply the givens rotation to the unit matrix Q = Q * G if (matrixQ) { - Map > q(matrixQ,n,n); + // FIXME if StorageOrder == RowMajor this operation is not very efficient + Map > q(matrixQ,n,n); q.applyOnTheRight(k,k+1,rot); } } diff --git a/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h b/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h new file mode 100644 index 0000000..3891cf8 --- /dev/null +++ b/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h @@ -0,0 +1,90 @@ +/* + Copyright (c) 2011, Intel Corporation. All rights reserved. + + Redistribution and use in source and binary forms, with or without modification, + are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of Intel Corporation nor the names of its contributors may + be used to endorse or promote products derived from this software without + specific prior written permission. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + + ******************************************************************************** + * Content : Eigen bindings to LAPACKe + * Self-adjoint eigenvalues/eigenvectors. + ******************************************************************************** +*/ + +#ifndef EIGEN_SAEIGENSOLVER_LAPACKE_H +#define EIGEN_SAEIGENSOLVER_LAPACKE_H + +namespace Eigen { + +/** \internal Specialization for the data types supported by LAPACKe */ + +#define EIGEN_LAPACKE_EIG_SELFADJ(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_NAME, EIGCOLROW, LAPACKE_COLROW ) \ +template<> template inline \ +SelfAdjointEigenSolver >& \ +SelfAdjointEigenSolver >::compute(const EigenBase& 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 = internal::convert_index(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(m_eivec.coeff(0,0)); \ + if(computeEigenvectors) m_eivec.setOnes(n,n); \ + m_info = Success; \ + m_isInitialized = true; \ + m_eigenvectorsOk = computeEigenvectors; \ + return *this; \ + } \ +\ + lda = internal::convert_index(m_eivec.outerStride()); \ + matrix_order=LAPACKE_COLROW; \ + char jobz, uplo='L'/*, range='A'*/; \ + jobz = computeEigenvectors ? 'V' : 'N'; \ +\ + 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; \ + return *this; \ +} + + +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_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 + +#endif // EIGEN_SAEIGENSOLVER_H diff --git a/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h b/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h deleted file mode 100644 index 17c0dad..0000000 --- a/eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h +++ /dev/null @@ -1,92 +0,0 @@ -/* - Copyright (c) 2011, Intel Corporation. All rights reserved. - - Redistribution and use in source and binary forms, with or without modification, - are permitted provided that the following conditions are met: - - * Redistributions of source code must retain the above copyright notice, this - list of conditions and the following disclaimer. - * Redistributions in binary form must reproduce the above copyright notice, - this list of conditions and the following disclaimer in the documentation - and/or other materials provided with the distribution. - * Neither the name of Intel Corporation nor the names of its contributors may - be used to endorse or promote products derived from this software without - specific prior written permission. - - THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND - ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED - WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE - DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR - ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES - (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; - LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON - ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT - (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS - SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - - ******************************************************************************** - * Content : Eigen bindings to Intel(R) MKL - * Self-adjoint eigenvalues/eigenvectors. - ******************************************************************************** -*/ - -#ifndef EIGEN_SAEIGENSOLVER_MKL_H -#define EIGEN_SAEIGENSOLVER_MKL_H - -#include "Eigen/src/Core/util/MKL_support.h" - -namespace Eigen { - -/** \internal Specialization for the data types supported by MKL */ - -#define EIGEN_MKL_EIG_SELFADJ(EIGTYPE, MKLTYPE, MKLRTYPE, MKLNAME, EIGCOLROW, MKLCOLROW ) \ -template<> inline \ -SelfAdjointEigenSolver >& \ -SelfAdjointEigenSolver >::compute(const Matrix& 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; \ - 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)); \ - if(computeEigenvectors) m_eivec.setOnes(n,n); \ - m_info = Success; \ - m_isInitialized = true; \ - m_eigenvectorsOk = computeEigenvectors; \ - return *this; \ - } \ -\ - lda = matrix.outerStride(); \ - matrix_order=MKLCOLROW; \ - 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() ); \ - m_info = (info==0) ? Success : NoConvergence; \ - m_isInitialized = true; \ - m_eigenvectorsOk = computeEigenvectors; \ - return *this; \ -} - - -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_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) - -} // end namespace Eigen - -#endif // EIGEN_SAEIGENSOLVER_H diff --git a/eigen/Eigen/src/Eigenvalues/Tridiagonalization.h b/eigen/Eigen/src/Eigenvalues/Tridiagonalization.h index a63c08a..1d102c1 100644 --- a/eigen/Eigen/src/Eigenvalues/Tridiagonalization.h +++ b/eigen/Eigen/src/Eigenvalues/Tridiagonalization.h @@ -18,8 +18,10 @@ namespace internal { template struct TridiagonalizationMatrixTReturnType; template struct traits > + : public traits { - typedef typename MatrixType::PlainObject ReturnType; + typedef typename MatrixType::PlainObject ReturnType; // FIXME shall it be a BandMatrix? + enum { Flags = 0 }; }; template @@ -67,7 +69,7 @@ template class Tridiagonalization typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::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 class Tridiagonalization >::type DiagonalReturnType; typedef typename internal::conditional::IsComplex, - typename internal::add_const_on_value_type >::RealReturnType>::type, - const Diagonal< - Block > + typename internal::add_const_on_value_type::RealReturnType>::type, + const Diagonal >::type SubDiagonalReturnType; /** \brief Return type of matrixQ() */ @@ -110,7 +110,7 @@ template 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 class Tridiagonalization * Example: \include Tridiagonalization_Tridiagonalization_MatrixType.cpp * Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out */ - Tridiagonalization(const MatrixType& matrix) - : m_matrix(matrix), + template + explicit Tridiagonalization(const EigenBase& matrix) + : m_matrix(matrix.derived()), m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1), m_isInitialized(false) { @@ -152,9 +153,10 @@ template class Tridiagonalization * Example: \include Tridiagonalization_compute.cpp * Output: \verbinclude Tridiagonalization_compute.out */ - Tridiagonalization& compute(const MatrixType& matrix) + template + Tridiagonalization& compute(const EigenBase& 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::DiagonalReturnType Tridiagonalization::diagonal() const { eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); - return m_matrix.diagonal(); + return m_matrix.diagonal().real(); } template @@ -313,8 +315,7 @@ typename Tridiagonalization::SubDiagonalReturnType Tridiagonalization::subDiagonal() const { eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); - Index n = m_matrix.rows(); - return Block(m_matrix, 1, 0, n-1,n-1).diagonal(); + return m_matrix.template diagonal<-1>().real(); } namespace internal { @@ -346,7 +347,6 @@ template 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(); @@ -438,7 +438,6 @@ struct tridiagonalization_inplace_selector { typedef typename Tridiagonalization::CoeffVectorType CoeffVectorType; typedef typename Tridiagonalization::HouseholderSequenceType HouseholderSequenceType; - typedef typename MatrixType::Index Index; template static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) { @@ -467,9 +466,10 @@ struct tridiagonalization_inplace_selector static void run(MatrixType& mat, DiagonalType& diag, SubDiagonalType& subdiag, bool extractQ) { using std::sqrt; + const RealScalar tol = (std::numeric_limits::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 template struct TridiagonalizationMatrixTReturnType : public ReturnByValue > { - typedef typename MatrixType::Index Index; public: /** \brief Constructor. * -- cgit v1.2.3