fix linking issue with manage_caching_sizes_second_if_negative

(transplanted from 8dd3ae282d
)

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

@@ -51,16 +51,16 @@
       #define EIGEN_SSE2_ON_MSVC_2008_OR_LATER
     #endif
   #endif
-#endif
-
-// Remember that usage of defined() in a #define is undefined by the standard
-#if (defined __SSE2__) && ( (!defined __GNUC__) || EIGEN_GNUC_AT_LEAST(4,2) )
-  #define EIGEN_SSE2_BUT_NOT_OLD_GCC
+#else
+  // Remember that usage of defined() in a #define is undefined by the standard
+  #if (defined __SSE2__) && ( (!defined __GNUC__) || EIGEN_GNUC_AT_LEAST(4,2) )
+    #define EIGEN_SSE2_ON_NON_MSVC_BUT_NOT_OLD_GCC
+  #endif
 #endif
 
 #ifndef EIGEN_DONT_VECTORIZE
 
-  #if defined (EIGEN_SSE2_BUT_NOT_OLD_GCC) || defined(EIGEN_SSE2_ON_MSVC_2008_OR_LATER)
+  #if defined (EIGEN_SSE2_ON_NON_MSVC_BUT_NOT_OLD_GCC) || defined(EIGEN_SSE2_ON_MSVC_2008_OR_LATER)
 
     // Defines symbols for compile-time detection of which instructions are
     // used.
@@ -143,6 +143,7 @@
 #ifdef EIGEN_HAS_ERRNO
 #include <cerrno>
 #endif
+#include <cstddef>
 #include <cstdlib>
 #include <cmath>
 #include <complex>
@@ -166,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
 
@@ -174,16 +175,10 @@
   #include <new>
 #endif
 
-// this needs to be done after all possible windows C header includes and before any Eigen source includes
-// (system C++ includes are supposed to be able to deal with this already):
-// windows.h defines min and max macros which would make Eigen fail to compile.
-#if defined(min) || defined(max)
-#error The preprocessor symbols 'min' or 'max' are defined. If you are compiling on Windows, do #define NOMINMAX to prevent windows.h from defining these symbols.
-#endif
-
 // defined in bits/termios.h
 #undef B0
 
+/** \brief Namespace containing all symbols from the %Eigen library. */
 namespace Eigen {
 
 inline static const char *SimdInstructionSetsInUse(void) {
@@ -239,6 +234,8 @@ inline static const char *SimdInstructionSetsInUse(void) {
 // we use size_t frequently and we'll never remember to prepend it with std:: everytime just to
 // ensure QNX/QCC support
 using std::size_t;
+// gcc 4.6.0 wants std:: for ptrdiff_t 
+using std::ptrdiff_t;
 
 /** \defgroup Core_Module Core module
   * This is the main module of Eigen providing dense matrix and vector support

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>

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/Cholesky/LLT.h

@@ -233,7 +233,7 @@ template<> struct llt_inplace<Lower>
 
     Index blockSize = size/8;
     blockSize = (blockSize/16)*16;
-    blockSize = std::min(std::max(blockSize,Index(8)), Index(128));
+    blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128));
 
     for (Index k=0; k<size; k+=blockSize)
     {
@@ -241,7 +241,7 @@ template<> struct llt_inplace<Lower>
       //       A00 |  -  |  -
       // lu  = A10 | A11 |  -
       //       A20 | A21 | A22
-      Index bs = std::min(blockSize, size-k);
+      Index bs = (std::min)(blockSize, size-k);
       Index rs = size - k - bs;
       Block<MatrixType,Dynamic,Dynamic> A11(m,k,   k,   bs,bs);
       Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k,   rs,bs);

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/ArrayWrapper.h

@@ -53,6 +53,12 @@ class ArrayWrapper : public ArrayBase<ArrayWrapper<ExpressionType> >
     EIGEN_DENSE_PUBLIC_INTERFACE(ArrayWrapper)
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(ArrayWrapper)
 
+    typedef typename internal::conditional<
+                       internal::is_lvalue<ExpressionType>::value,
+                       Scalar,
+                       const Scalar
+                     >::type ScalarWithConstIfNotLvalue;
+
     typedef typename internal::nested<ExpressionType>::type NestedExpressionType;
 
     inline ArrayWrapper(const ExpressionType& matrix) : m_expression(matrix) {}
@@ -62,6 +68,9 @@ class ArrayWrapper : public ArrayBase<ArrayWrapper<ExpressionType> >
     inline Index outerStride() const { return m_expression.outerStride(); }
     inline Index innerStride() const { return m_expression.innerStride(); }
 
+    inline ScalarWithConstIfNotLvalue* data() { return m_expression.data(); }
+    inline const Scalar* data() const { return m_expression.data(); }
+
     inline const CoeffReturnType coeff(Index row, Index col) const
     {
       return m_expression.coeff(row, col);
@@ -151,6 +160,12 @@ class MatrixWrapper : public MatrixBase<MatrixWrapper<ExpressionType> >
     EIGEN_DENSE_PUBLIC_INTERFACE(MatrixWrapper)
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(MatrixWrapper)
 
+    typedef typename internal::conditional<
+                       internal::is_lvalue<ExpressionType>::value,
+                       Scalar,
+                       const Scalar
+                     >::type ScalarWithConstIfNotLvalue;
+
     typedef typename internal::nested<ExpressionType>::type NestedExpressionType;
 
     inline MatrixWrapper(const ExpressionType& matrix) : m_expression(matrix) {}
@@ -160,6 +175,9 @@ class MatrixWrapper : public MatrixBase<MatrixWrapper<ExpressionType> >
     inline Index outerStride() const { return m_expression.outerStride(); }
     inline Index innerStride() const { return m_expression.innerStride(); }
 
+    inline ScalarWithConstIfNotLvalue* data() { return m_expression.data(); }
+    inline const Scalar* data() const { return m_expression.data(); }
+
     inline const CoeffReturnType coeff(Index row, Index col) const
     {
       return m_expression.coeff(row, col);

Eigen/src/Core/BandMatrix.h

@@ -87,7 +87,7 @@ class BandMatrixBase : public EigenBase<Derived>
       if (i<=supers())
       {
         start = supers()-i;
-        len = std::min(rows(),std::max<Index>(0,coeffs().rows() - (supers()-i)));
+        len = (std::min)(rows(),std::max<Index>(0,coeffs().rows() - (supers()-i)));
       }
       else if (i>=rows()-subs())
         len = std::max<Index>(0,coeffs().rows() - (i + 1 - rows() + subs()));
@@ -96,11 +96,11 @@ class BandMatrixBase : public EigenBase<Derived>
 
     /** \returns a vector expression of the main diagonal */
     inline Block<CoefficientsType,1,SizeAtCompileTime> diagonal()
-    { return Block<CoefficientsType,1,SizeAtCompileTime>(coeffs(),supers(),0,1,std::min(rows(),cols())); }
+    { return Block<CoefficientsType,1,SizeAtCompileTime>(coeffs(),supers(),0,1,(std::min)(rows(),cols())); }
 
     /** \returns a vector expression of the main diagonal (const version) */
     inline const Block<const CoefficientsType,1,SizeAtCompileTime> diagonal() const
-    { return Block<const CoefficientsType,1,SizeAtCompileTime>(coeffs(),supers(),0,1,std::min(rows(),cols())); }
+    { return Block<const CoefficientsType,1,SizeAtCompileTime>(coeffs(),supers(),0,1,(std::min)(rows(),cols())); }
 
     template<int Index> struct DiagonalIntReturnType {
       enum {
@@ -122,13 +122,13 @@ class BandMatrixBase : public EigenBase<Derived>
     /** \returns a vector expression of the \a N -th sub or super diagonal */
     template<int N> inline typename DiagonalIntReturnType<N>::Type diagonal()
     {
-      return typename DiagonalIntReturnType<N>::BuildType(coeffs(), supers()-N, std::max(0,N), 1, diagonalLength(N));
+      return typename DiagonalIntReturnType<N>::BuildType(coeffs(), supers()-N, (std::max)(0,N), 1, diagonalLength(N));
     }
 
     /** \returns a vector expression of the \a N -th sub or super diagonal */
     template<int N> inline const typename DiagonalIntReturnType<N>::Type diagonal() const
     {
-      return typename DiagonalIntReturnType<N>::BuildType(coeffs(), supers()-N, std::max(0,N), 1, diagonalLength(N));
+      return typename DiagonalIntReturnType<N>::BuildType(coeffs(), supers()-N, (std::max)(0,N), 1, diagonalLength(N));
     }
 
     /** \returns a vector expression of the \a i -th sub or super diagonal */
@@ -166,7 +166,7 @@ class BandMatrixBase : public EigenBase<Derived>
   protected:
 
     inline Index diagonalLength(Index i) const
-    { return i<0 ? std::min(cols(),rows()+i) : std::min(rows(),cols()-i); }
+    { return i<0 ? (std::min)(cols(),rows()+i) : (std::min)(rows(),cols()-i); }
 };
 
 /**
@@ -180,7 +180,7 @@ class BandMatrixBase : public EigenBase<Derived>
   * \param Cols Number of columns, or \b Dynamic
   * \param Supers Number of super diagonal
   * \param Subs Number of sub diagonal
-  * \param _Options A combination of either \b RowMajor or \b ColMajor, and of \b SelfAdjoint
+  * \param _Options A combination of either \b #RowMajor or \b #ColMajor, and of \b #SelfAdjoint
   *                 The former controls \ref TopicStorageOrders "storage order", and defaults to
   *                 column-major. The latter controls whether the matrix represents a selfadjoint 
   *                 matrix in which case either Supers of Subs have to be null.
@@ -284,6 +284,7 @@ class BandMatrixWrapper : public BandMatrixBase<BandMatrixWrapper<_CoefficientsT
       : m_coeffs(coeffs),
         m_rows(rows), m_supers(supers), m_subs(subs)
     {
+      EIGEN_UNUSED_VARIABLE(cols);
       //internal::assert(coeffs.cols()==cols() && (supers()+subs()+1)==coeffs.rows());
     }
 

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,

Eigen/src/Core/CwiseNullaryOp.h

@@ -742,7 +742,7 @@ struct setIdentity_impl<Derived, true>
   static EIGEN_STRONG_INLINE Derived& run(Derived& m)
   {
     m.setZero();
-    const Index size = std::min(m.rows(), m.cols());
+    const Index size = (std::min)(m.rows(), m.cols());
     for(Index i = 0; i < size; ++i) m.coeffRef(i,i) = typename Derived::Scalar(1);
     return m;
   }

Eigen/src/Core/DenseCoeffsBase.h

@@ -35,7 +35,7 @@ template<typename T> struct add_const_on_value_type_if_arithmetic
 /** \brief Base class providing read-only coefficient access to matrices and arrays.
   * \ingroup Core_Module
   * \tparam Derived Type of the derived class
-  * \tparam ReadOnlyAccessors Constant indicating read-only access
+  * \tparam #ReadOnlyAccessors Constant indicating read-only access
   *
   * This class defines the \c operator() \c const function and friends, which can be used to read specific
   * entries of a matrix or array.
@@ -212,7 +212,7 @@ class DenseCoeffsBase<Derived,ReadOnlyAccessors> : public EigenBase<Derived>
       * to ensure that a packet really starts there. This method is only available on expressions having the
       * PacketAccessBit.
       *
-      * The \a LoadMode parameter may have the value \a Aligned or \a Unaligned. Its effect is to select
+      * The \a LoadMode parameter may have the value \a #Aligned or \a #Unaligned. Its effect is to select
       * the appropriate vectorization instruction. Aligned access is faster, but is only possible for packets
       * starting at an address which is a multiple of the packet size.
       */
@@ -239,7 +239,7 @@ class DenseCoeffsBase<Derived,ReadOnlyAccessors> : public EigenBase<Derived>
       * to ensure that a packet really starts there. This method is only available on expressions having the
       * PacketAccessBit and the LinearAccessBit.
       *
-      * The \a LoadMode parameter may have the value \a Aligned or \a Unaligned. Its effect is to select
+      * The \a LoadMode parameter may have the value \a #Aligned or \a #Unaligned. Its effect is to select
       * the appropriate vectorization instruction. Aligned access is faster, but is only possible for packets
       * starting at an address which is a multiple of the packet size.
       */
@@ -275,7 +275,7 @@ class DenseCoeffsBase<Derived,ReadOnlyAccessors> : public EigenBase<Derived>
 /** \brief Base class providing read/write coefficient access to matrices and arrays.
   * \ingroup Core_Module
   * \tparam Derived Type of the derived class
-  * \tparam WriteAccessors Constant indicating read/write access
+  * \tparam #WriteAccessors Constant indicating read/write access
   *
   * This class defines the non-const \c operator() function and friends, which can be used to write specific
   * entries of a matrix or array. This class inherits DenseCoeffsBase<Derived, ReadOnlyAccessors> which
@@ -433,7 +433,7 @@ class DenseCoeffsBase<Derived, WriteAccessors> : public DenseCoeffsBase<Derived,
       * to ensure that a packet really starts there. This method is only available on expressions having the
       * PacketAccessBit.
       *
-      * The \a LoadMode parameter may have the value \a Aligned or \a Unaligned. Its effect is to select
+      * The \a LoadMode parameter may have the value \a #Aligned or \a #Unaligned. Its effect is to select
       * the appropriate vectorization instruction. Aligned access is faster, but is only possible for packets
       * starting at an address which is a multiple of the packet size.
       */
@@ -567,7 +567,7 @@ class DenseCoeffsBase<Derived, WriteAccessors> : public DenseCoeffsBase<Derived,
 /** \brief Base class providing direct read-only coefficient access to matrices and arrays.
   * \ingroup Core_Module
   * \tparam Derived Type of the derived class
-  * \tparam DirectAccessors Constant indicating direct access
+  * \tparam #DirectAccessors Constant indicating direct access
   *
   * This class defines functions to work with strides which can be used to access entries directly. This class
   * inherits DenseCoeffsBase<Derived, ReadOnlyAccessors> which defines functions to access entries read-only using
@@ -637,7 +637,7 @@ class DenseCoeffsBase<Derived, DirectAccessors> : public DenseCoeffsBase<Derived
 /** \brief Base class providing direct read/write coefficient access to matrices and arrays.
   * \ingroup Core_Module
   * \tparam Derived Type of the derived class
-  * \tparam DirectAccessors Constant indicating direct access
+  * \tparam #DirectWriteAccessors Constant indicating direct access
   *
   * This class defines functions to work with strides which can be used to access entries directly. This class
   * inherits DenseCoeffsBase<Derived, WriteAccessors> which defines functions to access entries read/write using

Eigen/src/Core/DenseStorage.h

@@ -58,7 +58,7 @@ struct plain_array
   #define EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(sizemask) \
     eigen_assert((reinterpret_cast<size_t>(array) & sizemask) == 0 \
               && "this assertion is explained here: " \
-              "http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html" \
+              "http://eigen.tuxfamily.org/dox-devel/TopicUnalignedArrayAssert.html" \
               " **** READ THIS WEB PAGE !!! ****");
 #endif
 

Eigen/src/Core/Diagonal.h

@@ -87,7 +87,7 @@ template<typename MatrixType, int DiagIndex> class Diagonal
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Diagonal)
 
     inline Index rows() const
-    { return m_index.value()<0 ? std::min(m_matrix.cols(),m_matrix.rows()+m_index.value()) : std::min(m_matrix.rows(),m_matrix.cols()-m_index.value()); }
+    { return m_index.value()<0 ? (std::min)(m_matrix.cols(),m_matrix.rows()+m_index.value()) : (std::min)(m_matrix.rows(),m_matrix.cols()-m_index.value()); }
 
     inline Index cols() const { return 1; }
 

Eigen/src/Core/Dot.h

@@ -116,7 +116,9 @@ MatrixBase<Derived>::eigen2_dot(const MatrixBase<OtherDerived>& other) const
 
 //---------- implementation of L2 norm and related functions ----------
 
-/** \returns the squared \em l2 norm of *this, i.e., for vectors, the dot product of *this with itself.
+/** \returns, for vectors, the squared \em l2 norm of \c *this, and for matrices the Frobenius norm.
+  * In both cases, it consists in the sum of the square of all the matrix entries.
+  * For vectors, this is also equals to the dot product of \c *this with itself.
   *
   * \sa dot(), norm()
   */
@@ -126,7 +128,9 @@ EIGEN_STRONG_INLINE typename NumTraits<typename internal::traits<Derived>::Scala
   return internal::real((*this).cwiseAbs2().sum());
 }
 
-/** \returns the \em l2 norm of *this, i.e., for vectors, the square root of the dot product of *this with itself.
+/** \returns, for vectors, the \em l2 norm of \c *this, and for matrices the Frobenius norm.
+  * In both cases, it consists in the square root of the sum of the square of all the matrix entries.
+  * For vectors, this is also equals to the square root of the dot product of \c *this with itself.
   *
   * \sa dot(), squaredNorm()
   */

Eigen/src/Core/Functors.h

@@ -116,7 +116,7 @@ struct functor_traits<scalar_conj_product_op<LhsScalar,RhsScalar> > {
   */
 template<typename Scalar> struct scalar_min_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_min_op)
-  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return std::min(a, b); }
+  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { using std::min; return (min)(a, b); }
   template<typename Packet>
   EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const
   { return internal::pmin(a,b); }
@@ -139,7 +139,7 @@ struct functor_traits<scalar_min_op<Scalar> > {
   */
 template<typename Scalar> struct scalar_max_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_max_op)
-  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return std::max(a, b); }
+  EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { using std::max; return (max)(a, b); }
   template<typename Packet>
   EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const
   { return internal::pmax(a,b); }
@@ -165,8 +165,10 @@ template<typename Scalar> struct scalar_hypot_op {
 //   typedef typename NumTraits<Scalar>::Real result_type;
   EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& _x, const Scalar& _y) const
   {
-    Scalar p = std::max(_x, _y);
-    Scalar q = std::min(_x, _y);
+    using std::max;
+    using std::min;
+    Scalar p = (max)(_x, _y);
+    Scalar q = (min)(_x, _y);
     Scalar qp = q/p;
     return p * sqrt(Scalar(1) + qp*qp);
   }
@@ -605,7 +607,7 @@ template <typename Scalar, bool RandomAccess> struct linspaced_op
   EIGEN_STRONG_INLINE const Packet packetOp(Index row, Index col) const
   {
     eigen_assert(col==0 || row==0);
-    return impl(col + row);
+    return impl.packetOp(col + row);
   }
 
   // This proxy object handles the actual required temporaries, the different
@@ -750,9 +752,9 @@ struct functor_traits<scalar_acos_op<Scalar> >
   */
 template<typename Scalar> struct scalar_asin_op {
   EIGEN_EMPTY_STRUCT_CTOR(scalar_asin_op)
-  inline const Scalar operator() (const Scalar& a) const { return acos(a); }
+  inline const Scalar operator() (const Scalar& a) const { return asin(a); }
   typedef typename packet_traits<Scalar>::type Packet;
-  inline Packet packetOp(const Packet& a) const { return internal::pacos(a); }
+  inline Packet packetOp(const Packet& a) const { return internal::pasin(a); }
 };
 template<typename Scalar>
 struct functor_traits<scalar_asin_op<Scalar> >

Eigen/src/Core/Fuzzy.h

@@ -34,9 +34,10 @@ struct isApprox_selector
 {
   static bool run(const Derived& x, const OtherDerived& y, typename Derived::RealScalar prec)
   {
+    using std::min;
     const typename internal::nested<Derived,2>::type nested(x);
     const typename internal::nested<OtherDerived,2>::type otherNested(y);
-    return (nested - otherNested).cwiseAbs2().sum() <= prec * prec * std::min(nested.cwiseAbs2().sum(), otherNested.cwiseAbs2().sum());
+    return (nested - otherNested).cwiseAbs2().sum() <= prec * prec * (min)(nested.cwiseAbs2().sum(), otherNested.cwiseAbs2().sum());
   }
 };
 

Eigen/src/Core/GenericPacketMath.h

@@ -134,12 +134,12 @@ pdiv(const Packet& a,
 /** \internal \returns the min of \a a and \a b  (coeff-wise) */
 template<typename Packet> inline Packet
 pmin(const Packet& a,
-        const Packet& b) { return std::min(a, b); }
+        const Packet& b) { using std::min; return (min)(a, b); }
 
 /** \internal \returns the max of \a a and \a b  (coeff-wise) */
 template<typename Packet> inline Packet
 pmax(const Packet& a,
-        const Packet& b) { return std::max(a, b); }
+        const Packet& b) { using std::max; return (max)(a, b); }
 
 /** \internal \returns the absolute value of \a a */
 template<typename Packet> inline Packet
@@ -286,7 +286,7 @@ pmadd(const Packet&  a,
 { return padd(pmul(a, b),c); }
 
 /** \internal \returns a packet version of \a *from.
-  * \If LoadMode equals Aligned, \a from must be 16 bytes aligned */
+  * If LoadMode equals #Aligned, \a from must be 16 bytes aligned */
 template<typename Packet, int LoadMode>
 inline Packet ploadt(const typename unpacket_traits<Packet>::type* from)
 {
@@ -297,7 +297,7 @@ inline Packet ploadt(const typename unpacket_traits<Packet>::type* from)
 }
 
 /** \internal copy the packet \a from to \a *to.
-  * If StoreMode equals Aligned, \a to must be 16 bytes aligned */
+  * If StoreMode equals #Aligned, \a to must be 16 bytes aligned */
 template<typename Scalar, typename Packet, int LoadMode>
 inline void pstoret(Scalar* to, const Packet& from)
 {

Eigen/src/Core/IO.h

@@ -141,7 +141,8 @@ struct significant_decimals_default_impl
   typedef typename NumTraits<Scalar>::Real RealScalar;
   static inline int run()
   {
-    return cast<RealScalar,int>(std::ceil(-log(NumTraits<RealScalar>::epsilon())/log(RealScalar(10))));
+    using std::ceil;
+    return cast<RealScalar,int>(ceil(-log(NumTraits<RealScalar>::epsilon())/log(RealScalar(10))));
   }
 };
 

Eigen/src/Core/Map.h

@@ -31,10 +31,10 @@
   *
   * \brief A matrix or vector expression mapping an existing array of data.
   *
-  * \param PlainObjectType the equivalent matrix type of the mapped data
-  * \param MapOptions specifies whether the pointer is \c Aligned, or \c Unaligned.
-  *                The default is \c Unaligned.
-  * \param StrideType optionnally specifies strides. By default, Map assumes the memory layout
+  * \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 optionally specifies strides. By default, Map assumes the memory layout
   *                   of an ordinary, contiguous array. This can be overridden by specifying strides.
   *                   The type passed here must be a specialization of the Stride template, see examples below.
   *
@@ -72,9 +72,9 @@
   * Example: \include Map_placement_new.cpp
   * Output: \verbinclude Map_placement_new.out
   *
-  * This class is the return type of Matrix::Map() but can also be used directly.
+  * This class is the return type of PlainObjectBase::Map() but can also be used directly.
   *
-  * \sa Matrix::Map(), \ref TopicStorageOrders
+  * \sa PlainObjectBase::Map(), \ref TopicStorageOrders
   */
 
 namespace internal {
@@ -102,7 +102,7 @@ struct traits<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:

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;
@@ -238,7 +238,7 @@ template<typename Derived> class MapBase<Derived, WriteAccessors>
                 (this->m_data + index * innerStride(), x);
     }
 
-    inline MapBase(PointerType data) : Base(data) {}
+    explicit inline MapBase(PointerType data) : Base(data) {}
     inline MapBase(PointerType data, Index size) : Base(data, size) {}
     inline MapBase(PointerType data, Index rows, Index cols) : Base(data, rows, cols) {}
 

Eigen/src/Core/MathFunctions.h

@@ -87,7 +87,8 @@ struct real_impl<std::complex<RealScalar> >
 {
   static inline RealScalar run(const std::complex<RealScalar>& x)
   {
-    return std::real(x);
+    using std::real;
+    return real(x);
   }
 };
 
@@ -122,7 +123,8 @@ struct imag_impl<std::complex<RealScalar> >
 {
   static inline RealScalar run(const std::complex<RealScalar>& x)
   {
-    return std::imag(x);
+    using std::imag;
+    return imag(x);
   }
 };
 
@@ -244,7 +246,8 @@ struct conj_impl<std::complex<RealScalar> >
 {
   static inline std::complex<RealScalar> run(const std::complex<RealScalar>& x)
   {
-    return std::conj(x);
+    using std::conj;
+    return conj(x);
   }
 };
 
@@ -270,7 +273,8 @@ struct abs_impl
   typedef typename NumTraits<Scalar>::Real RealScalar;
   static inline RealScalar run(const Scalar& x)
   {
-    return std::abs(x);
+    using std::abs;
+    return abs(x);
   }
 };
 
