bump

It did not work with clang, and I'm unsure how it worked for gcc/msvc since QuaternionBase was introduced
(transplanted from 84cf1b5b1d
)

CMakeLists.txt

@@ -279,9 +279,21 @@ install(FILES
   )
 
 if(EIGEN_BUILD_PKGCONFIG)
+    SET(path_separator ":")
+    STRING(REPLACE ${path_separator} ";" pkg_config_libdir_search "$ENV{PKG_CONFIG_LIBDIR}")
+    message(STATUS "searching for 'pkgconfig' directory in PKG_CONFIG_LIBDIR ( $ENV{PKG_CONFIG_LIBDIR} ), ${CMAKE_INSTALL_PREFIX}/share, and ${CMAKE_INSTALL_PREFIX}/lib")
+    FIND_PATH(pkg_config_libdir pkgconfig ${pkg_config_libdir_search} ${CMAKE_INSTALL_PREFIX}/share ${CMAKE_INSTALL_PREFIX}/lib ${pkg_config_libdir_search})
+    if(pkg_config_libdir)
+        SET(pkg_config_install_dir ${pkg_config_libdir})
+        message(STATUS "found ${pkg_config_libdir}/pkgconfig" )
+    else(pkg_config_libdir)
+        SET(pkg_config_install_dir ${CMAKE_INSTALL_PREFIX}/share)
+        message(STATUS "pkgconfig not found; installing in ${pkg_config_install_dir}" )
+    endif(pkg_config_libdir)
+
     configure_file(eigen3.pc.in eigen3.pc)
     install(FILES ${CMAKE_CURRENT_BINARY_DIR}/eigen3.pc
-        DESTINATION share/pkgconfig
+        DESTINATION ${pkg_config_install_dir}/pkgconfig
         )
 endif(EIGEN_BUILD_PKGCONFIG)
 
@@ -321,9 +333,9 @@ endif()
 configure_file(${CMAKE_BINARY_DIR}/DartConfiguration.tcl ${CMAKE_BINARY_DIR}/DartConfiguration.tcl)
 # restore default CMAKE_MAKE_PROGRAM
 set(CMAKE_MAKE_PROGRAM ${CMAKE_MAKE_PROGRAM_SAVE})
-# un-set temporary variables so that it is like they never existed. 
+# un-set temporary variables so that it is like they never existed.
 # CMake 2.6.3 introduces the more logical unset() syntax for this.
-set(CMAKE_MAKE_PROGRAM_SAVE) 
+set(CMAKE_MAKE_PROGRAM_SAVE)
 set(EIGEN_MAKECOMMAND_PLACEHOLDER)
 
 configure_file(${CMAKE_SOURCE_DIR}/CTestCustom.cmake.in ${CMAKE_BINARY_DIR}/CTestCustom.cmake)

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
 

Eigen/src/Cholesky/LDLT.h

@@ -158,10 +158,19 @@ template<typename _MatrixType, int _UpLo> class LDLT
     }
 
     /** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A.
+      *
+      * This function also supports in-place solves using the syntax <tt>x = decompositionObject.solve(x)</tt> .
       *
       * \note_about_checking_solutions
       *
-      * \sa solveInPlace(), MatrixBase::ldlt()
+      * More precisely, this method solves \f$ A x = b \f$ using the decomposition \f$ A = P^T L D L^* P \f$
+      * by solving the systems \f$ P^T y_1 = b \f$, \f$ L y_2 = y_1 \f$, \f$ D y_3 = y_2 \f$, 
+      * \f$ L^* y_4 = y_3 \f$ and \f$ P x = y_4 \f$ in succession. If the matrix \f$ A \f$ is singular, then
+      * \f$ D \f$ will also be singular (all the other matrices are invertible). In that case, the
+      * least-square solution of \f$ D y_3 = y_2 \f$ is computed. This does not mean that this function
+      * computes the least-square solution of \f$ A x = b \f$ is \f$ A \f$ is singular.
+      *
+      * \sa MatrixBase::ldlt()
       */
     template<typename Rhs>
     inline const internal::solve_retval<LDLT, Rhs>
@@ -376,7 +385,21 @@ struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
     dec().matrixL().solveInPlace(dst);
 
     // dst = D^-1 (L^-1 P b)
-    dst = dec().vectorD().asDiagonal().inverse() * dst;
+    // more precisely, use pseudo-inverse of D (see bug 241)
+    using std::abs;
+    using std::max;
+    typedef typename LDLTType::MatrixType MatrixType;
+    typedef typename LDLTType::Scalar Scalar;
+    typedef typename LDLTType::RealScalar RealScalar;
+    const Diagonal<const MatrixType> vectorD = dec().vectorD();
+    RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() * NumTraits<Scalar>::epsilon(),
+				 RealScalar(1) / NumTraits<RealScalar>::highest()); // motivated by LAPACK's xGELSS
+    for (Index i = 0; i < vectorD.size(); ++i) {
+      if(abs(vectorD(i)) > tolerance)
+	dst.row(i) /= vectorD(i);
+      else
+	dst.row(i).setZero();
+    }
 
     // dst = L^-T (D^-1 L^-1 P b)
     dec().matrixU().solveInPlace(dst);

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;

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 * sizeof(Scalar)) % 16) == 0)) ? AlignedBit : 0,
     FlagsLinearAccessBit = (RowsAtCompileTime == 1 || ColsAtCompileTime == 1) ? LinearAccessBit : 0,
     FlagsLvalueBit = is_lvalue<XprType>::value ? LvalueBit : 0,
     FlagsRowMajorBit = IsRowMajor ? RowMajorBit : 0,

Eigen/src/Core/Fuzzy.h

@@ -94,7 +94,7 @@ struct isMuchSmallerThan_scalar_selector<Derived, true>
   *
   * \note The fuzzy compares are done multiplicatively. Two vectors \f$ v \f$ and \f$ w \f$
   * are considered to be approximately equal within precision \f$ p \f$ if
-  * \f[ \Vert v - w \Vert \leqslant p\,\(min)(\Vert v\Vert, \Vert w\Vert). \f]
+  * \f[ \Vert v - w \Vert \leqslant p\,\min(\Vert v\Vert, \Vert w\Vert). \f]
   * For matrices, the comparison is done using the Hilbert-Schmidt norm (aka Frobenius norm
   * L2 norm).
   *

Eigen/src/Core/Map.h

@@ -34,7 +34,7 @@
   * \tparam PlainObjectType the equivalent matrix type of the mapped data
   * \tparam MapOptions specifies whether the pointer is \c #Aligned, or \c #Unaligned.
   *                The default is \c #Unaligned.
