diff options
author | Sergey Sharybin <sergey.vfx@gmail.com> | 2012-05-09 12:33:05 +0400 |
---|---|---|
committer | Sergey Sharybin <sergey.vfx@gmail.com> | 2012-05-09 12:33:05 +0400 |
commit | 91227ff6d0d49b67f4b332deeb29052502280981 (patch) | |
tree | e7dc3d4cc7baeb759be0f3d8907f6a93107b2fd7 /extern/Eigen3 | |
parent | e132098599a8562bec7bff9aeaa69a7891d50708 (diff) |
Update Eigen3 library, would be needed for some further integraiton.
Diffstat (limited to 'extern/Eigen3')
32 files changed, 245 insertions, 147 deletions
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 <intrin.h> #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 <Eigen/SVD> 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<Upper> template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower> { - typedef TriangularView<MatrixType, UnitLower> MatrixL; - typedef TriangularView<typename MatrixType::AdjointReturnType, UnitUpper> MatrixU; + typedef const TriangularView<const MatrixType, UnitLower> MatrixL; + typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU; inline static MatrixL getL(const MatrixType& m) { return m; } inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); } }; template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper> { - typedef TriangularView<typename MatrixType::AdjointReturnType, UnitLower> MatrixL; - typedef TriangularView<MatrixType, UnitUpper> MatrixU; + typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL; + typedef const TriangularView<const MatrixType, UnitUpper> 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<Upper> template<typename MatrixType> struct LLT_Traits<MatrixType,Lower> { - typedef TriangularView<MatrixType, Lower> MatrixL; - typedef TriangularView<typename MatrixType::AdjointReturnType, Upper> MatrixU; + typedef const TriangularView<const MatrixType, Lower> MatrixL; + typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> 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<typename MatrixType> struct LLT_Traits<MatrixType,Lower> template<typename MatrixType> struct LLT_Traits<MatrixType,Upper> { - typedef TriangularView<typename MatrixType::AdjointReturnType, Lower> MatrixL; - typedef TriangularView<MatrixType, Upper> MatrixU; + typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL; + typedef const TriangularView<const MatrixType, Upper> 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<int>(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<Block<XprType, BlockRows, BlockCols, InnerPanel, HasDirectAccess> MaskPacketAccessBit = (InnerSize == Dynamic || (InnerSize % packet_traits<Scalar>::size) == 0) && (InnerStrideAtCompileTime == 1) ? PacketAccessBit : 0, - MaskAlignedBit = (InnerPanel && (OuterStrideAtCompileTime!=Dynamic) && ((OuterStrideAtCompileTime % packet_traits<Scalar>::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<XprType>::value ? LvalueBit : 0, FlagsRowMajorBit = IsRowMajor ? RowMajorBit : 0, @@ -342,7 +342,7 @@ class Block<XprType,BlockRows,BlockCols, InnerPanel,true> } 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<typename Derived> 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<Derived>::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<Map<PlainObjectType, MapOptions, StrideType> > || HasNoOuterStride || ( OuterStrideAtCompileTime!=Dynamic && ((static_cast<int>(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<typename PlainObjectType, int MapOptions, typename StrideType> class Ma public: typedef MapBase<Map> Base; - EIGEN_DENSE_PUBLIC_INTERFACE(Map) typedef typename Base::PointerType PointerType; @@ -181,7 +180,6 @@ template<typename PlainObjectType, int MapOptions, typename StrideType> 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<typename Derived> class MapBase<Derived, ReadOnlyAccessors> EIGEN_STATIC_ASSERT(EIGEN_IMPLIES(internal::traits<Derived>::Flags&PacketAccessBit, internal::inner_stride_at_compile_time<Derived>::ret==1), PACKET_ACCESS_REQUIRES_TO_HAVE_INNER_STRIDE_FIXED_TO_1); - eigen_assert(EIGEN_IMPLIES(internal::traits<Derived>::Flags&AlignedBit, (size_t(m_data) % (sizeof(Scalar)*internal::packet_traits<Scalar>::size)) == 0) - && "data is not aligned"); + eigen_assert(EIGEN_IMPLIES(internal::traits<Derived>::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<int>(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<typename Derived> class MatrixBase // huuuge hack. make Eigen2's matrix.part<Diagonal>() 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<template<typename T, int n> class U> + // Note: replacing next line by "template<template<typename T, int n> class U>" produces a mysterious error C2082 in MSVC. + template<template<typename, int> class U> const DiagonalWrapper<ConstDiagonalReturnType> 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<typename Index> +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 <typename Derived, typename OtherDerived = Derived, bool IsVector = static_cast<bool>(Derived::IsVectorAtCompileTime)> struct conservative_resize_like_impl; template<typename MatrixTypeA, typename MatrixTypeB, bool SwapPointers> struct matrix_swap_impl; @@ -84,14 +97,12 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::type template<typename StrideType> struct StridedConstMapType { typedef Eigen::Map<const Derived, Unaligned, StrideType> type; }; template<typename StrideType> struct StridedAlignedMapType { typedef Eigen::Map<Derived, Aligned, StrideType> type; }; template<typename StrideType> struct StridedConstAlignedMapType { typedef Eigen::Map<const Derived, Aligned, StrideType> type; }; - protected: DenseStorage<Scalar, Base::MaxSizeAtCompileTime, Base::RowsAtCompileTime, Base::ColsAtCompileTime, Options> m_storage; public: - enum { NeedsToAlign = (!(Options&DontAlign)) - && SizeAtCompileTime!=Dynamic && ((static_cast<int>(sizeof(Scalar))*SizeAtCompileTime)%16)==0 }; + enum { NeedsToAlign = SizeAtCompileTime != Dynamic && (internal::traits<Derived>::Flags & AlignedBit) != 0 }; EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign) Base& base() { return *static_cast<Base*>(this); } @@ -200,11 +211,13 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::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<Derived>::type EIGEN_STRONG_INLINE void resizeLike(const EigenBase<OtherDerived>& _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<Derived>::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<Derived>::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<Derived> #else EIGEN_STATIC_ASSERT_SIZE_1x1(Derived) eigen_assert(this->rows() == 1 && this->cols() == 1); - return derived().coeff(row,col); + Matrix<Scalar,1,1> result = *this; + return result.coeff(row,col); #endif } @@ -160,7 +161,8 @@ class ProductBase : public MatrixBase<Derived> { EIGEN_STATIC_ASSERT_SIZE_1x1(Derived) eigen_assert(this->rows() == 1 && this->cols() == 1); - return derived().coeff(i); + Matrix<Scalar,1,1> 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<typename Dest> - 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<typename Dest> - inline void addTo(Dest& dst) const { scaleAndAddTo(dst,m_alpha); } + inline void addTo(Dest& dst) const { scaleAndAddTo(dst, Scalar(1)); } template<typename Dest> - inline void subTo(Dest& dst) const { scaleAndAddTo(dst,-m_alpha); } + inline void subTo(Dest& dst) const { scaleAndAddTo(dst, Scalar(-1)); } template<typename Dest> - 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<Replicate<MatrixType,RowFactor,ColFactor> > typedef typename MatrixType::Scalar Scalar; typedef typename traits<MatrixType>::StorageKind StorageKind; typedef typename traits<MatrixType>::XprKind XprKind; - typedef typename nested<MatrixType>::type MatrixTypeNested; + enum { + Factor = (RowFactor==Dynamic || ColFactor==Dynamic) ? Dynamic : RowFactor*ColFactor + }; + typedef typename nested<MatrixType,Factor>::type MatrixTypeNested; typedef typename remove_reference<MatrixTypeNested>::type _MatrixTypeNested; enum { RowsAtCompileTime = RowFactor==Dynamic || int(MatrixType::RowsAtCompileTime)==Dynamic @@ -72,6 +75,8 @@ struct traits<Replicate<MatrixType,RowFactor,ColFactor> > template<typename MatrixType,int RowFactor,int ColFactor> class Replicate : public internal::dense_xpr_base< Replicate<MatrixType,RowFactor,ColFactor> >::type { + typedef typename internal::traits<Replicate>::MatrixTypeNested MatrixTypeNested; + typedef typename internal::traits<Replicate>::_MatrixTypeNested _MatrixTypeNested; public: typedef typename internal::dense_xpr_base<Replicate>::type Base; @@ -124,7 +129,7 @@ template<typename MatrixType,int RowFactor,int ColFactor> class Replicate protected: - const typename MatrixType::Nested m_matrix; + const MatrixTypeNested m_matrix; const internal::variable_if_dynamic<Index, RowFactor> m_rowFactor; const internal::variable_if_dynamic<Index, ColFactor> 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<MatrixType,Mode>::solveInPlace(const MatrixBase<OtherDerived eigen_assert(cols() == rows()); eigen_assert( (Side==OnTheLeft && cols() == other.rows()) || (Side==OnTheRight && cols() == other.cols()) ); eigen_assert(!(Mode & ZeroDiag)); - eigen_assert(Mode & (Upper|Lower)); + eigen_assert((Mode & (Upper|Lower)) != 0); enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit && OtherDerived::IsVectorAtCompileTime }; typedef typename internal::conditional<copy, diff --git a/extern/Eigen3/Eigen/src/Core/arch/NEON/Complex.h b/extern/Eigen3/Eigen/src/Core/arch/NEON/Complex.h index 8e55548c946..212887184c2 100644 --- a/extern/Eigen3/Eigen/src/Core/arch/NEON/Complex.h +++ b/extern/Eigen3/Eigen/src/Core/arch/NEON/Complex.h @@ -27,8 +27,8 @@ namespace internal { -static uint32x4_t p4ui_CONJ_XOR = { 0x00000000, 0x80000000, 0x00000000, 0x80000000 }; -static uint32x2_t p2ui_CONJ_XOR = { 0x00000000, 0x80000000 }; +static uint32x4_t p4ui_CONJ_XOR = EIGEN_INIT_NEON_PACKET4(0x00000000, 0x80000000, 0x00000000, 0x80000000); +static uint32x2_t p2ui_CONJ_XOR = EIGEN_INIT_NEON_PACKET2(0x00000000, 0x80000000); //---------- float ---------- struct Packet2cf diff --git a/extern/Eigen3/Eigen/src/Core/arch/NEON/PacketMath.h b/extern/Eigen3/Eigen/src/Core/arch/NEON/PacketMath.h index 478ef8038c0..6c7cd159097 100644 --- a/extern/Eigen3/Eigen/src/Core/arch/NEON/PacketMath.h +++ b/extern/Eigen3/Eigen/src/Core/arch/NEON/PacketMath.h @@ -52,6 +52,16 @@ typedef uint32x4_t Packet4ui; #define _EIGEN_DECLARE_CONST_Packet4i(NAME,X) \ const Packet4i p4i_##NAME = pset1<Packet4i>(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<int> : 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<Packet4i>(const int& from) { template<> EIGEN_STRONG_INLINE Packet4f plset<float>(const float& a) { - Packet4f countdown = { 0, 1, 2, 3 }; + Packet4f countdown = EIGEN_INIT_NEON_PACKET4(0, 1, 2, 3); return vaddq_f32(pset1<Packet4f>(a), countdown); } template<> EIGEN_STRONG_INLINE Packet4i plset<int>(const int& a) { - Packet4i countdown = { 0, 1, 2, 3 }; + Packet4i countdown = EIGEN_INIT_NEON_PACKET4(0, 1, 2, 3); return vaddq_s32(pset1<Packet4i>(a), countdown); } @@ -395,25 +405,29 @@ template<> EIGEN_STRONG_INLINE int predux_max<Packet4i>(const Packet4i& a) return s[0]; } -template<int Offset> -struct palign_impl<Offset,Packet4f> -{ - EIGEN_STRONG_INLINE static void run(Packet4f& first, const Packet4f& second) - { - if (Offset!=0) - first = vextq_f32(first, second, Offset); - } -}; - -template<int Offset> -struct palign_impl<Offset,Packet4i> -{ - 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<Offset,Type>\ +{\ + 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<typename _LhsScalar, typename _RhsScalar, bool _ConjLhs=false, bool _ConjRhs=false> 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<typename CJ, typename A, typename B, typename C, typename T> 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<typename CJ, typename T> struct gebp_madd_selector<CJ,T,T,T,T> { - 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<Derived> Eigen::MatrixBase<Scalar, Derived>::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<Derived> Eigen::MatrixBase<Scalar, Derived>::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<bool Align> inline void* conditional_aligned_malloc(size_t size) { @@ -293,9 +301,8 @@ template<> inline void* conditional_aligned_malloc<false>(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<typename T> inline void destruct_elements_of_array(T *ptr, size_t size) *** Implementation of aligned new/delete-like functions *** *****************************************************************************/ +template<typename T> +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<typename T> inline T* aligned_new(size_t size) { + check_size_for_overflow<T>(size); T *result = reinterpret_cast<T*>(aligned_malloc(sizeof(T)*size)); return construct_elements_of_array(result, size); } template<typename T, bool Align> inline T* conditional_aligned_new(size_t size) { + check_size_for_overflow<T>(size); T *result = reinterpret_cast<T*>(conditional_aligned_malloc<Align>(sizeof(T)*size)); return construct_elements_of_array(result, size); } @@ -383,6 +399,8 @@ template<typename T, bool Align> inline void conditional_aligned_delete(T *ptr, template<typename T, bool Align> inline T* conditional_aligned_realloc_new(T* pts, size_t new_size, size_t old_size) { + check_size_for_overflow<T>(new_size); + check_size_for_overflow<T>(old_size); if(new_size < old_size) destruct_elements_of_array(pts+new_size, old_size-new_size); T *result = reinterpret_cast<T*>(conditional_aligned_realloc<Align>(reinterpret_cast<void*>(pts), sizeof(T)*new_size, sizeof(T)*old_size)); @@ -394,6 +412,7 @@ template<typename T, bool Align> inline T* conditional_aligned_realloc_new(T* pt template<typename T, bool Align> inline T* conditional_aligned_new_auto(size_t size) { + check_size_for_overflow<T>(size); T *result = reinterpret_cast<T*>(conditional_aligned_malloc<Align>(sizeof(T)*size)); if(NumTraits<T>::RequireInitialization) construct_elements_of_array(result, size); @@ -402,6 +421,8 @@ template<typename T, bool Align> inline T* conditional_aligned_new_auto(size_t s template<typename T, bool Align> inline T* conditional_aligned_realloc_new_auto(T* pts, size_t new_size, size_t old_size) { + check_size_for_overflow<T>(new_size); + check_size_for_overflow<T>(old_size); if(NumTraits<T>::RequireInitialization && (new_size < old_size)) destruct_elements_of_array(pts+new_size, old_size-new_size); T *result = reinterpret_cast<T*>(conditional_aligned_realloc<Align>(reinterpret_cast<void*>(pts), sizeof(T)*new_size, sizeof(T)*old_size)); @@ -536,6 +557,7 @@ template<typename T> class aligned_stack_memory_handler #endif #define ei_declare_aligned_stack_constructed_variable(TYPE,NAME,SIZE,BUFFER) \ + Eigen::internal::check_size_for_overflow<TYPE>(SIZE); \ TYPE* NAME = (BUFFER)!=0 ? (BUFFER) \ : reinterpret_cast<TYPE*>( \ (sizeof(TYPE)*SIZE<=EIGEN_STACK_ALLOCATION_LIMIT) ? EIGEN_ALIGNED_ALLOCA(sizeof(TYPE)*SIZE) \ @@ -545,6 +567,7 @@ template<typename T> class aligned_stack_memory_handler #else #define ei_declare_aligned_stack_constructed_variable(TYPE,NAME,SIZE,BUFFER) \ + Eigen::internal::check_size_for_overflow<TYPE>(SIZE); \ TYPE* NAME = (BUFFER)!=0 ? BUFFER : reinterpret_cast<TYPE*>(Eigen::internal::aligned_malloc(sizeof(TYPE)*SIZE)); \ Eigen::internal::aligned_stack_memory_handler<TYPE> 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<T>(num); return static_cast<pointer>( 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<Scalar>::Vectorizable && ( #if EIGEN_ALIGN_STATICALLY - ((!is_dynamic_size_storage) && (((MaxCols*MaxRows) % packet_traits<Scalar>::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<typename _MatrixType> 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<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eige EigenvectorsType matV(n,n); for (Index j=0; j<n; ++j) { - if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(j)), internal::real(m_eivalues.coeff(j)))) + if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(j)), internal::real(m_eivalues.coeff(j))) || j+1==n) { // we have a real eigen value matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>(); @@ -570,10 +570,13 @@ void EigenSolver<MatrixType>::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<typename _MatrixType> 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<MatrixType>& SelfAdjointEigenSolver<MatrixType> 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<MatrixType>& SelfAdjointEigenSolver<MatrixType> // 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<MatrixType>& SelfAdjointEigenSolver<MatrixType> internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(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<typename NewScalarType> inline typename internal::cast_return_type<Derived,Quaternion<NewScalarType> >::type cast() const { - return typename internal::cast_return_type<Derived,Quaternion<NewScalarType> >::type( - coeffs().template cast<NewScalarType>()); + return typename internal::cast_return_type<Derived,Quaternion<NewScalarType> >::type(derived()); } - + #ifdef EIGEN_QUATERNIONBASE_PLUGIN # include EIGEN_QUATERNIONBASE_PLUGIN #endif @@ -225,22 +224,25 @@ struct traits<Quaternion<_Scalar,_Options> > typedef _Scalar Scalar; typedef Matrix<_Scalar,4,1,_Options> Coefficients; enum{ - IsAligned = bool(EIGEN_ALIGN) && ((int(_Options)&Aligned)==Aligned), + IsAligned = internal::traits<Coefficients>::Flags & AlignedBit, Flags = IsAligned ? (AlignedBit | LvalueBit) : LvalueBit }; }; } template<typename _Scalar, int _Options> -class Quaternion : public QuaternionBase<Quaternion<_Scalar,_Options> >{ +class Quaternion : public QuaternionBase<Quaternion<_Scalar,_Options> > +{ typedef QuaternionBase<Quaternion<_Scalar,_Options> > Base; + enum { IsAligned = internal::traits<Quaternion>::IsAligned }; + public: typedef _Scalar Scalar; EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Quaternion) using Base::operator*=; - typedef typename internal::traits<Quaternion<Scalar,_Options> >::Coefficients Coefficients; + typedef typename internal::traits<Quaternion>::Coefficients Coefficients; typedef typename Base::AngleAxisType AngleAxisType; /** Default constructor leaving the quaternion uninitialized. */ @@ -271,9 +273,16 @@ public: template<typename Derived> explicit inline Quaternion(const MatrixBase<Derived>& other) { *this = other; } + /** Explicit copy constructor with scalar conversion */ + template<typename OtherScalar, int OtherOptions> + explicit inline Quaternion(const Quaternion<OtherScalar, OtherOptions>& other) + { m_coeffs = other.coeffs().template cast<Scalar>(); } + 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<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& oth Scalar scale0; Scalar scale1; - if (absD>=one) + if(absD>=one) { scale0 = Scalar(1) - t; scale1 = t; @@ -686,9 +695,8 @@ QuaternionBase<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& 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<Scalar>(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<Lhs>::IsProjective || - transform_traits<Lhs>::IsProjective> + transform_traits<Rhs>::IsProjective> struct transform_transform_product_impl; template< typename Other, @@ -1391,6 +1391,35 @@ struct transform_transform_product_impl<Transform<Scalar,Dim,LhsMode,LhsOptions> } }; +template<typename Scalar, int Dim, int LhsOptions, int RhsOptions> +struct transform_transform_product_impl<Transform<Scalar,Dim,AffineCompact,LhsOptions>,Transform<Scalar,Dim,Projective,RhsOptions>,true > +{ + typedef Transform<Scalar,Dim,AffineCompact,LhsOptions> Lhs; + typedef Transform<Scalar,Dim,Projective,RhsOptions> Rhs; + typedef Transform<Scalar,Dim,Projective> ResultType; + static ResultType run(const Lhs& lhs, const Rhs& rhs) + { + ResultType res; + res.matrix().template topRows<Dim>() = lhs.matrix() * rhs.matrix(); + res.matrix().row(Dim) = rhs.matrix().row(Dim); + return res; + } +}; + +template<typename Scalar, int Dim, int LhsOptions, int RhsOptions> +struct transform_transform_product_impl<Transform<Scalar,Dim,Projective,LhsOptions>,Transform<Scalar,Dim,AffineCompact,RhsOptions>,true > +{ + typedef Transform<Scalar,Dim,Projective,LhsOptions> Lhs; + typedef Transform<Scalar,Dim,AffineCompact,RhsOptions> Rhs; + typedef Transform<Scalar,Dim,Projective> ResultType; + static ResultType run(const Lhs& lhs, const Rhs& rhs) + { + ResultType res(lhs.matrix().template leftCols<Dim>() * 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<Architecture::SSE, Derived, OtherDerived, double, Aligned> */ 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<Architecture::SSE, Derived, OtherDerived, double, Aligned> */ 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<MatrixType>& FullPivLU<MatrixType>::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<MatrixType>& FullPivLU<MatrixType>::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<Scalar>::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<Architecture::SSE, float, MatrixType, ResultType> 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<MatrixAlignment>( 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<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig // only worsening the precision of U and V as we accumulate more rotations const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::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<RealScalar>::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<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows>::run(*this, matrix) @@ -617,10 +620,11 @@ JacobiSVD<MatrixType, QRPreconditioner>::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<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs> }; } // end namespace internal +/** \svd_module + * + * \return the singular value decomposition of \c *this computed by two-sided + * Jacobi transformations. + * + * \sa class JacobiSVD + */ template<typename Derived> JacobiSVD<typename MatrixBase<Derived>::PlainObject> MatrixBase<Derived>::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<ExpressionType,Mode>::solveInPlace(MatrixBase<OtherDer eigen_assert(m_matrix.cols() == m_matrix.rows()); eigen_assert(m_matrix.cols() == other.rows()); eigen_assert(!(Mode & ZeroDiag)); - eigen_assert(Mode & (Upper|Lower)); + eigen_assert((Mode & (Upper|Lower)) != 0); enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit }; @@ -298,7 +298,7 @@ void SparseTriangularView<ExpressionType,Mode>::solveInPlace(SparseMatrixBase<Ot eigen_assert(m_matrix.cols() == m_matrix.rows()); eigen_assert(m_matrix.cols() == other.rows()); eigen_assert(!(Mode & ZeroDiag)); - eigen_assert(Mode & (Upper|Lower)); + eigen_assert((Mode & (Upper|Lower)) != 0); // enum { copy = internal::traits<OtherDerived>::Flags & RowMajorBit }; |