@@ -305,7 +309,8 @@ struct abs2_impl<std::complex<RealScalar> >
 {
   static inline RealScalar run(const std::complex<RealScalar>& x)
   {
-    return std::norm(x);
+    using std::norm;
+    return norm(x);
   }
 };
 
@@ -369,10 +374,12 @@ struct hypot_impl
   typedef typename NumTraits<Scalar>::Real RealScalar;
   static inline RealScalar run(const Scalar& x, const Scalar& y)
   {
+    using std::max;
+    using std::min;
     RealScalar _x = abs(x);
     RealScalar _y = abs(y);
-    RealScalar p = std::max(_x, _y);
-    RealScalar q = std::min(_x, _y);
+    RealScalar p = (max)(_x, _y);
+    RealScalar q = (min)(_x, _y);
     RealScalar qp = q/p;
     return p * sqrt(RealScalar(1) + qp*qp);
   }
@@ -420,7 +427,8 @@ struct sqrt_default_impl
 {
   static inline Scalar run(const Scalar& x)
   {
-    return std::sqrt(x);
+    using std::sqrt;
+    return sqrt(x);
   }
 };
 
@@ -460,7 +468,7 @@ inline EIGEN_MATHFUNC_RETVAL(sqrt, Scalar) sqrt(const Scalar& x)
 // This macro instanciate all the necessary template mechanism which is common to all unary real functions.
 #define EIGEN_MATHFUNC_STANDARD_REAL_UNARY(NAME) \
   template<typename Scalar, bool IsInteger> struct NAME##_default_impl {            \
-    static inline Scalar run(const Scalar& x) { return std::NAME(x); }              \
+    static inline Scalar run(const Scalar& x) { using std::NAME; return NAME(x); }  \
   };                                                                                \
   template<typename Scalar> struct NAME##_default_impl<Scalar, true> {              \
     static inline Scalar run(const Scalar&) {                                       \
@@ -495,7 +503,8 @@ struct atan2_default_impl
   typedef Scalar retval;
   static inline Scalar run(const Scalar& x, const Scalar& y)
   {
-    return std::atan2(x, y);
+    using std::atan2;
+    return atan2(x, y);
   }
 };
 
@@ -534,7 +543,8 @@ struct pow_default_impl
   typedef Scalar retval;
   static inline Scalar run(const Scalar& x, const Scalar& y)
   {
-    return std::pow(x, y);
+    using std::pow;
+    return pow(x, y);
   }
 };
 
@@ -726,7 +736,8 @@ struct scalar_fuzzy_default_impl<Scalar, false, false>
   }
   static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec)
   {
-    return abs(x - y) <= std::min(abs(x), abs(y)) * prec;
+    using std::min;
+    return abs(x - y) <= (min)(abs(x), abs(y)) * prec;
   }
   static inline bool isApproxOrLessThan(const Scalar& x, const Scalar& y, const RealScalar& prec)
   {
@@ -764,7 +775,8 @@ struct scalar_fuzzy_default_impl<Scalar, true, false>
   }
   static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec)
   {
-    return abs2(x - y) <= std::min(abs2(x), abs2(y)) * prec * prec;
+    using std::min;
+    return abs2(x - y) <= (min)(abs2(x), abs2(y)) * prec * prec;
   }
 };
 

Eigen/src/Core/Matrix.h

@@ -43,8 +43,8 @@
   * \tparam _Cols Number of columns, or \b Dynamic
   *
   * The remaining template parameters are optional -- in most cases you don't have to worry about them.
-  * \tparam _Options \anchor matrix_tparam_options A combination of either \b RowMajor or \b ColMajor, and of either
-  *                 \b AutoAlign or \b DontAlign.
+  * \tparam _Options \anchor matrix_tparam_options A combination of either \b #RowMajor or \b #ColMajor, and of either
+  *                 \b #AutoAlign or \b #DontAlign.
   *                 The former controls \ref TopicStorageOrders "storage order", and defaults to column-major. The latter controls alignment, which is required
   *                 for vectorization. It defaults to aligning matrices except for fixed sizes that aren't a multiple of the packet size.
   * \tparam _MaxRows Maximum number of rows. Defaults to \a _Rows (\ref maxrows "note").
@@ -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

@@ -111,7 +111,7 @@ template<typename Derived> class MatrixBase
 
     /** \returns the size of the main diagonal, which is min(rows(),cols()).
       * \sa rows(), cols(), SizeAtCompileTime. */
-    inline Index diagonalSize() const { return std::min(rows(),cols()); }
+    inline Index diagonalSize() const { return (std::min)(rows(),cols()); }
 
     /** \brief The plain matrix type corresponding to this expression.
       *
@@ -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/NumTraits.h

@@ -87,8 +87,8 @@ template<typename T> struct GenericNumTraits
     // make sure to override this for floating-point types
     return Real(0);
   }
-  inline static T highest() { return std::numeric_limits<T>::max(); }
-  inline static T lowest()  { return IsInteger ? std::numeric_limits<T>::min() : (-std::numeric_limits<T>::max()); }
+  inline static T highest() { return (std::numeric_limits<T>::max)(); }
+  inline static T lowest()  { return IsInteger ? (std::numeric_limits<T>::min)() : (-(std::numeric_limits<T>::max)()); }
   
 #ifdef EIGEN2_SUPPORT
   enum {

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,14 +654,15 @@ 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
     {
       // The storage order does not allow us to use reallocation.
       typename Derived::PlainObject tmp(rows,cols);
-      const Index common_rows = std::min(rows, _this.rows());
-      const Index common_cols = std::min(cols, _this.cols());
+      const Index common_rows = (std::min)(rows, _this.rows());
+      const Index common_cols = (std::min)(cols, _this.cols());
       tmp.block(0,0,common_rows,common_cols) = _this.block(0,0,common_rows,common_cols);
       _this.derived().swap(tmp);
     }
@@ -681,8 +695,8 @@ struct internal::conservative_resize_like_impl
     {
       // The storage order does not allow us to use reallocation.
       typename Derived::PlainObject tmp(other);
-      const Index common_rows = std::min(tmp.rows(), _this.rows());
-      const Index common_cols = std::min(tmp.cols(), _this.cols());
+      const Index common_rows = (std::min)(tmp.rows(), _this.rows());
+      const Index common_cols = (std::min)(tmp.cols(), _this.cols());
       tmp.block(0,0,common_rows,common_cols) = _this.block(0,0,common_rows,common_cols);
       _this.derived().swap(tmp);
     }

Eigen/src/Core/Product.h

@@ -375,8 +375,23 @@ struct gemv_static_vector_if<Scalar,Size,Dynamic,true>
 template<typename Scalar,int Size,int MaxSize>
 struct gemv_static_vector_if<Scalar,Size,MaxSize,true>
 {
+  #if EIGEN_ALIGN_STATICALLY
   internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize),0> m_data;
   EIGEN_STRONG_INLINE Scalar* data() { return m_data.array; }
+  #else
+  // Some architectures cannot align on the stack,
+  // => let's manually enforce alignment by allocating more data and return the address of the first aligned element.
+  enum {
+    ForceAlignment  = internal::packet_traits<Scalar>::Vectorizable,
+    PacketSize      = internal::packet_traits<Scalar>::size
+  };
+  internal::plain_array<Scalar,EIGEN_SIZE_MIN_PREFER_FIXED(Size,MaxSize)+(ForceAlignment?PacketSize:0),0> m_data;
+  EIGEN_STRONG_INLINE Scalar* data() {
+    return ForceAlignment
+            ? reinterpret_cast<Scalar*>((reinterpret_cast<size_t>(m_data.array) & ~(size_t(15))) + 16)
+            : m_data.array;
+  }
+  #endif
 };
 
 template<> struct gemv_selector<OnTheRight,ColMajor,true>
@@ -411,28 +426,21 @@ template<> struct gemv_selector<OnTheRight,ColMajor,true>
 
     gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,MightCannotUseDest> static_dest;
 
-    bool alphaIsCompatible = (!ComplexByReal) || (imag(actualAlpha)==RealScalar(0));
+    // this is written like this (i.e., with a ?:) to workaround an ICE with ICC 12
+    bool alphaIsCompatible = (!ComplexByReal) ? true : (imag(actualAlpha)==RealScalar(0));
     bool evalToDest = EvalToDestAtCompileTime && alphaIsCompatible;
     
     RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
 
-    ResScalar* actualDestPtr;
-    bool freeDestPtr = false;
-    if (evalToDest)
-    {
-      actualDestPtr = &dest.coeffRef(0);
-    }
-    else
+    ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
+                                                  evalToDest ? dest.data() : static_dest.data());
+    
+    if(!evalToDest)
     {
       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       int size = dest.size();
       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       #endif
-      if((actualDestPtr = static_dest.data())==0)
-      {
-        freeDestPtr = true;
-        actualDestPtr = ei_aligned_stack_new(ResScalar,dest.size());
-      }
       if(!alphaIsCompatible)
       {
         MappedDest(actualDestPtr, dest.size()).setZero();
@@ -456,7 +464,6 @@ template<> struct gemv_selector<OnTheRight,ColMajor,true>
         dest += actualAlpha * MappedDest(actualDestPtr, dest.size());
       else
         dest = MappedDest(actualDestPtr, dest.size());
-      if(freeDestPtr) ei_aligned_stack_delete(ResScalar, actualDestPtr, dest.size());
     }
   }
 };
@@ -490,23 +497,15 @@ template<> struct gemv_selector<OnTheRight,RowMajor,true>
 
     gemv_static_vector_if<RhsScalar,_ActualRhsType::SizeAtCompileTime,_ActualRhsType::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs;
 
-    RhsScalar* actualRhsPtr;
-    bool freeRhsPtr = false;
-    if (DirectlyUseRhs)
-    {
-      actualRhsPtr = const_cast<RhsScalar*>(&actualRhs.coeffRef(0));
-    }
-    else
+    ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,actualRhs.size(),
+        DirectlyUseRhs ? const_cast<RhsScalar*>(actualRhs.data()) : static_rhs.data());
+
+    if(!DirectlyUseRhs)
     {
       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       int size = actualRhs.size();
       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       #endif
-      if((actualRhsPtr = static_rhs.data())==0)
-      {
-        freeRhsPtr = true;
-        actualRhsPtr = ei_aligned_stack_new(RhsScalar, actualRhs.size());
-      }
       Map<typename _ActualRhsType::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs;
     }
 
@@ -517,8 +516,6 @@ template<> struct gemv_selector<OnTheRight,RowMajor,true>
         actualRhsPtr, 1,
         &dest.coeffRef(0,0), dest.innerStride(),
         actualAlpha);
-
-    if((!DirectlyUseRhs) && freeRhsPtr) ei_aligned_stack_delete(RhsScalar, actualRhsPtr, prod.rhs().size());
   }
 };
 

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/SelfAdjointView.h

@@ -32,13 +32,13 @@
   * \brief Expression of a selfadjoint matrix from a triangular part of a dense matrix
   *
   * \param MatrixType the type of the dense matrix storing the coefficients
-  * \param TriangularPart can be either \c Lower or \c Upper
+  * \param TriangularPart can be either \c #Lower or \c #Upper
   *
   * This class is an expression of a sefladjoint matrix from a triangular part of a matrix
   * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView()
   * and most of the time this is the only way that it is used.
   *
-  * \sa class TriangularBase, MatrixBase::selfAdjointView()
+  * \sa class TriangularBase, MatrixBase::selfadjointView()
   */
 
 namespace internal {

Eigen/src/Core/SolveTriangular.h

@@ -74,26 +74,19 @@ struct triangular_solver_selector<Lhs,Rhs,Side,Mode,NoUnrolling,1>
     // FIXME find a way to allow an inner stride if packet_traits<Scalar>::size==1
 
     bool useRhsDirectly = Rhs::InnerStrideAtCompileTime==1 || rhs.innerStride()==1;
-    RhsScalar* actualRhs;
-    if(useRhsDirectly)
-    {
-      actualRhs = &rhs.coeffRef(0);
-    }
-    else
-    {
-      actualRhs = ei_aligned_stack_new(RhsScalar,rhs.size());
+
+    ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhs,rhs.size(),
+                                                  (useRhsDirectly ? rhs.data() : 0));
+                                                  
+    if(!useRhsDirectly)
       MappedRhs(actualRhs,rhs.size()) = rhs;
-    }
 
     triangular_solve_vector<LhsScalar, RhsScalar, typename Lhs::Index, Side, Mode, LhsProductTraits::NeedToConjugate,
                             (int(Lhs::Flags) & RowMajorBit) ? RowMajor : ColMajor>
       ::run(actualLhs.cols(), actualLhs.data(), actualLhs.outerStride(), actualRhs);
 
     if(!useRhsDirectly)
-    {
       rhs = MappedRhs(actualRhs, rhs.size());
-      ei_aligned_stack_delete(RhsScalar, actualRhs, rhs.size());
-    }
   }
 };
 
@@ -187,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/StableNorm.h

@@ -56,6 +56,7 @@ template<typename Derived>
 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
 MatrixBase<Derived>::stableNorm() const
 {
+  using std::min;
   const Index blockSize = 4096;
   RealScalar scale = 0;
   RealScalar invScale = 1;
@@ -68,7 +69,7 @@ MatrixBase<Derived>::stableNorm() const
   if (bi>0)
     internal::stable_norm_kernel(this->head(bi), ssq, scale, invScale);
   for (; bi<n; bi+=blockSize)
-    internal::stable_norm_kernel(this->segment(bi,std::min(blockSize, n - bi)).template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale);
+    internal::stable_norm_kernel(this->segment(bi,(min)(blockSize, n - bi)).template forceAlignedAccessIf<Alignment>(), ssq, scale, invScale);
   return scale * internal::sqrt(ssq);
 }
 
@@ -85,6 +86,9 @@ template<typename Derived>
 inline typename NumTraits<typename internal::traits<Derived>::Scalar>::Real
 MatrixBase<Derived>::blueNorm() const
 {
+  using std::pow;
+  using std::min;
+  using std::max;
   static Index nmax = -1;
   static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr;
   if(nmax <= 0)
@@ -99,25 +103,25 @@ MatrixBase<Derived>::blueNorm() const
     // For portability, the PORT subprograms "ilmaeh" and "rlmach"
     // are used. For any specific computer, each of the assignment
     // statements can be replaced
-    nbig  = std::numeric_limits<Index>::max();            // largest integer
+    nbig  = (std::numeric_limits<Index>::max)();            // largest integer
     ibeta = std::numeric_limits<RealScalar>::radix;         // base for floating-point numbers
     it    = std::numeric_limits<RealScalar>::digits;        // number of base-beta digits in mantissa
     iemin = std::numeric_limits<RealScalar>::min_exponent;  // minimum exponent
     iemax = std::numeric_limits<RealScalar>::max_exponent;  // maximum exponent
-    rbig  = std::numeric_limits<RealScalar>::max();         // largest floating-point number
+    rbig  = (std::numeric_limits<RealScalar>::max)();         // largest floating-point number
 
     iexp  = -((1-iemin)/2);
-    b1    = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp)));  // lower boundary of midrange
+    b1    = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));  // lower boundary of midrange
     iexp  = (iemax + 1 - it)/2;
-    b2    = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp)));   // upper boundary of midrange
+    b2    = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));   // upper boundary of midrange
 
     iexp  = (2-iemin)/2;
-    s1m   = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp)));   // scaling factor for lower range
+    s1m   = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));   // scaling factor for lower range
     iexp  = - ((iemax+it)/2);
-    s2m   = RealScalar(std::pow(RealScalar(ibeta),RealScalar(iexp)));   // scaling factor for upper range
+    s2m   = RealScalar(pow(RealScalar(ibeta),RealScalar(iexp)));   // scaling factor for upper range
 
     overfl  = rbig*s2m;             // overflow boundary for abig
-    eps     = RealScalar(std::pow(double(ibeta), 1-it));
+    eps     = RealScalar(pow(double(ibeta), 1-it));
     relerr  = internal::sqrt(eps);         // tolerance for neglecting asml
     abig    = RealScalar(1.0/eps - 1.0);
     if (RealScalar(nbig)>abig)  nmax = int(abig);  // largest safe n
@@ -163,8 +167,8 @@ MatrixBase<Derived>::blueNorm() const
   }
   else
     return internal::sqrt(amed);
-  asml = std::min(abig, amed);
-  abig = std::max(abig, amed);
+  asml = (min)(abig, amed);
+  abig = (max)(abig, amed);
   if(asml <= abig*relerr)
     return abig;
   else

Eigen/src/Core/Transpose.h

@@ -350,15 +350,14 @@ struct blas_traits<SelfCwiseBinaryOp<BinOp,NestedXpr,Rhs> >
 template<bool DestIsTransposed, typename OtherDerived>
 struct check_transpose_aliasing_compile_time_selector
 {
-  enum { ret = blas_traits<OtherDerived>::IsTransposed != DestIsTransposed
-  };
+  enum { ret = bool(blas_traits<OtherDerived>::IsTransposed) != DestIsTransposed };
 };
 
 template<bool DestIsTransposed, typename BinOp, typename DerivedA, typename DerivedB>
 struct check_transpose_aliasing_compile_time_selector<DestIsTransposed,CwiseBinaryOp<BinOp,DerivedA,DerivedB> >
 {
-  enum { ret =    blas_traits<DerivedA>::IsTransposed != DestIsTransposed
-               || blas_traits<DerivedB>::IsTransposed != DestIsTransposed
+  enum { ret =    bool(blas_traits<DerivedA>::IsTransposed) != DestIsTransposed
+               || bool(blas_traits<DerivedB>::IsTransposed) != DestIsTransposed
   };
 };
 
@@ -367,7 +366,7 @@ struct check_transpose_aliasing_run_time_selector
 {
   static bool run(const Scalar* dest, const OtherDerived& src)
   {
-    return (blas_traits<OtherDerived>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(Scalar*)extract_data(src));
+    return (bool(blas_traits<OtherDerived>::IsTransposed) != DestIsTransposed) && (dest!=0 && dest==(Scalar*)extract_data(src));
   }
 };
 

Eigen/src/Core/TriangularMatrix.h

@@ -111,6 +111,7 @@ template<typename Derived> class TriangularBase : public EigenBase<Derived>
       EIGEN_ONLY_USED_FOR_DEBUG(col);
       eigen_assert(col>=0 && col<cols() && row>=0 && row<rows());
       const int mode = int(Mode) & ~SelfAdjoint;
+      EIGEN_ONLY_USED_FOR_DEBUG(mode);
       eigen_assert((mode==Upper && col>=row)
                 || (mode==Lower && col<=row)
                 || ((mode==StrictlyUpper || mode==UnitUpper) && col>row)
@@ -134,13 +135,13 @@ template<typename Derived> class TriangularBase : public EigenBase<Derived>
   * \brief Base class for triangular part in a matrix
   *
   * \param MatrixType the type of the object in which we are taking the triangular part
-  * \param Mode the kind of triangular matrix expression to construct. Can be Upper,
-  *             Lower, UpperSelfadjoint, or LowerSelfadjoint. This is in fact a bit field;
-  *             it must have either Upper or Lower, and additionnaly it may have either
-  *             UnitDiag or Selfadjoint.
+  * \param Mode the kind of triangular matrix expression to construct. Can be #Upper,
+  *             #Lower, #UnitUpper, #UnitLower, #StrictlyUpper, or #StrictlyLower.
+  *             This is in fact a bit field; it must have either #Upper or #Lower, 
+  *             and additionnaly it may have #UnitDiag or #ZeroDiag or neither.
   *
   * This class represents a triangular part of a matrix, not necessarily square. Strictly speaking, for rectangular
-  * matrices one should speak ok "trapezoid" parts. This class is the return type
+  * matrices one should speak of "trapezoid" parts. This class is the return type
   * of MatrixBase::triangularView() and most of the time this is the only way it is used.
   *
   * \sa MatrixBase::triangularView()
@@ -448,6 +449,8 @@ struct triangular_assignment_selector
     col = (UnrollCount-1) / Derived1::RowsAtCompileTime,
     row = (UnrollCount-1) % Derived1::RowsAtCompileTime
   };
+  
+  typedef typename Derived1::Scalar Scalar;
 
   inline static void run(Derived1 &dst, const Derived2 &src)
   {
@@ -466,9 +469,9 @@ struct triangular_assignment_selector
     else if(ClearOpposite)
     {
       if (Mode&UnitDiag && row==col)
-        dst.coeffRef(row, col) = 1;
+        dst.coeffRef(row, col) = Scalar(1);
       else
-        dst.coeffRef(row, col) = 0;
+        dst.coeffRef(row, col) = Scalar(0);
     }
   }
 };
@@ -484,16 +487,17 @@ template<typename Derived1, typename Derived2, bool ClearOpposite>
 struct triangular_assignment_selector<Derived1, Derived2, Upper, Dynamic, ClearOpposite>
 {
   typedef typename Derived1::Index Index;
+  typedef typename Derived1::Scalar Scalar;
   inline static void run(Derived1 &dst, const Derived2 &src)
   {
     for(Index j = 0; j < dst.cols(); ++j)
     {
-      Index maxi = std::min(j, dst.rows()-1);
+      Index maxi = (std::min)(j, dst.rows()-1);
       for(Index i = 0; i <= maxi; ++i)
         dst.copyCoeff(i, j, src);
       if (ClearOpposite)
         for(Index i = maxi+1; i < dst.rows(); ++i)
-          dst.coeffRef(i, j) = 0;
+          dst.coeffRef(i, j) = Scalar(0);
     }
   }
 };
@@ -508,10 +512,10 @@ struct triangular_assignment_selector<Derived1, Derived2, Lower, Dynamic, ClearO
     {
       for(Index i = j; i < dst.rows(); ++i)
         dst.copyCoeff(i, j, src);
-      Index maxi = std::min(j, dst.rows());
+      Index maxi = (std::min)(j, dst.rows());
       if (ClearOpposite)
         for(Index i = 0; i < maxi; ++i)
-          dst.coeffRef(i, j) = 0;
+          dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
     }
   }
 };
@@ -524,7 +528,7 @@ struct triangular_assignment_selector<Derived1, Derived2, StrictlyUpper, Dynamic
   {
     for(Index j = 0; j < dst.cols(); ++j)
     {
-      Index maxi = std::min(j, dst.rows());
+      Index maxi = (std::min)(j, dst.rows());
       for(Index i = 0; i < maxi; ++i)
         dst.copyCoeff(i, j, src);
       if (ClearOpposite)
@@ -544,10 +548,10 @@ struct triangular_assignment_selector<Derived1, Derived2, StrictlyLower, Dynamic
     {
       for(Index i = j+1; i < dst.rows(); ++i)
         dst.copyCoeff(i, j, src);
-      Index maxi = std::min(j, dst.rows()-1);
+      Index maxi = (std::min)(j, dst.rows()-1);
       if (ClearOpposite)
         for(Index i = 0; i <= maxi; ++i)
-          dst.coeffRef(i, j) = 0;
+          dst.coeffRef(i, j) = static_cast<typename Derived1::Scalar>(0);
     }
   }
 };