-  * \tparam StrideType optionnally specifies strides. By default, Map assumes the memory layout
+  * \tparam StrideType optionally specifies strides. By default, Map assumes the memory layout
   *                   of an ordinary, contiguous array. This can be overridden by specifying strides.
   *                   The type passed here must be a specialization of the Stride template, see examples below.
   *
@@ -72,9 +72,9 @@
   * Example: \include Map_placement_new.cpp
   * Output: \verbinclude Map_placement_new.out
   *
-  * This class is the return type of Matrix::Map() but can also be used directly.
+  * This class is the return type of PlainObjectBase::Map() but can also be used directly.
   *
-  * \sa Matrix::Map(), \ref TopicStorageOrders
+  * \sa PlainObjectBase::Map(), \ref TopicStorageOrders
   */
 
 namespace internal {

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;

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;
 

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

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());
     }
 
@@ -425,9 +440,6 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::type
       * while the AlignedMap() functions return aligned Map objects and thus should be called only with 16-byte-aligned
       * \a data pointers.
       *
-      * These methods do not allow to specify strides. If you need to specify strides, you have to
-      * use the Map class directly.
-      *
       * \see class Map
       */
     //@{
@@ -584,6 +596,7 @@ class PlainObjectBase : public internal::dense_xpr_base<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
     }
@@ -641,6 +654,7 @@ struct internal::conservative_resize_like_impl
     if ( ( Derived::IsRowMajor && _this.cols() == cols) || // row-major and we change only the number of rows
          (!Derived::IsRowMajor && _this.rows() == rows) )  // column-major and we change only the number of columns
     {
+      internal::check_rows_cols_for_overflow(rows, cols);
       _this.derived().m_storage.conservativeResize(rows*cols,rows,cols);
     }
     else

Eigen/src/Core/ProductBase.h

@@ -256,16 +256,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; }
     

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,

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

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);
-  }
-};
+// 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);\
+    }\
+};\
 
-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);
-  }
-};
+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
 

Eigen/src/Core/products/GeneralBlockPanelKernel.h

@@ -118,14 +118,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);
     }

Eigen/src/Core/util/Macros.h

@@ -28,7 +28,7 @@
 
 #define EIGEN_WORLD_VERSION 3
 #define EIGEN_MAJOR_VERSION 0
-#define EIGEN_MINOR_VERSION 2
+#define EIGEN_MINOR_VERSION 4
 
 #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 +45,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 +130,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 +252,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

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) ) );
     }
 

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*sizeof(Scalar)) % 16) == 0))
 #else
              0
 #endif

Eigen/src/Eigen2Support/Geometry/AlignedBox.h

@@ -51,19 +51,19 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim==
   { if (AmbientDimAtCompileTime!=Dynamic) setNull(); }
 
   /** Constructs a null box with \a _dim the dimension of the ambient space. */
-  inline explicit AlignedBox(int _dim) : m_(min)(_dim), m_(max)(_dim)
+  inline explicit AlignedBox(int _dim) : m_min(_dim), m_max(_dim)
   { setNull(); }
 
   /** Constructs a box with extremities \a _min and \a _max. */
-  inline AlignedBox(const VectorType& _min, const VectorType& _max) : m_(min)(_min), m_(max)(_max) {}
+  inline AlignedBox(const VectorType& _min, const VectorType& _max) : m_min(_min), m_max(_max) {}
 
   /** Constructs a box containing a single point \a p. */
-  inline explicit AlignedBox(const VectorType& p) : m_(min)(p), m_(max)(p) {}
+  inline explicit AlignedBox(const VectorType& p) : m_min(p), m_max(p) {}
 
   ~AlignedBox() {}
 
   /** \returns the dimension in which the box holds */
-  inline int dim() const { return AmbientDimAtCompileTime==Dynamic ? m_min.size()-1 : AmbientDimAtCompileTime; }
+  inline int dim() const { return AmbientDimAtCompileTime==Dynamic ? m_min.size()-1 : int(AmbientDimAtCompileTime); }
 
   /** \returns true if the box is null, i.e, empty. */
   inline bool isNull() const { return (m_min.cwise() > m_max).any(); }
@@ -71,8 +71,8 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim==
   /** Makes \c *this a null/empty box. */
   inline void setNull()
   {
-    m_min.setConstant( std::numeric_limits<Scalar>::(max)());
-    m_max.setConstant(-std::numeric_limits<Scalar>::(max)());
+    m_min.setConstant( (std::numeric_limits<Scalar>::max)());
+    m_max.setConstant(-(std::numeric_limits<Scalar>::max)());
   }
 
   /** \returns the minimal corner */
@@ -90,19 +90,19 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim==
 
   /** \returns true if the box \a b is entirely inside the box \c *this. */
   inline bool contains(const AlignedBox& b) const
-  { return (m_min.cwise()<=b.(min)()).all() && (b.(max)().cwise()<=m_max).all(); }
+  { return (m_min.cwise()<=(b.min)()).all() && ((b.max)().cwise()<=m_max).all(); }
 
   /** Extends \c *this such that it contains the point \a p and returns a reference to \c *this. */
   inline AlignedBox& extend(const VectorType& p)
-  { m_min = m_min.cwise().(min)(p); m_max = m_max.cwise().(max)(p); return *this; }
+  { m_min = (m_min.cwise().min)(p); m_max = (m_max.cwise().max)(p); return *this; }
 
   /** Extends \c *this such that it contains the box \a b and returns a reference to \c *this. */
   inline AlignedBox& extend(const AlignedBox& b)
-  { m_min = m_min.cwise().(min)(b.m_min); m_max = m_max.cwise().(max)(b.m_max); return *this; }
+  { m_min = (m_min.cwise().min)(b.m_min); m_max = (m_max.cwise().max)(b.m_max); return *this; }
 
   /** Clamps \c *this by the box \a b and returns a reference to \c *this. */
   inline AlignedBox& clamp(const AlignedBox& b)
-  { m_min = m_min.cwise().(max)(b.m_min); m_max = m_max.cwise().(min)(b.m_max); return *this; }
+  { m_min = (m_min.cwise().max)(b.m_min); m_max = (m_max.cwise().min)(b.m_max); return *this; }
 
   /** Translate \c *this by the vector \a t and returns a reference to \c *this. */
   inline AlignedBox& translate(const VectorType& t)
