From 8044b00a7fa30af20cc184fa2991bd0acf0f9aa3 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 3 Apr 2014 23:41:47 +0200 Subject: [PATCH 01/16] bug #782: Workaround for gcc <= 4.4 compilation error on the NEON PacketMath code. --- Eigen/src/Core/arch/NEON/PacketMath.h | 1 + 1 file changed, 1 insertion(+) diff --git a/Eigen/src/Core/arch/NEON/PacketMath.h b/Eigen/src/Core/arch/NEON/PacketMath.h index 05e891df2..c2f28e25c 100644 --- a/Eigen/src/Core/arch/NEON/PacketMath.h +++ b/Eigen/src/Core/arch/NEON/PacketMath.h @@ -95,6 +95,7 @@ template<> struct packet_traits : default_packet_traits // workaround gcc 4.2, 4.3 and 4.4 compilatin issue EIGEN_STRONG_INLINE float32x4_t vld1q_f32(const float* x) { return ::vld1q_f32((const float32_t*)x); } EIGEN_STRONG_INLINE float32x2_t vld1_f32 (const float* x) { return ::vld1_f32 ((const float32_t*)x); } +EIGEN_STRONG_INLINE float32x2_t vld1_dup_f32 (const float* x) { return ::vld1_dup_f32 ((const float32_t*)x); } EIGEN_STRONG_INLINE void vst1q_f32(float* to, float32x4_t from) { ::vst1q_f32((float32_t*)to,from); } EIGEN_STRONG_INLINE void vst1_f32 (float* to, float32x2_t from) { ::vst1_f32 ((float32_t*)to,from); } #endif From 096af597992609be62ec1104fdee476e3065f2e4 Mon Sep 17 00:00:00 2001 From: Christoph Hertzberg Date: Fri, 4 Apr 2014 17:48:37 +0200 Subject: [PATCH 02/16] Fix bug #784: Assert if assigning a product to a triangularView does not match the size. --- Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h index ffa871cae..225b994d1 100644 --- a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h +++ b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h @@ -265,6 +265,8 @@ template template TriangularView& TriangularView::assignProduct(const ProductBase& prod, const Scalar& alpha) { + eigen_assert(m_matrix.rows() == prod.rows() && m_matrix.cols() == prod.cols()); + general_product_to_triangular_selector::run(m_matrix.const_cast_derived(), prod.derived(), alpha); return *this; From 5afcb4965c9ccdce05e1852dd22b552559696f61 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Fri, 4 Apr 2014 16:48:13 +0100 Subject: [PATCH 03/16] Remove out-dated comment in cholesky test. --- test/cholesky.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/test/cholesky.cpp b/test/cholesky.cpp index d4d90e467..569318f83 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -82,10 +82,6 @@ template void cholesky(const MatrixType& m) symm += a1 * a1.adjoint(); } - // to test if really Cholesky only uses the upper triangular part, uncomment the following - // FIXME: currently that fails !! - //symm.template part().setZero(); - { SquareMatrixType symmUp = symm.template triangularView(); SquareMatrixType symmLo = symm.template triangularView(); From a91a7a1964305311133858de96b845da49389922 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Mon, 7 Apr 2014 14:14:48 +0100 Subject: [PATCH 04/16] doc: Add references to Cholesky methods in SelfAdjointView. --- Eigen/Cholesky | 6 ++++-- Eigen/src/Cholesky/LDLT.h | 6 ++++-- Eigen/src/Cholesky/LLT.h | 6 ++++-- doc/Doxyfile.in | 3 ++- 4 files changed, 14 insertions(+), 7 deletions(-) diff --git a/Eigen/Cholesky b/Eigen/Cholesky index f727f5d89..7314d326c 100644 --- a/Eigen/Cholesky +++ b/Eigen/Cholesky @@ -10,9 +10,11 @@ * * * This module provides two variants of the Cholesky decomposition for selfadjoint (hermitian) matrices. - * Those decompositions are accessible via the following MatrixBase methods: - * - MatrixBase::llt(), + * Those decompositions are also accessible via the following methods: + * - MatrixBase::llt() * - MatrixBase::ldlt() + * - SelfAdjointView::llt() + * - SelfAdjointView::ldlt() * * \code * #include diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index b43e85e7f..efac7fe40 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -43,7 +43,7 @@ namespace internal { * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky * decomposition to determine whether a system of equations has a solution. * - * \sa MatrixBase::ldlt(), class LLT + * \sa MatrixBase::ldlt(), SelfAdjointView::ldlt(), class LLT */ template class LDLT { @@ -179,7 +179,7 @@ template class LDLT * least-square solution of \f$ D y_3 = y_2 \f$ is computed. This does not mean that this function * computes the least-square solution of \f$ A x = b \f$ is \f$ A \f$ is singular. * - * \sa MatrixBase::ldlt() + * \sa MatrixBase::ldlt(), SelfAdjointView::ldlt() */ template inline const internal::solve_retval @@ -582,6 +582,7 @@ MatrixType LDLT::reconstructedMatrix() const #ifndef __CUDACC__ /** \cholesky_module * \returns the Cholesky decomposition with full pivoting without square root of \c *this + * \sa MatrixBase::ldlt() */ template inline const LDLT::PlainObject, UpLo> @@ -592,6 +593,7 @@ SelfAdjointView::ldlt() const /** \cholesky_module * \returns the Cholesky decomposition with full pivoting without square root of \c *this + * \sa SelfAdjointView::ldlt() */ template inline const LDLT::PlainObject> diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 2201c641e..45ed8438f 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -41,7 +41,7 @@ template struct LLT_Traits; * Example: \include LLT_example.cpp * Output: \verbinclude LLT_example.out * - * \sa MatrixBase::llt(), class LDLT + * \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT */ /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH) * Note that during the decomposition, only the upper triangular part of A is considered. Therefore, @@ -115,7 +115,7 @@ template class LLT * Example: \include LLT_solve.cpp * Output: \verbinclude LLT_solve.out * - * \sa solveInPlace(), MatrixBase::llt() + * \sa solveInPlace(), MatrixBase::llt(), SelfAdjointView::llt() */ template inline const internal::solve_retval @@ -468,6 +468,7 @@ MatrixType LLT::reconstructedMatrix() const #ifndef __CUDACC__ /** \cholesky_module * \returns the LLT decomposition of \c *this + * \sa SelfAdjointView::llt() */ template inline const LLT::PlainObject> @@ -478,6 +479,7 @@ MatrixBase::llt() const /** \cholesky_module * \returns the LLT decomposition of \c *this + * \sa SelfAdjointView::llt() */ template inline const LLT::PlainObject, UpLo> diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in index 85af9f1d4..7bbf693a0 100644 --- a/doc/Doxyfile.in +++ b/doc/Doxyfile.in @@ -223,7 +223,7 @@ ALIASES = "only_for_vectors=This is only for vectors (either row- "note_about_using_kernel_to_study_multiple_solutions=If you need a complete analysis of the space of solutions, take the one solution obtained by this method and add to it elements of the kernel, as determined by kernel()." \ "note_about_checking_solutions=This method just tries to find as good a solution as possible. If you want to check whether a solution exists or if it is accurate, just call this function to get a result and then compute the error of this result, or use MatrixBase::isApprox() directly, for instance like this: \code bool a_solution_exists = (A*result).isApprox(b, precision); \endcode This method avoids dividing by zero, so that the non-existence of a solution doesn't by itself mean that you'll get \c inf or \c nan values." \ "note_try_to_help_rvo=This function returns the result by value. In order to make that efficient, it is implemented as just a return statement using a special constructor, hopefully allowing the compiler to perform a RVO (return value optimization)." \ - "nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\" + "nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\"" ALIASES += "eigenAutoToc= " @@ -1583,6 +1583,7 @@ PREDEFINED = EIGEN_EMPTY_STRUCT \ EIGEN_VECTORIZE \ EIGEN_QT_SUPPORT \ EIGEN_STRONG_INLINE=inline \ + EIGEN_DEVICE_FUNC= \ "EIGEN2_SUPPORT_STAGE=99" \ "EIGEN_MAKE_CWISE_BINARY_OP(METHOD,FUNCTOR)=template const CwiseBinaryOp, const Derived, const OtherDerived> METHOD(const EIGEN_CURRENT_STORAGE_BASE_CLASS &other) const;" \ "EIGEN_CWISE_PRODUCT_RETURN_TYPE(LHS,RHS)=CwiseBinaryOp, const LHS, const RHS>" From a1fcf599faa023d1752b37aef215d80286e76839 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 10 Apr 2014 11:19:37 -0700 Subject: [PATCH 05/16] Silenced a compilation warning produced by nvcc. --- Eigen/src/Core/util/Memory.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index 0f8ab065a..1a1d3e98d 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -89,7 +89,7 @@ inline void throw_std_bad_alloc() #ifdef EIGEN_EXCEPTIONS throw std::bad_alloc(); #else - std::size_t huge = -1; + std::size_t huge = -1ULL; new int[huge]; #endif } From 1b333c89c9762795fda69f29a1239541a712f4f1 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Thu, 10 Apr 2014 17:43:13 -0700 Subject: [PATCH 06/16] Updated my previous fix to avoid introducing a compilation warning on ARM platforms. --- Eigen/src/Core/util/Memory.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index 1a1d3e98d..4988be5d9 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -89,7 +89,7 @@ inline void throw_std_bad_alloc() #ifdef EIGEN_EXCEPTIONS throw std::bad_alloc(); #else - std::size_t huge = -1ULL; + std::size_t huge = static_cast(-1); new int[huge]; #endif } From 91288e9bf9d6a9c3d63a9152e863b4390d0a40c7 Mon Sep 17 00:00:00 2001 From: Freddie Witherden Date: Sat, 12 Apr 2014 12:53:09 +0100 Subject: [PATCH 07/16] Add include LevenbergMarquardt in CMakeLists.txt. This fixes bug #768. --- unsupported/Eigen/src/CMakeLists.txt | 1 + unsupported/Eigen/src/LevenbergMarquardt/CMakeLists.txt | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/unsupported/Eigen/src/CMakeLists.txt b/unsupported/Eigen/src/CMakeLists.txt index f3180b52b..8eb2808e3 100644 --- a/unsupported/Eigen/src/CMakeLists.txt +++ b/unsupported/Eigen/src/CMakeLists.txt @@ -2,6 +2,7 @@ ADD_SUBDIRECTORY(AutoDiff) ADD_SUBDIRECTORY(BVH) ADD_SUBDIRECTORY(FFT) ADD_SUBDIRECTORY(IterativeSolvers) +ADD_SUBDIRECTORY(LevenbergMarquardt) ADD_SUBDIRECTORY(MatrixFunctions) ADD_SUBDIRECTORY(MoreVectorization) ADD_SUBDIRECTORY(NonLinearOptimization) diff --git a/unsupported/Eigen/src/LevenbergMarquardt/CMakeLists.txt b/unsupported/Eigen/src/LevenbergMarquardt/CMakeLists.txt index 8513803ce..d9690854d 100644 --- a/unsupported/Eigen/src/LevenbergMarquardt/CMakeLists.txt +++ b/unsupported/Eigen/src/LevenbergMarquardt/CMakeLists.txt @@ -2,5 +2,5 @@ FILE(GLOB Eigen_LevenbergMarquardt_SRCS "*.h") INSTALL(FILES ${Eigen_LevenbergMarquardt_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/LevenbergMarquardt COMPONENT Devel + DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/LevenbergMarquardt COMPONENT Devel ) From a803ff18a963fc23ddb3dccf33ed7058af415f39 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Sat, 12 Apr 2014 20:24:05 -0700 Subject: [PATCH 08/16] Fixed a typo in cuda_basic.cu --- test/cuda_basic.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/cuda_basic.cu b/test/cuda_basic.cu index aa7f7a599..4c7e96c10 100644 --- a/test/cuda_basic.cu +++ b/test/cuda_basic.cu @@ -129,7 +129,7 @@ void test_cuda_basic() CALL_SUBTEST( run_and_compare_to_cuda(prod(), nthreads, in, out) ); CALL_SUBTEST( run_and_compare_to_cuda(diagonal(), nthreads, in, out) ); - CALL_SUBTEST( run_and_compare_to_c(), nthreads, in, out) ); + CALL_SUBTEST( run_and_compare_to_cuda(diagonal(), nthreads, in, out) ); CALL_SUBTEST( run_and_compare_to_cuda(eigenvalues(), nthreads, in, out) ); CALL_SUBTEST( run_and_compare_to_cuda(eigenvalues(), nthreads, in, out) ); From 7903d3f27b275040702ce30eac8d329d6f571205 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Sat, 12 Apr 2014 23:39:37 -0700 Subject: [PATCH 09/16] Updated the compiler flags to enable nvcc to work with clang. --- test/CMakeLists.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 62cbedae7..c2d827051 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -315,6 +315,9 @@ find_package(CUDA) if(CUDA_FOUND) set(CUDA_PROPAGATE_HOST_FLAGS OFF) + if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang") + set(CUDA_NVCC_FLAGS "-ccbin /usr/bin/clang" CACHE STRING "nvcc flags" FORCE) + endif() cuda_include_directories(${CMAKE_CURRENT_BINARY_DIR}) set(EIGEN_ADD_TEST_FILENAME_EXTENSION "cu") From 0587db8bf5ca5d5eb6fb8df1c02abddd5b5718ba Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 14 Apr 2014 11:43:08 +0200 Subject: [PATCH 10/16] bug #793: fix overflow in EigenSolver and add respective regression unit test --- Eigen/src/Eigenvalues/EigenSolver.h | 15 ++++++++++++++- test/eigensolver_generic.cpp | 13 +++++++++++++ 2 files changed, 27 insertions(+), 1 deletion(-) diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index bf20e03ef..739466949 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -366,6 +366,7 @@ EigenSolver::compute(const MatrixType& matrix, bool computeEigenvect { using std::sqrt; using std::abs; + using std::max; eigen_assert(matrix.cols() == matrix.rows()); // Reduce to real Schur form. @@ -390,7 +391,19 @@ EigenSolver::compute(const MatrixType& matrix, bool computeEigenvect 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 = (max)(abs(p),(max)(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); i += 2; diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp index 005af81eb..91383b5cf 100644 --- a/test/eigensolver_generic.cpp +++ b/test/eigensolver_generic.cpp @@ -121,5 +121,18 @@ void test_eigensolver_generic() } ); + // regression test for bug 793 +#ifdef EIGEN_TEST_PART_2 + { + MatrixXd a(3,3); + a << 0, 0, 1, + 1, 1, 1, + 1, 1e+200, 1; + Eigen::EigenSolver eig(a); + VERIFY_IS_APPROX(a * eig.pseudoEigenvectors(), eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()); + VERIFY_IS_APPROX(a * eig.eigenvectors(), eig.eigenvectors() * eig.eigenvalues().asDiagonal()); + } +#endif + TEST_SET_BUT_UNUSED_VARIABLE(s) } From 148acf8e4fb71294703d4d1deafaf52829535ab7 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 14 Apr 2014 13:52:16 +0200 Subject: [PATCH 11/16] bug #790: fix overflow in real_2x2_jacobi_svd --- Eigen/src/SVD/JacobiSVD.h | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index eee31ca97..439eb5d29 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -415,6 +415,7 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, JacobiRotation *j_right) { using std::sqrt; + using std::abs; Matrix m; m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); @@ -428,9 +429,11 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, } else { - RealScalar u = d / t; - rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u)); - rot1.s() = rot1.c() * u; + RealScalar t2d2 = numext::hypot(t,d); + rot1.c() = abs(t)/t2d2; + rot1.s() = d/t2d2; + if(tmakeJacobi(m,0,1); From 7098e6d976ee8d5b25776e749d3ef6e66a302829 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 14 Apr 2014 21:57:49 +0200 Subject: [PATCH 12/16] Add isfinite overload for complexes. --- Eigen/src/Core/MathFunctions.h | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 63fb92b75..20fc2be74 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -669,6 +669,15 @@ bool (isfinite)(const T& x) return x::highest() && x>NumTraits::lowest(); } +template +EIGEN_DEVICE_FUNC +bool (isfinite)(const std::complex& x) +{ + using std::real; + using std::imag; + return isfinite(real(x)) && isfinite(imag(x)); +} + } // end namespace numext namespace internal { From 3c66bb136bf2adcb9d73d3d66850a8b907bc9264 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 14 Apr 2014 22:00:27 +0200 Subject: [PATCH 13/16] bug #793: detect NaN and INF in EigenSolver instead of aborting with an assert. --- Eigen/src/Eigenvalues/EigenSolver.h | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index 739466949..d2563d470 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -275,10 +275,11 @@ template class EigenSolver */ EigenSolver& compute(const MatrixType& 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. */ @@ -302,6 +303,7 @@ template class EigenSolver EigenvalueType m_eivalues; bool m_isInitialized; bool m_eigenvectorsOk; + ComputationInfo m_info; RealSchur m_realSchur; MatrixType m_matT; @@ -367,12 +369,15 @@ EigenSolver::compute(const MatrixType& matrix, bool computeEigenvect using std::sqrt; using std::abs; using std::max; + using numext::isfinite; eigen_assert(matrix.cols() == matrix.rows()); // Reduce to real Schur form. m_realSchur.compute(matrix, computeEigenvectors); + + m_info = m_realSchur.info(); - if (m_realSchur.info() == Success) + if (m_info == Success) { m_matT = m_realSchur.matrixT(); if (computeEigenvectors) @@ -386,6 +391,13 @@ 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 @@ -406,6 +418,13 @@ EigenSolver::compute(const MatrixType& matrix, bool computeEigenvect 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; } } @@ -594,7 +613,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 } } From 1afd50e0f31d56f7ee228e915e15998422c3ea11 Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Mon, 14 Apr 2014 14:26:30 -0700 Subject: [PATCH 14/16] Fixed a typo in CXX11Meta.h --- unsupported/Eigen/CXX11/src/Core/util/CXX11Meta.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unsupported/Eigen/CXX11/src/Core/util/CXX11Meta.h b/unsupported/Eigen/CXX11/src/Core/util/CXX11Meta.h index d6b5d75d9..618e2eb7b 100644 --- a/unsupported/Eigen/CXX11/src/Core/util/CXX11Meta.h +++ b/unsupported/Eigen/CXX11/src/Core/util/CXX11Meta.h @@ -285,7 +285,7 @@ struct equal_op { template constexpr static inli struct not_equal_op { template constexpr static inline auto run(A a, B b) -> decltype(a != b) { return a != b; } }; struct lesser_op { template constexpr static inline auto run(A a, B b) -> decltype(a < b) { return a < b; } }; struct lesser_equal_op { template constexpr static inline auto run(A a, B b) -> decltype(a <= b) { return a <= b; } }; -struct greater_op { template constexpr static inline auto run(A a, B b) -> decltype(a < b) { return a < b; } }; +struct greater_op { template constexpr static inline auto run(A a, B b) -> decltype(a > b) { return a > b; } }; struct greater_equal_op { template constexpr static inline auto run(A a, B b) -> decltype(a >= b) { return a >= b; } }; /* generic unary operations */ From e0dbb68c2f17f3c8c6accc7dc0b2b8d544e2eebc Mon Sep 17 00:00:00 2001 From: Mark Borgerding Date: Tue, 15 Apr 2014 13:57:03 -0400 Subject: [PATCH 15/16] Check IMKL version for compatibility with Eigen --- Eigen/src/Core/util/MKL_support.h | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/Eigen/src/Core/util/MKL_support.h b/Eigen/src/Core/util/MKL_support.h index 1e6e355d6..8acca9c8c 100644 --- a/Eigen/src/Core/util/MKL_support.h +++ b/Eigen/src/Core/util/MKL_support.h @@ -54,8 +54,25 @@ #endif #if defined EIGEN_USE_MKL +# include +/*Check IMKL version for compatibility: < 10.3 is not usable with Eigen*/ +# ifndef INTEL_MKL_VERSION +# undef EIGEN_USE_MKL /* INTEL_MKL_VERSION is not even defined on older versions */ +# elif INTEL_MKL_VERSION < 100305 /* the intel-mkl-103-release-notes say this was when the lapacke.h interface was added*/ +# undef EIGEN_USE_MKL +# endif +# ifndef EIGEN_USE_MKL + /*If the MKL version is too old, undef everything*/ +# undef EIGEN_USE_MKL_ALL +# undef EIGEN_USE_BLAS +# undef EIGEN_USE_LAPACKE +# undef EIGEN_USE_MKL_VML +# undef EIGEN_USE_LAPACKE_STRICT +# undef EIGEN_USE_LAPACKE +# endif +#endif -#include +#if defined EIGEN_USE_MKL #include #define EIGEN_MKL_VML_THRESHOLD 128 From e5d0cb54a5f2a2200a4656d993c82a80f159a7c4 Mon Sep 17 00:00:00 2001 From: Benjamin Chretien Date: Thu, 17 Apr 2014 18:49:23 +0200 Subject: [PATCH 16/16] Fix typo in Reductions tutorial. --- doc/TutorialReductionsVisitorsBroadcasting.dox | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/TutorialReductionsVisitorsBroadcasting.dox b/doc/TutorialReductionsVisitorsBroadcasting.dox index 992cf6f34..eb6787dbc 100644 --- a/doc/TutorialReductionsVisitorsBroadcasting.dox +++ b/doc/TutorialReductionsVisitorsBroadcasting.dox @@ -32,7 +32,7 @@ Eigen also provides the \link MatrixBase::norm() norm() \endlink method, which r These operations can also operate on matrices; in that case, a n-by-p matrix is seen as a vector of size (n*p), so for example the \link MatrixBase::norm() norm() \endlink method returns the "Frobenius" or "Hilbert-Schmidt" norm. We refrain from speaking of the \f$\ell^2\f$ norm of a matrix because that can mean different things. -If you want other \f$\ell^p\f$ norms, use the \link MatrixBase::lpNorm() lpNnorm

() \endlink method. The template parameter \a p can take the special value \a Infinity if you want the \f$\ell^\infty\f$ norm, which is the maximum of the absolute values of the coefficients. +If you want other \f$\ell^p\f$ norms, use the \link MatrixBase::lpNorm() lpNorm

() \endlink method. The template parameter \a p can take the special value \a Infinity if you want the \f$\ell^\infty\f$ norm, which is the maximum of the absolute values of the coefficients. The following example demonstrates these methods.