@@ -560,7 +564,7 @@ struct triangular_assignment_selector<Derived1, Derived2, UnitUpper, Dynamic, Cl
   {
     for(Index j = 0; j < dst.cols(); ++j)
     {
-      Index maxi = std::min(j, dst.rows());
+      Index maxi = (std::min)(j, dst.rows());
       for(Index i = 0; i < maxi; ++i)
         dst.copyCoeff(i, j, src);
       if (ClearOpposite)
@@ -580,7 +584,7 @@ struct triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic, Cl
   {
     for(Index j = 0; j < dst.cols(); ++j)
     {
-      Index maxi = std::min(j, dst.rows());
+      Index maxi = (std::min)(j, dst.rows());
       for(Index i = maxi+1; i < dst.rows(); ++i)
         dst.copyCoeff(i, j, src);
       if (ClearOpposite)
@@ -756,8 +760,8 @@ typename internal::eigen2_part_return_type<Derived, Mode>::type MatrixBase<Deriv
 /**
   * \returns an expression of a triangular view extracted from the current matrix
   *
-  * The parameter \a Mode can have the following values: \c Upper, \c StrictlyUpper, \c UnitUpper,
-  * \c Lower, \c StrictlyLower, \c UnitLower.
+  * The parameter \a Mode can have the following values: \c #Upper, \c #StrictlyUpper, \c #UnitUpper,
+  * \c #Lower, \c #StrictlyLower, \c #UnitLower.
   *
   * Example: \include MatrixBase_extract.cpp
   * Output: \verbinclude MatrixBase_extract.out
@@ -792,7 +796,7 @@ bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const
   RealScalar maxAbsOnUpperPart = static_cast<RealScalar>(-1);
   for(Index j = 0; j < cols(); ++j)
   {
-    Index maxi = std::min(j, rows()-1);
+    Index maxi = (std::min)(j, rows()-1);
     for(Index i = 0; i <= maxi; ++i)
     {
       RealScalar absValue = internal::abs(coeff(i,j));
@@ -824,7 +828,7 @@ bool MatrixBase<Derived>::isLowerTriangular(RealScalar prec) const
   RealScalar threshold = maxAbsOnLowerPart * prec;
   for(Index j = 1; j < cols(); ++j)
   {
-    Index maxi = std::min(j, rows()-1);
+    Index maxi = (std::min)(j, rows()-1);
     for(Index i = 0; i < maxi; ++i)
       if(internal::abs(coeff(i, j)) > threshold) return false;
   }

Eigen/src/Core/VectorwiseOp.h

@@ -31,9 +31,9 @@
   *
   * \brief Generic expression of a partially reduxed matrix
   *
-  * \param MatrixType the type of the matrix we are applying the redux operation
-  * \param MemberOp type of the member functor
-  * \param Direction indicates the direction of the redux (Vertical or Horizontal)
+  * \tparam MatrixType the type of the matrix we are applying the redux operation
+  * \tparam MemberOp type of the member functor
+  * \tparam Direction indicates the direction of the redux (#Vertical or #Horizontal)
   *
   * This class represents an expression of a partial redux operator of a matrix.
   * It is the return type of some VectorwiseOp functions,
@@ -164,7 +164,7 @@ struct member_redux {
   * \brief Pseudo expression providing partial reduction operations
   *
   * \param ExpressionType the type of the object on which to do partial reductions
-  * \param Direction indicates the direction of the redux (Vertical or Horizontal)
+  * \param Direction indicates the direction of the redux (#Vertical or #Horizontal)
   *
   * This class represents a pseudo expression with partial reduction features.
   * It is the return type of DenseBase::colwise() and DenseBase::rowwise()

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
@@ -43,6 +43,7 @@ template<> struct packet_traits<std::complex<float> >  : default_packet_traits
   typedef Packet2cf type;
   enum {
     Vectorizable = 1,
+    AlignedOnScalar = 1,
     size = 2,
 
     HasAdd    = 1,

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 = { 3, 2, 1, 0 };
+  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 = { 3, 2, 1, 0 };
+  Packet4i countdown = EIGEN_INIT_NEON_PACKET4(0, 1, 2, 3);
   return vaddq_s32(pset1<Packet4i>(a), countdown);
 }
 
@@ -191,14 +201,14 @@ template<> EIGEN_STRONG_INLINE Packet4f ploaddup<Packet4f>(const float*   from)
 {
   float32x2_t lo, hi;
   lo = vdup_n_f32(*from);
-  hi = vdup_n_f32(*from);
+  hi = vdup_n_f32(*(from+1));
   return vcombine_f32(lo, hi);
 }
 template<> EIGEN_STRONG_INLINE Packet4i ploaddup<Packet4i>(const int*     from)
 {
   int32x2_t lo, hi;
   lo = vdup_n_s32(*from);
-  hi = vdup_n_s32(*from);
+  hi = vdup_n_s32(*(from+1));
   return vcombine_s32(lo, hi);
 }
 
@@ -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

@@ -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)
   {
@@ -81,6 +78,7 @@ inline void manage_caching_sizes(Action action, std::ptrdiff_t* l1=0, std::ptrdi
 template<typename LhsScalar, typename RhsScalar, int KcFactor>
 void computeProductBlockingSizes(std::ptrdiff_t& k, std::ptrdiff_t& m, std::ptrdiff_t& n)
 {
+  EIGEN_UNUSED_VARIABLE(n);
   // Explanations:
   // Let's recall the product algorithms form kc x nc horizontal panels B' on the rhs and
   // mc x kc blocks A' on the lhs. A' has to fit into L2 cache. Moreover, B' is processed
@@ -102,7 +100,6 @@ void computeProductBlockingSizes(std::ptrdiff_t& k, std::ptrdiff_t& m, std::ptrd
   k = std::min<std::ptrdiff_t>(k, l1/kdiv);
   std::ptrdiff_t _m = k>0 ? l2/(4 * sizeof(LhsScalar) * k) : 0;
   if(_m<m) m = _m & mr_mask;
-  n = n;
 }
 
 template<typename LhsScalar, typename RhsScalar>
@@ -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);
     }

Eigen/src/Core/products/GeneralMatrixMatrix.h

@@ -78,7 +78,7 @@ static void run(Index rows, Index cols, Index depth,
   typedef gebp_traits<LhsScalar,RhsScalar> Traits;
 
   Index kc = blocking.kc();                 // cache block size along the K direction
-  Index mc = std::min(rows,blocking.mc());  // cache block size along the M direction
+  Index mc = (std::min)(rows,blocking.mc());  // cache block size along the M direction
   //Index nc = blocking.nc(); // cache block size along the N direction
 
   gemm_pack_lhs<LhsScalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
@@ -94,15 +94,16 @@ static void run(Index rows, Index cols, Index depth,
     
     std::size_t sizeA = kc*mc;
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
-    LhsScalar* blockA = ei_aligned_stack_new(LhsScalar, sizeA);
-    RhsScalar* w = ei_aligned_stack_new(RhsScalar, sizeW);
+    ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, sizeA, 0);
+    ei_declare_aligned_stack_constructed_variable(RhsScalar, w, sizeW, 0);
+    
     RhsScalar* blockB = blocking.blockB();
     eigen_internal_assert(blockB!=0);
 
     // For each horizontal panel of the rhs, and corresponding vertical panel of the lhs...
     for(Index k=0; k<depth; k+=kc)
     {
-      const Index actual_kc = std::min(k+kc,depth)-k; // => rows of B', and cols of the A'
+      const Index actual_kc = (std::min)(k+kc,depth)-k; // => rows of B', and cols of the A'
 
       // In order to reduce the chance that a thread has to wait for the other,
       // let's start by packing A'.
@@ -139,7 +140,7 @@ static void run(Index rows, Index cols, Index depth,
       // Then keep going as usual with the remaining A'
       for(Index i=mc; i<rows; i+=mc)
       {
-        const Index actual_mc = std::min(i+mc,rows)-i;
+        const Index actual_mc = (std::min)(i+mc,rows)-i;
 
         // pack A_i,k to A'
         pack_lhs(blockA, &lhs(i,k), lhsStride, actual_kc, actual_mc);
@@ -154,9 +155,6 @@ static void run(Index rows, Index cols, Index depth,
         #pragma omp atomic
         --(info[j].users);
     }
-
-    ei_aligned_stack_delete(LhsScalar, blockA, kc*mc);
-    ei_aligned_stack_delete(RhsScalar, w, sizeW);
   }
   else
 #endif // EIGEN_HAS_OPENMP
@@ -167,15 +165,16 @@ static void run(Index rows, Index cols, Index depth,
     std::size_t sizeA = kc*mc;
     std::size_t sizeB = kc*cols;
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
-    LhsScalar *blockA = blocking.blockA()==0 ? ei_aligned_stack_new(LhsScalar, sizeA) : blocking.blockA();
-    RhsScalar *blockB = blocking.blockB()==0 ? ei_aligned_stack_new(RhsScalar, sizeB) : blocking.blockB();
-    RhsScalar *blockW = blocking.blockW()==0 ? ei_aligned_stack_new(RhsScalar, sizeW) : blocking.blockW();
+
+    ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, sizeA, blocking.blockA());
+    ei_declare_aligned_stack_constructed_variable(RhsScalar, blockB, sizeB, blocking.blockB());
+    ei_declare_aligned_stack_constructed_variable(RhsScalar, blockW, sizeW, blocking.blockW());
 
     // For each horizontal panel of the rhs, and corresponding panel of the lhs...
     // (==GEMM_VAR1)
     for(Index k2=0; k2<depth; k2+=kc)
     {
-      const Index actual_kc = std::min(k2+kc,depth)-k2;
+      const Index actual_kc = (std::min)(k2+kc,depth)-k2;
 
       // OK, here we have selected one horizontal panel of rhs and one vertical panel of lhs.
       // => Pack rhs's panel into a sequential chunk of memory (L2 caching)
@@ -188,7 +187,7 @@ static void run(Index rows, Index cols, Index depth,
       // (==GEPP_VAR1)
       for(Index i2=0; i2<rows; i2+=mc)
       {
-        const Index actual_mc = std::min(i2+mc,rows)-i2;
+        const Index actual_mc = (std::min)(i2+mc,rows)-i2;
 
         // We pack the lhs's block into a sequential chunk of memory (L1 caching)
         // Note that this block will be read a very high number of times, which is equal to the number of
@@ -200,10 +199,6 @@ static void run(Index rows, Index cols, Index depth,
 
       }
     }
-
-    if(blocking.blockA()==0) ei_aligned_stack_delete(LhsScalar, blockA, sizeA);
-    if(blocking.blockB()==0) ei_aligned_stack_delete(RhsScalar, blockB, sizeB);
-    if(blocking.blockW()==0) ei_aligned_stack_delete(RhsScalar, blockW, sizeW);
   }
 }
 

Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h

@@ -83,10 +83,10 @@ struct general_matrix_matrix_triangular_product<Index,LhsScalar,LhsStorageOrder,
     if(mc > Traits::nr)
       mc = (mc/Traits::nr)*Traits::nr;
 
-    LhsScalar* blockA = ei_aligned_stack_new(LhsScalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*size;
-    RhsScalar* allocatedBlockB = ei_aligned_stack_new(RhsScalar, sizeB);
+    ei_declare_aligned_stack_constructed_variable(LhsScalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(RhsScalar, allocatedBlockB, sizeB, 0);
     RhsScalar* blockB = allocatedBlockB + sizeW;
     
     gemm_pack_lhs<LhsScalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder> pack_lhs;
@@ -96,14 +96,14 @@ struct general_matrix_matrix_triangular_product<Index,LhsScalar,LhsStorageOrder,
 
     for(Index k2=0; k2<depth; k2+=kc)
     {
-      const Index actual_kc = std::min(k2+kc,depth)-k2;
+      const Index actual_kc = (std::min)(k2+kc,depth)-k2;
 
       // note that the actual rhs is the transpose/adjoint of mat
       pack_rhs(blockB, &rhs(k2,0), rhsStride, actual_kc, size);
 
       for(Index i2=0; i2<size; i2+=mc)
       {
-        const Index actual_mc = std::min(i2+mc,size)-i2;
+        const Index actual_mc = (std::min)(i2+mc,size)-i2;
 
         pack_lhs(blockA, &lhs(i2, k2), lhsStride, actual_kc, actual_mc);
 
@@ -112,7 +112,7 @@ struct general_matrix_matrix_triangular_product<Index,LhsScalar,LhsStorageOrder,
         //  2 - the actual_mc x actual_mc symmetric block => processed with a special kernel
         //  3 - after the diagonal => processed with gebp or skipped
         if (UpLo==Lower)
-          gebp(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, std::min(size,i2), alpha,
+          gebp(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, (std::min)(size,i2), alpha,
                -1, -1, 0, 0, allocatedBlockB);
 
         sybb(res+resStride*i2 + i2, resStride, blockA, blockB + actual_kc*i2, actual_mc, actual_kc, alpha, allocatedBlockB);
@@ -120,13 +120,11 @@ struct general_matrix_matrix_triangular_product<Index,LhsScalar,LhsStorageOrder,
         if (UpLo==Upper)
         {
           Index j2 = i2+actual_mc;
-          gebp(res+resStride*j2+i2, resStride, blockA, blockB+actual_kc*j2, actual_mc, actual_kc, std::max(Index(0), size-j2), alpha,
+          gebp(res+resStride*j2+i2, resStride, blockA, blockB+actual_kc*j2, actual_mc, actual_kc, (std::max)(Index(0), size-j2), alpha,
                -1, -1, 0, 0, allocatedBlockB);
         }
       }
     }
-    ei_aligned_stack_delete(LhsScalar, blockA, kc*mc);
-    ei_aligned_stack_delete(RhsScalar, allocatedBlockB, sizeB);
   }
 };
 

Eigen/src/Core/products/GeneralMatrixVector.h

@@ -134,7 +134,7 @@ EIGEN_DONT_INLINE static void run(
     }
     else
     {
-      skipColumns = std::min(skipColumns,cols);
+      skipColumns = (std::min)(skipColumns,cols);
       // note that the skiped columns are processed later.
     }
 
@@ -386,7 +386,7 @@ EIGEN_DONT_INLINE static void run(
     }
     else
     {
-      skipRows = std::min(skipRows,Index(rows));
+      skipRows = (std::min)(skipRows,Index(rows));
       // note that the skiped columns are processed later.
     }
     eigen_internal_assert(  alignmentPattern==NoneAligned

Eigen/src/Core/products/SelfadjointMatrixMatrix.h

@@ -114,7 +114,7 @@ struct symm_pack_rhs
     }
 
     // second part: diagonal block
-    for(Index j2=k2; j2<std::min(k2+rows,packet_cols); j2+=nr)
+    for(Index j2=k2; j2<(std::min)(k2+rows,packet_cols); j2+=nr)
     {
       // again we can split vertically in three different parts (transpose, symmetric, normal)
       // transpose
@@ -179,7 +179,7 @@ struct symm_pack_rhs
     for(Index j2=packet_cols; j2<cols; ++j2)
     {
       // transpose
-      Index half = std::min(end_k,j2);
+      Index half = (std::min)(end_k,j2);
       for(Index k=k2; k<half; k++)
       {
         blockB[count] = conj(rhs(j2,k));
@@ -261,12 +261,12 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,ConjugateLhs
     Index nc = cols;  // cache block size along the N direction
     computeProductBlockingSizes<Scalar,Scalar>(kc, mc, nc);
     // kc must smaller than mc
-    kc = std::min(kc,mc);
+    kc = (std::min)(kc,mc);
 
-    Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*cols;
-    Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar, sizeB);
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);
     Scalar* blockB = allocatedBlockB + sizeW;
 
     gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel;
@@ -276,7 +276,7 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,ConjugateLhs
 
     for(Index k2=0; k2<size; k2+=kc)
     {
-      const Index actual_kc = std::min(k2+kc,size)-k2;
+      const Index actual_kc = (std::min)(k2+kc,size)-k2;
 
       // we have selected one row panel of rhs and one column panel of lhs
       // pack rhs's panel into a sequential chunk of memory
@@ -289,7 +289,7 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,ConjugateLhs
       //  3 - the panel below the diagonal block => generic packed copy
       for(Index i2=0; i2<k2; i2+=mc)
       {
-        const Index actual_mc = std::min(i2+mc,k2)-i2;
+        const Index actual_mc = (std::min)(i2+mc,k2)-i2;
         // transposed packed copy
         pack_lhs_transposed(blockA, &lhs(k2, i2), lhsStride, actual_kc, actual_mc);
 
@@ -297,7 +297,7 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,ConjugateLhs
       }
       // the block diagonal
       {
-        const Index actual_mc = std::min(k2+kc,size)-k2;
+        const Index actual_mc = (std::min)(k2+kc,size)-k2;
         // symmetric packed copy
         pack_lhs(blockA, &lhs(k2,k2), lhsStride, actual_kc, actual_mc);
 
@@ -306,16 +306,13 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,ConjugateLhs
 
       for(Index i2=k2+kc; i2<size; i2+=mc)
       {
-        const Index actual_mc = std::min(i2+mc,size)-i2;
+        const Index actual_mc = (std::min)(i2+mc,size)-i2;
         gemm_pack_lhs<Scalar, Index, Traits::mr, Traits::LhsProgress, LhsStorageOrder,false>()
           (blockA, &lhs(i2, k2), lhsStride, actual_kc, actual_mc);
 
         gebp_kernel(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha);
       }
     }
-
-    ei_aligned_stack_delete(Scalar, blockA, kc*mc);
-    ei_aligned_stack_delete(Scalar, allocatedBlockB, sizeB);
   }
 };
 
@@ -343,11 +340,10 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,false,ConjugateLh
     Index mc = rows;  // cache block size along the M direction
     Index nc = cols;  // cache block size along the N direction
     computeProductBlockingSizes<Scalar,Scalar>(kc, mc, nc);
-
-    Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*cols;
-    Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar, sizeB);
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);
     Scalar* blockB = allocatedBlockB + sizeW;
 
     gebp_kernel<Scalar, Scalar, Index, Traits::mr, Traits::nr, ConjugateLhs, ConjugateRhs> gebp_kernel;
@@ -356,22 +352,19 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,false,ConjugateLh
 
     for(Index k2=0; k2<size; k2+=kc)
     {
-      const Index actual_kc = std::min(k2+kc,size)-k2;
+      const Index actual_kc = (std::min)(k2+kc,size)-k2;
 
       pack_rhs(blockB, _rhs, rhsStride, actual_kc, cols, k2);
 
       // => GEPP
       for(Index i2=0; i2<rows; i2+=mc)
       {
-        const Index actual_mc = std::min(i2+mc,rows)-i2;
+        const Index actual_mc = (std::min)(i2+mc,rows)-i2;
         pack_lhs(blockA, &lhs(i2, k2), lhsStride, actual_kc, actual_mc);
 
         gebp_kernel(res+i2, resStride, blockA, blockB, actual_mc, actual_kc, cols, alpha);
       }
     }
-
-    ei_aligned_stack_delete(Scalar, blockA, kc*mc);
-    ei_aligned_stack_delete(Scalar, allocatedBlockB, sizeB);
   }
 };
 

Eigen/src/Core/products/SelfadjointMatrixVector.h

@@ -62,17 +62,15 @@ static EIGEN_DONT_INLINE void product_selfadjoint_vector(
   // FIXME this copy is now handled outside product_selfadjoint_vector, so it could probably be removed.
   // if the rhs is not sequentially stored in memory we copy it to a temporary buffer,
   // this is because we need to extract packets
-  const Scalar* EIGEN_RESTRICT rhs = _rhs;
+  ei_declare_aligned_stack_constructed_variable(Scalar,rhs,size,rhsIncr==1 ? const_cast<Scalar*>(_rhs) : 0);  
   if (rhsIncr!=1)
   {
-    Scalar* r = ei_aligned_stack_new(Scalar, size);
     const Scalar* it = _rhs;
     for (Index i=0; i<size; ++i, it+=rhsIncr)
-      r[i] = *it;
-    rhs = r;
+      rhs[i] = *it;
   }
 
-  Index bound = std::max(Index(0),size-8) & 0xfffffffe;
+  Index bound = (std::max)(Index(0),size-8) & 0xfffffffe;
   if (FirstTriangular)
     bound = size - bound;
 
@@ -160,9 +158,6 @@ static EIGEN_DONT_INLINE void product_selfadjoint_vector(
     }
     res[j] += alpha * t2;
   }
-
-  if(rhsIncr!=1)
-    ei_aligned_stack_delete(Scalar, const_cast<Scalar*>(rhs), size);
 }
 
 } // end namespace internal 
@@ -211,40 +206,28 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true>
     
     internal::gemv_static_vector_if<ResScalar,Dest::SizeAtCompileTime,Dest::MaxSizeAtCompileTime,!EvalToDest> static_dest;
     internal::gemv_static_vector_if<RhsScalar,_ActualRhsType::SizeAtCompileTime,_ActualRhsType::MaxSizeAtCompileTime,!UseRhs> static_rhs;
+
+    ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
+                                                  EvalToDest ? dest.data() : static_dest.data());
+                                                  
+    ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,rhs.size(),
+        UseRhs ? const_cast<RhsScalar*>(rhs.data()) : static_rhs.data());
     
-    bool freeDestPtr = false;
-    ResScalar* actualDestPtr;
-    if(EvalToDest)
-      actualDestPtr = dest.data();
-    else
+    if(!EvalToDest)
     {
       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       int size = dest.size();
       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       #endif
-      if((actualDestPtr=static_dest.data())==0)
-      {
-        freeDestPtr = true;
-        actualDestPtr = ei_aligned_stack_new(ResScalar,dest.size());
-      }
       MappedDest(actualDestPtr, dest.size()) = dest;
     }
       
-    bool freeRhsPtr = false;
-    RhsScalar* actualRhsPtr;
-    if(UseRhs)
-      actualRhsPtr = const_cast<RhsScalar*>(rhs.data());
-    else
+    if(!UseRhs)
     {
       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       int size = rhs.size();
       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       #endif
-      if((actualRhsPtr=static_rhs.data())==0)
-      {
-        freeRhsPtr = true;
-        actualRhsPtr = ei_aligned_stack_new(RhsScalar,rhs.size());
-      }
       Map<typename _ActualRhsType::PlainObject>(actualRhsPtr, rhs.size()) = rhs;
     }
       
@@ -259,11 +242,7 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true>
       );
     
     if(!EvalToDest)
-    {
       dest = MappedDest(actualDestPtr, dest.size());
-      if(freeDestPtr) ei_aligned_stack_delete(ResScalar, actualDestPtr, dest.size());
-    }
-    if(freeRhsPtr) ei_aligned_stack_delete(RhsScalar, actualRhsPtr, rhs.size());
   }
 };
 

Eigen/src/Core/products/SelfadjointProduct.h

@@ -81,27 +81,17 @@ struct selfadjoint_product_selector<MatrixType,OtherType,UpLo,true>
       UseOtherDirectly = _ActualOtherType::InnerStrideAtCompileTime==1
     };
     internal::gemv_static_vector_if<Scalar,OtherType::SizeAtCompileTime,OtherType::MaxSizeAtCompileTime,!UseOtherDirectly> static_other;
-    
-    bool freeOtherPtr = false;
-    Scalar* actualOtherPtr;
-    if(UseOtherDirectly)
-      actualOtherPtr = const_cast<Scalar*>(actualOther.data());
-    else
-    {
-      if((actualOtherPtr=static_other.data())==0)
-      {
-        freeOtherPtr = true;
-        actualOtherPtr = ei_aligned_stack_new(Scalar,other.size());
-      }
+
+    ei_declare_aligned_stack_constructed_variable(Scalar, actualOtherPtr, other.size(),
+      (UseOtherDirectly ? const_cast<Scalar*>(actualOther.data()) : static_other.data()));
+      
+    if(!UseOtherDirectly)
       Map<typename _ActualOtherType::PlainObject>(actualOtherPtr, actualOther.size()) = actualOther;
-    }
     
     selfadjoint_rank1_update<Scalar,Index,StorageOrder,UpLo,
                               OtherBlasTraits::NeedToConjugate  && NumTraits<Scalar>::IsComplex,
                             (!OtherBlasTraits::NeedToConjugate) && NumTraits<Scalar>::IsComplex>
           ::run(other.size(), mat.data(), mat.outerStride(), actualOtherPtr, actualAlpha);
-    
-    if((!UseOtherDirectly) && freeOtherPtr) ei_aligned_stack_delete(Scalar, actualOtherPtr, other.size());
   }
 };
 

Eigen/src/Core/products/TriangularMatrixMatrix.h

@@ -96,33 +96,38 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,true,
                                            LhsStorageOrder,ConjugateLhs,
                                            RhsStorageOrder,ConjugateRhs,ColMajor>
 {
+  
+  typedef gebp_traits<Scalar,Scalar> Traits;
+  enum {
+    SmallPanelWidth   = EIGEN_PLAIN_ENUM_MAX(Traits::mr,Traits::nr),
+    IsLower = (Mode&Lower) == Lower,
+    SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1
+  };
 
   static EIGEN_DONT_INLINE void run(
-    Index rows, Index cols, Index depth,
+    Index _rows, Index _cols, Index _depth,
     const Scalar* _lhs, Index lhsStride,
     const Scalar* _rhs, Index rhsStride,
     Scalar* res,        Index resStride,
     Scalar alpha)
   {
+    // strip zeros
+    Index diagSize  = (std::min)(_rows,_depth);
+    Index rows      = IsLower ? _rows : diagSize;
+    Index depth     = IsLower ? diagSize : _depth;
+    Index cols      = _cols;
+    
     const_blas_data_mapper<Scalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
     const_blas_data_mapper<Scalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride);
 
-    typedef gebp_traits<Scalar,Scalar> Traits;
-    enum {
-      SmallPanelWidth   = EIGEN_PLAIN_ENUM_MAX(Traits::mr,Traits::nr),
-      IsLower = (Mode&Lower) == Lower,
-      SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1
-    };
-
     Index kc = depth; // cache block size along the K direction
     Index mc = rows;  // cache block size along the M direction
     Index nc = cols;  // cache block size along the N direction
     computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
-
-    Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*cols;
-    Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar, sizeB);
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);    
     Scalar* blockB = allocatedBlockB + sizeW;
 
     Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,LhsStorageOrder> triangularBuffer;
@@ -140,7 +145,7 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,true,
         IsLower ? k2>0 : k2<depth;
         IsLower ? k2-=kc : k2+=kc)
     {
-      Index actual_kc = std::min(IsLower ? k2 : depth-k2, kc);
+      Index actual_kc = (std::min)(IsLower ? k2 : depth-k2, kc);
       Index actual_k2 = IsLower ? k2-actual_kc : k2;
 
       // align blocks with the end of the triangular part for trapezoidal lhs
@@ -153,10 +158,11 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,true,
       pack_rhs(blockB, &rhs(actual_k2,0), rhsStride, actual_kc, cols);
 
       // the selected lhs's panel has to be split in three different parts:
-      //  1 - the part which is above the diagonal block => skip it
+      //  1 - the part which is zero => skip it
       //  2 - the diagonal block => special kernel
-      //  3 - the panel below the diagonal block => GEPP
-      // the block diagonal, if any
+      //  3 - the dense panel below (lower case) or above (upper case) the diagonal block => GEPP
+
+      // the block diagonal, if any:
       if(IsLower || actual_k2<rows)
       {
         // for each small vertical panels of lhs
@@ -194,13 +200,13 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,true,
           }
         }
       }
-      // the part below the diagonal => GEPP
+      // the part below (lower case) or above (upper case) the diagonal => GEPP
       {
         Index start = IsLower ? k2 : 0;
-        Index end   = IsLower ? rows : std::min(actual_k2,rows);
+        Index end   = IsLower ? rows : (std::min)(actual_k2,rows);
         for(Index i2=start; i2<end; i2+=mc)
         {
-          const Index actual_mc = std::min(i2+mc,end)-i2;
+          const Index actual_mc = (std::min)(i2+mc,end)-i2;
           gemm_pack_lhs<Scalar, Index, Traits::mr,Traits::LhsProgress, LhsStorageOrder,false>()
             (blockA, &lhs(i2, actual_k2), lhsStride, actual_kc, actual_mc);
 
@@ -208,10 +214,6 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,true,
         }
       }
     }
-
-    ei_aligned_stack_delete(Scalar, blockA, kc*mc);
-    ei_aligned_stack_delete(Scalar, allocatedBlockB, sizeB);
-//     delete[] allocatedBlockB;
   }
 };
 
@@ -223,33 +225,38 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,false,
                                            LhsStorageOrder,ConjugateLhs,
                                            RhsStorageOrder,ConjugateRhs,ColMajor>
 {
+  typedef gebp_traits<Scalar,Scalar> Traits;
+  enum {
+    SmallPanelWidth   = EIGEN_PLAIN_ENUM_MAX(Traits::mr,Traits::nr),
+    IsLower = (Mode&Lower) == Lower,
+    SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1
+  };
 
   static EIGEN_DONT_INLINE void run(
-    Index rows, Index cols, Index depth,
+    Index _rows, Index _cols, Index _depth,
     const Scalar* _lhs, Index lhsStride,
     const Scalar* _rhs, Index rhsStride,
     Scalar* res,        Index resStride,
     Scalar alpha)
   {
+    // strip zeros
+    Index diagSize  = (std::min)(_cols,_depth);
+    Index rows      = _rows;
+    Index depth     = IsLower ? _depth : diagSize;
+    Index cols      = IsLower ? diagSize : _cols;
+    
     const_blas_data_mapper<Scalar, Index, LhsStorageOrder> lhs(_lhs,lhsStride);
     const_blas_data_mapper<Scalar, Index, RhsStorageOrder> rhs(_rhs,rhsStride);
 
-    typedef gebp_traits<Scalar,Scalar> Traits;
-    enum {
-      SmallPanelWidth   = EIGEN_PLAIN_ENUM_MAX(Traits::mr,Traits::nr),
-      IsLower = (Mode&Lower) == Lower,
-      SetDiag = (Mode&(ZeroDiag|UnitDiag)) ? 0 : 1
-    };
-
     Index kc = depth; // cache block size along the K direction
     Index mc = rows;  // cache block size along the M direction
     Index nc = cols;  // cache block size along the N direction
     computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
 
-    Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*cols;
-    Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar,sizeB);
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);
     Scalar* blockB = allocatedBlockB + sizeW;
 
     Matrix<Scalar,SmallPanelWidth,SmallPanelWidth,RhsStorageOrder> triangularBuffer;
@@ -268,7 +275,7 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,false,
         IsLower ? k2<depth  : k2>0;
         IsLower ? k2+=kc   : k2-=kc)
     {
-      Index actual_kc = std::min(IsLower ? depth-k2 : k2, kc);
+      Index actual_kc = (std::min)(IsLower ? depth-k2 : k2, kc);
       Index actual_k2 = IsLower ? k2 : k2-actual_kc;
 
       // align blocks with the end of the triangular part for trapezoidal rhs
@@ -279,7 +286,7 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,false,
       }
 
       // remaining size
-      Index rs = IsLower ? std::min(cols,actual_k2) : cols - k2;
+      Index rs = IsLower ? (std::min)(cols,actual_k2) : cols - k2;
       // size of the triangular part
       Index ts = (IsLower && actual_k2>=cols) ? 0 : actual_kc;
 
@@ -320,7 +327,7 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,false,
 
       for (Index i2=0; i2<rows; i2+=mc)
       {
-        const Index actual_mc = std::min(mc,rows-i2);
+        const Index actual_mc = (std::min)(mc,rows-i2);
         pack_lhs(blockA, &lhs(i2, actual_k2), lhsStride, actual_kc, actual_mc);
 
         // triangular kernel
@@ -347,9 +354,6 @@ struct product_triangular_matrix_matrix<Scalar,Index,Mode,false,
                     -1, -1, 0, 0, allocatedBlockB);
       }
     }
-
-    ei_aligned_stack_delete(Scalar, blockA, kc*mc);
-    ei_aligned_stack_delete(Scalar, allocatedBlockB, sizeB);
   }
 };
 

Eigen/src/Core/products/TriangularMatrixVector.h

@@ -41,9 +41,6 @@ struct product_triangular_matrix_vector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,C
   static EIGEN_DONT_INLINE  void run(Index rows, Index cols, const LhsScalar* _lhs, Index lhsStride,
                                      const RhsScalar* _rhs, Index rhsIncr, ResScalar* _res, Index resIncr, ResScalar alpha)
   {
-    EIGEN_UNUSED_VARIABLE(resIncr);
-    eigen_assert(resIncr==1);
-    
     static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
 
     typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,ColMajor>, 0, OuterStride<> > LhsMap;
@@ -59,7 +56,7 @@ struct product_triangular_matrix_vector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,C
 
     for (Index pi=0; pi<cols; pi+=PanelWidth)
     {
-      Index actualPanelWidth = std::min(PanelWidth, cols-pi);
+      Index actualPanelWidth = (std::min)(PanelWidth, cols-pi);
       for (Index k=0; k<actualPanelWidth; ++k)
       {
         Index i = pi + k;
@@ -95,9 +92,6 @@ struct product_triangular_matrix_vector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,C
   static void run(Index rows, Index cols, const LhsScalar* _lhs, Index lhsStride,
                   const RhsScalar* _rhs, Index rhsIncr, ResScalar* _res, Index resIncr, ResScalar alpha)
   {
-    eigen_assert(rhsIncr==1);
-    EIGEN_UNUSED_VARIABLE(rhsIncr);
-    
     static const Index PanelWidth = EIGEN_TUNE_TRIANGULAR_PANEL_WIDTH;
 
     typedef Map<const Matrix<LhsScalar,Dynamic,Dynamic,RowMajor>, 0, OuterStride<> > LhsMap;
@@ -113,7 +107,7 @@ struct product_triangular_matrix_vector<Index,Mode,LhsScalar,ConjLhs,RhsScalar,C
     
     for (Index pi=0; pi<cols; pi+=PanelWidth)
     {
-      Index actualPanelWidth = std::min(PanelWidth, cols-pi);
+      Index actualPanelWidth = (std::min)(PanelWidth, cols-pi);
       for (Index k=0; k<actualPanelWidth; ++k)
       {
         Index i = pi + k;
@@ -185,7 +179,7 @@ struct TriangularProduct<Mode,false,Lhs,true,Rhs,false>
   template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const
   {
     eigen_assert(dst.rows()==m_lhs.rows() && dst.cols()==m_rhs.cols());
-    
+
     typedef TriangularProduct<(Mode & UnitDiag) | ((Mode & Lower) ? Upper : Lower),true,Transpose<const Rhs>,false,Transpose<const Lhs>,true> TriangularProductTranspose;
     Transpose<Dest> dstT(dst);
     internal::trmv_selector<(int(internal::traits<Rhs>::Flags)&RowMajorBit) ? ColMajor : RowMajor>::run(
@@ -235,23 +229,15 @@ template<> struct trmv_selector<ColMajor>
     
     RhsScalar compatibleAlpha = get_factor<ResScalar,RhsScalar>::run(actualAlpha);
 
-    ResScalar* actualDestPtr;
-    bool freeDestPtr = false;
-    if (evalToDest)
-    {
-      actualDestPtr = dest.data();
-    }
-    else
+    ei_declare_aligned_stack_constructed_variable(ResScalar,actualDestPtr,dest.size(),
+                                                  evalToDest ? dest.data() : static_dest.data());
+
+    if(!evalToDest)
     {
       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       int size = dest.size();
       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       #endif
-      if((actualDestPtr = static_dest.data())==0)
-      {
-        freeDestPtr = true;
-        actualDestPtr = ei_aligned_stack_new(ResScalar,dest.size());
-      }
       if(!alphaIsCompatible)
       {
         MappedDest(actualDestPtr, dest.size()).setZero();
@@ -277,7 +263,6 @@ template<> struct trmv_selector<ColMajor>
         dest += actualAlpha * MappedDest(actualDestPtr, dest.size());
       else
         dest = MappedDest(actualDestPtr, dest.size());
-      if(freeDestPtr) ei_aligned_stack_delete(ResScalar, actualDestPtr, dest.size());
     }
   }
 };
@@ -310,23 +295,15 @@ template<> struct trmv_selector<RowMajor>
 
     gemv_static_vector_if<RhsScalar,_ActualRhsType::SizeAtCompileTime,_ActualRhsType::MaxSizeAtCompileTime,!DirectlyUseRhs> static_rhs;
 
-    RhsScalar* actualRhsPtr;
-    bool freeRhsPtr = false;
-    if (DirectlyUseRhs)
-    {
-      actualRhsPtr = const_cast<RhsScalar*>(actualRhs.data());
-    }
-    else
+    ei_declare_aligned_stack_constructed_variable(RhsScalar,actualRhsPtr,actualRhs.size(),
+        DirectlyUseRhs ? const_cast<RhsScalar*>(actualRhs.data()) : static_rhs.data());
+
+    if(!DirectlyUseRhs)
     {
       #ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       int size = actualRhs.size();
       EIGEN_DENSE_STORAGE_CTOR_PLUGIN
       #endif
-      if((actualRhsPtr = static_rhs.data())==0)
-      {
-        freeRhsPtr = true;
-        actualRhsPtr = ei_aligned_stack_new(RhsScalar, actualRhs.size());
-      }
       Map<typename _ActualRhsType::PlainObject>(actualRhsPtr, actualRhs.size()) = actualRhs;
     }
     
@@ -340,8 +317,6 @@ template<> struct trmv_selector<RowMajor>
             actualRhsPtr,1,
             dest.data(),dest.innerStride(),
             actualAlpha);
-
-    if((!DirectlyUseRhs) && freeRhsPtr) ei_aligned_stack_delete(RhsScalar, actualRhsPtr, prod.rhs().size());
   }
 };
 

Eigen/src/Core/products/TriangularSolverMatrix.h

@@ -70,10 +70,10 @@ struct triangular_solve_matrix<Scalar,Index,OnTheLeft,Mode,Conjugate,TriStorageO
     Index nc = cols;  // cache block size along the N direction
     computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
 
-    Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*cols;
-    Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar, sizeB);
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);
     Scalar* blockB = allocatedBlockB + sizeW;
 
     conj_if<Conjugate> conj;
@@ -85,7 +85,7 @@ struct triangular_solve_matrix<Scalar,Index,OnTheLeft,Mode,Conjugate,TriStorageO
         IsLower ? k2<size : k2>0;
         IsLower ? k2+=kc : k2-=kc)
     {
-      const Index actual_kc = std::min(IsLower ? size-k2 : k2, kc);
+      const Index actual_kc = (std::min)(IsLower ? size-k2 : k2, kc);
 
       // We have selected and packed a big horizontal panel R1 of rhs. Let B be the packed copy of this panel,
       // and R2 the remaining part of rhs. The corresponding vertical panel of lhs is split into
@@ -164,7 +164,7 @@ struct triangular_solve_matrix<Scalar,Index,OnTheLeft,Mode,Conjugate,TriStorageO
         Index end   = IsLower ? size : k2-kc;
         for(Index i2=start; i2<end; i2+=mc)
         {
-          const Index actual_mc = std::min(mc,end-i2);
+          const Index actual_mc = (std::min)(mc,end-i2);
           if (actual_mc>0)
           {
             pack_lhs(blockA, &tri(i2, IsLower ? k2 : k2-kc), triStride, actual_kc, actual_mc);
@@ -174,9 +174,6 @@ struct triangular_solve_matrix<Scalar,Index,OnTheLeft,Mode,Conjugate,TriStorageO
         }
       }
     }
-
-    ei_aligned_stack_delete(Scalar, blockA, kc*mc);
-    ei_aligned_stack_delete(Scalar, allocatedBlockB, sizeB);
   }
 };
 
@@ -209,10 +206,10 @@ struct triangular_solve_matrix<Scalar,Index,OnTheRight,Mode,Conjugate,TriStorage
     Index nc = rows;  // cache block size along the N direction
     computeProductBlockingSizes<Scalar,Scalar,4>(kc, mc, nc);
 
-    Scalar* blockA = ei_aligned_stack_new(Scalar, kc*mc);
     std::size_t sizeW = kc*Traits::WorkSpaceFactor;
     std::size_t sizeB = sizeW + kc*size;
-    Scalar* allocatedBlockB = ei_aligned_stack_new(Scalar, sizeB);
+    ei_declare_aligned_stack_constructed_variable(Scalar, blockA, kc*mc, 0);
+    ei_declare_aligned_stack_constructed_variable(Scalar, allocatedBlockB, sizeB, 0);
     Scalar* blockB = allocatedBlockB + sizeW;
 
     conj_if<Conjugate> conj;
@@ -225,7 +222,7 @@ struct triangular_solve_matrix<Scalar,Index,OnTheRight,Mode,Conjugate,TriStorage
         IsLower ? k2>0 : k2<size;
         IsLower ? k2-=kc : k2+=kc)
     {
-      const Index actual_kc = std::min(IsLower ? k2 : size-k2, kc);
+      const Index actual_kc = (std::min)(IsLower ? k2 : size-k2, kc);
       Index actual_k2 = IsLower ? k2-actual_kc : k2 ;
 
       Index startPanel = IsLower ? 0 : k2+actual_kc;
@@ -254,7 +251,7 @@ struct triangular_solve_matrix<Scalar,Index,OnTheRight,Mode,Conjugate,TriStorage
 
       for(Index i2=0; i2<rows; i2+=mc)
       {
-        const Index actual_mc = std::min(mc,rows-i2);
+        const Index actual_mc = (std::min)(mc,rows-i2);
 
         // triangular solver kernel
         {
@@ -314,9 +311,6 @@ struct triangular_solve_matrix<Scalar,Index,OnTheRight,Mode,Conjugate,TriStorage
                       -1, -1, 0, 0, allocatedBlockB);
       }
     }
-
-    ei_aligned_stack_delete(Scalar, blockA, kc*mc);
-    ei_aligned_stack_delete(Scalar, allocatedBlockB, sizeB);
   }
 };
 

Eigen/src/Core/products/TriangularSolverVector.h

@@ -60,7 +60,7 @@ struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Con
         IsLower ? pi<size : pi>0;
         IsLower ? pi+=PanelWidth : pi-=PanelWidth)
     {
-      Index actualPanelWidth = std::min(IsLower ? size - pi : pi, PanelWidth);
+      Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth);
 
       Index r = IsLower ? pi : size - pi; // remaining size
       if (r > 0)
@@ -114,7 +114,7 @@ struct triangular_solve_vector<LhsScalar, RhsScalar, Index, OnTheLeft, Mode, Con
         IsLower ? pi<size : pi>0;
         IsLower ? pi+=PanelWidth : pi-=PanelWidth)
     {
-      Index actualPanelWidth = std::min(IsLower ? size - pi : pi, PanelWidth);
+      Index actualPanelWidth = (std::min)(IsLower ? size - pi : pi, PanelWidth);
       Index startBlock = IsLower ? pi : pi-actualPanelWidth;
       Index endBlock = IsLower ? pi + actualPanelWidth : 0;
 

Eigen/src/Core/util/Constants.h

@@ -161,23 +161,72 @@ const unsigned int HereditaryBits = RowMajorBit
                                   | EvalBeforeNestingBit
                                   | EvalBeforeAssigningBit;
 
-// Possible values for the Mode parameter of triangularView()
-enum {
-  Lower=0x1, Upper=0x2, UnitDiag=0x4, ZeroDiag=0x8,
-  UnitLower=UnitDiag|Lower, UnitUpper=UnitDiag|Upper,
-  StrictlyLower=ZeroDiag|Lower, StrictlyUpper=ZeroDiag|Upper,
-  SelfAdjoint=0x10};
+/** \defgroup enums Enumerations
+  * \ingroup Core_Module
+  *
+  * Various enumerations used in %Eigen. Many of these are used as template parameters.
+  */
+
+/** \ingroup enums
+  * Enum containing possible values for the \p Mode parameter of 
+  * MatrixBase::selfadjointView() and MatrixBase::triangularView(). */
+enum {
+  /** View matrix as a lower triangular matrix. */
+  Lower=0x1,                      
+  /** View matrix as an upper triangular matrix. */
+  Upper=0x2,                      
+  /** %Matrix has ones on the diagonal; to be used in combination with #Lower or #Upper. */
+  UnitDiag=0x4, 
+  /** %Matrix has zeros on the diagonal; to be used in combination with #Lower or #Upper. */
+  ZeroDiag=0x8,
+  /** View matrix as a lower triangular matrix with ones on the diagonal. */
+  UnitLower=UnitDiag|Lower, 
+  /** View matrix as an upper triangular matrix with ones on the diagonal. */
+  UnitUpper=UnitDiag|Upper,
+  /** View matrix as a lower triangular matrix with zeros on the diagonal. */
+  StrictlyLower=ZeroDiag|Lower, 
+  /** View matrix as an upper triangular matrix with zeros on the diagonal. */
+  StrictlyUpper=ZeroDiag|Upper,
+  /** Used in BandMatrix and SelfAdjointView to indicate that the matrix is self-adjoint. */
+  SelfAdjoint=0x10
+};
+
+/** \ingroup enums
+  * Enum for indicating whether an object is aligned or not. */
+enum { 
+  /** Object is not correctly aligned for vectorization. */
+  Unaligned=0, 
+  /** Object is aligned for vectorization. */
+  Aligned=1 
+};
 
-enum { Unaligned=0, Aligned=1 };
 enum { ConditionalJumpCost = 5 };
 
+/** \ingroup enums
+ * Enum used by DenseBase::corner() in Eigen2 compatibility mode. */
 // FIXME after the corner() API change, this was not needed anymore, except by AlignedBox
 // TODO: find out what to do with that. Adapt the AlignedBox API ?
 enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight };
 
-enum DirectionType { Vertical, Horizontal, BothDirections };
+/** \ingroup enums
+  * Enum containing possible values for the \p Direction parameter of
+  * Reverse, PartialReduxExpr and VectorwiseOp. */
+enum DirectionType { 
+  /** For Reverse, all columns are reversed; 
+    * for PartialReduxExpr and VectorwiseOp, act on columns. */
+  Vertical, 
+  /** For Reverse, all rows are reversed; 
+    * for PartialReduxExpr and VectorwiseOp, act on rows. */
+  Horizontal, 
+  /** For Reverse, both rows and columns are reversed; 
+    * not used for PartialReduxExpr and VectorwiseOp. */
+  BothDirections 
+};
+
 enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct };
 
+/** \internal \ingroup enums
+  * Enum to specify how to traverse the entries of a matrix. */
 enum {
   /** \internal Default traversal, no vectorization, no index-based access */
   DefaultTraversal,
@@ -196,14 +245,25 @@ enum {
   InvalidTraversal
 };
 
+/** \internal \ingroup enums
+  * Enum to specify whether to unroll loops when traversing over the entries of a matrix. */
 enum {
+  /** \internal Do not unroll loops. */
   NoUnrolling,
+  /** \internal Unroll only the inner loop, but not the outer loop. */
   InnerUnrolling,
+  /** \internal Unroll both the inner and the outer loop. If there is only one loop, 
+    * because linear traversal is used, then unroll that loop. */
   CompleteUnrolling
 };
 
+/** \ingroup enums
+  * Enum containing possible values for the \p _Options template parameter of
+  * Matrix, Array and BandMatrix. */
 enum {
+  /** Storage order is column major (see \ref TopicStorageOrders). */
   ColMajor = 0,
+  /** Storage order is row major (see \ref TopicStorageOrders). */
   RowMajor = 0x1,  // it is only a coincidence that this is equal to RowMajorBit -- don't rely on that
   /** \internal Align the matrix itself if it is vectorizable fixed-size */
   AutoAlign = 0,
@@ -211,11 +271,13 @@ enum {
   DontAlign = 0x2
 };
 
-/** \brief Enum for specifying whether to apply or solve on the left or right. 
-  */
+/** \ingroup enums
+  * Enum for specifying whether to apply or solve on the left or right. */
 enum {
-  OnTheLeft = 1,  /**< \brief Apply transformation on the left. */
-  OnTheRight = 2  /**< \brief Apply transformation on the right. */
+  /** Apply transformation on the left. */
+  OnTheLeft = 1,  
+  /** Apply transformation on the right. */
+  OnTheRight = 2  
 };
 
 /* the following could as well be written:
@@ -239,53 +301,108 @@ namespace {
   EIGEN_UNUSED Default_t Default;
 }
 
+/** \internal \ingroup enums
+  * Used in AmbiVector. */
 enum {
   IsDense         = 0,
   IsSparse
 };
 
+/** \ingroup enums
+  * Used as template parameter in DenseCoeffBase and MapBase to indicate 
+  * which accessors should be provided. */
 enum AccessorLevels {
-  ReadOnlyAccessors, WriteAccessors, DirectAccessors, DirectWriteAccessors
+  /** Read-only access via a member function. */
+  ReadOnlyAccessors, 
+  /** Read/write access via member functions. */
+  WriteAccessors, 
+  /** Direct read-only access to the coefficients. */
+  DirectAccessors, 
+  /** Direct read/write access to the coefficients. */
+  DirectWriteAccessors
 };
 
+/** \ingroup enums
+  * Enum with options to give to various decompositions. */
 enum DecompositionOptions {
-  Pivoting            = 0x01, // LDLT,
-  NoPivoting          = 0x02, // LDLT,
-  ComputeFullU        = 0x04, // SVD,
-  ComputeThinU        = 0x08, // SVD,
-  ComputeFullV        = 0x10, // SVD,
-  ComputeThinV        = 0x20, // SVD,
-  EigenvaluesOnly     = 0x40, // all eigen solvers
-  ComputeEigenvectors = 0x80, // all eigen solvers
+  /** \internal Not used (meant for LDLT?). */
+  Pivoting            = 0x01, 
+  /** \internal Not used (meant for LDLT?). */
+  NoPivoting          = 0x02, 
+  /** Used in JacobiSVD to indicate that the square matrix U is to be computed. */
+  ComputeFullU        = 0x04,
+  /** Used in JacobiSVD to indicate that the thin matrix U is to be computed. */
+  ComputeThinU        = 0x08,
+  /** Used in JacobiSVD to indicate that the square matrix V is to be computed. */
+  ComputeFullV        = 0x10,
+  /** Used in JacobiSVD to indicate that the thin matrix V is to be computed. */
+  ComputeThinV        = 0x20,
+  /** Used in SelfAdjointEigenSolver and GeneralizedSelfAdjointEigenSolver to specify
+    * that only the eigenvalues are to be computed and not the eigenvectors. */
+  EigenvaluesOnly     = 0x40,
+  /** Used in SelfAdjointEigenSolver and GeneralizedSelfAdjointEigenSolver to specify
+    * that both the eigenvalues and the eigenvectors are to be computed. */
+  ComputeEigenvectors = 0x80,
+  /** \internal */
   EigVecMask = EigenvaluesOnly | ComputeEigenvectors,
+  /** Used in GeneralizedSelfAdjointEigenSolver to indicate that it should
+    * solve the generalized eigenproblem \f$ Ax = \lambda B x \f$. */
   Ax_lBx              = 0x100,
+  /** Used in GeneralizedSelfAdjointEigenSolver to indicate that it should
+    * solve the generalized eigenproblem \f$ ABx = \lambda x \f$. */
   ABx_lx              = 0x200,
+  /** Used in GeneralizedSelfAdjointEigenSolver to indicate that it should
+    * solve the generalized eigenproblem \f$ BAx = \lambda x \f$. */
   BAx_lx              = 0x400,
+  /** \internal */
   GenEigMask = Ax_lBx | ABx_lx | BAx_lx
 };
 
+/** \ingroup enums
+  * Possible values for the \p QRPreconditioner template parameter of JacobiSVD. */
 enum QRPreconditioners {
+  /** Do not specify what is to be done if the SVD of a non-square matrix is asked for. */
   NoQRPreconditioner,
+  /** Use a QR decomposition without pivoting as the first step. */
   HouseholderQRPreconditioner,
+  /** Use a QR decomposition with column pivoting as the first step. */
   ColPivHouseholderQRPreconditioner,
+  /** Use a QR decomposition with full pivoting as the first step. */
   FullPivHouseholderQRPreconditioner
 };
 
-/** \brief Enum for reporting the status of a computation.
-  */
+#ifdef Success
+#error The preprocessor symbol 'Success' is defined, possibly by the X11 header file X.h
+#endif
+
+/** \ingroups enums
+  * Enum for reporting the status of a computation. */
 enum ComputationInfo {
-  Success = 0,        /**< \brief Computation was successful. */
-  NumericalIssue = 1, /**< \brief The provided data did not satisfy the prerequisites. */
-  NoConvergence = 2   /**< \brief Iterative procedure did not converge. */
+  /** Computation was successful. */
+  Success = 0,        
+  /** The provided data did not satisfy the prerequisites. */
+  NumericalIssue = 1, 
+  /** Iterative procedure did not converge. */
+  NoConvergence = 2
 };
 
+/** \ingroup enums
+  * Enum used to specify how a particular transformation is stored in a matrix.
+  * \sa Transform, Hyperplane::transform(). */
 enum TransformTraits {
+  /** Transformation is an isometry. */
   Isometry      = 0x1,
+  /** Transformation is an affine transformation stored as a (Dim+1)^2 matrix whose last row is 
+    * assumed to be [0 ... 0 1]. */
   Affine        = 0x2,
+  /** Transformation is an affine transformation stored as a (Dim) x (Dim+1) matrix. */
   AffineCompact = 0x10 | Affine,
+  /** Transformation is a general projective transformation stored as a (Dim+1)^2 matrix. */
   Projective    = 0x20
 };
 
+/** \internal \ingroup enums
+  * Enum used to choose between implementation depending on the computer architecture. */
 namespace Architecture
 {
   enum Type {
@@ -302,8 +419,12 @@ namespace Architecture
   };
 }
 
+/** \internal \ingroup enums
+  * Enum used as template parameter in GeneralProduct. */
 enum { CoeffBasedProductMode, LazyCoeffBasedProductMode, OuterProduct, InnerProduct, GemvProduct, GemmProduct };
 
+/** \internal \ingroup enums
+  * Enum used in experimental parallel implementation. */
 enum Action {GetAction, SetAction};
 
 /** The type used to identify a dense storage. */

Eigen/src/Core/util/Macros.h

@@ -1,3 +1,4 @@
+
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
@@ -26,9 +27,9 @@
 #ifndef EIGEN_MACROS_H
 #define EIGEN_MACROS_H
 
-#define EIGEN_WORLD_VERSION 2
-#define EIGEN_MAJOR_VERSION 95
-#define EIGEN_MINOR_VERSION 0
+#define EIGEN_WORLD_VERSION 3
+#define EIGEN_MAJOR_VERSION 0
+#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
@@ -399,7 +403,7 @@
 #define EIGEN_MAKE_CWISE_BINARY_OP(METHOD,FUNCTOR) \
   template<typename OtherDerived> \
   EIGEN_STRONG_INLINE const CwiseBinaryOp<FUNCTOR<Scalar>, const Derived, const OtherDerived> \
-  METHOD(const EIGEN_CURRENT_STORAGE_BASE_CLASS<OtherDerived> &other) const \
+  (METHOD)(const EIGEN_CURRENT_STORAGE_BASE_CLASS<OtherDerived> &other) const \
   { \
     return CwiseBinaryOp<FUNCTOR<Scalar>, const Derived, const OtherDerived>(derived(), other.derived()); \
   }

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                           ***
 *****************************************************************************/
@@ -156,7 +166,7 @@ inline void* generic_aligned_realloc(void* ptr, size_t size, size_t old_size)
 
   if (ptr != 0)
   {
-    std::memcpy(newptr, ptr, std::min(size,old_size));
+    std::memcpy(newptr, ptr, (std::min)(size,old_size));
     aligned_free(ptr);
   }
 
@@ -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));
@@ -468,36 +489,89 @@ inline static Index first_aligned(const Scalar* array, Index size)
 *** Implementation of runtime stack allocation (falling back to malloc)    ***
 *****************************************************************************/
 
-/** \internal
-  * Allocates an aligned buffer of SIZE bytes on the stack if SIZE is smaller than
-  * EIGEN_STACK_ALLOCATION_LIMIT, and if stack allocation is supported by the platform
-  * (currently, this is Linux only). Otherwise the memory is allocated on the heap.
-  * Data allocated with ei_aligned_stack_alloc \b must be freed by calling
-  * ei_aligned_stack_free(PTR,SIZE).
-  * \code
-  * float * data = ei_aligned_stack_alloc(float,array.size());
-  * // ...
-  * ei_aligned_stack_free(data,float,array.size());
-  * \endcode
-  */
-#if (defined __linux__)
-  #define ei_aligned_stack_alloc(SIZE) (SIZE<=EIGEN_STACK_ALLOCATION_LIMIT) \
-                                    ? alloca(SIZE) \
-                                    : Eigen::internal::aligned_malloc(SIZE)
-  #define ei_aligned_stack_free(PTR,SIZE) if(SIZE>EIGEN_STACK_ALLOCATION_LIMIT) Eigen::internal::aligned_free(PTR)
-#elif defined(_MSC_VER)
-  #define ei_aligned_stack_alloc(SIZE) (SIZE<=EIGEN_STACK_ALLOCATION_LIMIT) \
-                                    ? _alloca(SIZE) \
-                                    : Eigen::internal::aligned_malloc(SIZE)
-  #define ei_aligned_stack_free(PTR,SIZE) if(SIZE>EIGEN_STACK_ALLOCATION_LIMIT) Eigen::internal::aligned_free(PTR)
-#else
-  #define ei_aligned_stack_alloc(SIZE) Eigen::internal::aligned_malloc(SIZE)
-  #define ei_aligned_stack_free(PTR,SIZE) Eigen::internal::aligned_free(PTR)
+// you can overwrite Eigen's default behavior regarding alloca by defining EIGEN_ALLOCA
+// to the appropriate stack allocation function
+#ifndef EIGEN_ALLOCA
+  #if (defined __linux__)
+    #define EIGEN_ALLOCA alloca
+  #elif defined(_MSC_VER)
+    #define EIGEN_ALLOCA _alloca
+  #endif
 #endif
 
-#define ei_aligned_stack_new(TYPE,SIZE) Eigen::internal::construct_elements_of_array(reinterpret_cast<TYPE*>(ei_aligned_stack_alloc(sizeof(TYPE)*SIZE)), SIZE)
-#define ei_aligned_stack_delete(TYPE,PTR,SIZE) do {Eigen::internal::destruct_elements_of_array<TYPE>(PTR, SIZE); \
-                                                   ei_aligned_stack_free(PTR,sizeof(TYPE)*SIZE);} while(0)
+namespace internal {
+
+// This helper class construct the allocated memory, and takes care of destructing and freeing the handled data
+// at destruction time. In practice this helper class is mainly useful to avoid memory leak in case of exceptions.
+template<typename T> class aligned_stack_memory_handler
+{
+  public:
+    /* Creates a stack_memory_handler responsible for the buffer \a ptr of size \a size.
+     * Note that \a ptr can be 0 regardless of the other parameters.
+     * This constructor takes care of constructing/initializing the elements of the buffer if required by the scalar type T (see NumTraits<T>::RequireInitialization).
+     * In this case, the buffer elements will also be destructed when this handler will be destructed.
+     * Finally, if \a dealloc is true, then the pointer \a ptr is freed.
+     **/
+    aligned_stack_memory_handler(T* ptr, size_t size, bool dealloc)
+      : m_ptr(ptr), m_size(size), m_deallocate(dealloc)
+    {
+      if(NumTraits<T>::RequireInitialization && m_ptr)
+        Eigen::internal::construct_elements_of_array(m_ptr, size);
+    }
+    ~aligned_stack_memory_handler()
+    {
+      if(NumTraits<T>::RequireInitialization && m_ptr)
+        Eigen::internal::destruct_elements_of_array<T>(m_ptr, m_size);
+      if(m_deallocate)
+        Eigen::internal::aligned_free(m_ptr);
+    }
+  protected:
+    T* m_ptr;
+    size_t m_size;
+    bool m_deallocate;
+};
+
+}
+
+/** \internal
+  * Declares, allocates and construct an aligned buffer named NAME of SIZE elements of type TYPE on the stack
+  * if SIZE is smaller than EIGEN_STACK_ALLOCATION_LIMIT, and if stack allocation is supported by the platform
+  * (currently, this is Linux and Visual Studio only). Otherwise the memory is allocated on the heap.
+  * The allocated buffer is automatically deleted when exiting the scope of this declaration.
+  * If BUFFER is non nul, then the declared variable is simply an alias for BUFFER, and no allocation/deletion occurs.
+  * Here is an example:
+  * \code
+  * {
+  *   ei_declare_aligned_stack_constructed_variable(float,data,size,0);
+  *   // use data[0] to data[size-1]
+  * }
+  * \endcode
+  * The underlying stack allocation function can controlled with the EIGEN_ALLOCA preprocessor token.
+  */
+#ifdef EIGEN_ALLOCA
+
+  #ifdef __arm__
+    #define EIGEN_ALIGNED_ALLOCA(SIZE) reinterpret_cast<void*>((reinterpret_cast<size_t>(EIGEN_ALLOCA(SIZE+16)) & ~(size_t(15))) + 16)
+  #else
+    #define EIGEN_ALIGNED_ALLOCA EIGEN_ALLOCA
+  #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) \
+                    : 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,sizeof(TYPE)*SIZE>EIGEN_STACK_ALLOCATION_LIMIT)
+
+#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)
+    
+#endif
 
 
 /*****************************************************************************
@@ -612,12 +686,13 @@ public:
 
     size_type max_size() const throw()
     {
-        return std::numeric_limits<size_type>::max();
+        return (std::numeric_limits<size_type>::max)();
     }
 
-    pointer allocate( size_type num, const_pointer* hint = 0 )
+    pointer allocate( size_type num, const void* hint = 0 )
     {
-        static_cast<void>( hint ); // suppress unused variable warning
+        EIGEN_UNUSED_VARIABLE(hint);
+        internal::check_size_for_overflow<T>(num);
         return static_cast<pointer>( internal::aligned_malloc( num * sizeof(T) ) );
     }
 
@@ -852,7 +927,7 @@ inline int queryTopLevelCacheSize()
 {
   int l1, l2(-1), l3(-1);
   queryCacheSizes(l1,l2,l3);
-  return std::max(l2,l3);
+  return (std::max)(l2,l3);
 }
 
 } // end namespace internal

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

Eigen/src/Eigen2Support/Cwise.h

@@ -82,13 +82,17 @@ template<typename ExpressionType> class Cwise
     const EIGEN_CWISE_BINOP_RETURN_TYPE(internal::scalar_quotient_op)
     operator/(const MatrixBase<OtherDerived> &other) const;
 
+    /** \deprecated ArrayBase::min() */
     template<typename OtherDerived>
     const EIGEN_CWISE_BINOP_RETURN_TYPE(internal::scalar_min_op)