@@ -138,8 +138,8 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim==
   template<typename OtherScalarType>
   inline explicit AlignedBox(const AlignedBox<OtherScalarType,AmbientDimAtCompileTime>& other)
   {
-    m_min = other.(min)().template cast<Scalar>();
-    m_max = other.(max)().template cast<Scalar>();
+    m_min = (other.min)().template cast<Scalar>();
+    m_max = (other.max)().template cast<Scalar>();
   }
 
   /** \returns \c true if \c *this is approximately equal to \a other, within the precision

Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h

@@ -220,6 +220,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     const MatrixType& eigenvectors() const
     {
       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
+      eigen_assert(info() == Success && "Eigenvalue computation did not converge.");
       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
       return m_eivec;
     }
@@ -242,6 +243,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     const RealVectorType& eigenvalues() const
     {
       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
+      eigen_assert(info() == Success && "Eigenvalue computation did not converge.");
       return m_eivalues;
     }
 
@@ -266,6 +268,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     MatrixType operatorSqrt() const
     {
       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
+      eigen_assert(info() == Success && "Eigenvalue computation did not converge.");
       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
       return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
     }
@@ -291,6 +294,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     MatrixType operatorInverseSqrt() const
     {
       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
+      eigen_assert(info() == Success && "Eigenvalue computation did not converge.");
       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
       return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
     }
@@ -307,7 +311,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 +412,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 +423,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 +439,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;

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;
   

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.

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;
 

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 };
 

bench/btl/generic_bench/btl.hh

@@ -39,7 +39,7 @@
 #endif
 
 #if (defined __GNUC__)
-#define BTL_ASM_COMMENT(X)  asm("#"X)
+#define BTL_ASM_COMMENT(X)  asm("#" X)
 #else
 #define BTL_ASM_COMMENT(X)
 #endif

doc/C08_TutorialGeometry.dox

@@ -147,7 +147,7 @@ OpenGL compatibility \b 3D </td><td>\code
 glLoadMatrixf(t.data());\endcode</td></tr>
 <tr class="alt"><td>
 OpenGL compatibility \b 2D </td><td>\code
-Affine3f aux(Affine3f::Identity);
+Affine3f aux(Affine3f::Identity());
 aux.linear().topLeftCorner<2,2>() = t.linear();
 aux.translation().start<2>() = t.translation();
 glLoadMatrixf(aux.data());\endcode</td></tr>

doc/CMakeLists.txt

@@ -64,9 +64,14 @@ add_custom_target(
 
 add_dependencies(doc-eigen-prerequisites all_snippets all_examples)
 add_dependencies(doc-unsupported-prerequisites unsupported_snippets unsupported_examples)
+
 add_custom_target(doc ALL
   COMMAND doxygen Doxyfile-unsupported
   COMMAND doxygen
   COMMAND doxygen Doxyfile-unsupported # run doxygen twice to get proper eigen <=> unsupported cross references
+  COMMAND ${CMAKE_COMMAND} -E rename html eigen-doc
+  COMMAND ${CMAKE_COMMAND} -E tar cvfz eigen-doc/eigen-doc.tgz eigen-doc/*.html eigen-doc/*.map eigen-doc/*.png eigen-doc/*.css eigen-doc/*.js eigen-doc/*.txt eigen-doc/unsupported
+  COMMAND ${CMAKE_COMMAND} -E rename eigen-doc html
   WORKING_DIRECTORY ${Eigen_BINARY_DIR}/doc)
+
 add_dependencies(doc doc-eigen-prerequisites doc-unsupported-prerequisites)

doc/Overview.dox

@@ -8,7 +8,7 @@ o /** \mainpage Eigen
   | \ref QuickRefPage         "Short reference"
 </div>
 
-This is the API documentation for Eigen3.
+This is the API documentation for Eigen3. You can <a href="eigen-doc.tgz">download</a> it as a tgz archive for offline reading.
 
 Eigen2 users: here is a \ref Eigen2ToEigen3 guide to help porting your application.
 

doc/examples/TutorialLinAlgSelfAdjointEigenSolver.cpp

@@ -10,8 +10,9 @@ int main()
    A << 1, 2, 2, 3;
    cout << "Here is the matrix A:\n" << A << endl;
    SelfAdjointEigenSolver<Matrix2f> eigensolver(A);
+   if (eigensolver.info() != Success) abort();
    cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << endl;
-   cout << "Here's a matrix whose columns are eigenvectors of A "
+   cout << "Here's a matrix whose columns are eigenvectors of A \n"
         << "corresponding to these eigenvalues:\n"
         << eigensolver.eigenvectors() << endl;
 }

doc/unsupported_modules.dox

@@ -1,64 +0,0 @@
-
-namespace Eigen {
-
-/** \defgroup Unsupported_modules Unsupported modules
-  *
-  * The unsupported modules are contributions from various users. They are
-  * provided "as is", without any support. Nevertheless, some of them are
-  * subject to be included in Eigen in the future.
-  */
-
-// please list here all unsupported modules
-
-/** \ingroup Unsupported_modules
-  * \defgroup AdolcForward_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup AlignedVector3_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup AutoDiff_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup BVH_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup FFT_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup IterativeSolvers_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup MatrixFunctions_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup MoreVectorization */
-
-/** \ingroup Unsupported_modules
-  * \defgroup NonLinearOptimization_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup NumericalDiff_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup Polynomials_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup Skyline_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup CholmodSupport_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup SuperLUSupport_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup UmfPackSupport_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup SparseExtra_Module */
-
-/** \ingroup Unsupported_modules
-  * \defgroup MpfrcxxSupport_Module */
-
-} // namespace Eigen

test/CMakeLists.txt

@@ -121,7 +121,7 @@ ei_add_test(nullary)
 ei_add_test(nesting_ops "${CMAKE_CXX_FLAGS_DEBUG}")
 ei_add_test(zerosized)
 ei_add_test(dontalign)
-
+ei_add_test(sizeoverflow)
 ei_add_test(prec_inverse_4x4)
 
 string(TOLOWER "${CMAKE_CXX_COMPILER}" cmake_cxx_compiler_tolower)

test/cholesky.cpp

@@ -245,6 +245,22 @@ template<typename MatrixType> void cholesky_cplx(const MatrixType& m)
 
 }
 
