diff --git a/eigenlib/Eigen/Core b/eigenlib/Eigen/Core index 6e855427..a5025e37 100644 --- a/eigenlib/Eigen/Core +++ b/eigenlib/Eigen/Core @@ -167,7 +167,7 @@ #include #endif -#if (defined(_CPPUNWIND) || defined(__EXCEPTIONS)) && !defined(EIGEN_NO_EXCEPTIONS) +#if defined(_CPPUNWIND) || defined(__EXCEPTIONS) #define EIGEN_EXCEPTIONS #endif diff --git a/eigenlib/Eigen/SVD b/eigenlib/Eigen/SVD index d24471fd..7c987a9d 100644 --- a/eigenlib/Eigen/SVD +++ b/eigenlib/Eigen/SVD @@ -13,9 +13,9 @@ namespace Eigen { * * * - * This module provides SVD decomposition for (currently) real matrices. + * This module provides SVD decomposition for matrices (both real and complex). * This decomposition is accessible via the following MatrixBase method: - * - MatrixBase::svd() + * - MatrixBase::jacobiSvd() * * \code * #include diff --git a/eigenlib/Eigen/src/Cholesky/LDLT.h b/eigenlib/Eigen/src/Cholesky/LDLT.h index 5e2352ca..f47b2ea5 100644 --- a/eigenlib/Eigen/src/Cholesky/LDLT.h +++ b/eigenlib/Eigen/src/Cholesky/LDLT.h @@ -158,10 +158,19 @@ template class LDLT } /** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A. + * + * This function also supports in-place solves using the syntax x = decompositionObject.solve(x) . * * \note_about_checking_solutions * - * \sa solveInPlace(), MatrixBase::ldlt() + * More precisely, this method solves \f$ A x = b \f$ using the decomposition \f$ A = P^T L D L^* P \f$ + * by solving the systems \f$ P^T y_1 = b \f$, \f$ L y_2 = y_1 \f$, \f$ D y_3 = y_2 \f$, + * \f$ L^* y_4 = y_3 \f$ and \f$ P x = y_4 \f$ in succession. If the matrix \f$ A \f$ is singular, then + * \f$ D \f$ will also be singular (all the other matrices are invertible). In that case, the + * 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() */ template inline const internal::solve_retval @@ -376,7 +385,21 @@ struct solve_retval, Rhs> dec().matrixL().solveInPlace(dst); // dst = D^-1 (L^-1 P b) - dst = dec().vectorD().asDiagonal().inverse() * dst; + // more precisely, use pseudo-inverse of D (see bug 241) + using std::abs; + using std::max; + typedef typename LDLTType::MatrixType MatrixType; + typedef typename LDLTType::Scalar Scalar; + typedef typename LDLTType::RealScalar RealScalar; + const Diagonal vectorD = dec().vectorD(); + RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() * NumTraits::epsilon(), + RealScalar(1) / NumTraits::highest()); // motivated by LAPACK's xGELSS + for (Index i = 0; i < vectorD.size(); ++i) { + if(abs(vectorD(i)) > tolerance) + dst.row(i) /= vectorD(i); + else + dst.row(i).setZero(); + } // dst = L^-T (D^-1 L^-1 P b) dec().matrixU().solveInPlace(dst); diff --git a/eigenlib/Eigen/src/Core/Array.h b/eigenlib/Eigen/src/Core/Array.h index a3a2167a..a11fb1b5 100644 --- a/eigenlib/Eigen/src/Core/Array.h +++ b/eigenlib/Eigen/src/Core/Array.h @@ -68,10 +68,8 @@ class Array friend struct internal::conservative_resize_like_impl; using Base::m_storage; + public: - enum { NeedsToAlign = (!(Options&DontAlign)) - && SizeAtCompileTime!=Dynamic && ((static_cast(sizeof(Scalar))*SizeAtCompileTime)%16)==0 }; - EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign) using Base::base; using Base::coeff; diff --git a/eigenlib/Eigen/src/Core/Block.h b/eigenlib/Eigen/src/Core/Block.h index 2b251bc2..d508a60a 100644 --- a/eigenlib/Eigen/src/Core/Block.h +++ b/eigenlib/Eigen/src/Core/Block.h @@ -94,7 +94,7 @@ struct traits MaskPacketAccessBit = (InnerSize == Dynamic || (InnerSize % packet_traits::size) == 0) && (InnerStrideAtCompileTime == 1) ? PacketAccessBit : 0, - MaskAlignedBit = (InnerPanel && (OuterStrideAtCompileTime!=Dynamic) && ((OuterStrideAtCompileTime % packet_traits::size) == 0)) ? AlignedBit : 0, + MaskAlignedBit = (InnerPanel && (OuterStrideAtCompileTime!=Dynamic) && (((OuterStrideAtCompileTime * int(sizeof(Scalar))) % 16) == 0)) ? AlignedBit : 0, FlagsLinearAccessBit = (RowsAtCompileTime == 1 || ColsAtCompileTime == 1) ? LinearAccessBit : 0, FlagsLvalueBit = is_lvalue::value ? LvalueBit : 0, FlagsRowMajorBit = IsRowMajor ? RowMajorBit : 0, diff --git a/eigenlib/Eigen/src/Core/Fuzzy.h b/eigenlib/Eigen/src/Core/Fuzzy.h index 1926d6ab..d266eed0 100644 --- a/eigenlib/Eigen/src/Core/Fuzzy.h +++ b/eigenlib/Eigen/src/Core/Fuzzy.h @@ -94,7 +94,7 @@ struct isMuchSmallerThan_scalar_selector * * \note The fuzzy compares are done multiplicatively. Two vectors \f$ v \f$ and \f$ w \f$ * are considered to be approximately equal within precision \f$ p \f$ if - * \f[ \Vert v - w \Vert \leqslant p\,\(min)(\Vert v\Vert, \Vert w\Vert). \f] + * \f[ \Vert v - w \Vert \leqslant p\,\min(\Vert v\Vert, \Vert w\Vert). \f] * For matrices, the comparison is done using the Hilbert-Schmidt norm (aka Frobenius norm * L2 norm). * diff --git a/eigenlib/Eigen/src/Core/Map.h b/eigenlib/Eigen/src/Core/Map.h index 81e3979f..2bf80b3a 100644 --- a/eigenlib/Eigen/src/Core/Map.h +++ b/eigenlib/Eigen/src/Core/Map.h @@ -34,7 +34,7 @@ * \tparam PlainObjectType the equivalent matrix type of the mapped data * \tparam MapOptions specifies whether the pointer is \c #Aligned, or \c #Unaligned. * The default is \c #Unaligned. - * \tparam StrideType optionnally specifies strides. By default, Map assumes the memory layout + * \tparam StrideType optionally specifies strides. By default, Map assumes the memory layout * of an ordinary, contiguous array. This can be overridden by specifying strides. * The type passed here must be a specialization of the Stride template, see examples below. * @@ -72,9 +72,9 @@ * Example: \include Map_placement_new.cpp * Output: \verbinclude Map_placement_new.out * - * This class is the return type of Matrix::Map() but can also be used directly. + * This class is the return type of PlainObjectBase::Map() but can also be used directly. * - * \sa Matrix::Map(), \ref TopicStorageOrders + * \sa PlainObjectBase::Map(), \ref TopicStorageOrders */ namespace internal { @@ -102,7 +102,7 @@ struct traits > || HasNoOuterStride || ( OuterStrideAtCompileTime!=Dynamic && ((static_cast(sizeof(Scalar))*OuterStrideAtCompileTime)%16)==0 ) ), - Flags0 = TraitsBase::Flags, + Flags0 = TraitsBase::Flags & (~NestByRefBit), Flags1 = IsAligned ? (int(Flags0) | AlignedBit) : (int(Flags0) & ~AlignedBit), Flags2 = (bool(HasNoStride) || bool(PlainObjectType::IsVectorAtCompileTime)) ? int(Flags1) : int(Flags1 & ~LinearAccessBit), @@ -120,7 +120,6 @@ template class Ma public: typedef MapBase Base; - EIGEN_DENSE_PUBLIC_INTERFACE(Map) typedef typename Base::PointerType PointerType; @@ -181,7 +180,6 @@ template class Ma PlainObjectType::Base::_check_template_params(); } - EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Map) protected: diff --git a/eigenlib/Eigen/src/Core/MapBase.h b/eigenlib/Eigen/src/Core/MapBase.h index c23bcbfd..9426e2d2 100644 --- a/eigenlib/Eigen/src/Core/MapBase.h +++ b/eigenlib/Eigen/src/Core/MapBase.h @@ -170,8 +170,8 @@ template class MapBase EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(internal::traits::Flags&PacketAccessBit, internal::inner_stride_at_compile_time::ret==1), PACKET_ACCESS_REQUIRES_TO_HAVE_INNER_STRIDE_FIXED_TO_1); - eigen_assert(EIGEN_IMPLIES(internal::traits::Flags&AlignedBit, (size_t(m_data) % (sizeof(Scalar)*internal::packet_traits::size)) == 0) - && "data is not aligned"); + eigen_assert(EIGEN_IMPLIES(internal::traits::Flags&AlignedBit, (size_t(m_data) % 16) == 0) + && "data is not aligned"); } PointerType m_data; diff --git a/eigenlib/Eigen/src/Core/Matrix.h b/eigenlib/Eigen/src/Core/Matrix.h index 44de22cb..982c9256 100644 --- a/eigenlib/Eigen/src/Core/Matrix.h +++ b/eigenlib/Eigen/src/Core/Matrix.h @@ -153,10 +153,6 @@ class Matrix typedef typename Base::PlainObject PlainObject; - enum { NeedsToAlign = (!(Options&DontAlign)) - && SizeAtCompileTime!=Dynamic && ((static_cast(sizeof(Scalar))*SizeAtCompileTime)%16)==0 }; - EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign) - using Base::base; using Base::coeffRef; diff --git a/eigenlib/Eigen/src/Core/MatrixBase.h b/eigenlib/Eigen/src/Core/MatrixBase.h index db156f6e..62877bce 100644 --- a/eigenlib/Eigen/src/Core/MatrixBase.h +++ b/eigenlib/Eigen/src/Core/MatrixBase.h @@ -250,7 +250,8 @@ template class MatrixBase // huuuge hack. make Eigen2's matrix.part() work in eigen3. Problem: Diagonal is now a class template instead // of an integer constant. Solution: overload the part() method template wrt template parameters list. - template class U> + // Note: replacing next line by "template class U>" produces a mysterious error C2082 in MSVC. + template class U> const DiagonalWrapper part() const { return diagonal().asDiagonal(); } #endif // EIGEN2_SUPPORT diff --git a/eigenlib/Eigen/src/Core/PlainObjectBase.h b/eigenlib/Eigen/src/Core/PlainObjectBase.h index ed34b0ed..612254e9 100644 --- a/eigenlib/Eigen/src/Core/PlainObjectBase.h +++ b/eigenlib/Eigen/src/Core/PlainObjectBase.h @@ -34,6 +34,19 @@ namespace internal { +template +EIGEN_ALWAYS_INLINE void check_rows_cols_for_overflow(Index rows, Index cols) +{ + // http://hg.mozilla.org/mozilla-central/file/6c8a909977d3/xpcom/ds/CheckedInt.h#l242 + // we assume Index is signed + Index max_index = (size_t(1) << (8 * sizeof(Index) - 1)) - 1; // assume Index is signed + bool error = (rows < 0 || cols < 0) ? true + : (rows == 0 || cols == 0) ? false + : (rows > max_index / cols); + if (error) + throw_std_bad_alloc(); +} + template (Derived::IsVectorAtCompileTime)> struct conservative_resize_like_impl; template struct matrix_swap_impl; @@ -84,14 +97,12 @@ class PlainObjectBase : public internal::dense_xpr_base::type template struct StridedConstMapType { typedef Eigen::Map type; }; template struct StridedAlignedMapType { typedef Eigen::Map type; }; template struct StridedConstAlignedMapType { typedef Eigen::Map type; }; - protected: DenseStorage m_storage; public: - enum { NeedsToAlign = (!(Options&DontAlign)) - && SizeAtCompileTime!=Dynamic && ((static_cast(sizeof(Scalar))*SizeAtCompileTime)%16)==0 }; + enum { NeedsToAlign = SizeAtCompileTime != Dynamic && (internal::traits::Flags & AlignedBit) != 0 }; EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign) Base& base() { return *static_cast(this); } @@ -200,11 +211,13 @@ class PlainObjectBase : public internal::dense_xpr_base::type EIGEN_STRONG_INLINE void resize(Index rows, Index cols) { #ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO + internal::check_rows_cols_for_overflow(rows, cols); Index size = rows*cols; bool size_changed = size != this->size(); m_storage.resize(size, rows, cols); if(size_changed) EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED #else + internal::check_rows_cols_for_overflow(rows, cols); m_storage.resize(rows*cols, rows, cols); #endif } @@ -273,6 +286,7 @@ class PlainObjectBase : public internal::dense_xpr_base::type EIGEN_STRONG_INLINE void resizeLike(const EigenBase& _other) { const OtherDerived& other = _other.derived(); + internal::check_rows_cols_for_overflow(other.rows(), other.cols()); const Index othersize = other.rows()*other.cols(); if(RowsAtCompileTime == 1) { @@ -417,6 +431,7 @@ class PlainObjectBase : public internal::dense_xpr_base::type : m_storage(other.derived().rows() * other.derived().cols(), other.derived().rows(), other.derived().cols()) { _check_template_params(); + internal::check_rows_cols_for_overflow(other.derived().rows(), other.derived().cols()); Base::operator=(other.derived()); } @@ -425,9 +440,6 @@ class PlainObjectBase : public internal::dense_xpr_base::type * while the AlignedMap() functions return aligned Map objects and thus should be called only with 16-byte-aligned * \a data pointers. * - * These methods do not allow to specify strides. If you need to specify strides, you have to - * use the Map class directly. - * * \see class Map */ //@{ @@ -584,6 +596,7 @@ class PlainObjectBase : public internal::dense_xpr_base::type { eigen_assert(rows >= 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows) && cols >= 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols)); + internal::check_rows_cols_for_overflow(rows, cols); m_storage.resize(rows*cols,rows,cols); EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED } @@ -641,6 +654,7 @@ struct internal::conservative_resize_like_impl if ( ( Derived::IsRowMajor && _this.cols() == cols) || // row-major and we change only the number of rows (!Derived::IsRowMajor && _this.rows() == rows) ) // column-major and we change only the number of columns { + internal::check_rows_cols_for_overflow(rows, cols); _this.derived().m_storage.conservativeResize(rows*cols,rows,cols); } else diff --git a/eigenlib/Eigen/src/Core/ProductBase.h b/eigenlib/Eigen/src/Core/ProductBase.h index 3bd3487d..233ed646 100644 --- a/eigenlib/Eigen/src/Core/ProductBase.h +++ b/eigenlib/Eigen/src/Core/ProductBase.h @@ -256,16 +256,16 @@ class ScaledProduct : Base(prod.lhs(),prod.rhs()), m_prod(prod), m_alpha(x) {} template - inline void evalTo(Dest& dst) const { dst.setZero(); scaleAndAddTo(dst,m_alpha); } + inline void evalTo(Dest& dst) const { dst.setZero(); scaleAndAddTo(dst, Scalar(1)); } template - inline void addTo(Dest& dst) const { scaleAndAddTo(dst,m_alpha); } + inline void addTo(Dest& dst) const { scaleAndAddTo(dst, Scalar(1)); } template - inline void subTo(Dest& dst) const { scaleAndAddTo(dst,-m_alpha); } + inline void subTo(Dest& dst) const { scaleAndAddTo(dst, Scalar(-1)); } template - inline void scaleAndAddTo(Dest& dst,Scalar alpha) const { m_prod.derived().scaleAndAddTo(dst,alpha); } + inline void scaleAndAddTo(Dest& dst,Scalar alpha) const { m_prod.derived().scaleAndAddTo(dst,alpha * m_alpha); } const Scalar& alpha() const { return m_alpha; } diff --git a/eigenlib/Eigen/src/Core/SolveTriangular.h b/eigenlib/Eigen/src/Core/SolveTriangular.h index 71e129c7..a23014a3 100644 --- a/eigenlib/Eigen/src/Core/SolveTriangular.h +++ b/eigenlib/Eigen/src/Core/SolveTriangular.h @@ -180,7 +180,7 @@ void TriangularView::solveInPlace(const MatrixBase::Flags & RowMajorBit && OtherDerived::IsVectorAtCompileTime }; typedef typename internal::conditional(X) +#if defined(__llvm__) && !defined(__clang__) + //Special treatment for Apple's llvm-gcc, its NEON packet types are unions + #define EIGEN_INIT_NEON_PACKET2(X, Y) {{X, Y}} + #define EIGEN_INIT_NEON_PACKET4(X, Y, Z, W) {{X, Y, Z, W}} +#else + //Default initializer for packets + #define EIGEN_INIT_NEON_PACKET2(X, Y) {X, Y} + #define EIGEN_INIT_NEON_PACKET4(X, Y, Z, W) {X, Y, Z, W} +#endif + #ifndef __pld #define __pld(x) asm volatile ( " pld [%[addr]]\n" :: [addr] "r" (x) : "cc" ); #endif @@ -84,7 +94,7 @@ template<> struct packet_traits : default_packet_traits }; }; -#if EIGEN_GNUC_AT_MOST(4,4) +#if EIGEN_GNUC_AT_MOST(4,4) && !defined(__llvm__) // 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); } @@ -100,12 +110,12 @@ template<> EIGEN_STRONG_INLINE Packet4i pset1(const int& from) { template<> EIGEN_STRONG_INLINE Packet4f plset(const float& a) { - Packet4f countdown = { 0, 1, 2, 3 }; + Packet4f countdown = EIGEN_INIT_NEON_PACKET4(0, 1, 2, 3); return vaddq_f32(pset1(a), countdown); } template<> EIGEN_STRONG_INLINE Packet4i plset(const int& a) { - Packet4i countdown = { 0, 1, 2, 3 }; + Packet4i countdown = EIGEN_INIT_NEON_PACKET4(0, 1, 2, 3); return vaddq_s32(pset1(a), countdown); } @@ -395,25 +405,29 @@ template<> EIGEN_STRONG_INLINE int predux_max(const Packet4i& a) return s[0]; } -template -struct palign_impl -{ - EIGEN_STRONG_INLINE static void run(Packet4f& first, const Packet4f& second) - { - if (Offset!=0) - first = vextq_f32(first, second, Offset); - } -}; +// this PALIGN_NEON business is to work around a bug in LLVM Clang 3.0 causing incorrect compilation errors, +// see bug 347 and this LLVM bug: http://llvm.org/bugs/show_bug.cgi?id=11074 +#define PALIGN_NEON(Offset,Type,Command) \ +template<>\ +struct palign_impl\ +{\ + EIGEN_STRONG_INLINE static void run(Type& first, const Type& second)\ + {\ + if (Offset!=0)\ + first = Command(first, second, Offset);\ + }\ +};\ -template -struct palign_impl -{ - EIGEN_STRONG_INLINE static void run(Packet4i& first, const Packet4i& second) - { - if (Offset!=0) - first = vextq_s32(first, second, Offset); - } -}; +PALIGN_NEON(0,Packet4f,vextq_f32) +PALIGN_NEON(1,Packet4f,vextq_f32) +PALIGN_NEON(2,Packet4f,vextq_f32) +PALIGN_NEON(3,Packet4f,vextq_f32) +PALIGN_NEON(0,Packet4i,vextq_s32) +PALIGN_NEON(1,Packet4i,vextq_s32) +PALIGN_NEON(2,Packet4i,vextq_s32) +PALIGN_NEON(3,Packet4i,vextq_s32) + +#undef PALIGN_NEON } // end namespace internal diff --git a/eigenlib/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/eigenlib/Eigen/src/Core/products/GeneralBlockPanelKernel.h index 6f3f2717..cd1c37c7 100644 --- a/eigenlib/Eigen/src/Core/products/GeneralBlockPanelKernel.h +++ b/eigenlib/Eigen/src/Core/products/GeneralBlockPanelKernel.h @@ -30,19 +30,16 @@ namespace internal { template class gebp_traits; +inline std::ptrdiff_t manage_caching_sizes_second_if_negative(std::ptrdiff_t a, std::ptrdiff_t b) +{ + return a<=0 ? b : a; +} + /** \internal */ inline void manage_caching_sizes(Action action, std::ptrdiff_t* l1=0, std::ptrdiff_t* l2=0) { - static std::ptrdiff_t m_l1CacheSize = 0; - static std::ptrdiff_t m_l2CacheSize = 0; - if(m_l1CacheSize==0) - { - m_l1CacheSize = queryL1CacheSize(); - m_l2CacheSize = queryTopLevelCacheSize(); - - if(m_l1CacheSize<=0) m_l1CacheSize = 8 * 1024; - if(m_l2CacheSize<=0) m_l2CacheSize = 1 * 1024 * 1024; - } + static std::ptrdiff_t m_l1CacheSize = manage_caching_sizes_second_if_negative(queryL1CacheSize(),8 * 1024); + static std::ptrdiff_t m_l2CacheSize = manage_caching_sizes_second_if_negative(queryTopLevelCacheSize(),1*1024*1024); if(action==SetAction) { @@ -118,14 +115,14 @@ inline void computeProductBlockingSizes(std::ptrdiff_t& k, std::ptrdiff_t& m, st // FIXME (a bit overkill maybe ?) template struct gebp_madd_selector { - EIGEN_STRONG_INLINE EIGEN_ALWAYS_INLINE_ATTRIB static void run(const CJ& cj, A& a, B& b, C& c, T& /*t*/) + EIGEN_ALWAYS_INLINE static void run(const CJ& cj, A& a, B& b, C& c, T& /*t*/) { c = cj.pmadd(a,b,c); } }; template struct gebp_madd_selector { - EIGEN_STRONG_INLINE EIGEN_ALWAYS_INLINE_ATTRIB static void run(const CJ& cj, T& a, T& b, T& c, T& t) + EIGEN_ALWAYS_INLINE static void run(const CJ& cj, T& a, T& b, T& c, T& t) { t = b; t = cj.pmul(a,t); c = padd(c,t); } diff --git a/eigenlib/Eigen/src/Core/util/Macros.h b/eigenlib/Eigen/src/Core/util/Macros.h index 6c3f1e42..b7c2b79a 100644 --- a/eigenlib/Eigen/src/Core/util/Macros.h +++ b/eigenlib/Eigen/src/Core/util/Macros.h @@ -1,3 +1,4 @@ + // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // @@ -28,7 +29,7 @@ #define EIGEN_WORLD_VERSION 3 #define EIGEN_MAJOR_VERSION 0 -#define EIGEN_MINOR_VERSION 2 +#define EIGEN_MINOR_VERSION 5 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \ (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \ @@ -45,7 +46,7 @@ #define EIGEN_GNUC_AT_MOST(x,y) 0 #endif -#if EIGEN_GNUC_AT_MOST(4,3) +#if EIGEN_GNUC_AT_MOST(4,3) && !defined(__clang__) // see bug 89 #define EIGEN_SAFE_TO_USE_STANDARD_ASSERT_MACRO 0 #else @@ -130,31 +131,34 @@ #define EIGEN_MAKESTRING2(a) #a #define EIGEN_MAKESTRING(a) EIGEN_MAKESTRING2(a) -// EIGEN_ALWAYS_INLINE_ATTRIB should be use in the declaration of function -// which should be inlined even in debug mode. -// FIXME with the always_inline attribute, -// gcc 3.4.x reports the following compilation error: -// Eval.h:91: sorry, unimplemented: inlining failed in call to 'const Eigen::Eval Eigen::MatrixBase::eval() const' -// : function body not available -#if EIGEN_GNUC_AT_LEAST(4,0) -#define EIGEN_ALWAYS_INLINE_ATTRIB __attribute__((always_inline)) -#else -#define EIGEN_ALWAYS_INLINE_ATTRIB -#endif - #if EIGEN_GNUC_AT_LEAST(4,1) && !defined(__clang__) && !defined(__INTEL_COMPILER) #define EIGEN_FLATTEN_ATTRIB __attribute__((flatten)) #else #define EIGEN_FLATTEN_ATTRIB #endif -// EIGEN_FORCE_INLINE means "inline as much as possible" +// EIGEN_STRONG_INLINE is a stronger version of the inline, using __forceinline on MSVC, +// but it still doesn't use GCC's always_inline. This is useful in (common) situations where MSVC needs forceinline +// but GCC is still doing fine with just inline. #if (defined _MSC_VER) || (defined __INTEL_COMPILER) #define EIGEN_STRONG_INLINE __forceinline #else #define EIGEN_STRONG_INLINE inline #endif +// EIGEN_ALWAYS_INLINE is the stronget, it has the effect of making the function inline and adding every possible +// attribute to maximize inlining. This should only be used when really necessary: in particular, +// it uses __attribute__((always_inline)) on GCC, which most of the time is useless and can severely harm compile times. +// FIXME with the always_inline attribute, +// gcc 3.4.x reports the following compilation error: +// Eval.h:91: sorry, unimplemented: inlining failed in call to 'const Eigen::Eval Eigen::MatrixBase::eval() const' +// : function body not available +#if EIGEN_GNUC_AT_LEAST(4,0) +#define EIGEN_ALWAYS_INLINE __attribute__((always_inline)) inline +#else +#define EIGEN_ALWAYS_INLINE EIGEN_STRONG_INLINE +#endif + #if (defined __GNUC__) #define EIGEN_DONT_INLINE __attribute__((noinline)) #elif (defined _MSC_VER) @@ -249,7 +253,7 @@ #define EIGEN_UNUSED_VARIABLE(var) (void)var; #if (defined __GNUC__) -#define EIGEN_ASM_COMMENT(X) asm("#"X) +#define EIGEN_ASM_COMMENT(X) asm("#" X) #else #define EIGEN_ASM_COMMENT(X) #endif diff --git a/eigenlib/Eigen/src/Core/util/Memory.h b/eigenlib/Eigen/src/Core/util/Memory.h index a580b95a..023716dc 100644 --- a/eigenlib/Eigen/src/Core/util/Memory.h +++ b/eigenlib/Eigen/src/Core/util/Memory.h @@ -82,6 +82,16 @@ namespace internal { +inline void throw_std_bad_alloc() +{ + #ifdef EIGEN_EXCEPTIONS + throw std::bad_alloc(); + #else + std::size_t huge = -1; + new int[huge]; + #endif +} + /***************************************************************************** *** Implementation of handmade aligned functions *** *****************************************************************************/ @@ -192,7 +202,7 @@ inline void check_that_malloc_is_allowed() #endif /** \internal Allocates \a size bytes. The returned pointer is guaranteed to have 16 bytes alignment. - * On allocation error, the returned pointer is null, and if exceptions are enabled then a std::bad_alloc is thrown. + * On allocation error, the returned pointer is null, and std::bad_alloc is thrown. */ inline void* aligned_malloc(size_t size) { @@ -213,10 +223,9 @@ inline void* aligned_malloc(size_t size) result = handmade_aligned_malloc(size); #endif - #ifdef EIGEN_EXCEPTIONS - if(result == 0) - throw std::bad_alloc(); - #endif + if(!result && size) + throw_std_bad_alloc(); + return result; } @@ -241,7 +250,7 @@ inline void aligned_free(void *ptr) /** * \internal * \brief Reallocates an aligned block of memory. -* \throws std::bad_alloc if EIGEN_EXCEPTIONS are defined. +* \throws std::bad_alloc on allocation failure **/ inline void* aligned_realloc(void *ptr, size_t new_size, size_t old_size) { @@ -269,10 +278,9 @@ inline void* aligned_realloc(void *ptr, size_t new_size, size_t old_size) result = handmade_aligned_realloc(ptr,new_size,old_size); #endif -#ifdef EIGEN_EXCEPTIONS - if (result==0 && new_size!=0) - throw std::bad_alloc(); -#endif + if (!result && new_size) + throw_std_bad_alloc(); + return result; } @@ -281,7 +289,7 @@ inline void* aligned_realloc(void *ptr, size_t new_size, size_t old_size) *****************************************************************************/ /** \internal Allocates \a size bytes. If Align is true, then the returned ptr is 16-byte-aligned. - * On allocation error, the returned pointer is null, and if exceptions are enabled then a std::bad_alloc is thrown. + * On allocation error, the returned pointer is null, and a std::bad_alloc is thrown. */ template inline void* conditional_aligned_malloc(size_t size) { @@ -293,9 +301,8 @@ template<> inline void* conditional_aligned_malloc(size_t size) check_that_malloc_is_allowed(); void *result = std::malloc(size); - #ifdef EIGEN_EXCEPTIONS - if(!result) throw std::bad_alloc(); - #endif + if(!result && size) + throw_std_bad_alloc(); return result; } @@ -347,18 +354,27 @@ template inline void destruct_elements_of_array(T *ptr, size_t size) *** Implementation of aligned new/delete-like functions *** *****************************************************************************/ +template +EIGEN_ALWAYS_INLINE void check_size_for_overflow(size_t size) +{ + if(size > size_t(-1) / sizeof(T)) + throw_std_bad_alloc(); +} + /** \internal Allocates \a size objects of type T. The returned pointer is guaranteed to have 16 bytes alignment. - * On allocation error, the returned pointer is undefined, but if exceptions are enabled then a std::bad_alloc is thrown. + * On allocation error, the returned pointer is undefined, but a std::bad_alloc is thrown. * The default constructor of T is called. */ template inline T* aligned_new(size_t size) { + check_size_for_overflow(size); T *result = reinterpret_cast(aligned_malloc(sizeof(T)*size)); return construct_elements_of_array(result, size); } template inline T* conditional_aligned_new(size_t size) { + check_size_for_overflow(size); T *result = reinterpret_cast(conditional_aligned_malloc(sizeof(T)*size)); return construct_elements_of_array(result, size); } @@ -383,6 +399,8 @@ template inline void conditional_aligned_delete(T *ptr, template inline T* conditional_aligned_realloc_new(T* pts, size_t new_size, size_t old_size) { + check_size_for_overflow(new_size); + check_size_for_overflow(old_size); if(new_size < old_size) destruct_elements_of_array(pts+new_size, old_size-new_size); T *result = reinterpret_cast(conditional_aligned_realloc(reinterpret_cast(pts), sizeof(T)*new_size, sizeof(T)*old_size)); @@ -394,6 +412,7 @@ template inline T* conditional_aligned_realloc_new(T* pt template inline T* conditional_aligned_new_auto(size_t size) { + check_size_for_overflow(size); T *result = reinterpret_cast(conditional_aligned_malloc(sizeof(T)*size)); if(NumTraits::RequireInitialization) construct_elements_of_array(result, size); @@ -402,6 +421,8 @@ template inline T* conditional_aligned_new_auto(size_t s template inline T* conditional_aligned_realloc_new_auto(T* pts, size_t new_size, size_t old_size) { + check_size_for_overflow(new_size); + check_size_for_overflow(old_size); if(NumTraits::RequireInitialization && (new_size < old_size)) destruct_elements_of_array(pts+new_size, old_size-new_size); T *result = reinterpret_cast(conditional_aligned_realloc(reinterpret_cast(pts), sizeof(T)*new_size, sizeof(T)*old_size)); @@ -536,6 +557,7 @@ template class aligned_stack_memory_handler #endif #define ei_declare_aligned_stack_constructed_variable(TYPE,NAME,SIZE,BUFFER) \ + Eigen::internal::check_size_for_overflow(SIZE); \ TYPE* NAME = (BUFFER)!=0 ? (BUFFER) \ : reinterpret_cast( \ (sizeof(TYPE)*SIZE<=EIGEN_STACK_ALLOCATION_LIMIT) ? EIGEN_ALIGNED_ALLOCA(sizeof(TYPE)*SIZE) \ @@ -545,6 +567,7 @@ template class aligned_stack_memory_handler #else #define ei_declare_aligned_stack_constructed_variable(TYPE,NAME,SIZE,BUFFER) \ + Eigen::internal::check_size_for_overflow(SIZE); \ TYPE* NAME = (BUFFER)!=0 ? BUFFER : reinterpret_cast(Eigen::internal::aligned_malloc(sizeof(TYPE)*SIZE)); \ Eigen::internal::aligned_stack_memory_handler EIGEN_CAT(NAME,_stack_memory_destructor)((BUFFER)==0 ? NAME : 0,SIZE,true) @@ -669,6 +692,7 @@ public: pointer allocate( size_type num, const void* hint = 0 ) { EIGEN_UNUSED_VARIABLE(hint); + internal::check_size_for_overflow(num); return static_cast( internal::aligned_malloc( num * sizeof(T) ) ); } diff --git a/eigenlib/Eigen/src/Core/util/XprHelper.h b/eigenlib/Eigen/src/Core/util/XprHelper.h index 9047c5f8..c2078f13 100644 --- a/eigenlib/Eigen/src/Core/util/XprHelper.h +++ b/eigenlib/Eigen/src/Core/util/XprHelper.h @@ -125,10 +125,9 @@ class compute_matrix_flags aligned_bit = ( ((Options&DontAlign)==0) - && packet_traits::Vectorizable && ( #if EIGEN_ALIGN_STATICALLY - ((!is_dynamic_size_storage) && (((MaxCols*MaxRows) % packet_traits::size) == 0)) + ((!is_dynamic_size_storage) && (((MaxCols*MaxRows*int(sizeof(Scalar))) % 16) == 0)) #else 0 #endif diff --git a/eigenlib/Eigen/src/Eigen2Support/Geometry/AlignedBox.h b/eigenlib/Eigen/src/Eigen2Support/Geometry/AlignedBox.h index 6ad846d3..78df29d4 100644 --- a/eigenlib/Eigen/src/Eigen2Support/Geometry/AlignedBox.h +++ b/eigenlib/Eigen/src/Eigen2Support/Geometry/AlignedBox.h @@ -51,19 +51,19 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim== { if (AmbientDimAtCompileTime!=Dynamic) setNull(); } /** Constructs a null box with \a _dim the dimension of the ambient space. */ - inline explicit AlignedBox(int _dim) : m_(min)(_dim), m_(max)(_dim) + inline explicit AlignedBox(int _dim) : m_min(_dim), m_max(_dim) { setNull(); } /** Constructs a box with extremities \a _min and \a _max. */ - inline AlignedBox(const VectorType& _min, const VectorType& _max) : m_(min)(_min), m_(max)(_max) {} + inline AlignedBox(const VectorType& _min, const VectorType& _max) : m_min(_min), m_max(_max) {} /** Constructs a box containing a single point \a p. */ - inline explicit AlignedBox(const VectorType& p) : m_(min)(p), m_(max)(p) {} + inline explicit AlignedBox(const VectorType& p) : m_min(p), m_max(p) {} ~AlignedBox() {} /** \returns the dimension in which the box holds */ - inline int dim() const { return AmbientDimAtCompileTime==Dynamic ? m_min.size()-1 : AmbientDimAtCompileTime; } + inline int dim() const { return AmbientDimAtCompileTime==Dynamic ? m_min.size()-1 : int(AmbientDimAtCompileTime); } /** \returns true if the box is null, i.e, empty. */ inline bool isNull() const { return (m_min.cwise() > m_max).any(); } @@ -71,8 +71,8 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim== /** Makes \c *this a null/empty box. */ inline void setNull() { - m_min.setConstant( std::numeric_limits::(max)()); - m_max.setConstant(-std::numeric_limits::(max)()); + m_min.setConstant( (std::numeric_limits::max)()); + m_max.setConstant(-(std::numeric_limits::max)()); } /** \returns the minimal corner */ @@ -90,19 +90,19 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim== /** \returns true if the box \a b is entirely inside the box \c *this. */ inline bool contains(const AlignedBox& b) const - { return (m_min.cwise()<=b.(min)()).all() && (b.(max)().cwise()<=m_max).all(); } + { return (m_min.cwise()<=(b.min)()).all() && ((b.max)().cwise()<=m_max).all(); } /** Extends \c *this such that it contains the point \a p and returns a reference to \c *this. */ inline AlignedBox& extend(const VectorType& p) - { m_min = m_min.cwise().(min)(p); m_max = m_max.cwise().(max)(p); return *this; } + { m_min = (m_min.cwise().min)(p); m_max = (m_max.cwise().max)(p); return *this; } /** Extends \c *this such that it contains the box \a b and returns a reference to \c *this. */ inline AlignedBox& extend(const AlignedBox& b) - { m_min = m_min.cwise().(min)(b.m_min); m_max = m_max.cwise().(max)(b.m_max); return *this; } + { m_min = (m_min.cwise().min)(b.m_min); m_max = (m_max.cwise().max)(b.m_max); return *this; } /** Clamps \c *this by the box \a b and returns a reference to \c *this. */ inline AlignedBox& clamp(const AlignedBox& b) - { m_min = m_min.cwise().(max)(b.m_min); m_max = m_max.cwise().(min)(b.m_max); return *this; } + { m_min = (m_min.cwise().max)(b.m_min); m_max = (m_max.cwise().min)(b.m_max); return *this; } /** Translate \c *this by the vector \a t and returns a reference to \c *this. */ inline AlignedBox& translate(const VectorType& t) @@ -138,8 +138,8 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim== template inline explicit AlignedBox(const AlignedBox& other) { - m_min = other.(min)().template cast(); - m_max = other.(max)().template cast(); + m_min = (other.min)().template cast(); + m_max = (other.max)().template cast(); } /** \returns \c true if \c *this is approximately equal to \a other, within the precision diff --git a/eigenlib/Eigen/src/Eigenvalues/EigenSolver.h b/eigenlib/Eigen/src/Eigenvalues/EigenSolver.h index ac4c4242..f57353c0 100644 --- a/eigenlib/Eigen/src/Eigenvalues/EigenSolver.h +++ b/eigenlib/Eigen/src/Eigenvalues/EigenSolver.h @@ -291,7 +291,7 @@ template class EigenSolver ComputationInfo info() const { - eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); + eigen_assert(m_isInitialized && "EigenSolver is not initialized."); return m_realSchur.info(); } @@ -339,7 +339,7 @@ typename EigenSolver::EigenvectorsType EigenSolver::eige EigenvectorsType matV(n,n); for (Index j=0; j(); @@ -570,10 +570,13 @@ void EigenSolver::doComputeEigenvectors() } } + + // We handled a pair of complex conjugate eigenvalues, so need to skip them both + n--; } else { - eigen_assert("Internal bug in EigenSolver"); // this should not happen + eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen } } diff --git a/eigenlib/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/eigenlib/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 965dda88..ad107c63 100644 --- a/eigenlib/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/eigenlib/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -307,7 +307,8 @@ template class SelfAdjointEigenSolver /** \brief Maximum number of iterations. * - * Maximum number of iterations allowed for an eigenvalue to converge. + * The algorithm terminates if it does not converge within m_maxIterations * n iterations, where n + * denotes the size of the matrix. This value is currently set to 30 (copied from LAPACK). */ static const int m_maxIterations = 30; @@ -407,7 +408,7 @@ SelfAdjointEigenSolver& SelfAdjointEigenSolver Index end = n-1; Index start = 0; - Index iter = 0; // number of iterations we are working on one element + Index iter = 0; // total number of iterations while (end>0) { @@ -418,15 +419,14 @@ SelfAdjointEigenSolver& SelfAdjointEigenSolver // find the largest unreduced block while (end>0 && m_subdiag[end-1]==0) { - iter = 0; end--; } if (end<=0) break; - // if we spent too many iterations on the current element, we give up + // if we spent too many iterations, we give up iter++; - if(iter > m_maxIterations) break; + if(iter > m_maxIterations * n) break; start = end - 1; while (start>0 && m_subdiag[start-1]!=0) @@ -435,7 +435,7 @@ SelfAdjointEigenSolver& SelfAdjointEigenSolver internal::tridiagonal_qr_step(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n); } - if (iter <= m_maxIterations) + if (iter <= m_maxIterations * n) m_info = Success; else m_info = NoConvergence; diff --git a/eigenlib/Eigen/src/Geometry/Hyperplane.h b/eigenlib/Eigen/src/Geometry/Hyperplane.h index eb0a5877..d85d3e55 100644 --- a/eigenlib/Eigen/src/Geometry/Hyperplane.h +++ b/eigenlib/Eigen/src/Geometry/Hyperplane.h @@ -225,7 +225,7 @@ public: normal() = mat * normal(); else { - eigen_assert("invalid traits value in Hyperplane::transform()"); + eigen_assert(0 && "invalid traits value in Hyperplane::transform()"); } return *this; } diff --git a/eigenlib/Eigen/src/Geometry/Quaternion.h b/eigenlib/Eigen/src/Geometry/Quaternion.h index 2662d60f..9180db67 100644 --- a/eigenlib/Eigen/src/Geometry/Quaternion.h +++ b/eigenlib/Eigen/src/Geometry/Quaternion.h @@ -182,10 +182,9 @@ public: template inline typename internal::cast_return_type >::type cast() const { - return typename internal::cast_return_type >::type( - coeffs().template cast()); + return typename internal::cast_return_type >::type(derived()); } - + #ifdef EIGEN_QUATERNIONBASE_PLUGIN # include EIGEN_QUATERNIONBASE_PLUGIN #endif @@ -225,22 +224,25 @@ struct traits > typedef _Scalar Scalar; typedef Matrix<_Scalar,4,1,_Options> Coefficients; enum{ - IsAligned = bool(EIGEN_ALIGN) && ((int(_Options)&Aligned)==Aligned), + IsAligned = internal::traits::Flags & AlignedBit, Flags = IsAligned ? (AlignedBit | LvalueBit) : LvalueBit }; }; } template -class Quaternion : public QuaternionBase >{ +class Quaternion : public QuaternionBase > +{ typedef QuaternionBase > Base; + enum { IsAligned = internal::traits::IsAligned }; + public: typedef _Scalar Scalar; EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Quaternion) using Base::operator*=; - typedef typename internal::traits >::Coefficients Coefficients; + typedef typename internal::traits::Coefficients Coefficients; typedef typename Base::AngleAxisType AngleAxisType; /** Default constructor leaving the quaternion uninitialized. */ @@ -271,9 +273,16 @@ public: template explicit inline Quaternion(const MatrixBase& other) { *this = other; } + /** Explicit copy constructor with scalar conversion */ + template + explicit inline Quaternion(const Quaternion& other) + { m_coeffs = other.coeffs().template cast(); } + inline Coefficients& coeffs() { return m_coeffs;} inline const Coefficients& coeffs() const { return m_coeffs;} + EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(IsAligned) + protected: Coefficients m_coeffs; @@ -673,7 +682,7 @@ QuaternionBase::slerp(Scalar t, const QuaternionBase& oth Scalar scale0; Scalar scale1; - if (absD>=one) + if(absD>=one) { scale0 = Scalar(1) - t; scale1 = t; @@ -686,9 +695,8 @@ QuaternionBase::slerp(Scalar t, const QuaternionBase& oth scale0 = internal::sin( ( Scalar(1) - t ) * theta) / sinTheta; scale1 = internal::sin( ( t * theta) ) / sinTheta; - if (d<0) - scale1 = -scale1; } + if(d<0) scale1 = -scale1; return Quaternion(scale0 * coeffs() + scale1 * other.coeffs()); } diff --git a/eigenlib/Eigen/src/Geometry/Rotation2D.h b/eigenlib/Eigen/src/Geometry/Rotation2D.h index e1214bd3..cf36da1c 100644 --- a/eigenlib/Eigen/src/Geometry/Rotation2D.h +++ b/eigenlib/Eigen/src/Geometry/Rotation2D.h @@ -89,7 +89,7 @@ public: /** Concatenates two rotations */ inline Rotation2D& operator*=(const Rotation2D& other) - { return m_angle += other.m_angle; return *this; } + { m_angle += other.m_angle; return *this; } /** Applies the rotation to a 2D vector */ Vector2 operator* (const Vector2& vec) const diff --git a/eigenlib/Eigen/src/LU/FullPivLU.h b/eigenlib/Eigen/src/LU/FullPivLU.h index 633fb23f..46ae7d65 100644 --- a/eigenlib/Eigen/src/LU/FullPivLU.h +++ b/eigenlib/Eigen/src/LU/FullPivLU.h @@ -443,7 +443,6 @@ FullPivLU& FullPivLU::compute(const MatrixType& matrix) m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) m_maxpivot = RealScalar(0); - RealScalar cutoff(0); for(Index k = 0; k < size; ++k) { @@ -458,14 +457,7 @@ FullPivLU& FullPivLU::compute(const MatrixType& matrix) row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner, col_of_biggest_in_corner += k; // need to add k to them. - // when k==0, biggest_in_corner is the biggest coeff absolute value in the original matrix - if(k == 0) cutoff = biggest_in_corner * NumTraits::epsilon(); - - // if the pivot (hence the corner) is "zero", terminate to avoid generating nan/inf values. - // Notice that using an exact comparison (biggest_in_corner==0) here, as Golub-van Loan do in - // their pseudo-code, results in numerical instability! The cutoff here has been validated - // by running the unit test 'lu' with many repetitions. - if(biggest_in_corner < cutoff) + if(biggest_in_corner==RealScalar(0)) { // before exiting, make sure to initialize the still uninitialized transpositions // in a sane state without destroying what we already have. diff --git a/eigenlib/Eigen/src/LU/arch/Inverse_SSE.h b/eigenlib/Eigen/src/LU/arch/Inverse_SSE.h index 176c349c..4c6153f0 100644 --- a/eigenlib/Eigen/src/LU/arch/Inverse_SSE.h +++ b/eigenlib/Eigen/src/LU/arch/Inverse_SSE.h @@ -55,7 +55,7 @@ struct compute_inverse_size4 static void run(const MatrixType& matrix, ResultType& result) { - EIGEN_ALIGN16 const int _Sign_PNNP[4] = { 0x00000000, 0x80000000, 0x80000000, 0x00000000 }; + EIGEN_ALIGN16 const unsigned int _Sign_PNNP[4] = { 0x00000000, 0x80000000, 0x80000000, 0x00000000 }; // Load the full matrix into registers __m128 _L1 = matrix.template packet( 0); diff --git a/eigenlib/Eigen/src/SVD/JacobiSVD.h b/eigenlib/Eigen/src/SVD/JacobiSVD.h index 5f613999..3c423095 100644 --- a/eigenlib/Eigen/src/SVD/JacobiSVD.h +++ b/eigenlib/Eigen/src/SVD/JacobiSVD.h @@ -590,6 +590,9 @@ JacobiSVD::compute(const MatrixType& matrix, unsig // only worsening the precision of U and V as we accumulate more rotations const RealScalar precision = RealScalar(2) * NumTraits::epsilon(); + // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286) + const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits::denorm_min(); + /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ if(!internal::qr_preconditioner_impl::run(*this, matrix) @@ -617,10 +620,11 @@ JacobiSVD::compute(const MatrixType& matrix, unsig { // if this 2x2 sub-matrix is not diagonal already... // notice that this comparison will evaluate to false if any NaN is involved, ensuring that NaN's don't - // keep us iterating forever. + // keep us iterating forever. Similarly, small denormal numbers are considered zero. using std::max; - if((max)(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p))) - > (max)(internal::abs(m_workMatrix.coeff(p,p)),internal::abs(m_workMatrix.coeff(q,q)))*precision) + RealScalar threshold = (max)(considerAsZero, precision * (max)(internal::abs(m_workMatrix.coeff(p,p)), + internal::abs(m_workMatrix.coeff(q,q)))); + if((max)(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p))) > threshold) { finished = false; @@ -704,6 +708,13 @@ struct solve_retval, Rhs> }; } // end namespace internal +/** \svd_module + * + * \return the singular value decomposition of \c *this computed by two-sided + * Jacobi transformations. + * + * \sa class JacobiSVD + */ template JacobiSVD::PlainObject> MatrixBase::jacobiSvd(unsigned int computationOptions) const diff --git a/eigenlib/Eigen/src/Sparse/TriangularSolver.h b/eigenlib/Eigen/src/Sparse/TriangularSolver.h index 73468e04..62bb8bb4 100644 --- a/eigenlib/Eigen/src/Sparse/TriangularSolver.h +++ b/eigenlib/Eigen/src/Sparse/TriangularSolver.h @@ -171,7 +171,7 @@ void SparseTriangularView::solveInPlace(MatrixBase::Flags & RowMajorBit }; @@ -298,7 +298,7 @@ void SparseTriangularView::solveInPlace(SparseMatrixBase::Flags & RowMajorBit }; diff --git a/eigenlib/howto.txt b/eigenlib/howto.txt new file mode 100644 index 00000000..fd31fd3d --- /dev/null +++ b/eigenlib/howto.txt @@ -0,0 +1,7 @@ +Current Eigen Version 3.05 (10.02.2012) +To update the lib: +- download Eigen +- unzip it somewhere +- copy the two folders 'Eigen' and 'unsupported' here +- update this file +- commit \ No newline at end of file diff --git a/eigenlib/unsupported/Eigen/FFT b/eigenlib/unsupported/Eigen/FFT index 0b282305..c56bd63d 100644 --- a/eigenlib/unsupported/Eigen/FFT +++ b/eigenlib/unsupported/Eigen/FFT @@ -331,7 +331,7 @@ class FFT // if the vector is strided, then we need to copy it to a packed temporary Matrix tmp; if ( resize_input ) { - size_t ncopy = std::min(src.size(),src.size() + resize_input); + size_t ncopy = (std::min)(src.size(),src.size() + resize_input); tmp.setZero(src.size() + resize_input); if ( realfft && HasFlag(HalfSpectrum) ) { // pad at the Nyquist bin diff --git a/eigenlib/unsupported/Eigen/MPRealSupport b/eigenlib/unsupported/Eigen/MPRealSupport index 10fa23b3..8f239635 100644 --- a/eigenlib/unsupported/Eigen/MPRealSupport +++ b/eigenlib/unsupported/Eigen/MPRealSupport @@ -35,7 +35,8 @@ namespace Eigen { - /** \defgroup MPRealSupport_Module MPFRC++ Support module + /** \ingroup Unsupported_modules + * \defgroup MPRealSupport_Module MPFRC++ Support module * * \code * #include diff --git a/eigenlib/unsupported/Eigen/src/BVH/BVAlgorithms.h b/eigenlib/unsupported/Eigen/src/BVH/BVAlgorithms.h index e6fdf473..d65a9774 100644 --- a/eigenlib/unsupported/Eigen/src/BVH/BVAlgorithms.h +++ b/eigenlib/unsupported/Eigen/src/BVH/BVAlgorithms.h @@ -231,7 +231,7 @@ private: template typename Minimizer::Scalar BVMinimize(const BVH &tree, Minimizer &minimizer) { - return internal::minimize_helper(tree, minimizer, tree.getRootIndex(), std::numeric_limits::max()); + return internal::minimize_helper(tree, minimizer, tree.getRootIndex(), (std::numeric_limits::max)()); } /** Given two BVH's, runs the query on their cartesian product encapsulated by \a minimizer. @@ -264,7 +264,7 @@ typename Minimizer::Scalar BVMinimize(const BVH1 &tree1, const BVH2 &tree2, Mini ObjIter2 oBegin2 = ObjIter2(), oEnd2 = ObjIter2(), oCur2 = ObjIter2(); std::priority_queue, std::greater > todo; //smallest is at the top - Scalar minimum = std::numeric_limits::max(); + Scalar minimum = (std::numeric_limits::max)(); todo.push(std::make_pair(Scalar(), std::make_pair(tree1.getRootIndex(), tree2.getRootIndex()))); while(!todo.empty()) { diff --git a/eigenlib/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/eigenlib/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index cedb1d55..50c0ca84 100644 --- a/eigenlib/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/eigenlib/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -259,7 +259,7 @@ void MatrixExponential::computeUV(float) pade5(m_M); } else { const float maxnorm = 3.925724783138660f; - m_squarings = max(0, (int)ceil(log2(m_l1norm / maxnorm))); + m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm))); MatrixType A = m_M / pow(Scalar(2), Scalar(static_cast(m_squarings))); pade7(A); } @@ -281,7 +281,7 @@ void MatrixExponential::computeUV(double) pade9(m_M); } else { const double maxnorm = 5.371920351148152; - m_squarings = max(0, (int)ceil(log2(m_l1norm / maxnorm))); + m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm))); MatrixType A = m_M / pow(Scalar(2), Scalar(m_squarings)); pade13(A); } diff --git a/eigenlib/unsupported/doc/Overview.dox b/eigenlib/unsupported/doc/Overview.dox index c09f0685..458b507b 100644 --- a/eigenlib/unsupported/doc/Overview.dox +++ b/eigenlib/unsupported/doc/Overview.dox @@ -1,13 +1,22 @@ namespace Eigen { -o /** \mainpage Eigen's unsupported modules +/** \mainpage Eigen's unsupported modules This is the API documentation for Eigen's unsupported modules. These modules are contributions from various users. They are provided "as is", without any support. +Click on the \e Modules tab at the top of this page to get a list of all unsupported modules. + Don't miss the official Eigen documentation. + +\defgroup Unsupported_modules Unsupported modules + +The unsupported modules are contributions from various users. They are +provided "as is", without any support. Nevertheless, some of them are +subject to be included in Eigen in the future. + */ } diff --git a/eigenlib/unsupported/test/BVH.cpp b/eigenlib/unsupported/test/BVH.cpp index d773afb7..e77e84b6 100644 --- a/eigenlib/unsupported/test/BVH.cpp +++ b/eigenlib/unsupported/test/BVH.cpp @@ -90,13 +90,13 @@ struct BallPointStuff //this class provides functions to be both an intersector } double minimumOnVolume(const BoxType &r) { ++calls; return r.squaredExteriorDistance(p); } - double minimumOnObject(const BallType &b) { ++calls; return std::max(0., (b.center - p).squaredNorm() - SQR(b.radius)); } + double minimumOnObject(const BallType &b) { ++calls; return (std::max)(0., (b.center - p).squaredNorm() - SQR(b.radius)); } double minimumOnVolumeVolume(const BoxType &r1, const BoxType &r2) { ++calls; return r1.squaredExteriorDistance(r2); } - double minimumOnVolumeObject(const BoxType &r, const BallType &b) { ++calls; return SQR(std::max(0., r.exteriorDistance(b.center) - b.radius)); } - double minimumOnObjectVolume(const BallType &b, const BoxType &r) { ++calls; return SQR(std::max(0., r.exteriorDistance(b.center) - b.radius)); } - double minimumOnObjectObject(const BallType &b1, const BallType &b2){ ++calls; return SQR(std::max(0., (b1.center - b2.center).norm() - b1.radius - b2.radius)); } + double minimumOnVolumeObject(const BoxType &r, const BallType &b) { ++calls; return SQR((std::max)(0., r.exteriorDistance(b.center) - b.radius)); } + double minimumOnObjectVolume(const BallType &b, const BoxType &r) { ++calls; return SQR((std::max)(0., r.exteriorDistance(b.center) - b.radius)); } + double minimumOnObjectObject(const BallType &b1, const BallType &b2){ ++calls; return SQR((std::max)(0., (b1.center - b2.center).norm() - b1.radius - b2.radius)); } double minimumOnVolumeObject(const BoxType &r, const VectorType &v) { ++calls; return r.squaredExteriorDistance(v); } - double minimumOnObjectObject(const BallType &b, const VectorType &v){ ++calls; return SQR(std::max(0., (b.center - v).norm() - b.radius)); } + double minimumOnObjectObject(const BallType &b, const VectorType &v){ ++calls; return SQR((std::max)(0., (b.center - v).norm() - b.radius)); } VectorType p; int calls; diff --git a/eigenlib/unsupported/test/matrix_exponential.cpp b/eigenlib/unsupported/test/matrix_exponential.cpp index 5ea438c2..996b42a7 100644 --- a/eigenlib/unsupported/test/matrix_exponential.cpp +++ b/eigenlib/unsupported/test/matrix_exponential.cpp @@ -36,7 +36,7 @@ double binom(int n, int k) template double relerr(const MatrixBase& A, const MatrixBase& B) { - return std::sqrt((A - B).cwiseAbs2().sum() / std::min(A.cwiseAbs2().sum(), B.cwiseAbs2().sum())); + return std::sqrt((A - B).cwiseAbs2().sum() / (std::min)(A.cwiseAbs2().sum(), B.cwiseAbs2().sum())); } template diff --git a/eigenlib/unsupported/test/sparse_extra.cpp b/eigenlib/unsupported/test/sparse_extra.cpp index a004995f..b1fd481e 100644 --- a/eigenlib/unsupported/test/sparse_extra.cpp +++ b/eigenlib/unsupported/test/sparse_extra.cpp @@ -67,7 +67,7 @@ template void sparse_extra(const SparseMatrixType& re typedef typename SparseMatrixType::Scalar Scalar; enum { Flags = SparseMatrixType::Flags }; - double density = std::max(8./(rows*cols), 0.01); + double density = (std::max)(8./(rows*cols), 0.01); typedef Matrix DenseMatrix; typedef Matrix DenseVector; Scalar eps = 1e-6; diff --git a/eigenlib/unsupported/test/sparse_ldlt.cpp b/eigenlib/unsupported/test/sparse_ldlt.cpp index 4ceda318..03a26bcd 100644 --- a/eigenlib/unsupported/test/sparse_ldlt.cpp +++ b/eigenlib/unsupported/test/sparse_ldlt.cpp @@ -33,7 +33,7 @@ template void sparse_ldlt(int rows, int cols) { static bool odd = true; odd = !odd; - double density = std::max(8./(rows*cols), 0.01); + double density = (std::max)(8./(rows*cols), 0.01); typedef Matrix DenseMatrix; typedef Matrix DenseVector; diff --git a/eigenlib/unsupported/test/sparse_llt.cpp b/eigenlib/unsupported/test/sparse_llt.cpp index df198cd5..5f8a7ce3 100644 --- a/eigenlib/unsupported/test/sparse_llt.cpp +++ b/eigenlib/unsupported/test/sparse_llt.cpp @@ -31,7 +31,7 @@ template void sparse_llt(int rows, int cols) { - double density = std::max(8./(rows*cols), 0.01); + double density = (std::max)(8./(rows*cols), 0.01); typedef Matrix DenseMatrix; typedef Matrix DenseVector; diff --git a/eigenlib/unsupported/test/sparse_lu.cpp b/eigenlib/unsupported/test/sparse_lu.cpp index 188d291c..d58e85a0 100644 --- a/eigenlib/unsupported/test/sparse_lu.cpp +++ b/eigenlib/unsupported/test/sparse_lu.cpp @@ -35,7 +35,7 @@ template void sparse_lu(int rows, int cols) { - double density = std::max(8./(rows*cols), 0.01); + double density = (std::max)(8./(rows*cols), 0.01); typedef Matrix DenseMatrix; typedef Matrix DenseVector;