-    min(const MatrixBase<OtherDerived> &other) const;
+    (min)(const MatrixBase<OtherDerived> &other) const
+    { return EIGEN_CWISE_BINOP_RETURN_TYPE(internal::scalar_min_op)(_expression(), other.derived()); }
 
+    /** \deprecated ArrayBase::max() */
     template<typename OtherDerived>
     const EIGEN_CWISE_BINOP_RETURN_TYPE(internal::scalar_max_op)
-    max(const MatrixBase<OtherDerived> &other) const;
+    (max)(const MatrixBase<OtherDerived> &other) const
+    { return EIGEN_CWISE_BINOP_RETURN_TYPE(internal::scalar_max_op)(_expression(), other.derived()); }
 
     const EIGEN_CWISE_UNOP_RETURN_TYPE(internal::scalar_abs_op)      abs() const;
     const EIGEN_CWISE_UNOP_RETURN_TYPE(internal::scalar_abs2_op)     abs2() const;

Eigen/src/Eigen2Support/CwiseOperators.h

@@ -96,24 +96,6 @@ inline ExpressionType& Cwise<ExpressionType>::operator/=(const MatrixBase<OtherD
   return m_matrix.const_cast_derived() = *this / other;
 }
 
-/** \deprecated ArrayBase::min() */
-template<typename ExpressionType>
-template<typename OtherDerived>
-EIGEN_STRONG_INLINE const EIGEN_CWISE_BINOP_RETURN_TYPE(internal::scalar_min_op)
-Cwise<ExpressionType>::min(const MatrixBase<OtherDerived> &other) const
-{
-  return EIGEN_CWISE_BINOP_RETURN_TYPE(internal::scalar_min_op)(_expression(), other.derived());
-}
-
-/** \deprecated ArrayBase::max() */
-template<typename ExpressionType>
-template<typename OtherDerived>
-EIGEN_STRONG_INLINE const EIGEN_CWISE_BINOP_RETURN_TYPE(internal::scalar_max_op)
-Cwise<ExpressionType>::max(const MatrixBase<OtherDerived> &other) const
-{
-  return EIGEN_CWISE_BINOP_RETURN_TYPE(internal::scalar_max_op)(_expression(), other.derived());
-}
-
 /***************************************************************************
 * The following functions were defined in Array
 ***************************************************************************/