+// regression test for bug 241
+template<typename MatrixType> void cholesky_bug241(const MatrixType& m)
+{
+  eigen_assert(m.rows() == 2 && m.cols() == 2);
+
+  typedef typename MatrixType::Scalar Scalar;
+  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
+
+  MatrixType matA;
+  matA << 1, 1, 1, 1;
+  VectorType vecB;
+  vecB << 1, 1;
+  VectorType vecX = matA.ldlt().solve(vecB);
+  VERIFY_IS_APPROX(matA * vecX, vecB);
+}
+
 template<typename MatrixType> void cholesky_verify_assert()
 {
   MatrixType tmp;
@@ -271,6 +287,7 @@ void test_cholesky()
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1( cholesky(Matrix<double,1,1>()) );
     CALL_SUBTEST_3( cholesky(Matrix2d()) );
+    CALL_SUBTEST_3( cholesky_bug241(Matrix2d()) );
     CALL_SUBTEST_4( cholesky(Matrix3f()) );
     CALL_SUBTEST_5( cholesky(Matrix4d()) );
     s = internal::random<int>(1,200);

test/eigen2/main.h

@@ -29,6 +29,10 @@
 #include <string>
 #include <vector>
 
+#ifdef NDEBUG
+#undef NDEBUG
+#endif
+
 #ifndef EIGEN_TEST_FUNC
 #error EIGEN_TEST_FUNC must be defined
 #endif

test/geo_quaternion.cpp

@@ -120,6 +120,10 @@ template<typename Scalar, int Options> void quaternion(void)
   VERIFY_IS_APPROX(q1f.template cast<Scalar>(),q1);
   Quaternion<double> q1d = q1.template cast<double>();
   VERIFY_IS_APPROX(q1d.template cast<Scalar>(),q1);
+
+  // test bug 369 - improper alignment.
+  Quaternionx *q = new Quaternionx;
+  delete q;
 }
 
 template<typename Scalar> void mapQuaternion(void){
@@ -191,7 +195,6 @@ template<typename PlainObjectType> void check_const_correctness(const PlainObjec
   VERIFY( !(Map<ConstPlainObjectType, Aligned>::Flags & LvalueBit) );
 }
 
-
 void test_geo_quaternion()
 {
   for(int i = 0; i < g_repeat; i++) {

test/householder.cpp

@@ -48,7 +48,7 @@ template<typename MatrixType> void householder(const MatrixType& m)
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, Dynamic> VBlockMatrixType;
   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime> TMatrixType;
   
-  Matrix<Scalar, EIGEN_SIZE_MAX(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), 1> _tmp(std::max(rows,cols));
+  Matrix<Scalar, EIGEN_SIZE_MAX(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), 1> _tmp((std::max)(rows,cols));
   Scalar* tmp = &_tmp.coeffRef(0,0);
 
   Scalar beta;

test/jacobisvd.cpp

@@ -66,7 +66,7 @@ void jacobisvd_compare_to_full(const MatrixType& m,
   typedef typename MatrixType::Index Index;
   Index rows = m.rows();
   Index cols = m.cols();
-  Index diagSize = std::min(rows, cols);
+  Index diagSize = (std::min)(rows, cols);
 
   JacobiSVD<MatrixType, QRPreconditioner> svd(m, computationOptions);
 
@@ -244,6 +244,17 @@ void jacobisvd_inf_nan()
   svd.compute(m, ComputeFullU | ComputeFullV);
 }
 
+// Regression test for bug 286: JacobiSVD loops indefinitely with some
+// matrices containing denormal numbers.
+void jacobisvd_bug286()
+{
+  Matrix2d M;
+  M << -7.90884e-313, -4.94e-324,
+                 0, 5.60844e-313;
+  JacobiSVD<Matrix2d> svd;
+  svd.compute(M); // just check we don't loop indefinitely
+}
+
 void jacobisvd_preallocate()
 {
   Vector3f v(3.f, 2.f, 1.f);
@@ -280,8 +291,6 @@ void jacobisvd_preallocate()
   internal::set_is_malloc_allowed(false);
   svd2.compute(m, ComputeFullU|ComputeFullV);
   internal::set_is_malloc_allowed(true);
-
-  
 }
 
 void test_jacobisvd()
@@ -336,4 +345,7 @@ void test_jacobisvd()
 
   // Check that preallocation avoids subsequent mallocs
   CALL_SUBTEST_9( jacobisvd_preallocate() );
+
+  // Regression check for bug 286
+  CALL_SUBTEST_2( jacobisvd_bug286() );
 }

test/lu.cpp

@@ -64,7 +64,7 @@ template<typename MatrixType> void lu_non_invertible()
   typedef Matrix<typename MatrixType::Scalar, RowsAtCompileTime, RowsAtCompileTime>
           RMatrixType;
 
-  Index rank = internal::random<Index>(1, std::min(rows, cols)-1);
+  Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1);
 
   // The image of the zero matrix should consist of a single (zero) column vector
   VERIFY((MatrixType::Zero(rows,cols).fullPivLu().image(MatrixType::Zero(rows,cols)).cols() == 1));
@@ -84,8 +84,8 @@ template<typename MatrixType> void lu_non_invertible()
   MatrixType u(rows,cols);
   u = lu.matrixLU().template triangularView<Upper>();
   RMatrixType l = RMatrixType::Identity(rows,rows);
-  l.block(0,0,rows,std::min(rows,cols)).template triangularView<StrictlyLower>()
-    = lu.matrixLU().block(0,0,rows,std::min(rows,cols));
+  l.block(0,0,rows,(std::min)(rows,cols)).template triangularView<StrictlyLower>()
+    = lu.matrixLU().block(0,0,rows,(std::min)(rows,cols));
 
   VERIFY_IS_APPROX(lu.permutationP() * m1 * lu.permutationQ(), l*u);
 

test/main.h

@@ -23,9 +23,6 @@
 // License and a copy of the GNU General Public License along with
 // Eigen. If not, see <http://www.gnu.org/licenses/>.
 
-#define min(A,B) please_protect_your_min_with_parentheses
-#define max(A,B) please_protect_your_max_with_parentheses
-
 #include <cstdlib>
 #include <cerrno>
 #include <ctime>
@@ -33,6 +30,15 @@
 #include <string>
 #include <vector>
 #include <typeinfo>
+#include <limits>
+#include <algorithm>
+#include <sstream>
+#include <complex>
+#include <deque>
+#include <queue>
+
+#define min(A,B) please_protect_your_min_with_parentheses
+#define max(A,B) please_protect_your_max_with_parentheses
 
 #ifdef NDEBUG
 #undef NDEBUG
@@ -112,7 +118,7 @@ namespace Eigen
       }                                     \
       else if (Eigen::internal::push_assert)       \
       {                                     \
-        eigen_assert_list.push_back(std::string(EI_PP_MAKE_STRING(__FILE__)" ("EI_PP_MAKE_STRING(__LINE__)") : "#a) ); \
+        eigen_assert_list.push_back(std::string(EI_PP_MAKE_STRING(__FILE__) " (" EI_PP_MAKE_STRING(__LINE__) ") : " #a) ); \
       }
 
     #define VERIFY_RAISES_ASSERT(a)                                                   \

test/nullary.cpp

@@ -38,7 +38,7 @@ bool equalsIdentity(const MatrixType& A)
     }
   }
   for (Index i = 0; i < A.rows(); ++i) {
-    for (Index j = 0; j < std::min(i, A.cols()); ++j) {
+    for (Index j = 0; j < (std::min)(i, A.cols()); ++j) {
       offDiagOK = offDiagOK && (A(i,j) == zero);
     }
   }

test/packetmath.cpp

@@ -128,7 +128,7 @@ template<typename Scalar> void packetmath()
   {
     data1[i] = internal::random<Scalar>()/RealScalar(PacketSize);
     data2[i] = internal::random<Scalar>()/RealScalar(PacketSize);
-    refvalue = std::max(refvalue,internal::abs(data1[i]));
+    refvalue = (std::max)(refvalue,internal::abs(data1[i]));
   }
 
   internal::pstore(data2, internal::pload<Packet>(data1));
@@ -264,16 +264,16 @@ template<typename Scalar> void packetmath_real()
 
   ref[0] = data1[0];
   for (int i=0; i<PacketSize; ++i)
-    ref[0] = std::min(ref[0],data1[i]);
+    ref[0] = (std::min)(ref[0],data1[i]);
   VERIFY(internal::isApprox(ref[0], internal::predux_min(internal::pload<Packet>(data1))) && "internal::predux_min");
 
-  CHECK_CWISE2(std::min, internal::pmin);
-  CHECK_CWISE2(std::max, internal::pmax);
+  CHECK_CWISE2((std::min), internal::pmin);
+  CHECK_CWISE2((std::max), internal::pmax);
   CHECK_CWISE1(internal::abs, internal::pabs);
 
   ref[0] = data1[0];
   for (int i=0; i<PacketSize; ++i)
-    ref[0] = std::max(ref[0],data1[i]);
+    ref[0] = (std::max)(ref[0],data1[i]);
   VERIFY(internal::isApprox(ref[0], internal::predux_max(internal::pload<Packet>(data1))) && "internal::predux_max");
   
   for (int i=0; i<PacketSize; ++i)

test/prec_inverse_4x4.cpp

@@ -58,7 +58,7 @@ template<typename MatrixType> void inverse_general_4x4(int repeat)
     MatrixType inv = m.inverse();
     double error = double( (m*inv-MatrixType::Identity()).norm() * absdet / NumTraits<Scalar>::epsilon() );
     error_sum += error;
-    error_max = std::max(error_max, error);
+    error_max = (std::max)(error_max, error);
   }
   std::cerr << "inverse_general_4x4, Scalar = " << type_name<Scalar>() << std::endl;
   double error_avg = error_sum / repeat;

test/product.h

@@ -29,7 +29,7 @@ template<typename Derived1, typename Derived2>
 bool areNotApprox(const MatrixBase<Derived1>& m1, const MatrixBase<Derived2>& m2, typename Derived1::RealScalar epsilon = NumTraits<typename Derived1::RealScalar>::dummy_precision())
 {
   return !((m1-m2).cwiseAbs2().maxCoeff() < epsilon * epsilon
-                          * std::max(m1.cwiseAbs2().maxCoeff(), m2.cwiseAbs2().maxCoeff()));
+                          * (std::max)(m1.cwiseAbs2().maxCoeff(), m2.cwiseAbs2().maxCoeff()));
 }
 
 template<typename MatrixType> void product(const MatrixType& m)
@@ -102,7 +102,7 @@ template<typename MatrixType> void product(const MatrixType& m)
 
   // test the previous tests were not screwed up because operator* returns 0
   // (we use the more accurate default epsilon)
-  if (!NumTraits<Scalar>::IsInteger && std::min(rows,cols)>1)
+  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows,cols)>1)
   {
     VERIFY(areNotApprox(m1.transpose()*m2,m2.transpose()*m1));
   }
