diff options
Diffstat (limited to 'Eigen/src/Eigenvalues')
| -rw-r--r-- | Eigen/src/Eigenvalues/CMakeLists.txt | 6 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h (renamed from Eigen/src/Eigenvalues/ComplexSchur_MKL.h) | 34 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/EigenSolver.h | 36 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/GeneralizedEigenSolver.h | 193 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/RealQZ.h | 32 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/RealSchur.h | 12 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/RealSchur_LAPACKE.h (renamed from Eigen/src/Eigenvalues/RealSchur_MKL.h) | 34 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 40 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h (renamed from Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h) | 36 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/Tridiagonalization.h | 4 |
10 files changed, 254 insertions, 173 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/ComplexSchur_MKL.h b/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h index e20c3725b..4980a3ede 100644 --- a/Eigen/src/Eigenvalues/ComplexSchur_MKL.h +++ b/Eigen/src/Eigenvalues/ComplexSchur_LAPACKE.h @@ -25,21 +25,19 @@ 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) \ +#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 EigenBase<InputType>& matrix, bool computeU) \ @@ -60,18 +58,18 @@ ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Eigen m_matUisUptodate = computeU; \ return *this; \ } \ - lapack_int n = matrix.cols(), sdim, info; \ - 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 = m_matT.outerStride(); \ + 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 \ @@ -83,11 +81,11 @@ ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Eigen \ } -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 532ca7d63..f205b185d 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -324,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 { @@ -345,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>(); @@ -451,26 +453,6 @@ EigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool comput 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() @@ -503,7 +485,7 @@ 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; @@ -557,7 +539,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors() } else { - std::complex<Scalar> cc = cdiv<Scalar>(Scalar(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); } @@ -580,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); } @@ -594,7 +576,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors() 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))) @@ -604,7 +586,7 @@ 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); } diff --git a/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h b/Eigen/src/Eigenvalues/GeneralizedEigenSolver.h index a9d6790d5..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 @@ -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 * @@ -126,10 +134,9 @@ template<typename _MatrixType> class GeneralizedEigenSolver : 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/RealQZ.h b/Eigen/src/Eigenvalues/RealQZ.h index a62071d42..b3a910dd9 100644 --- a/Eigen/src/Eigenvalues/RealQZ.h +++ b/Eigen/src/Eigenvalues/RealQZ.h @@ -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 f4ded69b6..d6a339f07 100644 --- a/Eigen/src/Eigenvalues/RealSchur.h +++ b/Eigen/src/Eigenvalues/RealSchur.h @@ -253,19 +253,25 @@ RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const EigenBase<InputType> if (maxIters == -1) maxIters = m_maxIterationsPerRow * matrix.rows(); + Scalar scale = matrix.derived().cwiseAbs().maxCoeff(); + // Step 1. Reduce to Hessenberg form - m_hess.compute(matrix.derived()); + 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; diff --git a/Eigen/src/Eigenvalues/RealSchur_MKL.h b/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h index e80926400..2c2251715 100644 --- a/Eigen/src/Eigenvalues/RealSchur_MKL.h +++ b/Eigen/src/Eigenvalues/RealSchur_LAPACKE.h @@ -25,39 +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) \ +#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 EigenBase<InputType>& matrix, bool computeU) \ { \ eigen_assert(matrix.cols() == matrix.rows()); \ \ - lapack_int n = matrix.cols(), sdim, info; \ - 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 = m_matT.outerStride(); \ + 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 \ @@ -69,11 +67,11 @@ RealSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const EigenBas \ } -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 469ea5e4e..a9f56c4f5 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -414,7 +414,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> if(n==1) { - m_eivalues.coeffRef(0,0) = numext::real(matrix(0,0)); + m_eivalues.coeffRef(0,0) = numext::real(matrix.diagonal()[0]); if(computeEigenvectors) m_eivec.setOnes(n,n); m_info = Success; @@ -458,7 +458,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> { m_eivec.setIdentity(diag.size(), diag.size()); } - m_info = computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec); + m_info = internal::computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec); m_isInitialized = true; m_eigenvectorsOk = computeEigenvectors; @@ -492,15 +492,16 @@ ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag 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(subdiag[i]),(abs(diag[i])+abs(diag[i+1]))) || abs(subdiag[i]) <= considerAsZero) + 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 && subdiag[end-1]==0) + while (end>0 && subdiag[end-1]==RealScalar(0)) { end--; } @@ -568,8 +569,8 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 EIGEN_USING_STD_MATH(atan2) EIGEN_USING_STD_MATH(cos) EIGEN_USING_STD_MATH(sin) - const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0); - const Scalar s_sqrt3 = sqrt(Scalar(3.0)); + 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 @@ -739,14 +740,18 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false> 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 = numext::maxi(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) { @@ -774,10 +779,11 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false> 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; @@ -809,14 +815,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 3499dc78a..3891cf883 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_LAPACKE.h @@ -25,21 +25,19 @@ 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 ) \ +#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 EigenBase<InputType>& matrix, int options) \ @@ -49,7 +47,7 @@ SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(c && (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; \ @@ -64,12 +62,12 @@ SelfAdjointEigenSolver<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(c return *this; \ } \ \ - lda = m_eivec.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 2030b5be1..1d102c17b 100644 --- a/Eigen/src/Eigenvalues/Tridiagonalization.h +++ b/Eigen/src/Eigenvalues/Tridiagonalization.h @@ -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; |