From 91227ff6d0d49b67f4b332deeb29052502280981 Mon Sep 17 00:00:00 2001 From: Sergey Sharybin Date: Wed, 9 May 2012 08:33:05 +0000 Subject: Update Eigen3 library, would be needed for some further integraiton. --- extern/Eigen3/Eigen/Core | 2 +- extern/Eigen3/Eigen/SVD | 4 +- extern/Eigen3/Eigen/src/Cholesky/LDLT.h | 8 +-- extern/Eigen3/Eigen/src/Cholesky/LLT.h | 8 +-- extern/Eigen3/Eigen/src/Core/Array.h | 4 +- extern/Eigen3/Eigen/src/Core/Block.h | 4 +- extern/Eigen3/Eigen/src/Core/DenseBase.h | 4 +- extern/Eigen3/Eigen/src/Core/Map.h | 4 +- extern/Eigen3/Eigen/src/Core/MapBase.h | 4 +- extern/Eigen3/Eigen/src/Core/Matrix.h | 4 -- extern/Eigen3/Eigen/src/Core/MatrixBase.h | 3 +- extern/Eigen3/Eigen/src/Core/PlainObjectBase.h | 23 +++++++-- extern/Eigen3/Eigen/src/Core/ProductBase.h | 14 +++--- extern/Eigen3/Eigen/src/Core/Replicate.h | 9 +++- extern/Eigen3/Eigen/src/Core/SolveTriangular.h | 2 +- extern/Eigen3/Eigen/src/Core/arch/NEON/Complex.h | 4 +- .../Eigen3/Eigen/src/Core/arch/NEON/PacketMath.h | 58 ++++++++++++++-------- .../src/Core/products/GeneralBlockPanelKernel.h | 21 ++++---- extern/Eigen3/Eigen/src/Core/util/Macros.h | 36 ++++++++------ extern/Eigen3/Eigen/src/Core/util/Memory.h | 54 ++++++++++++++------ extern/Eigen3/Eigen/src/Core/util/XprHelper.h | 3 +- extern/Eigen3/Eigen/src/Eigenvalues/EigenSolver.h | 9 ++-- .../Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 12 ++--- extern/Eigen3/Eigen/src/Geometry/Hyperplane.h | 2 +- extern/Eigen3/Eigen/src/Geometry/Quaternion.h | 26 ++++++---- extern/Eigen3/Eigen/src/Geometry/Rotation2D.h | 2 +- extern/Eigen3/Eigen/src/Geometry/Transform.h | 31 +++++++++++- .../Eigen3/Eigen/src/Geometry/arch/Geometry_SSE.h | 4 +- extern/Eigen3/Eigen/src/LU/FullPivLU.h | 10 +--- extern/Eigen3/Eigen/src/LU/arch/Inverse_SSE.h | 2 +- extern/Eigen3/Eigen/src/SVD/JacobiSVD.h | 17 +++++-- extern/Eigen3/Eigen/src/Sparse/TriangularSolver.h | 4 +- 32 files changed, 245 insertions(+), 147 deletions(-) (limited to 'extern/Eigen3') diff --git a/extern/Eigen3/Eigen/Core b/extern/Eigen3/Eigen/Core index 6e855427c33..a5025e37ead 100644 --- a/extern/Eigen3/Eigen/Core +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/SVD b/extern/Eigen3/Eigen/SVD index d24471fd724..7c987a9dd36 100644 --- a/extern/Eigen3/Eigen/SVD +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Cholesky/LDLT.h b/extern/Eigen3/Eigen/src/Cholesky/LDLT.h index f47b2ea5669..a19e947a4c6 100644 --- a/extern/Eigen3/Eigen/src/Cholesky/LDLT.h +++ b/extern/Eigen3/Eigen/src/Cholesky/LDLT.h @@ -331,16 +331,16 @@ template<> struct ldlt_inplace template struct LDLT_Traits { - typedef TriangularView MatrixL; - typedef TriangularView MatrixU; + typedef const TriangularView MatrixL; + typedef const TriangularView MatrixU; inline static MatrixL getL(const MatrixType& m) { return m; } inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); } }; template struct LDLT_Traits { - typedef TriangularView MatrixL; - typedef TriangularView MatrixU; + typedef const TriangularView MatrixL; + typedef const TriangularView MatrixU; inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); } inline static MatrixU getU(const MatrixType& m) { return m; } }; diff --git a/extern/Eigen3/Eigen/src/Cholesky/LLT.h b/extern/Eigen3/Eigen/src/Cholesky/LLT.h index a4ee5b11cb9..3bb76b5787f 100644 --- a/extern/Eigen3/Eigen/src/Cholesky/LLT.h +++ b/extern/Eigen3/Eigen/src/Cholesky/LLT.h @@ -274,8 +274,8 @@ template<> struct llt_inplace template struct LLT_Traits { - typedef TriangularView MatrixL; - typedef TriangularView MatrixU; + typedef const TriangularView MatrixL; + typedef const TriangularView MatrixU; inline static MatrixL getL(const MatrixType& m) { return m; } inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); } static bool inplace_decomposition(MatrixType& m) @@ -284,8 +284,8 @@ template struct LLT_Traits template struct LLT_Traits { - typedef TriangularView MatrixL; - typedef TriangularView MatrixU; + typedef const TriangularView MatrixL; + typedef const TriangularView MatrixU; inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); } inline static MatrixU getU(const MatrixType& m) { return m; } static bool inplace_decomposition(MatrixType& m) diff --git a/extern/Eigen3/Eigen/src/Core/Array.h b/extern/Eigen3/Eigen/src/Core/Array.h index a3a2167ad3e..a11fb1b53d5 100644 --- a/extern/Eigen3/Eigen/src/Core/Array.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Core/Block.h b/extern/Eigen3/Eigen/src/Core/Block.h index 2b251bc2ca9..d470bc13400 100644 --- a/extern/Eigen3/Eigen/src/Core/Block.h +++ b/extern/Eigen3/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, @@ -342,7 +342,7 @@ class Block } const typename XprType::Nested m_xpr; - int m_outerStride; + Index m_outerStride; }; diff --git a/extern/Eigen3/Eigen/src/Core/DenseBase.h b/extern/Eigen3/Eigen/src/Core/DenseBase.h index 838fa40307a..920904f243a 100644 --- a/extern/Eigen3/Eigen/src/Core/DenseBase.h +++ b/extern/Eigen3/Eigen/src/Core/DenseBase.h @@ -169,8 +169,8 @@ template class DenseBase IsRowMajor = int(Flags) & RowMajorBit, /**< True if this expression has row-major storage order. */ - InnerSizeAtCompileTime = int(IsVectorAtCompileTime) ? SizeAtCompileTime - : int(IsRowMajor) ? ColsAtCompileTime : RowsAtCompileTime, + InnerSizeAtCompileTime = int(IsVectorAtCompileTime) ? int(SizeAtCompileTime) + : int(IsRowMajor) ? int(ColsAtCompileTime) : int(RowsAtCompileTime), CoeffReadCost = internal::traits::CoeffReadCost, /**< This is a rough measure of how expensive it is to read one coefficient from diff --git a/extern/Eigen3/Eigen/src/Core/Map.h b/extern/Eigen3/Eigen/src/Core/Map.h index dd0673609c5..2bf80b3af3d 100644 --- a/extern/Eigen3/Eigen/src/Core/Map.h +++ b/extern/Eigen3/Eigen/src/Core/Map.h @@ -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/extern/Eigen3/Eigen/src/Core/MapBase.h b/extern/Eigen3/Eigen/src/Core/MapBase.h index c23bcbfdcca..9426e2d24dd 100644 --- a/extern/Eigen3/Eigen/src/Core/MapBase.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Core/Matrix.h b/extern/Eigen3/Eigen/src/Core/Matrix.h index 44de22cb4d5..982c9256af0 100644 --- a/extern/Eigen3/Eigen/src/Core/Matrix.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Core/MatrixBase.h b/extern/Eigen3/Eigen/src/Core/MatrixBase.h index db156f6e9d0..62877bce09e 100644 --- a/extern/Eigen3/Eigen/src/Core/MatrixBase.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Core/PlainObjectBase.h b/extern/Eigen3/Eigen/src/Core/PlainObjectBase.h index c70db92479a..612254e9da9 100644 --- a/extern/Eigen3/Eigen/src/Core/PlainObjectBase.h +++ b/extern/Eigen3/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()); } @@ -581,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 } @@ -638,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/extern/Eigen3/Eigen/src/Core/ProductBase.h b/extern/Eigen3/Eigen/src/Core/ProductBase.h index 3bd3487d6a2..91975880fdc 100644 --- a/extern/Eigen3/Eigen/src/Core/ProductBase.h +++ b/extern/Eigen3/Eigen/src/Core/ProductBase.h @@ -152,7 +152,8 @@ class ProductBase : public MatrixBase #else EIGEN_STATIC_ASSERT_SIZE_1x1(Derived) eigen_assert(this->rows() == 1 && this->cols() == 1); - return derived().coeff(row,col); + Matrix result = *this; + return result.coeff(row,col); #endif } @@ -160,7 +161,8 @@ class ProductBase : public MatrixBase { EIGEN_STATIC_ASSERT_SIZE_1x1(Derived) eigen_assert(this->rows() == 1 && this->cols() == 1); - return derived().coeff(i); + Matrix result = *this; + return result.coeff(i); } const Scalar& coeffRef(Index row, Index col) const @@ -256,16 +258,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/extern/Eigen3/Eigen/src/Core/Replicate.h b/extern/Eigen3/Eigen/src/Core/Replicate.h index d2f9712db0f..4c171f8d580 100644 --- a/extern/Eigen3/Eigen/src/Core/Replicate.h +++ b/extern/Eigen3/Eigen/src/Core/Replicate.h @@ -48,7 +48,10 @@ struct traits > typedef typename MatrixType::Scalar Scalar; typedef typename traits::StorageKind StorageKind; typedef typename traits::XprKind XprKind; - typedef typename nested::type MatrixTypeNested; + enum { + Factor = (RowFactor==Dynamic || ColFactor==Dynamic) ? Dynamic : RowFactor*ColFactor + }; + typedef typename nested::type MatrixTypeNested; typedef typename remove_reference::type _MatrixTypeNested; enum { RowsAtCompileTime = RowFactor==Dynamic || int(MatrixType::RowsAtCompileTime)==Dynamic @@ -72,6 +75,8 @@ struct traits > template class Replicate : public internal::dense_xpr_base< Replicate >::type { + typedef typename internal::traits::MatrixTypeNested MatrixTypeNested; + typedef typename internal::traits::_MatrixTypeNested _MatrixTypeNested; public: typedef typename internal::dense_xpr_base::type Base; @@ -124,7 +129,7 @@ template class Replicate protected: - const typename MatrixType::Nested m_matrix; + const MatrixTypeNested m_matrix; const internal::variable_if_dynamic m_rowFactor; const internal::variable_if_dynamic m_colFactor; }; diff --git a/extern/Eigen3/Eigen/src/Core/SolveTriangular.h b/extern/Eigen3/Eigen/src/Core/SolveTriangular.h index 71e129c7f12..a23014a343f 100644 --- a/extern/Eigen3/Eigen/src/Core/SolveTriangular.h +++ b/extern/Eigen3/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); - } -}; - -template -struct palign_impl -{ - EIGEN_STRONG_INLINE static void run(Packet4i& first, const Packet4i& second) - { - if (Offset!=0) - first = vextq_s32(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);\ + }\ +};\ + +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/extern/Eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/extern/Eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h index 6f3f2717007..cd1c37c780e 100644 --- a/extern/Eigen3/Eigen/src/Core/products/GeneralBlockPanelKernel.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Core/util/Macros.h b/extern/Eigen3/Eigen/src/Core/util/Macros.h index 6c3f1e421f0..b7c2b79af92 100644 --- a/extern/Eigen3/Eigen/src/Core/util/Macros.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Core/util/Memory.h b/extern/Eigen3/Eigen/src/Core/util/Memory.h index a580b95ad0d..023716dc9e0 100644 --- a/extern/Eigen3/Eigen/src/Core/util/Memory.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Core/util/XprHelper.h b/extern/Eigen3/Eigen/src/Core/util/XprHelper.h index 9047c5f8350..c2078f13786 100644 --- a/extern/Eigen3/Eigen/src/Core/util/XprHelper.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Eigenvalues/EigenSolver.h b/extern/Eigen3/Eigen/src/Eigenvalues/EigenSolver.h index ac4c4242dd4..f57353c065f 100644 --- a/extern/Eigen3/Eigen/src/Eigenvalues/EigenSolver.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/extern/Eigen3/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 965dda88bda..ad107c63282 100644 --- a/extern/Eigen3/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Geometry/Hyperplane.h b/extern/Eigen3/Eigen/src/Geometry/Hyperplane.h index eb0a5877168..d85d3e553f8 100644 --- a/extern/Eigen3/Eigen/src/Geometry/Hyperplane.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Geometry/Quaternion.h b/extern/Eigen3/Eigen/src/Geometry/Quaternion.h index 2662d60fed1..9180db67d84 100644 --- a/extern/Eigen3/Eigen/src/Geometry/Quaternion.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Geometry/Rotation2D.h b/extern/Eigen3/Eigen/src/Geometry/Rotation2D.h index e1214bd3ebb..cf36da1c50c 100644 --- a/extern/Eigen3/Eigen/src/Geometry/Rotation2D.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Geometry/Transform.h b/extern/Eigen3/Eigen/src/Geometry/Transform.h index 19d012572d4..a694673ebed 100644 --- a/extern/Eigen3/Eigen/src/Geometry/Transform.h +++ b/extern/Eigen3/Eigen/src/Geometry/Transform.h @@ -61,7 +61,7 @@ template< typename Lhs, typename Rhs, bool AnyProjective = transform_traits::IsProjective || - transform_traits::IsProjective> + transform_traits::IsProjective> struct transform_transform_product_impl; template< typename Other, @@ -1391,6 +1391,35 @@ struct transform_transform_product_impl } }; +template +struct transform_transform_product_impl,Transform,true > +{ + typedef Transform Lhs; + typedef Transform Rhs; + typedef Transform ResultType; + static ResultType run(const Lhs& lhs, const Rhs& rhs) + { + ResultType res; + res.matrix().template topRows() = lhs.matrix() * rhs.matrix(); + res.matrix().row(Dim) = rhs.matrix().row(Dim); + return res; + } +}; + +template +struct transform_transform_product_impl,Transform,true > +{ + typedef Transform Lhs; + typedef Transform Rhs; + typedef Transform ResultType; + static ResultType run(const Lhs& lhs, const Rhs& rhs) + { + ResultType res(lhs.matrix().template leftCols() * rhs.matrix()); + res.matrix().col(Dim) += lhs.matrix().col(Dim); + return res; + } +}; + } // end namespace internal #endif // EIGEN_TRANSFORM_H diff --git a/extern/Eigen3/Eigen/src/Geometry/arch/Geometry_SSE.h b/extern/Eigen3/Eigen/src/Geometry/arch/Geometry_SSE.h index cbe695c7259..2af32678d1c 100644 --- a/extern/Eigen3/Eigen/src/Geometry/arch/Geometry_SSE.h +++ b/extern/Eigen3/Eigen/src/Geometry/arch/Geometry_SSE.h @@ -96,7 +96,7 @@ struct quat_product */ t1 = padd(pmul(a_ww, b_xy), pmul(a_yy, b_zw)); t2 = psub(pmul(a_zz, b_xy), pmul(a_xx, b_zw)); -#ifdef __SSE3__ +#ifdef EIGEN_VECTORIZE_SSE3 EIGEN_UNUSED_VARIABLE(mask) pstore(&res.x(), _mm_addsub_pd(t1, preverse(t2))); #else @@ -110,7 +110,7 @@ struct quat_product */ t1 = psub(pmul(a_ww, b_zw), pmul(a_yy, b_xy)); t2 = padd(pmul(a_zz, b_zw), pmul(a_xx, b_xy)); -#ifdef __SSE3__ +#ifdef EIGEN_VECTORIZE_SSE3 EIGEN_UNUSED_VARIABLE(mask) pstore(&res.z(), preverse(_mm_addsub_pd(preverse(t1), t2))); #else diff --git a/extern/Eigen3/Eigen/src/LU/FullPivLU.h b/extern/Eigen3/Eigen/src/LU/FullPivLU.h index 633fb23fdbe..46ae7d651c8 100644 --- a/extern/Eigen3/Eigen/src/LU/FullPivLU.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/LU/arch/Inverse_SSE.h b/extern/Eigen3/Eigen/src/LU/arch/Inverse_SSE.h index 176c349ce44..4c6153f0aff 100644 --- a/extern/Eigen3/Eigen/src/LU/arch/Inverse_SSE.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/SVD/JacobiSVD.h b/extern/Eigen3/Eigen/src/SVD/JacobiSVD.h index 5f61399988c..3c423095c31 100644 --- a/extern/Eigen3/Eigen/src/SVD/JacobiSVD.h +++ b/extern/Eigen3/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/extern/Eigen3/Eigen/src/Sparse/TriangularSolver.h b/extern/Eigen3/Eigen/src/Sparse/TriangularSolver.h index 73468e0446c..62bb8bb44c9 100644 --- a/extern/Eigen3/Eigen/src/Sparse/TriangularSolver.h +++ b/extern/Eigen3/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 }; -- cgit v1.2.3