@@ -111,7 +111,7 @@ template<typename MatrixType> void product(const MatrixType& m)
   res = square;
   res.noalias() += m1 * m2.transpose();
   VERIFY_IS_APPROX(res, square + m1 * m2.transpose());
-  if (!NumTraits<Scalar>::IsInteger && std::min(rows,cols)>1)
+  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows,cols)>1)
   {
     VERIFY(areNotApprox(res,square + m2 * m1.transpose()));
   }
@@ -123,7 +123,7 @@ template<typename MatrixType> void product(const MatrixType& m)
   res = square;
   res.noalias() -= m1 * m2.transpose();
   VERIFY_IS_APPROX(res, square - (m1 * m2.transpose()));
-  if (!NumTraits<Scalar>::IsInteger && std::min(rows,cols)>1)
+  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows,cols)>1)
   {
     VERIFY(areNotApprox(res,square - m2 * m1.transpose()));
   }
@@ -147,7 +147,7 @@ template<typename MatrixType> void product(const MatrixType& m)
   res2 = square2;
   res2.noalias() += m1.transpose() * m2;
   VERIFY_IS_APPROX(res2, square2 + m1.transpose() * m2);
-  if (!NumTraits<Scalar>::IsInteger && std::min(rows,cols)>1)
+  if (!NumTraits<Scalar>::IsInteger && (std::min)(rows,cols)>1)
   {
     VERIFY(areNotApprox(res2,square2 + m2.transpose() * m1));
   }

test/product_extra.cpp

@@ -116,6 +116,16 @@ template<typename MatrixType> void product_extra(const MatrixType& m)
   VERIFY_IS_APPROX(tmp, m1 * m1.adjoint() * s1);
 }
 