Eigen/src/Eigen2Support/Geometry/AlignedBox.h

@@ -63,7 +63,7 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim==
   ~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,18 +71,18 @@ 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 */
-  inline const VectorType& min() const { return m_min; }
+  inline const VectorType& (min)() const { return m_min; }
   /** \returns a non const reference to the minimal corner */
-  inline VectorType& min() { return m_min; }
+  inline VectorType& (min)() { return m_min; }
   /** \returns the maximal corner */
-  inline const VectorType& max() const { return m_max; }
+  inline const VectorType& (max)() const { return m_max; }
   /** \returns a non const reference to the maximal corner */
-  inline VectorType& max() { return m_max; }
+  inline VectorType& (max)() { return m_max; }
 
   /** \returns true if the point \a p is inside the box \c *this. */
   inline bool contains(const VectorType& p) const
@@ -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/Eigen2Support/SVD.h

@@ -64,9 +64,9 @@ template<typename MatrixType> class SVD
     SVD() {} // a user who relied on compiler-generated default compiler reported problems with MSVC in 2.0.7
     
     SVD(const MatrixType& matrix)
-      : m_matU(matrix.rows(), std::min(matrix.rows(), matrix.cols())),
+      : m_matU(matrix.rows(), (std::min)(matrix.rows(), matrix.cols())),
         m_matV(matrix.cols(),matrix.cols()),
-        m_sigma(std::min(matrix.rows(),matrix.cols()))
+        m_sigma((std::min)(matrix.rows(),matrix.cols()))
     {
       compute(matrix);
     }