+// Regression test for bug reported at http://forum.kde.org/viewtopic.php?f=74&t=96947
+void mat_mat_scalar_scalar_product()
+{
+  Eigen::Matrix2Xd dNdxy(2, 3);
+  dNdxy << -0.5, 0.5, 0,
+           -0.3, 0, 0.3;
+  double det = 6.0, wt = 0.5;
+  VERIFY_IS_APPROX(dNdxy.transpose()*dNdxy*det*wt, det*wt*dNdxy.transpose()*dNdxy);
+}
+  
 void zero_sized_objects()
 {
   // Bug 127
@@ -145,6 +155,7 @@ void test_product_extra()
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1( product_extra(MatrixXf(internal::random<int>(1,320), internal::random<int>(1,320))) );
     CALL_SUBTEST_2( product_extra(MatrixXd(internal::random<int>(1,320), internal::random<int>(1,320))) );
+    CALL_SUBTEST_2( mat_mat_scalar_scalar_product() );
     CALL_SUBTEST_3( product_extra(MatrixXcf(internal::random<int>(1,150), internal::random<int>(1,150))) );
     CALL_SUBTEST_4( product_extra(MatrixXcd(internal::random<int>(1,150), internal::random<int>(1,150))) );
     CALL_SUBTEST_5( zero_sized_objects() );

test/qr_colpivoting.cpp

@@ -31,7 +31,7 @@ template<typename MatrixType> void qr()
   typedef typename MatrixType::Index Index;
 
   Index rows = internal::random<Index>(2,200), cols = internal::random<Index>(2,200), cols2 = internal::random<Index>(2,200);
-  Index rank = internal::random<Index>(1, std::min(rows, cols)-1);
+  Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1);
 
   typedef typename MatrixType::Scalar Scalar;
   typedef typename MatrixType::RealScalar RealScalar;
@@ -64,7 +64,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
 {
   enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
   typedef typename MatrixType::Scalar Scalar;
-  int rank = internal::random<int>(1, std::min(int(Rows), int(Cols))-1);
+  int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols))-1);
   Matrix<Scalar,Rows,Cols> m1;
   createRandomPIMatrixOfRank(rank,Rows,Cols,m1);
   ColPivHouseholderQR<Matrix<Scalar,Rows,Cols> > qr(m1);

test/qr_fullpivoting.cpp

@@ -31,7 +31,7 @@ template<typename MatrixType> void qr()
   typedef typename MatrixType::Index Index;
 
   Index rows = internal::random<Index>(20,200), cols = internal::random<int>(20,200), cols2 = internal::random<int>(20,200);
-  Index rank = internal::random<Index>(1, std::min(rows, cols)-1);
+  Index rank = internal::random<Index>(1, (std::min)(rows, cols)-1);
 
   typedef typename MatrixType::Scalar Scalar;
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixQType;

test/redux.cpp

@@ -43,8 +43,8 @@ template<typename MatrixType> void matrixRedux(const MatrixType& m)
   {
     s += m1(i,j);
     p *= m1(i,j);
-    minc = std::min(internal::real(minc), internal::real(m1(i,j)));
-    maxc = std::max(internal::real(maxc), internal::real(m1(i,j)));
+    minc = (std::min)(internal::real(minc), internal::real(m1(i,j)));
+    maxc = (std::max)(internal::real(maxc), internal::real(m1(i,j)));
   }
   const Scalar mean = s/Scalar(RealScalar(rows*cols));
 
@@ -86,8 +86,8 @@ template<typename VectorType> void vectorRedux(const VectorType& w)
     {
       s += v[j];
       p *= v[j];
-      minc = std::min(minc, internal::real(v[j]));
-      maxc = std::max(maxc, internal::real(v[j]));
+      minc = (std::min)(minc, internal::real(v[j]));
+      maxc = (std::max)(maxc, internal::real(v[j]));
     }
     VERIFY_IS_MUCH_SMALLER_THAN(internal::abs(s - v.head(i).sum()), Scalar(1));
     VERIFY_IS_APPROX(p, v.head(i).prod());
@@ -103,8 +103,8 @@ template<typename VectorType> void vectorRedux(const VectorType& w)
     {
       s += v[j];
       p *= v[j];
-      minc = std::min(minc, internal::real(v[j]));
-      maxc = std::max(maxc, internal::real(v[j]));
+      minc = (std::min)(minc, internal::real(v[j]));
+      maxc = (std::max)(maxc, internal::real(v[j]));
     }
     VERIFY_IS_MUCH_SMALLER_THAN(internal::abs(s - v.tail(size-i).sum()), Scalar(1));
     VERIFY_IS_APPROX(p, v.tail(size-i).prod());
@@ -120,8 +120,8 @@ template<typename VectorType> void vectorRedux(const VectorType& w)
     {
       s += v[j];
       p *= v[j];
-      minc = std::min(minc, internal::real(v[j]));
-      maxc = std::max(maxc, internal::real(v[j]));
+      minc = (std::min)(minc, internal::real(v[j]));
+      maxc = (std::max)(maxc, internal::real(v[j]));
     }
     VERIFY_IS_MUCH_SMALLER_THAN(internal::abs(s - v.segment(i, size-2*i).sum()), Scalar(1));
     VERIFY_IS_APPROX(p, v.segment(i, size-2*i).prod());

test/sizeoverflow.cpp

@@ -0,0 +1,81 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2011 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+
+#define VERIFY_THROWS_BADALLOC(a) {                           \
+    bool threw = false;                                       \
+    try {                                                     \
+      a;                                                      \
+    }                                                         \
+    catch (std::bad_alloc&) { threw = true; }                 \
+    VERIFY(threw && "should have thrown bad_alloc: " #a);     \
+  }
+
+typedef DenseIndex Index;
+
+template<typename MatrixType>
+void triggerMatrixBadAlloc(Index rows, Index cols)
+{
+  VERIFY_THROWS_BADALLOC( MatrixType m(rows, cols) );
+  VERIFY_THROWS_BADALLOC( MatrixType m; m.resize(rows, cols) );
+  VERIFY_THROWS_BADALLOC( MatrixType m; m.conservativeResize(rows, cols) );
+}
+
+template<typename VectorType>
+void triggerVectorBadAlloc(Index size)
+{
+  VERIFY_THROWS_BADALLOC( VectorType v(size) );
+  VERIFY_THROWS_BADALLOC( VectorType v; v.resize(size) );
+  VERIFY_THROWS_BADALLOC( VectorType v; v.conservativeResize(size) );
+}
+
+void test_sizeoverflow()
+{
+  // there are 2 levels of overflow checking. first in PlainObjectBase.h we check for overflow in rows*cols computations.
+  // this is tested in tests of the form times_itself_gives_0 * times_itself_gives_0
+  // Then in Memory.h we check for overflow in size * sizeof(T) computations.
+  // this is tested in tests of the form times_4_gives_0 * sizeof(float)
+  
+  size_t times_itself_gives_0 = size_t(1) << (8 * sizeof(Index) / 2);
+  VERIFY(times_itself_gives_0 * times_itself_gives_0 == 0);
+
+  size_t times_4_gives_0 = size_t(1) << (8 * sizeof(Index) - 2);
+  VERIFY(times_4_gives_0 * 4 == 0);
+
+  size_t times_8_gives_0 = size_t(1) << (8 * sizeof(Index) - 3);
+  VERIFY(times_8_gives_0 * 8 == 0);
+
+  triggerMatrixBadAlloc<MatrixXf>(times_itself_gives_0, times_itself_gives_0);
+  triggerMatrixBadAlloc<MatrixXf>(times_itself_gives_0 / 4, times_itself_gives_0);
+  triggerMatrixBadAlloc<MatrixXf>(times_4_gives_0, 1);
+
+  triggerMatrixBadAlloc<MatrixXd>(times_itself_gives_0, times_itself_gives_0);
+  triggerMatrixBadAlloc<MatrixXd>(times_itself_gives_0 / 8, times_itself_gives_0);
+  triggerMatrixBadAlloc<MatrixXd>(times_8_gives_0, 1);
+  
+  triggerVectorBadAlloc<VectorXf>(times_4_gives_0);
+  
+  triggerVectorBadAlloc<VectorXd>(times_8_gives_0);
+}