@@ -108,13 +108,13 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
 {
   const int m = matrix.rows();
   const int n = matrix.cols();
-  const int nu = std::min(m,n);
+  const int nu = (std::min)(m,n);
   ei_assert(m>=n && "In Eigen 2.0, SVD only works for MxN matrices with M>=N. Sorry!");
   ei_assert(m>1 && "In Eigen 2.0, SVD doesn't work on 1x1 matrices");
 
   m_matU.resize(m, nu);
   m_matU.setZero();
-  m_sigma.resize(std::min(m,n));
+  m_sigma.resize((std::min)(m,n));
   m_matV.resize(n,n);
 
   RowVector e(n);
@@ -126,9 +126,9 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
 
   // Reduce A to bidiagonal form, storing the diagonal elements
   // in s and the super-diagonal elements in e.
-  int nct = std::min(m-1,n);
-  int nrt = std::max(0,std::min(n-2,m));
-  for (k = 0; k < std::max(nct,nrt); ++k)
+  int nct = (std::min)(m-1,n);
+  int nrt = (std::max)(0,(std::min)(n-2,m));
+  for (k = 0; k < (std::max)(nct,nrt); ++k)
   {
     if (k < nct)
     {
@@ -193,7 +193,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
 
 
   // Set up the final bidiagonal matrix or order p.
-  int p = std::min(n,m+1);
+  int p = (std::min)(n,m+1);
   if (nct < n)
     m_sigma[nct] = matA(nct,nct);
   if (m < p)
@@ -380,7 +380,7 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
       case 3:
       {
         // Calculate the shift.
-        Scalar scale = std::max(std::max(std::max(std::max(
+        Scalar scale = (std::max)((std::max)((std::max)((std::max)(
                         ei_abs(m_sigma[p-1]),ei_abs(m_sigma[p-2])),ei_abs(e[p-2])),
                         ei_abs(m_sigma[k])),ei_abs(e[k]));
         Scalar sp = m_sigma[p-1]/scale;

Eigen/src/Eigenvalues/ComplexSchur.h

@@ -423,7 +423,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
     JacobiRotation<ComplexScalar> rot;
     rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
     m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
-    m_matT.topRows(std::min(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
+    m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
     if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
 
     for(Index i=il+1 ; i<iu ; i++)
@@ -431,7 +431,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
       rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
       m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
       m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
-      m_matT.topRows(std::min(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
+      m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
       if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
     }
   }

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,10 +339,11 @@ 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>();
+      matV.col(j).normalize();
     }
     else
     {
@@ -434,7 +435,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
   Scalar norm = 0.0;
   for (Index j = 0; j < size; ++j)
   {
-    norm += m_matT.row(j).segment(std::max(j-1,Index(0)), size-std::max(j-1,Index(0))).cwiseAbs().sum();
+    norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
   }
   
   // Backsubstitute to find vectors of upper triangular form
@@ -449,7 +450,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
     Scalar q = m_eivalues.coeff(n).imag();
 
     // Scalar vector
-    if (q == 0)
+    if (q == Scalar(0))
     {
       Scalar lastr=0, lastw=0;
       Index l = n;
@@ -490,12 +491,12 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
 
           // Overflow control
           Scalar t = internal::abs(m_matT.coeff(i,n));
-          if ((eps * t) * t > 1)
+          if ((eps * t) * t > Scalar(1))
             m_matT.col(n).tail(size-i) /= t;
         }
       }
     }
-    else if (q < 0) // Complex vector
+    else if (q < Scalar(0) && n > 0) // Complex vector
     {
       Scalar lastra=0, lastsa=0, lastw=0;
       Index l = n-1;
@@ -529,7 +530,7 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
         else
         {
           l = i;
-          if (m_eivalues.coeff(i).imag() == 0)
+          if (m_eivalues.coeff(i).imag() == RealScalar(0))
           {
             std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
             m_matT.coeffRef(i,n-1) = internal::real(cc);
@@ -562,12 +563,20 @@ void EigenSolver<MatrixType>::doComputeEigenvectors()
           }
 
           // Overflow control
-          Scalar t = std::max(internal::abs(m_matT.coeff(i,n-1)),internal::abs(m_matT.coeff(i,n)));
-          if ((eps * t) * t > 1)
+          using std::max;
+          Scalar t = (max)(internal::abs(m_matT.coeff(i,n-1)),internal::abs(m_matT.coeff(i,n)));
+          if ((eps * t) * t > Scalar(1))
             m_matT.block(i, n-1, size-i, 2) /= t;
 
         }
       }
+      
+      // We handled a pair of complex conjugate eigenvalues, so need to skip them both
+      n--;
+    }
+    else
+    {
+      eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen
     }
   }
 

Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h

@@ -98,8 +98,8 @@ class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT
       *                   Only the lower triangular part of the matrix is referenced.
       * \param[in]  matB  Positive-definite matrix in matrix pencil.
       *                   Only the lower triangular part of the matrix is referenced.
-      * \param[in]  options A or-ed set of flags {ComputeEigenvectors,EigenvaluesOnly} | {Ax_lBx,ABx_lx,BAx_lx}.
-      *                     Default is ComputeEigenvectors|Ax_lBx.
+      * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
+      *                     Default is #ComputeEigenvectors|#Ax_lBx.
       *
       * This constructor calls compute(const MatrixType&, const MatrixType&, int)
       * to compute the eigenvalues and (if requested) the eigenvectors of the
@@ -131,8 +131,8 @@ class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixT
       *                   Only the lower triangular part of the matrix is referenced.
       * \param[in]  matB  Positive-definite matrix in matrix pencil.
       *                   Only the lower triangular part of the matrix is referenced.
-      * \param[in]  options A or-ed set of flags {ComputeEigenvectors,EigenvaluesOnly} | {Ax_lBx,ABx_lx,BAx_lx}.
-      *                     Default is ComputeEigenvectors|Ax_lBx.
+      * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
+      *                     Default is #ComputeEigenvectors|#Ax_lBx.
       *
       * \returns    Reference to \c *this
       *

Eigen/src/Eigenvalues/RealSchur.h

@@ -290,7 +290,7 @@ inline typename MatrixType::Scalar RealSchur<MatrixType>::computeNormOfT()
   //               + m_matT.bottomLeftCorner(size-1,size-1).diagonal().cwiseAbs().sum();
   Scalar norm = 0.0;
   for (Index j = 0; j < size; ++j)
-    norm += m_matT.row(j).segment(std::max(j-1,Index(0)), size-std::max(j-1,Index(0))).cwiseAbs().sum();
+    norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
   return norm;
 }
 
@@ -324,11 +324,11 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu, bool computeU, Scal
   m_matT.coeffRef(iu,iu) += exshift;
   m_matT.coeffRef(iu-1,iu-1) += exshift;
 
-  if (q >= 0) // Two real eigenvalues
+  if (q >= Scalar(0)) // Two real eigenvalues
   {
     Scalar z = internal::sqrt(internal::abs(q));
     JacobiRotation<Scalar> rot;
-    if (p >= 0)
+    if (p >= Scalar(0))
       rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
     else
       rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
@@ -369,7 +369,7 @@ inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& ex
   {
     Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
     s = s * s + shiftInfo.coeff(2);
-    if (s > 0)
+    if (s > Scalar(0))
     {
       s = internal::sqrt(s);
       if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
@@ -442,7 +442,7 @@ inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Inde
 
       // These Householder transformations form the O(n^3) part of the algorithm
       m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
-      m_matT.block(0, k, std::min(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
+      m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
       if (computeU)
         m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
     }

Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h

@@ -147,11 +147,11 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       *
       * \param[in]  matrix  Selfadjoint matrix whose eigendecomposition is to
       *    be computed. Only the lower triangular part of the matrix is referenced.
-      * \param[in]  options Can be ComputeEigenvectors (default) or EigenvaluesOnly.
+      * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
       *
       * This constructor calls compute(const MatrixType&, int) to compute the
       * eigenvalues of the matrix \p matrix. The eigenvectors are computed if
-      * \p options equals ComputeEigenvectors.
+      * \p options equals #ComputeEigenvectors.
       *
       * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.cpp
       * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.out
@@ -171,11 +171,11 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       *
       * \param[in]  matrix  Selfadjoint matrix whose eigendecomposition is to
       *    be computed. Only the lower triangular part of the matrix is referenced.
-      * \param[in]  options Can be ComputeEigenvectors (default) or EigenvaluesOnly.
+      * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
       * \returns    Reference to \c *this
       *
       * This function computes the eigenvalues of \p matrix.  The eigenvalues()
-      * function can be used to retrieve them.  If \p options equals ComputeEigenvectors,
+      * function can be used to retrieve them.  If \p options equals #ComputeEigenvectors,
       * then the eigenvectors are also computed and can be retrieved by
       * calling eigenvectors().
       *
@@ -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;
 
@@ -387,7 +388,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
   {
     m_eivalues.coeffRef(0,0) = internal::real(matrix.coeff(0,0));
     if(computeEigenvectors)
-      m_eivec.setOnes();
+      m_eivec.setOnes(n,n);
     m_info = Success;
     m_isInitialized = true;
     m_eigenvectorsOk = computeEigenvectors;
@@ -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;

Eigen/src/Geometry/AlignedBox.h

@@ -111,13 +111,13 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim)
   }
 
   /** \returns the minimal corner */
-  inline const VectorType& min() const { return m_min; }
+  inline const VectorType& (min)() const { return m_min; }
   /** \returns a non const reference to the minimal corner */
-  inline VectorType& min() { return m_min; }
+  inline VectorType& (min)() { return m_min; }
   /** \returns the maximal corner */
-  inline const VectorType& max() const { return m_max; }
+  inline const VectorType& (max)() const { return m_max; }
   /** \returns a non const reference to the maximal corner */
-  inline VectorType& max() { return m_max; }
+  inline VectorType& (max)() { return m_max; }
 
   /** \returns the center of the box */
   inline const CwiseUnaryOp<internal::scalar_quotient1_op<Scalar>,
@@ -196,7 +196,7 @@ 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.array()<=b.min().array()).all() && (b.max().array()<=m_max.array()).all(); }
+  { return (m_min.array()<=(b.min)().array()).all() && ((b.max)().array()<=m_max.array()).all(); }
 
   /** Extends \c *this such that it contains the point \a p and returns a reference to \c *this. */
   template<typename Derived>
@@ -287,8 +287,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/Geometry/AngleAxis.h

@@ -171,6 +171,9 @@ template<typename Scalar>
 template<typename QuatDerived>
 AngleAxis<Scalar>& AngleAxis<Scalar>::operator=(const QuaternionBase<QuatDerived>& q)
 {
+  using std::acos;
+  using std::min;
+  using std::max;
   Scalar n2 = q.vec().squaredNorm();
   if (n2 < NumTraits<Scalar>::dummy_precision()*NumTraits<Scalar>::dummy_precision())
   {
@@ -179,7 +182,7 @@ AngleAxis<Scalar>& AngleAxis<Scalar>::operator=(const QuaternionBase<QuatDerived
   }
   else
   {
-    m_angle = Scalar(2)*std::acos(std::min(std::max(Scalar(-1),q.w()),Scalar(1)));
+    m_angle = Scalar(2)*acos((min)((max)(Scalar(-1),q.w()),Scalar(1)));
     m_axis = q.vec() / internal::sqrt(n2);
   }
   return *this;

Eigen/src/Geometry/Hyperplane.h

@@ -189,7 +189,7 @@ public:
     *
     * \note If \a other is approximately parallel to *this, this method will return any point on *this.
     */
-  VectorType intersection(const Hyperplane& other)
+  VectorType intersection(const Hyperplane& other) const
   {
     EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(VectorType, 2)
     Scalar det = coeffs().coeff(0) * other.coeffs().coeff(1) - coeffs().coeff(1) * other.coeffs().coeff(0);
@@ -213,8 +213,8 @@ public:
   /** Applies the transformation matrix \a mat to \c *this and returns a reference to \c *this.
     *
     * \param mat the Dim x Dim transformation matrix
-    * \param traits specifies whether the matrix \a mat represents an Isometry
-    *               or a more generic Affine transformation. The default is Affine.
+    * \param traits specifies whether the matrix \a mat represents an #Isometry
+    *               or a more generic #Affine transformation. The default is #Affine.
     */
   template<typename XprType>
   inline Hyperplane& transform(const MatrixBase<XprType>& mat, TransformTraits traits = Affine)
@@ -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;
   }
@@ -233,8 +233,8 @@ public:
   /** Applies the transformation \a t to \c *this and returns a reference to \c *this.
     *
     * \param t the transformation of dimension Dim
-    * \param traits specifies whether the transformation \a t represents an Isometry
-    *               or a more generic Affine transformation. The default is Affine.
+    * \param traits specifies whether the transformation \a t represents an #Isometry
+    *               or a more generic #Affine transformation. The default is #Affine.
     *               Other kind of transformations are not supported.
     */
   template<int TrOptions>

Eigen/src/Geometry/ParametrizedLine.h

@@ -34,7 +34,7 @@
   *
   * A parametrized line is defined by an origin point \f$ \mathbf{o} \f$ and a unit
   * direction vector \f$ \mathbf{d} \f$ such that the line corresponds to
-  * the set \f$ l(t) = \mathbf{o} + t \mathbf{d} \f$, \f$ l \in \mathbf{R} \f$.
+  * the set \f$ l(t) = \mathbf{o} + t \mathbf{d} \f$, \f$ t \in \mathbf{R} \f$.
   *
   * \param _Scalar the scalar type, i.e., the type of the coefficients
   * \param _AmbientDim the dimension of the ambient space, can be a compile time value or Dynamic.
@@ -107,7 +107,7 @@ public:
   { return origin() + direction().dot(p-origin()) * direction(); }
 
   template <int OtherOptions>
-  Scalar intersection(const Hyperplane<_Scalar, _AmbientDim, OtherOptions>& hyperplane);
+  Scalar intersection(const Hyperplane<_Scalar, _AmbientDim, OtherOptions>& hyperplane) const;
 
   /** \returns \c *this with scalar type casted to \a NewScalarType
     *
@@ -159,7 +159,7 @@ inline ParametrizedLine<_Scalar, _AmbientDim,_Options>::ParametrizedLine(const H
   */
 template <typename _Scalar, int _AmbientDim, int _Options>
 template <int OtherOptions>
-inline _Scalar ParametrizedLine<_Scalar, _AmbientDim,_Options>::intersection(const Hyperplane<_Scalar, _AmbientDim, OtherOptions>& hyperplane)
+inline _Scalar ParametrizedLine<_Scalar, _AmbientDim,_Options>::intersection(const Hyperplane<_Scalar, _AmbientDim, OtherOptions>& hyperplane) const
 {
   return -(hyperplane.offset()+hyperplane.normal().dot(origin()))
           / hyperplane.normal().dot(direction());

Eigen/src/Geometry/Quaternion.h

@@ -49,6 +49,9 @@ public:
   typedef typename internal::traits<Derived>::Scalar Scalar;
   typedef typename NumTraits<Scalar>::Real RealScalar;
   typedef typename internal::traits<Derived>::Coefficients Coefficients;
+  enum {
+    Flags = Eigen::internal::traits<Derived>::Flags
+  };
 
  // typedef typename Matrix<Scalar,4,1> Coefficients;
   /** the type of a 3D vector */
@@ -179,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
@@ -222,21 +224,25 @@ struct traits<Quaternion<_Scalar,_Options> >
   typedef _Scalar Scalar;
   typedef Matrix<_Scalar,4,1,_Options> Coefficients;
   enum{
-    PacketAccess = _Options & DontAlign ? Unaligned : 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. */
@@ -267,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;
   
@@ -294,33 +307,53 @@ typedef Quaternion<double> Quaterniond;
 ***************************************************************************/
 
 namespace internal {
-template<typename _Scalar, int _PacketAccess>
-struct traits<Map<Quaternion<_Scalar>, _PacketAccess> >:
-traits<Quaternion<_Scalar> >
-{
-  typedef _Scalar Scalar;
-  typedef Map<Matrix<_Scalar,4,1>, _PacketAccess> Coefficients;
-  enum {
-    PacketAccess = _PacketAccess
+  template<typename _Scalar, int _Options>
+  struct traits<Map<Quaternion<_Scalar>, _Options> >:
+  traits<Quaternion<_Scalar, _Options> >
+  {
+    typedef _Scalar Scalar;
+    typedef Map<Matrix<_Scalar,4,1>, _Options> Coefficients;
+
+    typedef traits<Quaternion<_Scalar, _Options> > TraitsBase;
+    enum {
+      IsAligned = TraitsBase::IsAligned,
+
+      Flags = TraitsBase::Flags
+    };
+  };
+}
+
+namespace internal {
+  template<typename _Scalar, int _Options>
+  struct traits<Map<const Quaternion<_Scalar>, _Options> >:
+  traits<Quaternion<_Scalar> >
+  {
+    typedef _Scalar Scalar;
+    typedef Map<const Matrix<_Scalar,4,1>, _Options> Coefficients;
+
+    typedef traits<Quaternion<_Scalar, _Options> > TraitsBase;
+    enum {
+      IsAligned = TraitsBase::IsAligned,
+      Flags = TraitsBase::Flags & ~LvalueBit
+    };
   };
-};
 }
 
 /** \brief Quaternion expression mapping a constant memory buffer
   *
   * \param _Scalar the type of the Quaternion coefficients
-  * \param PacketAccess see class Map
+  * \param _Options see class Map
   *
   * This is a specialization of class Map for Quaternion. This class allows to view
   * a 4 scalar memory buffer as an Eigen's Quaternion object.
   *
   * \sa class Map, class Quaternion, class QuaternionBase
   */
-template<typename _Scalar, int PacketAccess>
-class Map<const Quaternion<_Scalar>, PacketAccess >
-  : public QuaternionBase<Map<const Quaternion<_Scalar>, PacketAccess> >
+template<typename _Scalar, int _Options>
+class Map<const Quaternion<_Scalar>, _Options >
+  : public QuaternionBase<Map<const Quaternion<_Scalar>, _Options> >
 {
-    typedef QuaternionBase<Map<Quaternion<_Scalar>, PacketAccess> > Base;
+    typedef QuaternionBase<Map<const Quaternion<_Scalar>, _Options> > Base;
 
   public:
     typedef _Scalar Scalar;
@@ -333,7 +366,7 @@ class Map<const Quaternion<_Scalar>, PacketAccess >
       * The pointer \a coeffs must reference the four coeffecients of Quaternion in the following order:
       * \code *coeffs == {x, y, z, w} \endcode
       *
-      * If the template parameter PacketAccess is set to Aligned, then the pointer coeffs must be aligned. */
+      * If the template parameter _Options is set to #Aligned, then the pointer coeffs must be aligned. */
     EIGEN_STRONG_INLINE Map(const Scalar* coeffs) : m_coeffs(coeffs) {}
 
     inline const Coefficients& coeffs() const { return m_coeffs;}
@@ -345,18 +378,18 @@ class Map<const Quaternion<_Scalar>, PacketAccess >
 /** \brief Expression of a quaternion from a memory buffer
   *
   * \param _Scalar the type of the Quaternion coefficients
-  * \param PacketAccess see class Map
+  * \param _Options see class Map
   *
   * This is a specialization of class Map for Quaternion. This class allows to view
   * a 4 scalar memory buffer as an Eigen's  Quaternion object.
   *
   * \sa class Map, class Quaternion, class QuaternionBase
   */
-template<typename _Scalar, int PacketAccess>
-class Map<Quaternion<_Scalar>, PacketAccess >
-  : public QuaternionBase<Map<Quaternion<_Scalar>, PacketAccess> >
+template<typename _Scalar, int _Options>
+class Map<Quaternion<_Scalar>, _Options >
+  : public QuaternionBase<Map<Quaternion<_Scalar>, _Options> >
 {
-    typedef QuaternionBase<Map<Quaternion<_Scalar>, PacketAccess> > Base;
+    typedef QuaternionBase<Map<Quaternion<_Scalar>, _Options> > Base;
 
   public:
     typedef _Scalar Scalar;
@@ -369,7 +402,7 @@ class Map<Quaternion<_Scalar>, PacketAccess >
       * The pointer \a coeffs must reference the four coeffecients of Quaternion in the following order:
       * \code *coeffs == {x, y, z, w} \endcode
       *
-      * If the template parameter PacketAccess is set to Aligned, then the pointer coeffs must be aligned. */
+      * If the template parameter _Options is set to #Aligned, then the pointer coeffs must be aligned. */
     EIGEN_STRONG_INLINE Map(Scalar* coeffs) : m_coeffs(coeffs) {}
 
     inline Coefficients& coeffs() { return m_coeffs; }
@@ -399,7 +432,7 @@ typedef Map<Quaternion<double>, Aligned>  QuaternionMapAlignedd;
 // Generic Quaternion * Quaternion product
 // This product can be specialized for a given architecture via the Arch template argument.
 namespace internal {
-template<int Arch, class Derived1, class Derived2, typename Scalar, int PacketAccess> struct quat_product
+template<int Arch, class Derived1, class Derived2, typename Scalar, int _Options> struct quat_product
 {
   EIGEN_STRONG_INLINE static Quaternion<Scalar> run(const QuaternionBase<Derived1>& a, const QuaternionBase<Derived2>& b){
     return Quaternion<Scalar>
@@ -423,7 +456,7 @@ QuaternionBase<Derived>::operator* (const QuaternionBase<OtherDerived>& other) c
    YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
   return internal::quat_product<Architecture::Target, Derived, OtherDerived,
                          typename internal::traits<Derived>::Scalar,
-                         internal::traits<Derived>::PacketAccess && internal::traits<OtherDerived>::PacketAccess>::run(*this, other);
+                         internal::traits<Derived>::IsAligned && internal::traits<OtherDerived>::IsAligned>::run(*this, other);
 }
 
 /** \sa operator*(Quaternion) */
@@ -551,6 +584,7 @@ template<class Derived>
 template<typename Derived1, typename Derived2>
 inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b)
 {
+  using std::max;
   Vector3 v0 = a.normalized();
   Vector3 v1 = b.normalized();
   Scalar c = v1.dot(v0);
@@ -565,7 +599,7 @@ inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Deri
   //    which yields a singular value problem
   if (c < Scalar(-1)+NumTraits<Scalar>::dummy_precision())
   {
-    c = std::max<Scalar>(c,-1);
+    c = max<Scalar>(c,-1);
     Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose();
     JacobiSVD<Matrix<Scalar,2,3> > svd(m, ComputeFullV);
     Vector3 axis = svd.matrixV().col(2);
@@ -625,10 +659,11 @@ template <class OtherDerived>
 inline typename internal::traits<Derived>::Scalar
 QuaternionBase<Derived>::angularDistance(const QuaternionBase<OtherDerived>& other) const
 {
+  using std::acos;
   double d = internal::abs(this->dot(other));
   if (d>=1.0)
     return Scalar(0);
-  return static_cast<Scalar>(2 * std::acos(d));
+  return static_cast<Scalar>(2 * acos(d));
 }
 
 /** \returns the spherical linear interpolation between the two quaternions
@@ -639,6 +674,7 @@ template <class OtherDerived>
 Quaternion<typename internal::traits<Derived>::Scalar>
 QuaternionBase<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& other) const
 {
+  using std::acos;
   static const Scalar one = Scalar(1) - NumTraits<Scalar>::epsilon();
   Scalar d = this->dot(other);
   Scalar absD = internal::abs(d);
@@ -646,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;
@@ -654,14 +690,13 @@ QuaternionBase<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& oth
   else
   {
     // theta is the angle between the 2 quaternions
-    Scalar theta = std::acos(absD);
+    Scalar theta = acos(absD);
     Scalar sinTheta = internal::sin(theta);
 
     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());
 }

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

Eigen/src/Geometry/Transform.h

@@ -86,11 +86,11 @@ template<typename TransformType> struct transform_take_affine_part;
   * \tparam _Scalar the scalar type, i.e., the type of the coefficients
   * \tparam _Dim the dimension of the space
   * \tparam _Mode the type of the transformation. Can be:
-  *              - Affine: the transformation is stored as a (Dim+1)^2 matrix,
-  *                        where the last row is assumed to be [0 ... 0 1].
-  *              - AffineCompact: the transformation is stored as a (Dim)x(Dim+1) matrix.
-  *              - Projective: the transformation is stored as a (Dim+1)^2 matrix
-  *                            without any assumption.
+  *              - #Affine: the transformation is stored as a (Dim+1)^2 matrix,
+  *                         where the last row is assumed to be [0 ... 0 1].
+  *              - #AffineCompact: the transformation is stored as a (Dim)x(Dim+1) matrix.
+  *              - #Projective: the transformation is stored as a (Dim+1)^2 matrix
+  *                             without any assumption.
   * \tparam _Options has the same meaning as in class Matrix. It allows to specify DontAlign and/or RowMajor.
   *                  These Options are passed directly to the underlying matrix type.
   *
@@ -672,7 +672,7 @@ Transform<Scalar,Dim,Mode,Options>::Transform(const QMatrix& other)
   *
   * This function is available only if the token EIGEN_QT_SUPPORT is defined.
   */
-template<typename Scalar, int Dim, int Mode,int Otpions>
+template<typename Scalar, int Dim, int Mode,int Options>
 Transform<Scalar,Dim,Mode,Options>& Transform<Scalar,Dim,Mode,Options>::operator=(const QMatrix& other)
 {
   EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE)
@@ -718,9 +718,13 @@ Transform<Scalar,Dim,Mode,Options>& Transform<Scalar,Dim,Mode,Options>::operator
 {
   check_template_params();
   EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE)
-  m_matrix << other.m11(), other.m21(), other.dx(),
-              other.m12(), other.m22(), other.dy(),
-              other.m13(), other.m23(), other.m33();
+  if (Mode == int(AffineCompact))
+    m_matrix << other.m11(), other.m21(), other.dx(),
+                other.m12(), other.m22(), other.dy();
+  else
+    m_matrix << other.m11(), other.m21(), other.dx(),
+                other.m12(), other.m22(), other.dy(),
+                other.m13(), other.m23(), other.m33();
   return *this;
 }
 
@@ -732,9 +736,14 @@ template<typename Scalar, int Dim, int Mode, int Options>
 QTransform Transform<Scalar,Dim,Mode,Options>::toQTransform(void) const
 {
   EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE)
-  return QTransform(matrix.coeff(0,0), matrix.coeff(1,0), matrix.coeff(2,0)
-                    matrix.coeff(0,1), matrix.coeff(1,1), matrix.coeff(2,1)
-                    matrix.coeff(0,2), matrix.coeff(1,2), matrix.coeff(2,2));
+  if (Mode == int(AffineCompact))
+    return QTransform(m_matrix.coeff(0,0), m_matrix.coeff(1,0),
+                      m_matrix.coeff(0,1), m_matrix.coeff(1,1),
+                      m_matrix.coeff(0,2), m_matrix.coeff(1,2));
+  else
+    return QTransform(m_matrix.coeff(0,0), m_matrix.coeff(1,0), m_matrix.coeff(2,0),
+                      m_matrix.coeff(0,1), m_matrix.coeff(1,1), m_matrix.coeff(2,1),
+                      m_matrix.coeff(0,2), m_matrix.coeff(1,2), m_matrix.coeff(2,2));
 }
 #endif
 
@@ -1082,10 +1091,10 @@ struct projective_transform_inverse<TransformType, Projective>
   *
   * \param hint allows to optimize the inversion process when the transformation
   * is known to be not a general transformation (optional). The possible values are:
-  *  - Projective if the transformation is not necessarily affine, i.e., if the
+  *  - #Projective if the transformation is not necessarily affine, i.e., if the
   *    last row is not guaranteed to be [0 ... 0 1]
-  *  - Affine if the last row can be assumed to be [0 ... 0 1]
-  *  - Isometry if the transformation is only a concatenations of translations
+  *  - #Affine if the last row can be assumed to be [0 ... 0 1]
+  *  - #Isometry if the transformation is only a concatenations of translations
   *    and rotations.
   *  The default is the template class parameter \c Mode.
   *

Eigen/src/Householder/Householder.h

@@ -72,8 +72,9 @@ void MatrixBase<Derived>::makeHouseholder(
 
   if(tailSqNorm == RealScalar(0) && internal::imag(c0)==RealScalar(0))
   {
-    tau = 0;
+    tau = RealScalar(0);
     beta = internal::real(c0);
+    essential.setZero();
   }
   else
   {

Eigen/src/Jacobi/Jacobi.h

@@ -104,9 +104,9 @@ bool JacobiRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z)
   else
   {
     RealScalar tau = (x-z)/(RealScalar(2)*internal::abs(y));
-    RealScalar w = internal::sqrt(internal::abs2(tau) + 1);
+    RealScalar w = internal::sqrt(internal::abs2(tau) + RealScalar(1));
     RealScalar t;
-    if(tau>0)
+    if(tau>RealScalar(0))
     {
       t = RealScalar(1) / (tau + w);
     }
@@ -114,8 +114,8 @@ bool JacobiRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z)
     {
       t = RealScalar(1) / (tau - w);
     }
-    RealScalar sign_t = t > 0 ? 1 : -1;
-    RealScalar n = RealScalar(1) / internal::sqrt(internal::abs2(t)+1);
+    RealScalar sign_t = t > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
+    RealScalar n = RealScalar(1) / internal::sqrt(internal::abs2(t)+RealScalar(1));
     m_s = - sign_t * (internal::conj(y) / internal::abs(y)) * internal::abs(t) * n;
     m_c = n;
     return true;
@@ -221,15 +221,15 @@ template<typename Scalar>
 void JacobiRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, internal::false_type)
 {
 
-  if(q==0)
+  if(q==Scalar(0))
   {
     m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1);
-    m_s = 0;
+    m_s = Scalar(0);
     if(r) *r = internal::abs(p);
   }
-  else if(p==0)
+  else if(p==Scalar(0))
   {
-    m_c = 0;
+    m_c = Scalar(0);
     m_s = q<Scalar(0) ? Scalar(1) : Scalar(-1);
     if(r) *r = internal::abs(q);
   }

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.
@@ -533,7 +525,7 @@ template<typename MatrixType>
 MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
 {
   eigen_assert(m_isInitialized && "LU is not initialized.");
-  const Index smalldim = std::min(m_lu.rows(), m_lu.cols());
+  const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
   // LU
   MatrixType res(m_lu.rows(),m_lu.cols());
   // FIXME the .toDenseMatrix() should not be needed...
@@ -695,7 +687,7 @@ struct solve_retval<FullPivLU<_MatrixType>, Rhs>
     const Index rows = dec().rows(), cols = dec().cols(),
               nonzero_pivots = dec().nonzeroPivots();
     eigen_assert(rhs().rows() == rows);
-    const Index smalldim = std::min(rows, cols);
+    const Index smalldim = (std::min)(rows, cols);
 
     if(nonzero_pivots == 0)
     {

Eigen/src/LU/PartialPivLU.h

@@ -253,7 +253,7 @@ struct partial_lu_impl
   {
     const Index rows = lu.rows();
     const Index cols = lu.cols();
-    const Index size = std::min(rows,cols);
+    const Index size = (std::min)(rows,cols);
     nb_transpositions = 0;
     int first_zero_pivot = -1;
     for(Index k = 0; k < size; ++k)
@@ -268,7 +268,7 @@ struct partial_lu_impl
 
       row_transpositions[k] = row_of_biggest_in_col;
 
-      if(biggest_in_corner != 0)
+      if(biggest_in_corner != RealScalar(0))
       {
         if(k != row_of_biggest_in_col)
         {
@@ -313,7 +313,7 @@ struct partial_lu_impl
     MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
     MatrixType lu(lu1,0,0,rows,cols);
 
-    const Index size = std::min(rows,cols);
+    const Index size = (std::min)(rows,cols);
 
     // if the matrix is too small, no blocking:
     if(size<=16)
@@ -327,14 +327,14 @@ struct partial_lu_impl
     {
       blockSize = size/8;
       blockSize = (blockSize/16)*16;
-      blockSize = std::min(std::max(blockSize,Index(8)), maxBlockSize);
+      blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
     }
 
     nb_transpositions = 0;
     int first_zero_pivot = -1;
     for(Index k = 0; k < size; k+=blockSize)
     {
-      Index bs = std::min(size-k,blockSize); // actual size of the block
+      Index bs = (std::min)(size-k,blockSize); // actual size of the block
       Index trows = rows - k - bs; // trailing rows
       Index tsize = size - k - bs; // trailing size
 

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);
@@ -182,8 +182,8 @@ struct compute_inverse_size4<Architecture::SSE, double, MatrixType, ResultType>
   };
   static void run(const MatrixType& matrix, ResultType& result)
   {
-    const EIGEN_ALIGN16 long long int _Sign_NP[2] = { 0x8000000000000000ll, 0x0000000000000000ll };
-    const EIGEN_ALIGN16 long long int _Sign_PN[2] = { 0x0000000000000000ll, 0x8000000000000000ll };
+    const __m128d _Sign_NP = _mm_castsi128_pd(_mm_set_epi32(0x0,0x0,0x80000000,0x0));
+    const __m128d _Sign_PN = _mm_castsi128_pd(_mm_set_epi32(0x80000000,0x0,0x0,0x0));
 
     // The inverse is calculated using "Divide and Conquer" technique. The
     // original matrix is divide into four 2x2 sub-matrices. Since each
@@ -316,8 +316,8 @@ struct compute_inverse_size4<Architecture::SSE, double, MatrixType, ResultType>
     iB1 = _mm_sub_pd(_mm_mul_pd(C1, dB), iB1);
     iB2 = _mm_sub_pd(_mm_mul_pd(C2, dB), iB2);
 
-    d1 = _mm_xor_pd(rd, _mm_load_pd((double*)_Sign_PN));
-    d2 = _mm_xor_pd(rd, _mm_load_pd((double*)_Sign_NP));
+    d1 = _mm_xor_pd(rd, _Sign_PN);
+    d2 = _mm_xor_pd(rd, _Sign_NP);
 
     //  iC = B*|C| - A*C#*D;
     dC = _mm_shuffle_pd(dC,dC,0);

Eigen/src/QR/ColPivHouseholderQR.h

@@ -93,7 +93,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
       */
     ColPivHouseholderQR(Index rows, Index cols)
       : m_qr(rows, cols),
-        m_hCoeffs(std::min(rows,cols)),
+        m_hCoeffs((std::min)(rows,cols)),
         m_colsPermutation(cols),
         m_colsTranspositions(cols),
         m_temp(cols),
@@ -103,7 +103,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
 
     ColPivHouseholderQR(const MatrixType& matrix)
       : m_qr(matrix.rows(), matrix.cols()),
-        m_hCoeffs(std::min(matrix.rows(),matrix.cols())),
+        m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
         m_colsPermutation(matrix.cols()),
         m_colsTranspositions(matrix.cols()),
         m_temp(matrix.cols()),
@@ -330,12 +330,12 @@ template<typename _MatrixType> class ColPivHouseholderQR
       */
     inline Index nonzeroPivots() const
     {
-      eigen_assert(m_isInitialized && "LU is not initialized.");
+      eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
       return m_nonzero_pivots;
     }
 
     /** \returns the absolute value of the biggest pivot, i.e. the biggest
-      *          diagonal coefficient of U.
+      *          diagonal coefficient of R.
       */
     RealScalar maxPivot() const { return m_maxpivot; }
 
@@ -387,7 +387,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
   for(Index k = 0; k < cols; ++k)
     m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
 
-  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / rows;
+  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
 
   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
   m_maxpivot = RealScalar(0);
@@ -413,7 +413,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
     // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so)
     // repetitions of the unit test, with the result of solve() filled with large values of the order
     // of 1/(size*epsilon).
-    if(biggest_col_sq_norm < threshold_helper * (rows-k))
+    if(biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
     {
       m_nonzero_pivots = k;
       m_hCoeffs.tail(size-k).setZero();

Eigen/src/QR/FullPivHouseholderQR.h

@@ -93,21 +93,21 @@ template<typename _MatrixType> class FullPivHouseholderQR
       */
     FullPivHouseholderQR(Index rows, Index cols)
       : m_qr(rows, cols),
-        m_hCoeffs(std::min(rows,cols)),
+        m_hCoeffs((std::min)(rows,cols)),
         m_rows_transpositions(rows),
         m_cols_transpositions(cols),
         m_cols_permutation(cols),
-        m_temp(std::min(rows,cols)),
+        m_temp((std::min)(rows,cols)),
         m_isInitialized(false),
         m_usePrescribedThreshold(false) {}
 
     FullPivHouseholderQR(const MatrixType& matrix)
       : m_qr(matrix.rows(), matrix.cols()),
-        m_hCoeffs(std::min(matrix.rows(), matrix.cols())),
+        m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
         m_rows_transpositions(matrix.rows()),
         m_cols_transpositions(matrix.cols()),
         m_cols_permutation(matrix.cols()),
-        m_temp(std::min(matrix.rows(), matrix.cols())),
+        m_temp((std::min)(matrix.rows(), matrix.cols())),
         m_isInitialized(false),
         m_usePrescribedThreshold(false)
     {
@@ -379,7 +379,7 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
 {
   Index rows = matrix.rows();
   Index cols = matrix.cols();
-  Index size = std::min(rows,cols);
+  Index size = (std::min)(rows,cols);
 
   m_qr = matrix;
   m_hCoeffs.resize(size);
@@ -493,7 +493,7 @@ struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
       RealScalar biggest_in_upper_part_of_c = c.topRows(   dec().rank()     ).cwiseAbs().maxCoeff();
       RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff();
       // FIXME brain dead
-      const RealScalar m_precision = NumTraits<Scalar>::epsilon() * std::min(rows,cols);
+      const RealScalar m_precision = NumTraits<Scalar>::epsilon() * (std::min)(rows,cols);
       // this internal:: prefix is needed by at least gcc 3.4 and ICC
       if(!internal::isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
         return;
@@ -520,7 +520,7 @@ typename FullPivHouseholderQR<MatrixType>::MatrixQType FullPivHouseholderQR<Matr
   // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
   Index rows = m_qr.rows();
   Index cols = m_qr.cols();
-  Index size = std::min(rows,cols);
+  Index size = (std::min)(rows,cols);
   MatrixQType res = MatrixQType::Identity(rows, rows);
   Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
   for (Index k = size-1; k >= 0; k--)

Eigen/src/QR/HouseholderQR.h

@@ -88,13 +88,13 @@ template<typename _MatrixType> class HouseholderQR
       */
     HouseholderQR(Index rows, Index cols)
       : m_qr(rows, cols),
-        m_hCoeffs(std::min(rows,cols)),
+        m_hCoeffs((std::min)(rows,cols)),
         m_temp(cols),
         m_isInitialized(false) {}
 
     HouseholderQR(const MatrixType& matrix)
       : m_qr(matrix.rows(), matrix.cols()),
-        m_hCoeffs(std::min(matrix.rows(),matrix.cols())),
+        m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
         m_temp(matrix.cols()),
         m_isInitialized(false)
     {
@@ -210,7 +210,7 @@ void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename
   typedef typename MatrixQR::RealScalar RealScalar;
   Index rows = mat.rows();
   Index cols = mat.cols();
-  Index size = std::min(rows,cols);
+  Index size = (std::min)(rows,cols);
 
   eigen_assert(hCoeffs.size() == size);
 
@@ -250,7 +250,7 @@ void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
 
   Index rows = mat.rows();
   Index cols = mat.cols();
-  Index size = std::min(rows, cols);
+  Index size = (std::min)(rows, cols);
 
   typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
   TempType tempVector;
@@ -260,12 +260,12 @@ void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
     tempData = tempVector.data();
   }
 
-  Index blockSize = std::min(maxBlockSize,size);
+  Index blockSize = (std::min)(maxBlockSize,size);
 
-  int k = 0;
+  Index k = 0;
   for (k = 0; k < size; k += blockSize)
   {
-    Index bs = std::min(size-k,blockSize);  // actual size of the block
+    Index bs = (std::min)(size-k,blockSize);  // actual size of the block
     Index tcols = cols - k - bs;            // trailing columns
     Index brows = rows-k;                   // rows of the block
 
@@ -299,7 +299,7 @@ struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
   template<typename Dest> void evalTo(Dest& dst) const
   {
     const Index rows = dec().rows(), cols = dec().cols();
-    const Index rank = std::min(rows, cols);
+    const Index rank = (std::min)(rows, cols);
     eigen_assert(rhs().rows() == rows);
 
     typename Rhs::PlainObject c(rhs());
@@ -327,7 +327,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
 {
   Index rows = matrix.rows();
   Index cols = matrix.cols();
-  Index size = std::min(rows,cols);
+  Index size = (std::min)(rows,cols);
 
   m_qr = matrix;
   m_hCoeffs.resize(size);

Eigen/src/SVD/JacobiSVD.h

@@ -271,13 +271,13 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
   RealScalar d = m.coeff(1,0) - m.coeff(0,1);
   if(t == RealScalar(0))
   {
-    rot1.c() = 0;
-    rot1.s() = d > 0 ? 1 : -1;
+    rot1.c() = RealScalar(0);
+    rot1.s() = d > RealScalar(0) ? RealScalar(1) : RealScalar(-1);
   }
   else
   {
     RealScalar u = d / t;
-    rot1.c() = RealScalar(1) / sqrt(1 + abs2(u));
+    rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + abs2(u));
     rot1.s() = rot1.c() * u;
   }
   m.applyOnTheLeft(0,1,rot1);
@@ -292,7 +292,7 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
   *
   * \class JacobiSVD
   *
-  * \brief Two-sided Jacobi SVD decomposition of a square matrix
+  * \brief Two-sided Jacobi SVD decomposition of a rectangular matrix
   *
   * \param MatrixType the type of the matrix of which we are computing the SVD decomposition
   * \param QRPreconditioner this optional parameter allows to specify the type of QR decomposition that will be used internally
@@ -403,8 +403,8 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
      *
      * \param matrix the matrix to decompose
      * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
-     *                           By default, none is computed. This is a bit-field, the possible bits are ComputeFullU, ComputeThinU,
-     *                           ComputeFullV, ComputeThinV.
+     *                           By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU,
+     *                           #ComputeFullV, #ComputeThinV.
      *
      * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
      * available with the (non-default) FullPivHouseholderQR preconditioner.
@@ -422,8 +422,8 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
      *
      * \param matrix the matrix to decompose
      * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed.
-     *                           By default, none is computed. This is a bit-field, the possible bits are ComputeFullU, ComputeThinU,
-     *                           ComputeFullV, ComputeThinV.
+     *                           By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU,
+     *                           #ComputeFullV, #ComputeThinV.
      *
      * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not
      * available with the (non-default) FullPivHouseholderQR preconditioner.
@@ -444,7 +444,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
     /** \returns the \a U matrix.
      *
      * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p,
-     * the U matrix is n-by-n if you asked for ComputeFullU, and is n-by-m if you asked for ComputeThinU.
+     * the U matrix is n-by-n if you asked for #ComputeFullU, and is n-by-m if you asked for #ComputeThinU.
      *
      * The \a m first columns of \a U are the left singular vectors of the matrix being decomposed.
      *
@@ -460,7 +460,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
     /** \returns the \a V matrix.
      *
      * For the SVD decomposition of a n-by-p matrix, letting \a m be the minimum of \a n and \a p,
-     * the V matrix is p-by-p if you asked for ComputeFullV, and is p-by-m if you asked for ComputeThinV.
+     * the V matrix is p-by-p if you asked for #ComputeFullV, and is p-by-m if you asked for ComputeThinV.
      *
      * The \a m first columns of \a V are the right singular vectors of the matrix being decomposed.
      *
@@ -569,7 +569,7 @@ void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, u
               "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
               "Use the ColPivHouseholderQR preconditioner instead.");
   }
-  m_diagSize = std::min(m_rows, m_cols);
+  m_diagSize = (std::min)(m_rows, m_cols);
   m_singularValues.resize(m_diagSize);
   m_matrixU.resize(m_rows, m_computeFullU ? m_rows
                           : m_computeThinU ? m_diagSize
@@ -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,9 +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.
-        if(std::max(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p)))
-            > std::max(internal::abs(m_workMatrix.coeff(p,p)),internal::abs(m_workMatrix.coeff(q,q)))*precision)
+        // keep us iterating forever. Similarly, small denormal numbers are considered zero.
+        using std::max;
+        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;
 
@@ -688,7 +693,7 @@ struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
     // A = U S V^*
     // So A^{-1} = V S^{-1} U^*
 
-    Index diagSize = std::min(dec().rows(), dec().cols());
+    Index diagSize = (std::min)(dec().rows(), dec().cols());
     typename JacobiSVDType::SingularValuesType invertedSingVals(diagSize);
 
     Index nonzeroSingVals = dec().nonzeroSingularValues();
@@ -703,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

Eigen/src/Sparse/AmbiVector.h

@@ -97,7 +97,7 @@ class AmbiVector
     void reallocateSparse()
     {
       Index copyElements = m_allocatedElements;
-      m_allocatedElements = std::min(Index(m_allocatedElements*1.5),m_size);
+      m_allocatedElements = (std::min)(Index(m_allocatedElements*1.5),m_size);
       Index allocSize = m_allocatedElements * sizeof(ListEl);
       allocSize = allocSize/sizeof(Scalar) + (allocSize%sizeof(Scalar)>0?1:0);
       Scalar* newBuffer = new Scalar[allocSize];

Eigen/src/Sparse/CompressedStorage.h

@@ -216,7 +216,7 @@ class CompressedStorage
     {
       Scalar* newValues  = new Scalar[size];
       Index* newIndices = new Index[size];
-      size_t copySize = std::min(size, m_size);
+      size_t copySize = (std::min)(size, m_size);
       // copy
       memcpy(newValues,  m_values,  copySize * sizeof(Scalar));
       memcpy(newIndices, m_indices, copySize * sizeof(Index));

Eigen/src/Sparse/DynamicSparseMatrix.h

@@ -141,7 +141,7 @@ class DynamicSparseMatrix
     {
       if (outerSize()>0)
       {
-        Index reserveSizePerVector = std::max(reserveSize/outerSize(),Index(4));
+        Index reserveSizePerVector = (std::max)(reserveSize/outerSize(),Index(4));
         for (Index j=0; j<outerSize(); ++j)
         {
           m_data[j].reserve(reserveSizePerVector);

Eigen/src/Sparse/SparseFuzzy.h

@@ -35,7 +35,7 @@
 //   const typename internal::nested<Derived,2>::type nested(derived());
 //   const typename internal::nested<OtherDerived,2>::type otherNested(other.derived());
 //   return    (nested - otherNested).cwise().abs2().sum()
-//          <= prec * prec * std::min(nested.cwise().abs2().sum(), otherNested.cwise().abs2().sum());
+//          <= prec * prec * (std::min)(nested.cwise().abs2().sum(), otherNested.cwise().abs2().sum());
 // }
 
 #endif // EIGEN_SPARSE_FUZZY_H

Eigen/src/Sparse/SparseMatrix.h

@@ -257,7 +257,7 @@ class SparseMatrix
           // furthermore we bound the realloc ratio to:
           //   1) reduce multiple minor realloc when the matrix is almost filled
           //   2) avoid to allocate too much memory when the matrix is almost empty
-          reallocRatio = std::min(std::max(reallocRatio,1.5f),8.f);
+          reallocRatio = (std::min)((std::max)(reallocRatio,1.5f),8.f);
         }
       }
       m_data.resize(m_data.size()+1,reallocRatio);

Eigen/src/Sparse/SparseMatrixBase.h

@@ -223,7 +223,7 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived>
       // thanks to shallow copies, we always eval to a tempary
       Derived temp(other.rows(), other.cols());
 
-      temp.reserve(std::max(this->rows(),this->cols())*2);
+      temp.reserve((std::max)(this->rows(),this->cols())*2);
       for (Index j=0; j<outerSize; ++j)
       {
         temp.startVec(j);
@@ -253,7 +253,7 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived>
         // eval without temporary
         derived().resize(other.rows(), other.cols());
         derived().setZero();
-        derived().reserve(std::max(this->rows(),this->cols())*2);
+        derived().reserve((std::max)(this->rows(),this->cols())*2);
         for (Index j=0; j<outerSize; ++j)
         {
           derived().startVec(j);

Eigen/src/Sparse/SparseSelfAdjointView.h

@@ -31,13 +31,13 @@
   * \brief Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix.
   *
   * \param MatrixType the type of the dense matrix storing the coefficients
-  * \param UpLo can be either \c Lower or \c Upper
+  * \param UpLo can be either \c #Lower or \c #Upper
   *
   * This class is an expression of a sefladjoint matrix from a triangular part of a matrix
   * with given dense storage of the coefficients. It is the return type of MatrixBase::selfadjointView()
   * and most of the time this is the only way that it is used.
   *
-  * \sa SparseMatrixBase::selfAdjointView()
+  * \sa SparseMatrixBase::selfadjointView()
   */
 template<typename Lhs, typename Rhs, int UpLo>
 class SparseSelfAdjointTimeDenseProduct;
@@ -383,7 +383,7 @@ void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixTyp
         continue;
                   
       Index ip = perm ? perm[i] : i;
-      count[DstUpLo==Lower ? std::min(ip,jp) : std::max(ip,jp)]++;
+      count[DstUpLo==Lower ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
     }
   }
   dest._outerIndexPtr()[0] = 0;
@@ -403,8 +403,8 @@ void permute_symm_to_symm(const MatrixType& mat, SparseMatrix<typename MatrixTyp
         continue;
                   
       Index ip = perm? perm[i] : i;
-      Index k = count[DstUpLo==Lower ? std::min(ip,jp) : std::max(ip,jp)]++;
-      dest._innerIndexPtr()[k] = DstUpLo==Lower ? std::max(ip,jp) : std::min(ip,jp);
+      Index k = count[DstUpLo==Lower ? (std::min)(ip,jp) : (std::max)(ip,jp)]++;
+      dest._innerIndexPtr()[k] = DstUpLo==Lower ? (std::max)(ip,jp) : (std::min)(ip,jp);
       
       if((DstUpLo==Lower && ip<jp) || (DstUpLo==Upper && ip>jp))
         dest._valuePtr()[k] = conj(it.value());

Eigen/src/Sparse/SparseSparseProduct.h

@@ -45,7 +45,7 @@ static void sparse_product_impl2(const Lhs& lhs, const Rhs& rhs, ResultType& res
   // estimate the number of non zero entries
   float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols()));
   float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
-  float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f);
+  float ratioRes = (std::min)(ratioLhs * avgNnzPerRhsColumn, 1.f);
 
 //  int t200 = rows/(log2(200)*1.39);
 //  int t = (rows*100)/139;
@@ -131,7 +131,7 @@ static void sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
   // estimate the number of non zero entries
   float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols()));
   float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
-  float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f);
+  float ratioRes = (std::min)(ratioLhs * avgNnzPerRhsColumn, 1.f);
 
   // mimics a resizeByInnerOuter:
   if(ResultType::IsRowMajor)
@@ -143,7 +143,7 @@ static void sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
   for (Index j=0; j<cols; ++j)
   {
     // let's do a more accurate determination of the nnz ratio for the current column j of res
-    //float ratioColRes = std::min(ratioLhs * rhs.innerNonZeros(j), 1.f);
+    //float ratioColRes = (std::min)(ratioLhs * rhs.innerNonZeros(j), 1.f);
     // FIXME find a nice way to get the number of nonzeros of a sub matrix (here an inner vector)
     float ratioColRes = ratioRes;
     tempVector.init(ratioColRes);

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

Eigen/src/StlSupport/StdVector.h

@@ -28,32 +28,24 @@
 
 #include "Eigen/src/StlSupport/details.h"
 
-// Define the explicit instantiation (e.g. necessary for the Intel compiler)
-#if defined(__INTEL_COMPILER) || defined(__GNUC__)
-  #define EIGEN_EXPLICIT_STL_VECTOR_INSTANTIATION(...) template class std::vector<__VA_ARGS__, EIGEN_ALIGNED_ALLOCATOR<__VA_ARGS__> >;
-#else
-  #define EIGEN_EXPLICIT_STL_VECTOR_INSTANTIATION(...)
-#endif
-
 /**
  * This section contains a convenience MACRO which allows an easy specialization of
  * std::vector such that for data types with alignment issues the correct allocator
  * is used automatically.
  */
 #define EIGEN_DEFINE_STL_VECTOR_SPECIALIZATION(...) \
-EIGEN_EXPLICIT_STL_VECTOR_INSTANTIATION(__VA_ARGS__) \
 namespace std \
 { \
-  template<typename _Ay> \
-  class vector<__VA_ARGS__, _Ay>  \
+  template<> \
+  class vector<__VA_ARGS__, std::allocator<__VA_ARGS__> >  \
     : public vector<__VA_ARGS__, EIGEN_ALIGNED_ALLOCATOR<__VA_ARGS__> > \
   { \
     typedef vector<__VA_ARGS__, EIGEN_ALIGNED_ALLOCATOR<__VA_ARGS__> > vector_base; \
   public: \
     typedef __VA_ARGS__ value_type; \
-    typedef typename vector_base::allocator_type allocator_type; \
-    typedef typename vector_base::size_type size_type;  \
-    typedef typename vector_base::iterator iterator;  \
+    typedef vector_base::allocator_type allocator_type; \
+    typedef vector_base::size_type size_type;  \
+    typedef vector_base::iterator iterator;  \
     explicit vector(const allocator_type& a = allocator_type()) : vector_base(a) {}  \
     template<typename InputIterator> \
     vector(InputIterator first, InputIterator last, const allocator_type& a = allocator_type()) : vector_base(first, last, a) {} \

bench/BenchTimer.h

@@ -119,7 +119,7 @@ public:
 
   inline double getCpuTime()
   {
-#ifdef WIN32
+#ifdef _WIN32
     LARGE_INTEGER query_ticks;
     QueryPerformanceCounter(&query_ticks);
     return query_ticks.QuadPart/m_frequency;
@@ -132,7 +132,7 @@ public:
 
   inline double getRealTime()
   {
-#ifdef WIN32
+#ifdef _WIN32
     SYSTEMTIME st;
     GetSystemTime(&st);
     return (double)st.wSecond + 1.e-3 * (double)st.wMilliseconds;

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

debug/gdb/printers.py

@@ -56,12 +56,12 @@ class EigenMatrixPrinter:
 		template_params = m.split(',')
 		template_params = map(lambda x:x.replace(" ", ""), template_params)
 
-		if template_params[1] == '-0x00000000000000001':
+		if template_params[1] == '-0x00000000000000001' or template_params[1] == '-0x000000001':
 			self.rows = val['m_storage']['m_rows']
 		else:
 			self.rows = int(template_params[1])
 		
-		if template_params[2] == '-0x00000000000000001':
+		if template_params[2] == '-0x00000000000000001' or template_params[2] == '-0x000000001':
 			self.cols = val['m_storage']['m_cols']
 		else:
 			self.cols = int(template_params[2])

doc/A05_PortingFrom2To3.dox

@@ -287,7 +287,7 @@ The EIGEN_DONT_ALIGN option still exists in Eigen 3, but it has a new cousin: EI
 
 A common issue with Eigen 2 was that when mapping an array with Map, there was no way to tell Eigen that your array was aligned. There was a ForceAligned option but it didn't mean that; it was just confusing and has been removed.
 
-New in Eigen3 is the Aligned option. See the documentation of class Map. Use it like this:
+New in Eigen3 is the #Aligned option. See the documentation of class Map. Use it like this:
 \code
 Map<Vector4f, Aligned> myMappedVector(some_aligned_array);
 \endcode

doc/C00_QuickStartGuide.dox

@@ -63,7 +63,7 @@ The output is as follows:
 
 The second example starts by declaring a 3-by-3 matrix \c m which is initialized using the \link DenseBase::Random(Index,Index) Random() \endlink method with random values between -1 and 1. The next line applies a linear mapping such that the values are between 10 and 110. The function call \link DenseBase::Constant(Index,Index,const Scalar&) MatrixXd::Constant\endlink(3,3,1.2) returns a 3-by-3 matrix expression having all coefficients equal to 1.2. The rest is standard arithmetics.
 
-The next line of the \c main function introduces a new type: \c VectorXd. This represents a (column) vector of arbitrary size. Here, the vector \c v is created to contains \c 3 coefficients which are left unitialized. The one but last line uses the so-called comma-initializer, explained in \ref TutorialAdvancedInitialization, to set all coefficients of the vector \c v to be as follows:
+The next line of the \c main function introduces a new type: \c VectorXd. This represents a (column) vector of arbitrary size. Here, the vector \c v is created to contain \c 3 coefficients which are left unitialized. The one but last line uses the so-called comma-initializer, explained in \ref TutorialAdvancedInitialization, to set all coefficients of the vector \c v to be as follows:
 
 \f[
 v =

doc/C03_TutorialArrayClass.dox

@@ -37,7 +37,7 @@ we won't explain it again here and just refer to \ref TutorialMatrixClass.
 
 Eigen also provides typedefs for some common cases, in a way that is similar to the Matrix typedefs
 but with some slight differences, as the word "array" is used for both 1-dimensional and 2-dimensional arrays.
-We adopt that convention that typedefs of the form ArrayNt stand for 1-dimensional arrays, where N and t are
+We adopt the convention that typedefs of the form ArrayNt stand for 1-dimensional arrays, where N and t are
 the size and the scalar type, as in the Matrix typedefs explained on \ref TutorialMatrixClass "this page". For 2-dimensional arrays, we
 use typedefs of the form ArrayNNt. Some examples are shown in the following table:
 
@@ -104,8 +104,8 @@ This provides a functionality that is not directly available for Matrix objects.
 
 First of all, of course you can multiply an array by a scalar, this works in the same way as matrices. Where arrays
 are fundamentally different from matrices, is when you multiply two together. Matrices interpret
-multiplication as the matrix product and arrays interpret multiplication as the coefficient-wise product. Thus, two 
-arrays can be multiplied if they have the same size.
+multiplication as matrix product and arrays interpret multiplication as coefficient-wise product. Thus, two 
+arrays can be multiplied if and only if they have the same dimensions.
 
 <table class="example">
 <tr><th>Example:</th><th>Output:</th></tr>
@@ -119,8 +119,8 @@ arrays can be multiplied if they have the same size.
 
 \section TutorialArrayClassCwiseOther Other coefficient-wise operations
 
-The Array class defined other coefficient-wise operations besides the addition, subtraction and multiplication
-operators described about. For example, the \link ArrayBase::abs() .abs() \endlink method takes the absolute
+The Array class defines other coefficient-wise operations besides the addition, subtraction and multiplication
+operators described above. For example, the \link ArrayBase::abs() .abs() \endlink method takes the absolute
 value of each coefficient, while \link ArrayBase::sqrt() .sqrt() \endlink computes the square root of the
 coefficients. If you have two arrays of the same size, you can call \link ArrayBase::min() .min() \endlink to
 construct the array whose coefficients are the minimum of the corresponding coefficients of the two given
@@ -155,7 +155,7 @@ this doesn't have any runtime cost (provided that you let your compiler optimize
 Both \link MatrixBase::array() .array() \endlink and \link ArrayBase::matrix() .matrix() \endlink 
 can be used as rvalues and as lvalues.
 
-Mixing matrices and arrays in an expression is forbidden with Eigen. For instance, you cannot add a amtrix and
+Mixing matrices and arrays in an expression is forbidden with Eigen. For instance, you cannot add a matrix and
 array directly; the operands of a \c + operator should either both be matrices or both be arrays. However,
 it is easy to convert from one to the other with \link MatrixBase::array() .array() \endlink and 
 \link ArrayBase::matrix() .matrix()\endlink. The exception to this rule is the assignment operator: it is

doc/C07_TutorialReductionsVisitorsBroadcasting.dox

@@ -93,7 +93,7 @@ Array.
 
 The arguments passed to a visitor are pointers to the variables where the
 row and column position are to be stored. These variables should be of type
-\link DenseBase::Index Index \endlink (FIXME: link ok?), as shown below:
+\link DenseBase::Index Index \endlink, as shown below:
 
 <table class="example">
 <tr><th>Example:</th><th>Output:</th></tr>
@@ -141,7 +141,7 @@ return a 'column-vector'</b>
 
 \subsection TutorialReductionsVisitorsBroadcastingPartialReductionsCombined Combining partial reductions with other operations
 It is also possible to use the result of a partial reduction to do further processing.
-Here is another example that aims to find the column whose sum of elements is the maximum
+Here is another example that finds the column whose sum of elements is the maximum
  within a matrix. With column-wise partial reductions this can be coded as:
 
 <table class="example">

doc/C08_TutorialGeometry.dox

@@ -6,7 +6,7 @@ namespace Eigen {
 \li \b Previous: \ref TutorialReductionsVisitorsBroadcasting
 \li \b Next: \ref TutorialSparse
 
-In this tutorial, we will shortly introduce the many possibilities offered by the \ref Geometry_Module "geometry module", namely 2D and 3D rotations and projective or affine transformations.
+In this tutorial, we will briefly introduce the many possibilities offered by the \ref Geometry_Module "geometry module", namely 2D and 3D rotations and projective or affine transformations.
 
 \b Table \b of \b contents
   - \ref TutorialGeoElementaryTransformations
@@ -48,7 +48,7 @@ Rotation2D<float> rot2(angle_in_radian);\endcode</td></tr>
 AngleAxis<float> aa(angle_in_radian, Vector3f(ax,ay,az));\endcode</td></tr>
 <tr><td>
 3D rotation as a \ref Quaternion "quaternion"</td><td>\code
-Quaternion<float> q = AngleAxis<float>(angle_in_radian, axis);\endcode</td></tr>
+Quaternion<float> q;  q = AngleAxis<float>(angle_in_radian, axis);\endcode</td></tr>
 <tr class="alt"><td>
 N-D Scaling</td><td>\code
 Scaling<float,2>(sx, sy)
@@ -78,7 +78,7 @@ representations are rotation matrices, while for other usages Quaternion is the
 representation of choice as they are compact, fast and stable. Finally Rotation2D and
 AngleAxis are mainly convenient types to create other rotation objects.
 
-<strong>Notes on Translation and Scaling</strong>\n Likewise AngleAxis, these classes were
+<strong>Notes on Translation and Scaling</strong>\n Like AngleAxis, these classes were
 designed to simplify the creation/initialization of linear (Matrix) and affine (Transform)
 transformations. Nevertheless, unlike AngleAxis which is inefficient to use, these classes
 might still be interesting to write generic and efficient algorithms taking as input any
@@ -88,13 +88,13 @@ Any of the above transformation types can be converted to any other types of the
 or to a more generic type. Here are some additional examples:
 <table class="manual">
 <tr><td>\code
-Rotation2Df r = Matrix2f(..);       // assumes a pure rotation matrix
-AngleAxisf aa = Quaternionf(..);
-AngleAxisf aa = Matrix3f(..);       // assumes a pure rotation matrix
-Matrix2f m = Rotation2Df(..);
-Matrix3f m = Quaternionf(..);       Matrix3f m = Scaling3f(..);
-Affine3f m = AngleAxis3f(..);       Affine3f m = Scaling3f(..);
-Affine3f m = Translation3f(..);     Affine3f m = Matrix3f(..);
+Rotation2Df r;  r  = Matrix2f(..);       // assumes a pure rotation matrix
+AngleAxisf aa;  aa = Quaternionf(..);
+AngleAxisf aa;  aa = Matrix3f(..);       // assumes a pure rotation matrix
+Matrix2f m;     m  = Rotation2Df(..);
+Matrix3f m;     m  = Quaternionf(..);       Matrix3f m;   m = Scaling3f(..);
+Affine3f m;     m  = AngleAxis3f(..);       Affine3f m;   m = Scaling3f(..);
+Affine3f m;     m  = Translation3f(..);     Affine3f m;   m = Matrix3f(..);
 \endcode</td></tr>
 </table>
 
@@ -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>
@@ -186,7 +186,7 @@ matNxN = t.extractRotation();
 While transformation objects can be created and updated concatenating elementary transformations,
 the Transform class also features a procedural API:
 <table class="manual">
-<tr><th></th><th>procedurale API</th><th>equivalent natural API </th></tr>
+<tr><th></th><th>procedural API</th><th>equivalent natural API </th></tr>
 <tr><td>Translation</td><td>\code
 t.translate(Vector_(tx,ty,..));
 t.pretranslate(Vector_(tx,ty,..));
@@ -234,7 +234,7 @@ t = Translation_(..) * t * RotationType(..) * Translation_(..) * Scaling_(..);
 <table class="manual">
 <tr><td style="max-width:30em;">
 Euler angles might be convenient to create rotation objects.
-On the other hand, since there exist 24 differents convension,they are pretty confusing to use. This example shows how
+On the other hand, since there exist 24 different conventions, they are pretty confusing to use. This example shows how
 to create a rotation matrix according to the 2-1-2 convention.</td><td>\code
 Matrix3f m;
 m = AngleAxisf(angle1, Vector3f::UnitZ())

doc/C09_TutorialSparse.dox

@@ -4,7 +4,7 @@ namespace Eigen {
     \ingroup Tutorial
 
 \li \b Previous: \ref TutorialGeometry
-\li \b Next: TODO
+\li \b Next: \ref TutorialMapClass
 
 \b Table \b of \b contents \n
   - \ref TutorialSparseIntro
@@ -55,17 +55,17 @@ and its internal representation using the Compressed Column Storage format:
 </table>
 Outer indices:<table class="manual"><tr><td>0</td><td>2</td><td>4</td><td>5</td><td>6</td><td>\em 7 </td></tr></table>
 
-As you can guess, here the storage order is even more important than with dense matrix. We will therefore often make a clear difference between the \em inner and \em outer dimensions. For instance, it is easy to loop over the coefficients of an \em inner \em vector (e.g., a column of a column-major matrix), but completely inefficient to do the same for an \em outer \em vector (e.g., a row of a col-major matrix).
+As you might guess, here the storage order is even more important than with dense matrices. We will therefore often make a clear difference between the \em inner and \em outer dimensions. For instance, it is efficient to loop over the coefficients of an \em inner \em vector (e.g., a column of a column-major matrix), but completely inefficient to do the same for an \em outer \em vector (e.g., a row of a column-major matrix).
 
 The SparseVector class implements the same compressed storage scheme but, of course, without any outer index buffer.
 
-Since all nonzero coefficients of such a matrix are sequentially stored in memory, random insertion of new nonzeros can be extremely costly. To overcome this limitation, Eigen's sparse module provides a DynamicSparseMatrix class which is basically implemented as an array of SparseVector. In other words, a DynamicSparseMatrix is a SparseMatrix where the values and inner-indices arrays have been splitted into multiple small and resizable arrays. Assuming the number of nonzeros per inner vector is relatively low, this slight modification allow for very fast random insertion at the cost of a slight memory overhead and a lost of compatibility with other sparse libraries used by some of our highlevel solvers. Note that the major memory overhead comes from the extra memory preallocated by each inner vector to avoid an expensive memory reallocation at every insertion.
+Since all nonzero coefficients of such a matrix are sequentially stored in memory, inserting a new nonzero near the "beginning" of the matrix can be extremely costly. As described below (\ref TutorialSparseFilling), one strategy is to fill nonzero coefficients in order. In cases where this is not possible, Eigen's sparse module also provides a DynamicSparseMatrix class which allows efficient random insertion. DynamicSparseMatrix is essentially implemented as an array of SparseVector, where the values and inner-indices arrays have been split into multiple small and resizable arrays. Assuming the number of nonzeros per inner vector is relatively small, this modification allows for very fast random insertion at the cost of a slight memory overhead (due to extra memory preallocated by each inner vector to avoid an expensive memory reallocation at every insertion) and a loss of compatibility with other sparse libraries used by some of our high-level solvers. Once complete, a DynamicSparseMatrix can be converted to a SparseMatrix to permit usage of these sparse libraries.
 
-To summarize, it is recommanded to use a SparseMatrix whenever this is possible, and reserve the use of DynamicSparseMatrix for matrix assembly purpose when a SparseMatrix is not flexible enough. The respective pro/cons of both representations are summarized in the following table:
+To summarize, it is recommended to use SparseMatrix whenever possible, and reserve the use of DynamicSparseMatrix to assemble a sparse matrix in cases when a SparseMatrix is not flexible enough. The respective pros/cons of both representations are summarized in the following table:
 
 <table class="manual">
 <tr><td></td> <td>SparseMatrix</td><td>DynamicSparseMatrix</td></tr>
-<tr><td>memory usage</td><td>***</td><td>**</td></tr>
+<tr><td>memory efficiency</td><td>***</td><td>**</td></tr>
 <tr><td>sorted insertion</td><td>***</td><td>***</td></tr>
 <tr><td>random insertion \n in sorted inner vector</td><td>**</td><td>**</td></tr>
 <tr><td>sorted insertion \n in random inner vector</td><td>-</td><td>***</td></tr>
@@ -82,7 +82,7 @@ To summarize, it is recommanded to use a SparseMatrix whenever this is possible,
 
 \b Matrix \b and \b vector \b properties \n
 
-Here mat and vec represents any sparse-matrix and sparse-vector types respectively.
+Here mat and vec represent any sparse-matrix and sparse-vector type, respectively.
 
 <table class="manual">
 <tr><td>Standard \n dimensions</td><td>\code
@@ -96,7 +96,7 @@ mat.innerSize()
 mat.outerSize()\endcode</td>
 <td></td>
 </tr>
-<tr><td>Number of non \n zero coefficiens</td><td>\code
+<tr><td>Number of non \n zero coefficients</td><td>\code
 mat.nonZeros() \endcode</td>
 <td>\code
 vec.nonZeros() \endcode</td></tr>
@@ -105,12 +105,12 @@ vec.nonZeros() \endcode</td></tr>
 
 \b Iterating \b over \b the \b nonzero \b coefficients \n
 
-Iterating over the coefficients of a sparse matrix can be done only in the same order than the storage order. Here is an example:
+Iterating over the coefficients of a sparse matrix can be done only in the same order as the storage order. Here is an example:
 <table class="manual">
 <tr><td>
 \code
 SparseMatrixType mat(rows,cols);
-for (int k=0; k\<m1.outerSize(); ++k)
+for (int k=0; k<m1.outerSize(); ++k)
   for (SparseMatrixType::InnerIterator it(mat,k); it; ++it)
   {
     it.value();
@@ -254,7 +254,7 @@ if(!lu_of_A.solve(b,&x)) {
 
 See also the class SparseLLT, class SparseLU, and class SparseLDLT.
 
-\li \b Next: TODO
+\li \b Next: \ref TutorialMapClass
 
 */
 

doc/C10_TutorialMapClass.dox

@@ -0,0 +1,96 @@
+namespace Eigen {
+
+/** \page TutorialMapClass Tutorial page 10 - Interfacing with C/C++ arrays and external libraries: the %Map class
+
+\ingroup Tutorial
+
+\li \b Previous: \ref TutorialSparse
+\li \b Next: \ref TODO
+
+This tutorial page explains how to work with "raw" C++ arrays.  This can be useful in a variety of contexts, particularly when "importing" vectors and matrices from other libraries into Eigen.
+
+\b Table \b of \b contents
+  - \ref TutorialMapIntroduction
+  - \ref TutorialMapTypes
+  - \ref TutorialMapUsing
+  - \ref TutorialMapPlacementNew
+
+\section TutorialMapIntroduction Introduction
+
+Occasionally you may have a pre-defined array of numbers that you want to use within Eigen as a vector or matrix. While one option is to make a copy of the data, most commonly you probably want to re-use this memory as an Eigen type. Fortunately, this is very easy with the Map class.
+
+\section TutorialMapTypes Map types and declaring Map variables
+
+A Map object has a type defined by its Eigen equivalent:
+\code
+Map<Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime> >
+\endcode
+Note that, in this default case, a Map requires just a single template parameter.  
+
+To construct a Map variable, you need two other pieces of information: a pointer to the region of memory defining the array of coefficients, and the desired shape of the matrix or vector.  For example, to define a matrix of \c float with sizes determined at compile time, you might do the following:
+\code
+Map<MatrixXf> mf(pf,rows,columns);
+\endcode
+where \c pf is a \c float \c * pointing to the array of memory.  A fixed-size read-only vector of integers might be declared as
+\code
+Map<const Vector4i> mi(pi);
+\endcode
+where \c pi is an \c int \c *. In this case the size does not have to be passed to the constructor, because it is already specified by the Matrix/Array type.
+
+Note that Map does not have a default constructor; you \em must pass a pointer to intialize the object. However, you can work around this requirement (see \ref TutorialMapPlacementNew).
+
+Map is flexible enough to accomodate a variety of different data representations.  There are two other (optional) template parameters:
+\code
+Map<typename MatrixType,
+    int MapOptions,
+    typename StrideType>
+\endcode
+\li \c MapOptions specifies whether the pointer is \c #Aligned, or \c #Unaligned.  The default is \c #Unaligned.
+\li \c StrideType allows you to specify a custom layout for the memory array, using the Stride class.  One example would be to specify that the data array is organized in row-major format:
+<table class="example">
+<tr><th>Example:</th><th>Output:</th></tr>
+<tr>
+<td>\include Tutorial_Map_rowmajor.cpp </td>
+<td>\verbinclude Tutorial_Map_rowmajor.out </td>
+</table>
+However, Stride is even more flexible than this; for details, see the documentation for the Map and Stride classes.
+
+\section TutorialMapUsing Using Map variables
+
+You can use a Map object just like any other Eigen type:
+<table class="example">
+<tr><th>Example:</th><th>Output:</th></tr>
+<tr>
+<td>\include Tutorial_Map_using.cpp </td>
+<td>\verbinclude Tutorial_Map_using.out </td>
+</table>
+
+However, when writing functions taking Eigen types, it is important to realize that a Map type is \em not identical to its Dense equivalent.  See \ref TopicFunctionTakingEigenTypesMultiarguments for details.
+
+\section TutorialMapPlacementNew Changing the mapped array
+
+It is possible to change the array of a Map object after declaration, using the C++ "placement new" syntax:
+<table class="example">
+<tr><th>Example:</th><th>Output:</th></tr>
+<tr>
+<td>\include Map_placement_new.cpp </td>
+<td>\verbinclude Map_placement_new.out </td>
+</table>
+Despite appearances, this does not invoke the memory allocator, because the syntax specifies the location for storing the result.
+
+This syntax makes it possible to declare a Map object without first knowing the mapped array's location in memory:
+\code
+Map<Matrix3f> A(NULL);  // don't try to use this matrix yet!
+VectorXf b(n_matrices);
+for (int i = 0; i < n_matrices; i++)
+{
+  new (&A) Map<Matrix3f>(get_matrix_pointer(i));
+  b(i) = A.trace();
+}
+\endcode
+
+\li \b Next: \ref TODO
+
+*/
+
+}

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/Doxyfile.in

@@ -488,7 +488,7 @@ SHOW_FILES             = YES
 # Namespaces page.  This will remove the Namespaces entry from the Quick Index
 # and from the Folder Tree View (if specified). The default is YES.
 
-SHOW_NAMESPACES        = NO
+SHOW_NAMESPACES        = YES
 
 # The FILE_VERSION_FILTER tag can be used to specify a program or script that
 # doxygen should invoke to get the current version for each file (typically from

doc/I03_InsideEigenExample.dox

@@ -400,7 +400,7 @@ inline void writePacket(int index, const PacketScalar& x)
   internal::pstoret<Scalar, PacketScalar, StoreMode>(m_storage.data() + index, x);
 }
 \endcode
-Here, \a StoreMode is \a Aligned, indicating that we are doing a 128-bit-aligned write access, \a PacketScalar is a type representing a "SSE packet of 4 floats" and internal::pstoret is a function writing such a packet in memory. Their definitions are architecture-specific, we find them in src/Core/arch/SSE/PacketMath.h:
+Here, \a StoreMode is \a #Aligned, indicating that we are doing a 128-bit-aligned write access, \a PacketScalar is a type representing a "SSE packet of 4 floats" and internal::pstoret is a function writing such a packet in memory. Their definitions are architecture-specific, we find them in src/Core/arch/SSE/PacketMath.h:
 
 The line in src/Core/arch/SSE/PacketMath.h that determines the PacketScalar type (via a typedef in Matrix.h) is:
 \code
@@ -442,7 +442,7 @@ class CwiseBinaryOp
     }
 };
 \endcode
-Here, \a m_lhs is the vector \a v, and \a m_rhs is the vector \a w. So the packet() function here is Matrix::packet(). The template parameter \a LoadMode is \a Aligned. So we're looking at
+Here, \a m_lhs is the vector \a v, and \a m_rhs is the vector \a w. So the packet() function here is Matrix::packet(). The template parameter \a LoadMode is \a #Aligned. So we're looking at
 \code
 class Matrix
 {

doc/I13_FunctionsTakingEigenTypes.dox

@@ -47,9 +47,9 @@ void print_block(const DenseBase<Derived>& b, int x, int y, int r, int c)
 Prints the maximum coefficient of the array or array-expression.
 \code
 template <typename Derived>
-void print_max(const ArrayBase<Derived>& a, const ArrayBase<Derived>& b)
+void print_max_coeff(const ArrayBase<Derived> &a)
 {
-  std::cout << "max: " << (a.max(b)).maxCoeff() << std::endl;
+  std::cout << "max: " << a.maxCoeff() << std::endl;
 }
 \endcode
 <b> %MatrixBase Example </b><br/><br/>
@@ -63,6 +63,21 @@ void print_inv_cond(const MatrixBase<Derived>& a)
   std::cout << "inv cond: " << sing_vals(sing_vals.size()-1) / sing_vals(0) << std::endl;
 }
 \endcode
+<b> Multiple templated arguments example </b><br/><br/>
+Calculate the Euclidean distance between two points.
+\code
+template <typename DerivedA,typename DerivedB>
+typename DerivedA::Scalar squaredist(const MatrixBase<DerivedA>& p1,const MatrixBase<DerivedB>& p2)
+{
+  return (p1-p2).squaredNorm();
+}
+\endcode
+Notice that we used two template parameters, one per argument. This permits the function to handle inputs of different types, e.g.,
+\code
+squaredist(v1,2*v2)
+\endcode
+where the first argument \c v1 is a vector and the second argument \c 2*v2 is an expression.
+<br/><br/>
 
 These examples are just intended to give the reader a first impression of how functions can be written which take a plain and constant Matrix or Array argument. They are also intended to give the reader an idea about the most common base classes being the optimal candidates for functions. In the next section we will look in more detail at an example and the different ways it can be implemented, while discussing each implementation's problems and advantages. For the discussion below, Matrix and Array as well as MatrixBase and ArrayBase can be exchanged and all arguments still hold.
 
@@ -128,6 +143,8 @@ The implementation above does now not only work with temporary expressions but i
 
 \b Note: The const cast hack will only work with templated functions. It will not work with the MatrixXf implementation because it is not possible to cast a Block expression to a Matrix reference!
 
+
+
 \section TopicResizingInGenericImplementations How to resize matrices in generic implementations?
 
 One might think we are done now, right? This is not completely true because in order for our covariance function to be generically applicable, we want the follwing code to work

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.
 
@@ -27,6 +27,7 @@ For a first contact with Eigen, the best place is to have a look at the \ref Get
     - \ref TutorialReductionsVisitorsBroadcasting
     - \ref TutorialGeometry
     - \ref TutorialSparse
+    - \ref TutorialMapClass
   - \ref QuickRefPage
   - <b>Advanced topics</b>
     - \ref TopicAliasing