test/sparse.h

@@ -29,6 +29,15 @@
 #include "main.h"
 
 #if EIGEN_GNUC_AT_LEAST(4,0) && !defined __ICC && !defined(__clang__)
+
+#ifdef min
+#undef min
+#endif
+
+#ifdef max
+#undef max
+#endif
+
 #include <tr1/unordered_map>
 #define EIGEN_UNORDERED_MAP_SUPPORT
 namespace std {

test/sparse_basic.cpp

@@ -33,7 +33,7 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
   typedef typename SparseMatrixType::Scalar Scalar;
   enum { Flags = SparseMatrixType::Flags };
 
-  double density = std::max(8./(rows*cols), 0.01);
+  double density = (std::max)(8./(rows*cols), 0.01);
   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
   typedef Matrix<Scalar,Dynamic,1> DenseVector;
   Scalar eps = 1e-6;
@@ -206,7 +206,7 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
     initSparse<Scalar>(density, refMat2, m2);
     int j0 = internal::random(0,rows-2);
     int j1 = internal::random(0,rows-2);
-    int n0 = internal::random<int>(1,rows-std::max(j0,j1));
+    int n0 = internal::random<int>(1,rows-(std::max)(j0,j1));
     VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0));
     VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0),
                      refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0));

test/sparse_product.cpp

@@ -32,7 +32,7 @@ template<typename SparseMatrixType> void sparse_product(const SparseMatrixType&
   typedef typename SparseMatrixType::Scalar Scalar;
   enum { Flags = SparseMatrixType::Flags };
 
-  double density = std::max(8./(rows*cols), 0.01);
+  double density = (std::max)(8./(rows*cols), 0.01);
   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
   typedef Matrix<Scalar,Dynamic,1> DenseVector;
 

test/sparse_solvers.cpp

@@ -47,7 +47,7 @@ initSPD(double density,
 
 template<typename Scalar> void sparse_solvers(int rows, int cols)
 {
-  double density = std::max(8./(rows*cols), 0.01);
+  double density = (std::max)(8./(rows*cols), 0.01);
   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
   typedef Matrix<Scalar,Dynamic,1> DenseVector;
   // Scalar eps = 1e-6;

test/sparse_vector.cpp

@@ -26,8 +26,8 @@
 
 template<typename Scalar> void sparse_vector(int rows, int cols)
 {
-  double densityMat = std::max(8./(rows*cols), 0.01);
-  double densityVec = std::max(8./float(rows), 0.1);
+  double densityMat = (std::max)(8./(rows*cols), 0.01);
+  double densityVec = (std::max)(8./float(rows), 0.1);
   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
   typedef Matrix<Scalar,Dynamic,1> DenseVector;
   typedef SparseVector<Scalar> SparseVectorType;

test/stable_norm.cpp

@@ -68,8 +68,8 @@ template<typename MatrixType> void stable_norm(const MatrixType& m)
   Index rows = m.rows();
   Index cols = m.cols();
 
-  Scalar big = internal::random<Scalar>() * (std::numeric_limits<RealScalar>::max() * RealScalar(1e-4));
-  Scalar small = internal::random<Scalar>() * (std::numeric_limits<RealScalar>::min() * RealScalar(1e4));
+  Scalar big = internal::random<Scalar>() * ((std::numeric_limits<RealScalar>::max)() * RealScalar(1e-4));
+  Scalar small = internal::random<Scalar>() * ((std::numeric_limits<RealScalar>::min)() * RealScalar(1e4));
 
   MatrixType  vzero = MatrixType::Zero(rows, cols),
               vrand = MatrixType::Random(rows, cols),

unsupported/Eigen/FFT

@@ -331,7 +331,7 @@ class FFT
         // if the vector is strided, then we need to copy it to a packed temporary
         Matrix<src_type,1,Dynamic> tmp;
         if ( resize_input ) {
-          size_t ncopy = std::min(src.size(),src.size() + resize_input);
+          size_t ncopy = (std::min)(src.size(),src.size() + resize_input);
           tmp.setZero(src.size() + resize_input);
           if ( realfft && HasFlag(HalfSpectrum) ) {
             // pad at the Nyquist bin

unsupported/Eigen/MPRealSupport

@@ -35,7 +35,8 @@
 
 namespace Eigen {
 
-  /** \defgroup MPRealSupport_Module MPFRC++ Support module
+  /** \ingroup Unsupported_modules
+    * \defgroup MPRealSupport_Module MPFRC++ Support module
     *
     * \code
     * #include <Eigen/MPRealSupport>

unsupported/Eigen/src/BVH/BVAlgorithms.h

@@ -231,7 +231,7 @@ private:
 template<typename BVH, typename Minimizer>
 typename Minimizer::Scalar BVMinimize(const BVH &tree, Minimizer &minimizer)
 {
-  return internal::minimize_helper(tree, minimizer, tree.getRootIndex(), std::numeric_limits<typename Minimizer::Scalar>::max());
+  return internal::minimize_helper(tree, minimizer, tree.getRootIndex(), (std::numeric_limits<typename Minimizer::Scalar>::max)());
 }
 
 /**  Given two BVH's, runs the query on their cartesian product encapsulated by \a minimizer.
@@ -264,7 +264,7 @@ typename Minimizer::Scalar BVMinimize(const BVH1 &tree1, const BVH2 &tree2, Mini
   ObjIter2 oBegin2 = ObjIter2(), oEnd2 = ObjIter2(), oCur2 = ObjIter2();
   std::priority_queue<QueueElement, std::vector<QueueElement>, std::greater<QueueElement> > todo; //smallest is at the top
 
-  Scalar minimum = std::numeric_limits<Scalar>::max();
+  Scalar minimum = (std::numeric_limits<Scalar>::max)();
   todo.push(std::make_pair(Scalar(), std::make_pair(tree1.getRootIndex(), tree2.getRootIndex())));
 
   while(!todo.empty()) {

unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h

@@ -259,7 +259,7 @@ void MatrixExponential<MatrixType>::computeUV(float)
     pade5(m_M);
   } else {
     const float maxnorm = 3.925724783138660f;
-    m_squarings = max(0, (int)ceil(log2(m_l1norm / maxnorm)));
+    m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm)));
     MatrixType A = m_M / pow(Scalar(2), Scalar(static_cast<RealScalar>(m_squarings)));
     pade7(A);
   }
@@ -281,7 +281,7 @@ void MatrixExponential<MatrixType>::computeUV(double)
     pade9(m_M);
   } else {
     const double maxnorm = 5.371920351148152;
-    m_squarings = max(0, (int)ceil(log2(m_l1norm / maxnorm)));
+    m_squarings = (max)(0, (int)ceil(log2(m_l1norm / maxnorm)));
     MatrixType A = m_M / pow(Scalar(2), Scalar(m_squarings));
     pade13(A);
   }

unsupported/doc/Overview.dox

@@ -1,13 +1,22 @@
 namespace Eigen {
 
-o /** \mainpage Eigen's unsupported modules
+/** \mainpage Eigen's unsupported modules
 
 This is the API documentation for Eigen's unsupported modules.
 
 These modules are contributions from various users. They are provided "as is", without any support.
 
+Click on the \e Modules tab at the top of this page to get a list of all unsupported modules.
+
 Don't miss the <a href="..//index.html">official Eigen documentation</a>.
 
+
+\defgroup Unsupported_modules Unsupported modules
+
+The unsupported modules are contributions from various users. They are
+provided "as is", without any support. Nevertheless, some of them are
+subject to be included in Eigen in the future.
+
 */
 
 }

unsupported/test/BVH.cpp

@@ -90,13 +90,13 @@ struct BallPointStuff //this class provides functions to be both an intersector
   }
 
   double minimumOnVolume(const BoxType &r) { ++calls; return r.squaredExteriorDistance(p); }
-  double minimumOnObject(const BallType &b) { ++calls; return std::max(0., (b.center - p).squaredNorm() - SQR(b.radius)); }
+  double minimumOnObject(const BallType &b) { ++calls; return (std::max)(0., (b.center - p).squaredNorm() - SQR(b.radius)); }
   double minimumOnVolumeVolume(const BoxType &r1, const BoxType &r2) { ++calls; return r1.squaredExteriorDistance(r2); }
-  double minimumOnVolumeObject(const BoxType &r, const BallType &b) { ++calls; return SQR(std::max(0., r.exteriorDistance(b.center) - b.radius)); }
-  double minimumOnObjectVolume(const BallType &b, const BoxType &r) { ++calls; return SQR(std::max(0., r.exteriorDistance(b.center) - b.radius)); }
-  double minimumOnObjectObject(const BallType &b1, const BallType &b2){ ++calls; return SQR(std::max(0., (b1.center - b2.center).norm() - b1.radius - b2.radius)); }
+  double minimumOnVolumeObject(const BoxType &r, const BallType &b) { ++calls; return SQR((std::max)(0., r.exteriorDistance(b.center) - b.radius)); }
+  double minimumOnObjectVolume(const BallType &b, const BoxType &r) { ++calls; return SQR((std::max)(0., r.exteriorDistance(b.center) - b.radius)); }
+  double minimumOnObjectObject(const BallType &b1, const BallType &b2){ ++calls; return SQR((std::max)(0., (b1.center - b2.center).norm() - b1.radius - b2.radius)); }
   double minimumOnVolumeObject(const BoxType &r, const VectorType &v) { ++calls; return r.squaredExteriorDistance(v); }
-  double minimumOnObjectObject(const BallType &b, const VectorType &v){ ++calls; return SQR(std::max(0., (b.center - v).norm() - b.radius)); }
+  double minimumOnObjectObject(const BallType &b, const VectorType &v){ ++calls; return SQR((std::max)(0., (b.center - v).norm() - b.radius)); }
 
   VectorType p;
   int calls;

unsupported/test/matrix_exponential.cpp

@@ -36,7 +36,7 @@ double binom(int n, int k)
 template <typename Derived, typename OtherDerived>
 double relerr(const MatrixBase<Derived>& A, const MatrixBase<OtherDerived>& B)
 {
-  return std::sqrt((A - B).cwiseAbs2().sum() / std::min(A.cwiseAbs2().sum(), B.cwiseAbs2().sum()));
+  return std::sqrt((A - B).cwiseAbs2().sum() / (std::min)(A.cwiseAbs2().sum(), B.cwiseAbs2().sum()));
 }
 
 template <typename T>

unsupported/test/sparse_extra.cpp

@@ -67,7 +67,7 @@ template<typename SparseMatrixType> void sparse_extra(const SparseMatrixType& re
   typedef typename SparseMatrixType::Scalar Scalar;
   enum { Flags = SparseMatrixType::Flags };
 
-  double density = std::max(8./(rows*cols), 0.01);
+  double density = (std::max)(8./(rows*cols), 0.01);
   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
   typedef Matrix<Scalar,Dynamic,1> DenseVector;
   Scalar eps = 1e-6;

unsupported/test/sparse_ldlt.cpp

@@ -33,7 +33,7 @@ template<typename Scalar> void sparse_ldlt(int rows, int cols)
 {
   static bool odd = true;
   odd = !odd;
-  double density = std::max(8./(rows*cols), 0.01);
+  double density = (std::max)(8./(rows*cols), 0.01);
   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
   typedef Matrix<Scalar,Dynamic,1> DenseVector;
 

unsupported/test/sparse_llt.cpp

@@ -31,7 +31,7 @@
 
 template<typename Scalar> void sparse_llt(int rows, int cols)
 {
-  double density = std::max(8./(rows*cols), 0.01);
+  double density = (std::max)(8./(rows*cols), 0.01);
   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
   typedef Matrix<Scalar,Dynamic,1> DenseVector;
 

unsupported/test/sparse_lu.cpp

@@ -35,7 +35,7 @@
 
 template<typename Scalar> void sparse_lu(int rows, int cols)
 {
-  double density = std::max(8./(rows*cols), 0.01);
+  double density = (std::max)(8./(rows*cols), 0.01);
   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
   typedef Matrix<Scalar,Dynamic,1> DenseVector;