bump to 3.0.1

(transplanted from 9464745385
)

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>
@@ -184,6 +185,7 @@
 // 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 +241,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/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

@@ -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.

Eigen/src/Core/DenseBase.h

@@ -56,7 +56,13 @@ template<typename Derived> class DenseBase
     class InnerIterator;
 
     typedef typename internal::traits<Derived>::StorageKind StorageKind;
-    typedef typename internal::traits<Derived>::Index Index; /**< The type of indices */
+
+    /** \brief The type of indices 
+      * \details To change this, \c \#define the preprocessor symbol \c EIGEN_DEFAULT_DENSE_INDEX_TYPE.
+      * \sa \ref TopicPreprocessorDirectives.
+      */
+    typedef typename internal::traits<Derived>::Index Index; 
+
     typedef typename internal::traits<Derived>::Scalar Scalar;
     typedef typename internal::packet_traits<Scalar>::type PacketScalar;
     typedef typename NumTraits<Scalar>::Real RealScalar;

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,14 +58,14 @@ 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
 
 template <typename T, int Size, int MatrixOrArrayOptions>
 struct plain_array<T, Size, MatrixOrArrayOptions, 16>
 {
-  EIGEN_ALIGN16 T array[Size];
+  EIGEN_USER_ALIGN16 T array[Size];
   plain_array() { EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(0xf) }
   plain_array(constructor_without_unaligned_array_assert) {}
 };
@@ -73,7 +73,7 @@ struct plain_array<T, Size, MatrixOrArrayOptions, 16>
 template <typename T, int MatrixOrArrayOptions, int Alignment>
 struct plain_array<T, 0, MatrixOrArrayOptions, Alignment>
 {
-  EIGEN_ALIGN16 T array[1];
+  EIGEN_USER_ALIGN16 T array[1];
   plain_array() {}
   plain_array(constructor_without_unaligned_array_assert) {}
 };

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
@@ -270,6 +270,14 @@ Packet psqrt(const Packet& a) { return sqrt(a); }
 * The following functions might not have to be overwritten for vectorized types
 ***************************************************************************/
 
+/** \internal copy a packet with constant coeficient \a a (e.g., [a,a,a,a]) to \a *to. \a to must be 16 bytes aligned */
+// NOTE: this function must really be templated on the packet type (think about different packet types for the same scalar type)
+template<typename Packet>
+inline void pstore1(typename unpacket_traits<Packet>::type* to, const typename unpacket_traits<Packet>::type& a)
+{
+  pstore(to, pset1<Packet>(a));
+}
+
 /** \internal \returns a * b + c (coeff-wise) */
 template<typename Packet> inline Packet
 pmadd(const Packet&  a,
@@ -278,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)
 {
@@ -289,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 optionnally 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.
   *
@@ -95,7 +95,7 @@ struct traits<Map<PlainObjectType, MapOptions, StrideType> >
     HasNoInnerStride = InnerStrideAtCompileTime == 1,
     HasNoOuterStride = StrideType::OuterStrideAtCompileTime == 0,
     HasNoStride = HasNoInnerStride && HasNoOuterStride,
-    IsAligned = int(int(MapOptions)&Aligned)==Aligned,
+    IsAligned = bool(EIGEN_ALIGN) && ((int(MapOptions)&Aligned)==Aligned),
     IsDynamicSize = PlainObjectType::SizeAtCompileTime==Dynamic,
     KeepsPacketAccess = bool(HasNoInnerStride)
                         && ( bool(IsDynamicSize)

Eigen/src/Core/MapBase.h

@@ -85,6 +85,8 @@ template<typename Derived> class MapBase<Derived, ReadOnlyAccessors>
     using Base::rowStride;
     using Base::colStride;
 
+    // bug 217 - compile error on ICC 11.1
+    using Base::operator=;
 
     typedef typename Base::CoeffReturnType CoeffReturnType;
 
@@ -236,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").

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

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)
@@ -107,17 +111,17 @@ MatrixBase<Derived>::blueNorm() const
     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/TriangularMatrix.h

@@ -134,13 +134,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 +448,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 +468,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,6 +486,7 @@ 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)
@@ -493,7 +496,7 @@ struct triangular_assignment_selector<Derived1, Derived2, Upper, Dynamic, ClearO
         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);
     }
   }
 };
@@ -511,7 +514,7 @@ struct triangular_assignment_selector<Derived1, Derived2, Lower, Dynamic, ClearO
       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);
     }
   }
 };
@@ -547,7 +550,7 @@ struct triangular_assignment_selector<Derived1, Derived2, StrictlyLower, Dynamic
       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);
     }
   }
 };
@@ -756,8 +759,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

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

@@ -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

@@ -100,12 +100,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 = { 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 = { 0, 1, 2, 3 };
   return vaddq_s32(pset1<Packet4i>(a), countdown);
 }
 
@@ -191,14 +191,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);
 }
 

Eigen/src/Core/arch/SSE/PacketMath.h

@@ -305,6 +305,19 @@ template<> EIGEN_STRONG_INLINE void pstoreu<double>(double* to, const Packet2d&
 template<> EIGEN_STRONG_INLINE void pstoreu<float>(float*  to, const Packet4f& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((double*)to, _mm_castps_pd(from)); }
 template<> EIGEN_STRONG_INLINE void pstoreu<int>(int*      to, const Packet4i& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu((double*)to, _mm_castsi128_pd(from)); }
 
+// some compilers might be tempted to perform multiple moves instead of using a vector path.
+template<> EIGEN_STRONG_INLINE void pstore1<Packet4f>(float* to, const float& a)
+{
+  Packet4f pa = _mm_set_ss(a);
+  pstore(to, vec4f_swizzle1(pa,0,0,0,0));
+}
+// some compilers might be tempted to perform multiple moves instead of using a vector path.
+template<> EIGEN_STRONG_INLINE void pstore1<Packet2d>(double* to, const double& a)
+{
+  Packet2d pa = _mm_set_sd(a);
+  pstore(to, vec2d_swizzle1(pa,0,0));
+}
+
 template<> EIGEN_STRONG_INLINE void prefetch<float>(const float*   addr) { _mm_prefetch((const char*)(addr), _MM_HINT_T0); }
 template<> EIGEN_STRONG_INLINE void prefetch<double>(const double* addr) { _mm_prefetch((const char*)(addr), _MM_HINT_T0); }
 template<> EIGEN_STRONG_INLINE void prefetch<int>(const int*       addr) { _mm_prefetch((const char*)(addr), _MM_HINT_T0); }

Eigen/src/Core/products/GeneralBlockPanelKernel.h

@@ -199,7 +199,7 @@ public:
   EIGEN_STRONG_INLINE void unpackRhs(DenseIndex n, const RhsScalar* rhs, RhsScalar* b)
   {
     for(DenseIndex k=0; k<n; k++)
-      pstore(&b[k*RhsPacketSize], pset1<RhsPacket>(rhs[k]));
+      pstore1<RhsPacket>(&b[k*RhsPacketSize], rhs[k]);
   }
 
   EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacket& dest) const
@@ -270,7 +270,7 @@ public:
   EIGEN_STRONG_INLINE void unpackRhs(DenseIndex n, const RhsScalar* rhs, RhsScalar* b)
   {
     for(DenseIndex k=0; k<n; k++)
-      pstore(&b[k*RhsPacketSize], pset1<RhsPacket>(rhs[k]));
+      pstore1<RhsPacket>(&b[k*RhsPacketSize], rhs[k]);
   }
 
   EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacket& dest) const
@@ -363,8 +363,8 @@ public:
     {
       if(Vectorizable)
       {
-        pstore((RealScalar*)&b[k*ResPacketSize*2+0], pset1<RealPacket>(real(rhs[k])));
-        pstore((RealScalar*)&b[k*ResPacketSize*2+ResPacketSize], pset1<RealPacket>(imag(rhs[k])));
+        pstore1<RealPacket>((RealScalar*)&b[k*ResPacketSize*2+0],             real(rhs[k]));
+        pstore1<RealPacket>((RealScalar*)&b[k*ResPacketSize*2+ResPacketSize], imag(rhs[k]));
       }
       else
         b[k] = rhs[k];
@@ -475,7 +475,7 @@ public:
   EIGEN_STRONG_INLINE void unpackRhs(DenseIndex n, const RhsScalar* rhs, RhsScalar* b)
   {
     for(DenseIndex k=0; k<n; k++)
-      pstore(&b[k*RhsPacketSize], pset1<RhsPacket>(rhs[k]));
+      pstore1<RhsPacket>(&b[k*RhsPacketSize], rhs[k]);
   }
 
   EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacket& dest) const
@@ -1009,12 +1009,7 @@ EIGEN_ASM_COMMENT("mybegin4");
     for(Index j2=packet_cols; j2<cols; j2++)
     {
       // unpack B
-      {
-        traits.unpackRhs(depth, &blockB[j2*strideB+offsetB], unpackedB);
-//         const RhsScalar* blB = &blockB[j2*strideB+offsetB];
-//         for(Index k=0; k<depth; k++)
-//           pstore(&unpackedB[k*RhsPacketSize], pset1<RhsPacket>(blB[k]));
-      }
+      traits.unpackRhs(depth, &blockB[j2*strideB+offsetB], unpackedB);
 
       for(Index i=0; i<peeled_mc; i+=mr)
       {

Eigen/src/Core/products/GeneralMatrixMatrix.h

@@ -94,8 +94,9 @@ 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);
 
@@ -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,9 +165,10 @@ 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)
@@ -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;
@@ -125,8 +125,6 @@ struct general_matrix_matrix_triangular_product<Index,LhsScalar,LhsStorageOrder,
         }
       }
     }
-    ei_aligned_stack_delete(LhsScalar, blockA, kc*mc);
-    ei_aligned_stack_delete(RhsScalar, allocatedBlockB, sizeB);
   }
 };
 

Eigen/src/Core/products/SelfadjointMatrixMatrix.h

@@ -263,10 +263,10 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,ConjugateLhs
     // kc must smaller than 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;
@@ -313,9 +313,6 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,true,ConjugateLhs
         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;
@@ -369,9 +365,6 @@ struct product_selfadjoint_matrix<Scalar,Index,LhsStorageOrder,false,ConjugateLh
         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,14 +62,12 @@ 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;
@@ -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;
@@ -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,7 +200,7 @@ 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);
@@ -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;
@@ -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;
@@ -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;
@@ -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;
@@ -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;
@@ -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/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/DisableStupidWarnings.h

@@ -24,10 +24,12 @@
   // 2536 - type qualifiers are meaningless here
   //        ICC 12 generates this warning when a function return type is const qualified, even if that type is a template-parameter-dependent
   //        typedef that may be a reference type.
+  // 279  - controlling expression is constant
+  //        ICC 12 generates this warning on assert(constant_expression_depending_on_template_params) and frankly this is a legitimate use case.
   #ifndef EIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
     #pragma warning push
   #endif
-  #pragma warning disable 2196 2536
+  #pragma warning disable 2196 2536 279
 #elif defined __clang__
   // -Wconstant-logical-operand - warning: use of logical && with constant operand; switch to bitwise & or remove constant
   //     this is really a stupid warning as it warns on compile-time expressions involving enums

Eigen/src/Core/util/Macros.h

@@ -26,9 +26,9 @@
 #ifndef EIGEN_MACROS_H
 #define EIGEN_MACROS_H
 
-#define EIGEN_WORLD_VERSION 2
-#define EIGEN_MAJOR_VERSION 94
-#define EIGEN_MINOR_VERSION 0
+#define EIGEN_WORLD_VERSION 3
+#define EIGEN_MAJOR_VERSION 0
+#define EIGEN_MINOR_VERSION 1
 
 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
                                       (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \
@@ -261,9 +261,7 @@
  * If we made alignment depend on whether or not EIGEN_VECTORIZE is defined, it would be impossible to link
  * vectorized and non-vectorized code.
  */
-#if !EIGEN_ALIGN_STATICALLY
-  #define EIGEN_ALIGN_TO_BOUNDARY(n)
-#elif (defined __GNUC__) || (defined __PGI) || (defined __IBMCPP__)
+#if (defined __GNUC__) || (defined __PGI) || (defined __IBMCPP__)
   #define EIGEN_ALIGN_TO_BOUNDARY(n) __attribute__((aligned(n)))
 #elif (defined _MSC_VER)
   #define EIGEN_ALIGN_TO_BOUNDARY(n) __declspec(align(n))
@@ -276,6 +274,14 @@
 
 #define EIGEN_ALIGN16 EIGEN_ALIGN_TO_BOUNDARY(16)
 
+#if EIGEN_ALIGN_STATICALLY
+#define EIGEN_USER_ALIGN_TO_BOUNDARY(n) EIGEN_ALIGN_TO_BOUNDARY(n)
+#define EIGEN_USER_ALIGN16 EIGEN_ALIGN16
+#else
+#define EIGEN_USER_ALIGN_TO_BOUNDARY(n)
+#define EIGEN_USER_ALIGN16
+#endif
+
 #ifdef EIGEN_DONT_USE_RESTRICT_KEYWORD
   #define EIGEN_RESTRICT
 #endif

Eigen/src/Core/util/Memory.h

@@ -167,14 +167,36 @@ inline void* generic_aligned_realloc(void* ptr, size_t size, size_t old_size)
 *** Implementation of portable aligned versions of malloc/free/realloc     ***
 *****************************************************************************/
 
+#ifdef EIGEN_NO_MALLOC
+inline void check_that_malloc_is_allowed()
+{
+  eigen_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
+}
+#elif defined EIGEN_RUNTIME_NO_MALLOC
+inline bool is_malloc_allowed_impl(bool update, bool new_value = false)
+{
+  static bool value = true;
+  if (update == 1)
+    value = new_value;
+  return value;
+}
+inline bool is_malloc_allowed() { return is_malloc_allowed_impl(false); }
+inline bool set_is_malloc_allowed(bool new_value) { return is_malloc_allowed_impl(true, new_value); }
+inline void check_that_malloc_is_allowed()
+{
+  eigen_assert(is_malloc_allowed() && "heap allocation is forbidden (EIGEN_RUNTIME_NO_MALLOC is defined and g_is_malloc_allowed is false)");
+}
+#else 
+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.
   */
 inline void* aligned_malloc(size_t size)
 {
-  #ifdef EIGEN_NO_MALLOC
-    eigen_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
-  #endif
+  check_that_malloc_is_allowed();
 
   void *result;
   #if !EIGEN_ALIGN
@@ -268,9 +290,7 @@ template<bool Align> inline void* conditional_aligned_malloc(size_t size)
 
 template<> inline void* conditional_aligned_malloc<false>(size_t size)
 {
-  #ifdef EIGEN_NO_MALLOC
-    eigen_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
-  #endif
+  check_that_malloc_is_allowed();
 
   void *result = std::malloc(size);
   #ifdef EIGEN_EXCEPTIONS
@@ -448,36 +468,87 @@ 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)
+        Eigen::internal::construct_elements_of_array(m_ptr, size);
+    }
+    ~aligned_stack_memory_handler()
+    {
+      if(NumTraits<T>::RequireInitialization)
+        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) \
+    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) \
+    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
 
 
 /*****************************************************************************

Eigen/src/Eigen2Support/CMakeLists.txt

@@ -4,3 +4,5 @@ INSTALL(FILES
   ${Eigen_Eigen2Support_SRCS}
   DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigen2Support COMPONENT Devel
   )
+
+ADD_SUBDIRECTORY(Geometry)

Eigen/src/Eigen2Support/Geometry/CMakeLists.txt

@@ -1,6 +1,6 @@
-FILE(GLOB Eigen_Geometry_SRCS "*.h")
+FILE(GLOB Eigen_Eigen2Support_Geometry_SRCS "*.h")
 
 INSTALL(FILES
-  ${Eigen_Geometry_SRCS}
+  ${Eigen_Eigen2Support_Geometry_SRCS}
   DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigen2Support/Geometry
   )

Eigen/src/Eigenvalues/EigenSolver.h

@@ -343,6 +343,7 @@ typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eige
     {
       // we have a real eigen value
       matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
+      matV.col(j).normalize();
     }
     else
     {
@@ -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,13 +563,18 @@ 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;
 
         }
       }
     }
+    else
+    {
+      eigen_assert("Internal bug in EigenSolver"); // this should not happen
+    }
   }
 
   // Back transformation to get eigenvectors of original matrix

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

@@ -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))

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().
       *

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

@@ -232,13 +232,15 @@ template<typename MatrixType,typename Lhs>
 struct traits<homogeneous_left_product_impl<Homogeneous<MatrixType,Vertical>,Lhs> >
 {
   typedef typename take_matrix_for_product<Lhs>::type LhsMatrixType;
+  typedef typename remove_all<MatrixType>::type MatrixTypeCleaned;
+  typedef typename remove_all<LhsMatrixType>::type LhsMatrixTypeCleaned;
   typedef typename make_proper_matrix_type<
-                 typename traits<MatrixType>::Scalar,
-                 LhsMatrixType::RowsAtCompileTime,
-                 MatrixType::ColsAtCompileTime,
-                 MatrixType::PlainObject::Options,
-                 LhsMatrixType::MaxRowsAtCompileTime,
-                 MatrixType::MaxColsAtCompileTime>::type ReturnType;
+                 typename traits<MatrixTypeCleaned>::Scalar,
+                 LhsMatrixTypeCleaned::RowsAtCompileTime,
+                 MatrixTypeCleaned::ColsAtCompileTime,
+                 MatrixTypeCleaned::PlainObject::Options,
+                 LhsMatrixTypeCleaned::MaxRowsAtCompileTime,
+                 MatrixTypeCleaned::MaxColsAtCompileTime>::type ReturnType;
 };
 
 template<typename MatrixType,typename Lhs>
@@ -246,7 +248,8 @@ struct homogeneous_left_product_impl<Homogeneous<MatrixType,Vertical>,Lhs>
   : public ReturnByValue<homogeneous_left_product_impl<Homogeneous<MatrixType,Vertical>,Lhs> >
 {
   typedef typename traits<homogeneous_left_product_impl>::LhsMatrixType LhsMatrixType;
-  typedef typename remove_all<typename LhsMatrixType::Nested>::type LhsMatrixTypeNested;
+  typedef typename remove_all<LhsMatrixType>::type LhsMatrixTypeCleaned;
+  typedef typename remove_all<typename LhsMatrixTypeCleaned::Nested>::type LhsMatrixTypeNested;
   typedef typename MatrixType::Index Index;
   homogeneous_left_product_impl(const Lhs& lhs, const MatrixType& rhs)
     : m_lhs(take_matrix_for_product<Lhs>::run(lhs)),
@@ -267,7 +270,7 @@ struct homogeneous_left_product_impl<Homogeneous<MatrixType,Vertical>,Lhs>
             .template replicate<MatrixType::ColsAtCompileTime>(m_rhs.cols());
   }
 
-  const typename LhsMatrixType::Nested m_lhs;
+  const typename LhsMatrixTypeCleaned::Nested m_lhs;
   const typename MatrixType::Nested m_rhs;
 };
 

Eigen/src/Geometry/Hyperplane.h

@@ -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)
@@ -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/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 */
@@ -222,7 +225,8 @@ struct traits<Quaternion<_Scalar,_Options> >
   typedef _Scalar Scalar;
   typedef Matrix<_Scalar,4,1,_Options> Coefficients;
   enum{
-    PacketAccess = _Options & DontAlign ? Unaligned : Aligned
+    IsAligned = bool(EIGEN_ALIGN) && ((int(_Options)&Aligned)==Aligned),
+    Flags = IsAligned ? (AlignedBit | LvalueBit) : LvalueBit
   };
 };
 }
@@ -294,33 +298,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 +357,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 +369,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 +393,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 +423,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 +447,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 +575,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 +590,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 +650,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 +665,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);
@@ -654,7 +681,7 @@ 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;

Eigen/src/Geometry/Transform.h

@@ -43,7 +43,9 @@ struct transform_traits
 
 template< typename TransformType,
           typename MatrixType,
-          bool IsProjective = transform_traits<TransformType>::IsProjective>
+          int Case = transform_traits<TransformType>::IsProjective ? 0
+                   : int(MatrixType::RowsAtCompileTime) == int(transform_traits<TransformType>::HDim) ? 1
+                   : 2>
 struct transform_right_product_impl;
 
 template< typename Other,
@@ -84,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.
   *
@@ -199,7 +201,7 @@ public:
   typedef _Scalar Scalar;
   typedef DenseIndex Index;
   /** type of the matrix used to represent the transformation */
-  typedef Matrix<Scalar,Rows,HDim,Options> MatrixType;
+  typedef typename internal::make_proper_matrix_type<Scalar,Rows,HDim,Options>::type MatrixType;
   /** constified MatrixType */
   typedef const MatrixType ConstMatrixType;
   /** type of the matrix used to represent the linear part of the transformation */
@@ -670,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)
@@ -716,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;
 }
 
@@ -730,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
 
@@ -1080,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.
   *
@@ -1204,7 +1215,7 @@ struct transform_product_result
 };
 
 template< typename TransformType, typename MatrixType >
-struct transform_right_product_impl< TransformType, MatrixType, true >
+struct transform_right_product_impl< TransformType, MatrixType, 0 >
 {
   typedef typename MatrixType::PlainObject ResultType;
 
@@ -1215,7 +1226,7 @@ struct transform_right_product_impl< TransformType, MatrixType, true >
 };
 
 template< typename TransformType, typename MatrixType >
-struct transform_right_product_impl< TransformType, MatrixType, false >
+struct transform_right_product_impl< TransformType, MatrixType, 1 >
 {
   enum { 
     Dim = TransformType::Dim, 
@@ -1228,20 +1239,39 @@ struct transform_right_product_impl< TransformType, MatrixType, false >
 
   EIGEN_STRONG_INLINE static ResultType run(const TransformType& T, const MatrixType& other)
   {
-    EIGEN_STATIC_ASSERT(OtherRows==Dim || OtherRows==HDim, YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES);
+    EIGEN_STATIC_ASSERT(OtherRows==HDim, YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES);
 
     typedef Block<ResultType, Dim, OtherCols> TopLeftLhs;
-    typedef Block<const MatrixType, Dim, OtherCols> TopLeftRhs;
 
     ResultType res(other.rows(),other.cols());
+    TopLeftLhs(res, 0, 0, Dim, other.cols()).noalias() = T.affine() * other;
+    res.row(OtherRows-1) = other.row(OtherRows-1);
+    
+    return res;
+  }
+};
 
-    TopLeftLhs(res, 0, 0, Dim, other.cols()) =
-      ( T.linear() * TopLeftRhs(other, 0, 0, Dim, other.cols()) ).colwise() +
-        T.translation();
+template< typename TransformType, typename MatrixType >
+struct transform_right_product_impl< TransformType, MatrixType, 2 >
+{
+  enum { 
+    Dim = TransformType::Dim, 
+    HDim = TransformType::HDim,
+    OtherRows = MatrixType::RowsAtCompileTime,
+    OtherCols = MatrixType::ColsAtCompileTime
+  };
 
-    // we need to take .rows() because OtherRows might be Dim or HDim
-    if (OtherRows==HDim)
-      res.row(other.rows()) = other.row(other.rows());
+  typedef typename MatrixType::PlainObject ResultType;
+
+  EIGEN_STRONG_INLINE static ResultType run(const TransformType& T, const MatrixType& other)
+  {
+    EIGEN_STATIC_ASSERT(OtherRows==Dim, YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES);
+
+    typedef Block<ResultType, Dim, OtherCols> TopLeftLhs;
+
+    ResultType res(other.rows(),other.cols());
+    TopLeftLhs(res, 0, 0, Dim, other.cols()).noalias() = T.linear() * other;
+    TopLeftLhs(res, 0, 0, Dim, other.cols()).colwise() += T.translation();
 
     return res;
   }

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

@@ -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)
         {

Eigen/src/QR/ColPivHouseholderQR.h

@@ -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/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
@@ -376,7 +376,12 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
       * The default constructor is useful in cases in which the user intends to
       * perform decompositions via JacobiSVD::compute(const MatrixType&).
       */
-    JacobiSVD() : m_isInitialized(false) {}
+    JacobiSVD()
+      : m_isInitialized(false),
+        m_isAllocated(false),
+        m_computationOptions(0),
+        m_rows(-1), m_cols(-1)
+    {}
 
 
     /** \brief Default Constructor with memory preallocation
@@ -386,6 +391,10 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
       * \sa JacobiSVD()
       */
     JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
+      : m_isInitialized(false),
+        m_isAllocated(false),
+        m_computationOptions(0),
+        m_rows(-1), m_cols(-1)
     {
       allocate(rows, cols, computationOptions);
     }
@@ -394,33 +403,48 @@ 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.
      */
     JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
+      : m_isInitialized(false),
+        m_isAllocated(false),
+        m_computationOptions(0),
+        m_rows(-1), m_cols(-1)
     {
       compute(matrix, computationOptions);
     }
 
-    /** \brief Method performing the decomposition of given matrix.
+    /** \brief Method performing the decomposition of given matrix using custom options.
      *
      * \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.
      */
-    JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions = 0);
+    JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions);
+
+    /** \brief Method performing the decomposition of given matrix using current options.
+     *
+     * \param matrix the matrix to decompose
+     *
+     * This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int).
+     */
+    JacobiSVD& compute(const MatrixType& matrix)
+    {
+      return compute(matrix, m_computationOptions);
+    }
 
     /** \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.
      *
@@ -436,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.
      *
@@ -494,16 +518,17 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
     inline Index cols() const { return m_cols; }
 
   private:
-    void allocate(Index rows, Index cols, unsigned int computationOptions = 0);
+    void allocate(Index rows, Index cols, unsigned int computationOptions);
 
   protected:
     MatrixUType m_matrixU;
     MatrixVType m_matrixV;
     SingularValuesType m_singularValues;
     WorkMatrixType m_workMatrix;
-    bool m_isInitialized;
+    bool m_isInitialized, m_isAllocated;
     bool m_computeFullU, m_computeThinU;
     bool m_computeFullV, m_computeThinV;
+    unsigned int m_computationOptions;
     Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
 
     template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
@@ -515,9 +540,21 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
 template<typename MatrixType, int QRPreconditioner>
 void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Index rows, Index cols, unsigned int computationOptions)
 {
+  eigen_assert(rows >= 0 && cols >= 0);
+
+  if (m_isAllocated &&
+      rows == m_rows &&
+      cols == m_cols &&
+      computationOptions == m_computationOptions)
+  {
+    return;
+  }
+
   m_rows = rows;
   m_cols = cols;
   m_isInitialized = false;
+  m_isAllocated = true;
+  m_computationOptions = computationOptions;
   m_computeFullU = (computationOptions & ComputeFullU) != 0;
   m_computeThinU = (computationOptions & ComputeThinU) != 0;
   m_computeFullV = (computationOptions & ComputeFullV) != 0;
@@ -581,8 +618,9 @@ 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)
+        using std::max;
+        if(max(internal::abs(m_workMatrix.coeff(p,q)),internal::abs(m_workMatrix.coeff(q,p)))
+            > max(internal::abs(m_workMatrix.coeff(p,p)),internal::abs(m_workMatrix.coeff(q,q)))*precision)
         {
           finished = false;
 

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;

Eigen/src/Sparse/SparseSparseProduct.h

@@ -133,7 +133,12 @@ static void sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
   float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
   float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f);
 
-  res.resize(rows, cols);
+  // mimics a resizeByInnerOuter:
+  if(ResultType::IsRowMajor)
+    res.resize(cols, rows);
+  else
+    res.resize(rows, cols);
+
   res.reserve(Index(ratioRes*rows*cols));
   for (Index j=0; j<cols; ++j)
   {

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;

debug/gdb/printers.py

@@ -56,20 +56,22 @@ class EigenMatrixPrinter:
 		template_params = m.split(',')
 		template_params = map(lambda x:x.replace(" ", ""), template_params)
 
-		self.rows = int(template_params[1])
-		self.cols = int(template_params[2])
+		if template_params[1] == '-0x00000000000000001':
+			self.rows = val['m_storage']['m_rows']
+		else:
+			self.rows = int(template_params[1])
+		
+		if template_params[2] == '-0x00000000000000001':
+			self.cols = val['m_storage']['m_cols']
+		else:
+			self.cols = int(template_params[2])
+		
 		self.options = 0 # default value
 		if len(template_params) > 3:
 			self.options = template_params[3];
 		
 		self.rowMajor = (int(self.options) & 0x1)
 
-		if self.rows == 10000:
-			self.rows = val['m_storage']['m_rows']
-
-		if self.cols == 10000:
-			self.cols = val['m_storage']['m_cols']
-
 		self.innerType = self.type.template_argument(0)
 
 		self.val = val

doc/A05_PortingFrom2To3.dox

@@ -21,6 +21,7 @@ and gives tips to help porting your application from Eigen2 to Eigen3.
   - \ref LazyVsNoalias
   - \ref AlignMacros
   - \ref AlignedMap
+  - \ref StdContainers
   - \ref eiPrefix
 
 \section CompatibilitySupport Eigen2 compatibility support
@@ -286,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
@@ -295,6 +296,18 @@ There also are related convenience static methods, which actually are the prefer
 result = Vector4f::MapAligned(some_aligned_array);
 \endcode
 
+\section StdContainers STL Containers
+
+In Eigen2, #include<Eigen/StdVector> tweaked std::vector to automatically align elements. The problem was that that was quite invasive. In Eigen3, we only override standard behavior if you use Eigen::aligned_allocator<T> as your allocator type. So for example, if you use std::vector<Matrix4f>, you need to do the following change (note that aligned_allocator is under namespace Eigen):
+
+<table class="manual">
+<tr><th>Eigen 2</th><th>Eigen 3</th></tr>
+<tr>
+  <td> \code std::vector<Matrix4f> \endcode </td>
+  <td> \code std::vector<Matrix4f, aligned_allocator<Matrix4f> > \endcode </td>
+</tr>
+</table>
+
 \section eiPrefix Internal ei_ prefix
 
 In Eigen2, global internal functions and structures were prefixed by \c ei_. In Eigen3, they all have been moved into the more explicit \c internal namespace. So, e.g., \c ei_sqrt(x) now becomes \c internal::sqrt(x). Of course it is not recommended to rely on Eigen's internal features.

doc/A10_Eigen2SupportModes.dox

@@ -47,6 +47,7 @@ The parts of the API that are still not 100% Eigen 2 compatible in this mode are
 \li Dot products over complex numbers. Eigen 2's dot product was linear in the first variable. Eigen 3's dot product is linear in the second variable. In other words, the Eigen 2 code \code x.dot(y) \endcode is equivalent to the Eigen 3 code \code y.dot(x) \endcode In yet other words, dot products are complex-conjugated in Eigen 3 compared to Eigen 2. The switch to the new convention was commanded by common usage, especially with the notation \f$ x^Ty \f$ for dot products of column-vectors.
 \li The Sparse module.
 \li Certain fine details of linear algebraic decompositions. For example, LDLT decomposition is now pivoting in Eigen 3 whereas it wasn't in Eigen 2, so code that was relying on its underlying matrix structure will break.
+\li Usage of Eigen types in STL containers, \ref Eigen2ToEigen3 "as explained on this page".
 
 \section Stage20 Stage 20: define EIGEN2_SUPPORT_STAGE20_RESOLVE_API_CONFLICTS
 

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

doc/C08_TutorialGeometry.dox

@@ -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)
@@ -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>
 

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/I14_PreprocessorDirectives.dox

@@ -2,7 +2,7 @@ namespace Eigen {
 
 /** \page TopicPreprocessorDirectives Preprocessor directives
 
-You can control some aspects of Eigen by defining the preprocessor tokens using \c #define. These macros
+You can control some aspects of Eigen by defining the preprocessor tokens using \c \#define. These macros
 should be defined before any Eigen headers are included. Often they are best set in the project options.
 
 This page lists the preprocesor tokens recognised by Eigen.
@@ -22,6 +22,8 @@ This page lists the preprocesor tokens recognised by Eigen.
    Eigen3; see \ref Eigen2SupportModes.
  - \b EIGEN_DEFAULT_TO_ROW_MAJOR - when defined, the default storage order for matrices becomes row-major
    instead of column-major. Not defined by default.
+ - \b EIGEN_DEFAULT_DENSE_INDEX_TYPE - the type for column and row indices in matrices, vectors and array
+   (DenseBase::Index). Set to \c std::ptrdiff_t by default.
  - \b EIGEN_DEFAULT_IO_FORMAT - the IOFormat to use when printing a matrix if no #IOFormat is specified.
    Defaults to the #IOFormat constructed by the default constructor IOFormat().
  - \b EIGEN_INITIALIZE_MATRICES_BY_ZERO - if defined, all entries of newly constructed matrices and arrays are
@@ -43,7 +45,8 @@ run time. However, these assertions do cost time and can thus be turned off.
 
 \section TopicPreprocessorDirectivesPerformance Alignment, vectorization and performance tweaking
 
- - \b EIGEN_DONT_ALIGN - disables alignment.
+ - \b EIGEN_DONT_ALIGN - disables alignment completely. Eigen will not try to align its objects and does not
+   expect that any objects passed to it are aligned. This will turn off vectorization. Not defined by default.
  - \b EIGEN_DONT_ALIGN_STATICALLY - disables alignment of arrays on the stack. Not defined by default, unless
    \c EIGEN_DONT_ALIGN is defined.
  - \b EIGEN_DONT_VECTORIZE - disables explicit vectorization when defined. Not defined by default, unless 

doc/QuickReference.dox

@@ -535,7 +535,7 @@ vec.reverse()           mat.colwise().reverse()   mat.rowwise().reverse()
 vec.reverseInPlace()
 \endcode
 
-\subsection QuickRef_Reverse Replicate
+\subsection QuickRef_Replicate Replicate
 Vectors, matrices, rows, and/or columns can be replicated in any direction (see DenseBase::replicate(), VectorwiseOp::replicate())
 \code
 vec.replicate(times)                                          vec.replicate<Times>
@@ -594,7 +594,7 @@ unit or null diagonal (read/write):
 </td><td>\code
 m.triangularView<Xxx>()
 \endcode \n
-\c Xxx = Upper, Lower, StrictlyUpper, StrictlyLower, UnitUpper, UnitLower
+\c Xxx = ::Upper, ::Lower, ::StrictlyUpper, ::StrictlyLower, ::UnitUpper, ::UnitLower
 </td></tr>
 <tr><td>
 Writing to a specific triangular part:\n (only the referenced triangular part is evaluated)
@@ -645,14 +645,14 @@ m3 -= s1 * m3.adjoint() * m1.selfadjointView<Eigen::Lower>();\endcode
 </td></tr>
 <tr><td>
 Rank 1 and rank K update: \n
-\f$ upper(M_1) += s1 M_2^* M_2 \f$ \n
-\f$ lower(M_1) -= M_2 M_2^* \f$
+\f$ upper(M_1) \mathrel{{+}{=}} s_1 M_2^* M_2 \f$ \n
+\f$ lower(M_1) \mathbin{{-}{=}} M_2 M_2^* \f$
 </td><td>\n \code
 M1.selfadjointView<Eigen::Upper>().rankUpdate(M2,s1);
 m1.selfadjointView<Eigen::Lower>().rankUpdate(m2.adjoint(),-1); \endcode
 </td></tr>
 <tr><td>
-Rank 2 update: (\f$ M += s u v^* + s v u^* \f$)
+Rank 2 update: (\f$ M \mathrel{{+}{=}} s u v^* + s v u^* \f$)
 </td><td>\code
 M.selfadjointView<Eigen::Upper>().rankUpdate(u,v,s);
 \endcode

scripts/eigen_gen_docs

@@ -5,15 +5,17 @@
 # will be used
 USER=${USER:-'orzel'}
 
+#ulimit -v 1024000
+
 # step 1 : build
 # todo if 'build is not there, create one:
-#mkdir build
-(cd build && cmake .. && make -j3 doc) || { echo "make failed"; exit 1; }
+mkdir build -p
+(cd build && cmake .. && make doc) || { echo "make failed"; exit 1; }
 #todo: n+1 where n = number of cpus
 
 #step 2 : upload
 # (the '/' at the end of path are very important, see rsync documentation)
-rsync -az build/doc/html/ $USER@ssh.tuxfamily.org:eigen/eigen.tuxfamily.org-web/htdocs/dox-devel/ || { echo "upload failed"; exit 1; }
+rsync -az build/doc/html/ $USER@ssh.tuxfamily.org:eigen/eigen.tuxfamily.org-web/htdocs/dox-3.0/ || { echo "upload failed"; exit 1; }
 
 echo "Uploaded successfully"
 

test/CMakeLists.txt

@@ -42,6 +42,7 @@ ei_add_test(linearstructure)
 ei_add_test(integer_types)
 ei_add_test(cwiseop)
 ei_add_test(unalignedcount)
+ei_add_test(exceptions)
 ei_add_test(redux)
 ei_add_test(visitor)
 ei_add_test(block)
@@ -119,6 +120,7 @@ ei_add_test(eigen2support)
 ei_add_test(nullary)
 ei_add_test(nesting_ops "${CMAKE_CXX_FLAGS_DEBUG}")
 ei_add_test(zerosized)
+ei_add_test(dontalign)
 
 ei_add_test(prec_inverse_4x4)
 

test/array_for_matrix.cpp

@@ -57,10 +57,10 @@ template<typename MatrixType> void array_for_matrix(const MatrixType& m)
   VERIFY_IS_APPROX(m3, (m1.array() - s1).matrix());
 
   // reductions
-  VERIFY_IS_APPROX(m1.colwise().sum().sum(), m1.sum());
-  VERIFY_IS_APPROX(m1.rowwise().sum().sum(), m1.sum());
-  if (!internal::isApprox(m1.sum(), (m1+m2).sum()))
-    VERIFY_IS_NOT_APPROX(((m1+m2).rowwise().sum()).sum(), m1.sum());
+  VERIFY_IS_MUCH_SMALLER_THAN(m1.colwise().sum().sum() - m1.sum(), m1.cwiseAbs().maxCoeff());
+  VERIFY_IS_MUCH_SMALLER_THAN(m1.rowwise().sum().sum() - m1.sum(), m1.cwiseAbs().maxCoeff());
+  VERIFY_IS_MUCH_SMALLER_THAN(m1.colwise().sum() + m2.colwise().sum() - (m1+m2).colwise().sum(), (m1+m2).cwiseAbs().maxCoeff());
+  VERIFY_IS_MUCH_SMALLER_THAN(m1.rowwise().sum() - m2.rowwise().sum() - (m1-m2).rowwise().sum(), (m1-m2).cwiseAbs().maxCoeff());
   VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op<Scalar>()));
 
   // vector-wise ops
@@ -158,27 +158,28 @@ template<typename VectorType> void lpNorm(const VectorType& v)
 
 void test_array_for_matrix()
 {
+  int maxsize = 40;
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1( array_for_matrix(Matrix<float, 1, 1>()) );
     CALL_SUBTEST_2( array_for_matrix(Matrix2f()) );
     CALL_SUBTEST_3( array_for_matrix(Matrix4d()) );
-    CALL_SUBTEST_4( array_for_matrix(MatrixXcf(3, 3)) );
-    CALL_SUBTEST_5( array_for_matrix(MatrixXf(8, 12)) );
-    CALL_SUBTEST_6( array_for_matrix(MatrixXi(8, 12)) );
+    CALL_SUBTEST_4( array_for_matrix(MatrixXcf(internal::random<int>(1,maxsize), internal::random<int>(1,maxsize))) );
+    CALL_SUBTEST_5( array_for_matrix(MatrixXf(internal::random<int>(1,maxsize), internal::random<int>(1,maxsize))) );
+    CALL_SUBTEST_6( array_for_matrix(MatrixXi(internal::random<int>(1,maxsize), internal::random<int>(1,maxsize))) );
   }
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1( comparisons(Matrix<float, 1, 1>()) );
     CALL_SUBTEST_2( comparisons(Matrix2f()) );
     CALL_SUBTEST_3( comparisons(Matrix4d()) );
-    CALL_SUBTEST_5( comparisons(MatrixXf(8, 12)) );
-    CALL_SUBTEST_6( comparisons(MatrixXi(8, 12)) );
+    CALL_SUBTEST_5( comparisons(MatrixXf(internal::random<int>(1,maxsize), internal::random<int>(1,maxsize))) );
+    CALL_SUBTEST_6( comparisons(MatrixXi(internal::random<int>(1,maxsize), internal::random<int>(1,maxsize))) );
   }
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1( lpNorm(Matrix<float, 1, 1>()) );
     CALL_SUBTEST_2( lpNorm(Vector2f()) );
     CALL_SUBTEST_7( lpNorm(Vector3d()) );
     CALL_SUBTEST_8( lpNorm(Vector4f()) );
-    CALL_SUBTEST_5( lpNorm(VectorXf(16)) );
-    CALL_SUBTEST_4( lpNorm(VectorXcf(10)) );
+    CALL_SUBTEST_5( lpNorm(VectorXf(internal::random<int>(1,maxsize))) );
+    CALL_SUBTEST_4( lpNorm(VectorXcf(internal::random<int>(1,maxsize))) );
   }
 }

test/dontalign.cpp

@@ -0,0 +1,78 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2011 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#if defined EIGEN_TEST_PART_1 || defined EIGEN_TEST_PART_2 || defined EIGEN_TEST_PART_3 || defined EIGEN_TEST_PART_4
+#define EIGEN_DONT_ALIGN
+#elif defined EIGEN_TEST_PART_5 || defined EIGEN_TEST_PART_6 || defined EIGEN_TEST_PART_7 || defined EIGEN_TEST_PART_8
+#define EIGEN_DONT_ALIGN_STATICALLY
+#endif
+
+#include "main.h"
+#include <Eigen/Dense>
+
+template<typename MatrixType>
+void dontalign(const MatrixType& m)
+{
+  typedef typename MatrixType::Index Index;
+  typedef typename MatrixType::Scalar Scalar;
+  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
+  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
+
+  Index rows = m.rows();
+  Index cols = m.cols();
+
+  MatrixType a = MatrixType::Random(rows,cols);
+  SquareMatrixType square = SquareMatrixType::Random(rows,rows);
+  VectorType v = VectorType::Random(rows);
+
+  VERIFY_IS_APPROX(v, square * square.colPivHouseholderQr().solve(v));
+  square = square.inverse().eval();
+  a = square * a;
+  square = square*square;
+  v = square * v;
+  v = a.adjoint() * v;
+  VERIFY(square.determinant() != Scalar(0));
+
+  // bug 219: MapAligned() was giving an assert with EIGEN_DONT_ALIGN, because Map Flags were miscomputed
+  Scalar* array = internal::aligned_new<Scalar>(rows);
+  v = VectorType::MapAligned(array, rows);
+  internal::aligned_delete(array, rows);
+}
+
+void test_dontalign()
+{
+#if defined EIGEN_TEST_PART_1 || defined EIGEN_TEST_PART_5
+  dontalign(Matrix3d());
+  dontalign(Matrix4f());
+#elif defined EIGEN_TEST_PART_2 || defined EIGEN_TEST_PART_6
+  dontalign(Matrix3cd());
+  dontalign(Matrix4cf());
+#elif defined EIGEN_TEST_PART_3 || defined EIGEN_TEST_PART_7
+  dontalign(Matrix<float, 32, 32>());
+  dontalign(Matrix<std::complex<float>, 32, 32>());
+#elif defined EIGEN_TEST_PART_4 || defined EIGEN_TEST_PART_8
+  dontalign(MatrixXd(32, 32));
+  dontalign(MatrixXcf(32, 32));
+#endif
+}

test/dynalloc.cpp

@@ -70,11 +70,10 @@ void check_aligned_stack_alloc()
 {
   for(int i = 1; i < 1000; i++)
   {
-    float *p = ei_aligned_stack_new(float,i);
+    ei_declare_aligned_stack_constructed_variable(float,p,i,0);
     VERIFY(size_t(p)%ALIGNMENT==0);
     // if the buffer is wrongly allocated this will give a bad write --> check with valgrind
     for(int j = 0; j < i; j++) p[j]=0;
-    ei_aligned_stack_delete(float,p,i);
   }
 }
 
@@ -120,11 +119,12 @@ void test_dynalloc()
   }
   
   // check static allocation, who knows ?
+  #if EIGEN_ALIGN_STATICALLY
   {
     MyStruct foo0;  VERIFY(size_t(foo0.avec.data())%ALIGNMENT==0);
     MyClassA fooA;  VERIFY(size_t(fooA.avec.data())%ALIGNMENT==0);
   }
-
+  
   // dynamic allocation, single object
   for (int i=0; i<g_repeat*100; ++i)
   {
@@ -143,5 +143,6 @@ void test_dynalloc()
     delete[] foo0;
     delete[] fooA;
   }
+  #endif
   
 }

test/eigen2/eigen2_adjoint.cpp

@@ -68,7 +68,7 @@ template<typename MatrixType> void adjoint(const MatrixType& m)
   VERIFY(ei_isApprox((s1 * v1 + s2 * v2).eigen2_dot(v3),   s1 * v1.eigen2_dot(v3) + s2 * v2.eigen2_dot(v3), largerEps));
   VERIFY(ei_isApprox(v3.eigen2_dot(s1 * v1 + s2 * v2),     ei_conj(s1)*v3.eigen2_dot(v1)+ei_conj(s2)*v3.eigen2_dot(v2), largerEps));
   VERIFY_IS_APPROX(ei_conj(v1.eigen2_dot(v2)),               v2.eigen2_dot(v1));
-  VERIFY_IS_APPROX(ei_abs(v1.eigen2_dot(v1)),                v1.squaredNorm());
+  VERIFY_IS_APPROX(ei_real(v1.eigen2_dot(v1)),               v1.squaredNorm());
   if(NumTraits<Scalar>::HasFloatingPoint)
     VERIFY_IS_APPROX(v1.squaredNorm(),                      v1.norm() * v1.norm());
   VERIFY_IS_MUCH_SMALLER_THAN(ei_abs(vzero.eigen2_dot(v1)),  static_cast<RealScalar>(1));

test/eigen2/eigen2_dynalloc.cpp

@@ -70,11 +70,10 @@ void check_aligned_stack_alloc()
 {
   for(int i = 1; i < 1000; i++)
   {
-    float *p = ei_aligned_stack_new(float,i);
+    ei_declare_aligned_stack_constructed_variable(float, p, i, 0);
     VERIFY(std::size_t(p)%ALIGNMENT==0);
     // if the buffer is wrongly allocated this will give a bad write --> check with valgrind
     for(int j = 0; j < i; j++) p[j]=0;
-    ei_aligned_stack_delete(float,p,i);
   }
 }
 

test/eigen2/eigen2_eigensolver.cpp

@@ -89,10 +89,9 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
     VERIFY((symmA * _evec).isApprox(symmB * (_evec * _eval.asDiagonal()), largerEps));
 
     // compare with eigen
-//     std::cerr << _eval.transpose() << "\n" << eiSymmGen.eigenvalues().transpose() << "\n\n";
-//     std::cerr << _evec.format(6) << "\n\n" << eiSymmGen.eigenvectors().format(6) << "\n\n\n";
+    MatrixType normalized_eivec = eiSymmGen.eigenvectors()*eiSymmGen.eigenvectors().colwise().norm().asDiagonal().inverse();
     VERIFY_IS_APPROX(_eval, eiSymmGen.eigenvalues());
-    VERIFY_IS_APPROX(_evec.cwise().abs(), eiSymmGen.eigenvectors().cwise().abs());
+    VERIFY_IS_APPROX(_evec.cwiseAbs(), normalized_eivec.cwiseAbs());
 
     Gsl::free(gSymmA);
     Gsl::free(gSymmB);

test/eigen2/eigen2_stdvector.cpp

@@ -32,7 +32,7 @@ void check_stdvector_matrix(const MatrixType& m)
   int rows = m.rows();
   int cols = m.cols();
   MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols);
-  std::vector<MatrixType> v(10, MatrixType(rows,cols)), w(20, y);
+  std::vector<MatrixType, aligned_allocator<MatrixType> > v(10, MatrixType(rows,cols)), w(20, y);
   v[5] = x;
   w[6] = v[5];
   VERIFY_IS_APPROX(w[6], v[5]);
@@ -67,7 +67,7 @@ void check_stdvector_transform(const TransformType&)
 {
   typedef typename TransformType::MatrixType MatrixType;
   TransformType x(MatrixType::Random()), y(MatrixType::Random());
-  std::vector<TransformType> v(10), w(20, y);
+  std::vector<TransformType, aligned_allocator<TransformType> > v(10), w(20, y);
   v[5] = x;
   w[6] = v[5];
   VERIFY_IS_APPROX(w[6], v[5]);
@@ -102,7 +102,7 @@ void check_stdvector_quaternion(const QuaternionType&)
 {
   typedef typename QuaternionType::Coefficients Coefficients;
   QuaternionType x(Coefficients::Random()), y(Coefficients::Random());
-  std::vector<QuaternionType> v(10), w(20, y);
+  std::vector<QuaternionType, aligned_allocator<QuaternionType> > v(10), w(20, y);
   v[5] = x;
   w[6] = v[5];
   VERIFY_IS_APPROX(w[6], v[5]);

test/eigen2/main.h

@@ -111,8 +111,10 @@ namespace Eigen
   #else // EIGEN_DEBUG_ASSERTS
 
     #undef eigen_assert
+
+    // see bug 89. The copy_bool here is working around a bug in gcc <= 4.3
     #define eigen_assert(a) \
-      if( (!(a)) && (!no_more_assert) )     \
+      if( (!Eigen::internal::copy_bool(a)) && (!no_more_assert) )	\
       {                                     \
         Eigen::no_more_assert = true;       \
         throw Eigen::eigen_assert_exception(); \

test/eigensolver_generic.cpp

@@ -61,6 +61,7 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
   VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
   VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(),
                    ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
+  VERIFY_IS_APPROX(ei1.eigenvectors().colwise().norm(), RealVectorType::Ones(rows).transpose());
   VERIFY_IS_APPROX(a.eigenvalues(), ei1.eigenvalues());
 
   EigenSolver<MatrixType> eiNoEivecs(a, false);

test/exceptions.cpp

@@ -0,0 +1,124 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+
+// Various sanity tests with exceptions:
+//  - no memory leak when a custom scalar type trow an exceptions
+//  - todo: complete the list of tests!
+
+#define EIGEN_STACK_ALLOCATION_LIMIT 100000000
+
+#include "main.h"
+
+struct my_exception
+{
+  my_exception() {}
+  ~my_exception() {}
+};
+    
+class ScalarWithExceptions
+{
+  public:
+    ScalarWithExceptions() { init(); }
+    ScalarWithExceptions(const float& _v) { init(); *v = _v; }
+    ScalarWithExceptions(const ScalarWithExceptions& other) { init(); *v = *(other.v); }
+    ~ScalarWithExceptions() {
+      delete v;
+      instances--;
+    }
+
+    void init() {
+      v = new float;
+      instances++;
+    }
+
+    ScalarWithExceptions operator+(const ScalarWithExceptions& other) const
+    {
+      countdown--;
+      if(countdown<=0)
+        throw my_exception();
+      return ScalarWithExceptions(*v+*other.v);
+    }
+    
+    ScalarWithExceptions operator-(const ScalarWithExceptions& other) const
+    { return ScalarWithExceptions(*v-*other.v); }
+    
+    ScalarWithExceptions operator*(const ScalarWithExceptions& other) const
+    { return ScalarWithExceptions((*v)*(*other.v)); }
+    
+    ScalarWithExceptions& operator+=(const ScalarWithExceptions& other)
+    { *v+=*other.v; return *this; }
+    ScalarWithExceptions& operator-=(const ScalarWithExceptions& other)
+    { *v-=*other.v; return *this; }
+    ScalarWithExceptions& operator=(const ScalarWithExceptions& other)
+    { *v = *(other.v); return *this; }
+  
+    bool operator==(const ScalarWithExceptions& other) const
+    { return *v==*other.v; }
+    bool operator!=(const ScalarWithExceptions& other) const
+    { return *v!=*other.v; }
+    
+    float* v;
+    static int instances;
+    static int countdown;
+};
+
+int ScalarWithExceptions::instances = 0;
+int ScalarWithExceptions::countdown = 0;
+
+
+#define CHECK_MEMLEAK(OP) {                                 \
+    ScalarWithExceptions::countdown = 100;                  \
+    int before = ScalarWithExceptions::instances;           \
+    bool exception_thrown = false;                         \
+    try { OP; }                              \
+    catch (my_exception) {                                  \
+      exception_thrown = true;                              \
+      VERIFY(ScalarWithExceptions::instances==before && "memory leak detected in " && EIGEN_MAKESTRING(OP)); \
+    } \
+    VERIFY(exception_thrown && " no exception thrown in " && EIGEN_MAKESTRING(OP)); \
+  }
+
+void memoryleak()
+{
+  typedef Eigen::Matrix<ScalarWithExceptions,Dynamic,1> VectorType;
+  typedef Eigen::Matrix<ScalarWithExceptions,Dynamic,Dynamic> MatrixType;
+  
+  {
+    int n = 50;
+    VectorType v0(n), v1(n);
+    MatrixType m0(n,n), m1(n,n), m2(n,n);
+    v0.setOnes(); v1.setOnes();
+    m0.setOnes(); m1.setOnes(); m2.setOnes();
+    CHECK_MEMLEAK(v0 = m0 * m1 * v1);
+    CHECK_MEMLEAK(m2 = m0 * m1 * m2);
+    CHECK_MEMLEAK((v0+v1).dot(v0+v1));
+  }
+  VERIFY(ScalarWithExceptions::instances==0 && "global memory leak detected in " && EIGEN_MAKESTRING(OP)); \
+}
+
+void test_exceptions()
+{
+  CALL_SUBTEST( memoryleak() );
+}

test/geo_homogeneous.cpp

@@ -87,15 +87,29 @@ template<typename Scalar,int Size> void homogeneous(void)
                     m0.transpose().rowwise().homogeneous() * t3);
 
   // test product with a Transform object
-  Transform<Scalar, Size, Affine> Rt;
-  Matrix<Scalar, Size, Dynamic> pts, Rt_pts1;
+  Transform<Scalar, Size, Affine> aff;
+  Transform<Scalar, Size, AffineCompact> caff;
+  Transform<Scalar, Size, Projective> proj;
+  Matrix<Scalar, Size, Dynamic>   pts;
+  Matrix<Scalar, Size+1, Dynamic> pts1, pts2;
 
-  Rt.setIdentity();
-  pts.setRandom(Size,5);
-
-  Rt_pts1 = Rt * pts.colwise().homogeneous();
-  // std::cerr << (Rt_pts1 - pts).sum() << "\n";
-  VERIFY_IS_MUCH_SMALLER_THAN( (Rt_pts1 - pts).sum(), Scalar(1));
+  aff.affine().setRandom();
+  proj = caff = aff;
+  pts.setRandom(Size,internal::random<int>(1,20));
+  
+  pts1 = pts.colwise().homogeneous();
+  VERIFY_IS_APPROX(aff  * pts.colwise().homogeneous(), (aff  * pts1).colwise().hnormalized());
+  VERIFY_IS_APPROX(caff * pts.colwise().homogeneous(), (caff * pts1).colwise().hnormalized());
+  VERIFY_IS_APPROX(proj * pts.colwise().homogeneous(), (proj * pts1));
+  
+  VERIFY_IS_APPROX((aff  * pts1).colwise().hnormalized(),  aff  * pts);
+  VERIFY_IS_APPROX((caff * pts1).colwise().hnormalized(), caff * pts);
+  
+  pts2 = pts1;
+  pts2.row(Size).setRandom();
+  VERIFY_IS_APPROX((aff  * pts2).colwise().hnormalized(), aff  * pts2.colwise().hnormalized());
+  VERIFY_IS_APPROX((caff * pts2).colwise().hnormalized(), caff * pts2.colwise().hnormalized());
+  VERIFY_IS_APPROX((proj * pts2).colwise().hnormalized(), (proj * pts2.colwise().hnormalized().colwise().homogeneous()).colwise().hnormalized());
 }
 
 void test_geo_homogeneous()

test/geo_hyperplane.cpp

@@ -150,8 +150,9 @@ template<typename Scalar> void hyperplane_alignment()
   VERIFY_IS_APPROX(p1->coeffs(), p2->coeffs());
   VERIFY_IS_APPROX(p1->coeffs(), p3->coeffs());
   
-  #ifdef EIGEN_VECTORIZE
-  VERIFY_RAISES_ASSERT((::new(reinterpret_cast<void*>(array3u)) Plane3a));
+  #if defined(EIGEN_VECTORIZE) && EIGEN_ALIGN_STATICALLY
+  if(internal::packet_traits<Scalar>::Vectorizable)
+    VERIFY_RAISES_ASSERT((::new(reinterpret_cast<void*>(array3u)) Plane3a));
   #endif
 }
 

test/geo_parametrizedline.cpp

@@ -90,8 +90,9 @@ template<typename Scalar> void parametrizedline_alignment()
   VERIFY_IS_APPROX(p1->direction(), p2->direction());
   VERIFY_IS_APPROX(p1->direction(), p3->direction());
   
-  #ifdef EIGEN_VECTORIZE
-  VERIFY_RAISES_ASSERT((::new(reinterpret_cast<void*>(array3u)) Line4a));
+  #if defined(EIGEN_VECTORIZE) && EIGEN_ALIGN_STATICALLY
+  if(internal::packet_traits<Scalar>::Vectorizable)
+    VERIFY_RAISES_ASSERT((::new(reinterpret_cast<void*>(array3u)) Line4a));
   #endif
 }
 
@@ -100,6 +101,7 @@ void test_geo_parametrizedline()
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1( parametrizedline(ParametrizedLine<float,2>()) );
     CALL_SUBTEST_2( parametrizedline(ParametrizedLine<float,3>()) );
+    CALL_SUBTEST_2( parametrizedline_alignment<float>() );
     CALL_SUBTEST_3( parametrizedline(ParametrizedLine<double,4>()) );
     CALL_SUBTEST_3( parametrizedline_alignment<double>() );
     CALL_SUBTEST_4( parametrizedline(ParametrizedLine<std::complex<double>,5>()) );

test/geo_quaternion.cpp

@@ -72,7 +72,7 @@ template<typename Scalar, int Options> void quaternion(void)
 
   if((q1.coeffs()-q2.coeffs()).norm() > 10*largeEps)
   {
-    VERIFY(internal::isApprox(q1.angularDistance(q2), refangle, largeEps));
+    VERIFY_IS_MUCH_SMALLER_THAN(internal::abs(q1.angularDistance(q2) - refangle), Scalar(1));
   }
 
   // rotation matrix conversion
@@ -90,7 +90,15 @@ template<typename Scalar, int Options> void quaternion(void)
   // angle-axis conversion
   AngleAxisx aa = AngleAxisx(q1);
   VERIFY_IS_APPROX(q1 * v1, Quaternionx(aa) * v1);
-  VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1);
+
+  // Do not execute the test if the rotation angle is almost zero, or
+  // the rotation axis and v1 are almost parallel.
+  if (internal::abs(aa.angle()) > 5*test_precision<Scalar>()
+      && (aa.axis() - v1.normalized()).norm() < 1.99
+      && (aa.axis() + v1.normalized()).norm() < 1.99) 
+  {
+    VERIFY_IS_NOT_APPROX(q1 * v1, Quaternionx(AngleAxisx(aa.angle()*2,aa.axis())) * v1);
+  }
 
   // from two vector creation
   VERIFY_IS_APPROX( v2.normalized(),(q2.setFromTwoVectors(v1, v2)*v1).normalized());
@@ -117,6 +125,7 @@ template<typename Scalar, int Options> void quaternion(void)
 template<typename Scalar> void mapQuaternion(void){
   typedef Map<Quaternion<Scalar>, Aligned> MQuaternionA;
   typedef Map<Quaternion<Scalar> > MQuaternionUA;
+  typedef Map<const Quaternion<Scalar> > MCQuaternionUA;
   typedef Quaternion<Scalar> Quaternionx;
 
   EIGEN_ALIGN16 Scalar array1[4];
@@ -124,6 +133,7 @@ template<typename Scalar> void mapQuaternion(void){
   EIGEN_ALIGN16 Scalar array3[4+1];
   Scalar* array3unaligned = array3+1;
 
+//  std::cerr << array1 << " " << array2 << " " << array3 << "\n";
   MQuaternionA(array1).coeffs().setRandom();
   (MQuaternionA(array2)) = MQuaternionA(array1);
   (MQuaternionUA(array3unaligned)) = MQuaternionA(array1);
@@ -131,11 +141,14 @@ template<typename Scalar> void mapQuaternion(void){
   Quaternionx q1 = MQuaternionA(array1);
   Quaternionx q2 = MQuaternionA(array2);
   Quaternionx q3 = MQuaternionUA(array3unaligned);
+  Quaternionx q4 = MCQuaternionUA(array3unaligned);
 
   VERIFY_IS_APPROX(q1.coeffs(), q2.coeffs());
   VERIFY_IS_APPROX(q1.coeffs(), q3.coeffs());
+  VERIFY_IS_APPROX(q4.coeffs(), q3.coeffs());
   #ifdef EIGEN_VECTORIZE
-  VERIFY_RAISES_ASSERT((MQuaternionA(array3unaligned)));
+  if(internal::packet_traits<Scalar>::Vectorizable)
+    VERIFY_RAISES_ASSERT((MQuaternionA(array3unaligned)));
   #endif
 }
 
@@ -158,21 +171,39 @@ template<typename Scalar> void quaternionAlignment(void){
 
   VERIFY_IS_APPROX(q1->coeffs(), q2->coeffs());
   VERIFY_IS_APPROX(q1->coeffs(), q3->coeffs());
-  #ifdef EIGEN_VECTORIZE
-  VERIFY_RAISES_ASSERT((::new(reinterpret_cast<void*>(arrayunaligned)) QuaternionA));
+  #if defined(EIGEN_VECTORIZE) && EIGEN_ALIGN_STATICALLY
+  if(internal::packet_traits<Scalar>::Vectorizable)
+    VERIFY_RAISES_ASSERT((::new(reinterpret_cast<void*>(arrayunaligned)) QuaternionA));
   #endif
 }
 
+template<typename PlainObjectType> void check_const_correctness(const PlainObjectType&)
+{
+  // there's a lot that we can't test here while still having this test compile!
+  // the only possible approach would be to run a script trying to compile stuff and checking that it fails.
+  // CMake can help with that.
+
+  // verify that map-to-const don't have LvalueBit
+  typedef typename internal::add_const<PlainObjectType>::type ConstPlainObjectType;
+  VERIFY( !(internal::traits<Map<ConstPlainObjectType> >::Flags & LvalueBit) );
+  VERIFY( !(internal::traits<Map<ConstPlainObjectType, Aligned> >::Flags & LvalueBit) );
+  VERIFY( !(Map<ConstPlainObjectType>::Flags & LvalueBit) );
+  VERIFY( !(Map<ConstPlainObjectType, Aligned>::Flags & LvalueBit) );
+}
+
+
 void test_geo_quaternion()
 {
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1(( quaternion<float,AutoAlign>() ));
+    CALL_SUBTEST_1( check_const_correctness(Quaternionf()) );
     CALL_SUBTEST_2(( quaternion<double,AutoAlign>() ));
+    CALL_SUBTEST_2( check_const_correctness(Quaterniond()) );
     CALL_SUBTEST_3(( quaternion<float,DontAlign>() ));
     CALL_SUBTEST_4(( quaternion<double,DontAlign>() ));
     CALL_SUBTEST_5(( quaternionAlignment<float>() ));
     CALL_SUBTEST_6(( quaternionAlignment<double>() ));
-    CALL_SUBTEST( mapQuaternion<float>() );
-    CALL_SUBTEST( mapQuaternion<double>() );
+    CALL_SUBTEST_1( mapQuaternion<float>() );
+    CALL_SUBTEST_2( mapQuaternion<double>() );
   }
 }

test/geo_transformations.cpp

@@ -442,8 +442,9 @@ template<typename Scalar> void transform_alignment()
   
   VERIFY_IS_APPROX( (*p1) * (*p1), (*p2)*(*p3));
   
-  #ifdef EIGEN_VECTORIZE
-  VERIFY_RAISES_ASSERT((::new(reinterpret_cast<void*>(array3u)) Projective3a));
+  #if defined(EIGEN_VECTORIZE) && EIGEN_ALIGN_STATICALLY
+  if(internal::packet_traits<Scalar>::Vectorizable)
+    VERIFY_RAISES_ASSERT((::new(reinterpret_cast<void*>(array3u)) Projective3a));
   #endif
 }
 
@@ -455,7 +456,8 @@ void test_geo_transformations()
     
     CALL_SUBTEST_2(( transformations<float,AffineCompact,AutoAlign>() ));
     CALL_SUBTEST_2(( non_projective_only<float,AffineCompact,AutoAlign>() ));
-
+    CALL_SUBTEST_2(( transform_alignment<float>() ));
+    
     CALL_SUBTEST_3(( transformations<double,Projective,AutoAlign>() ));
     CALL_SUBTEST_3(( transformations<double,Projective,DontAlign>() ));
     CALL_SUBTEST_3(( transform_alignment<double>() ));

test/jacobisvd.cpp

@@ -23,6 +23,9 @@
 // License and a copy of the GNU General Public License along with
 // Eigen. If not, see <http://www.gnu.org/licenses/>.
 
+// discard stack allocation as that too bypasses malloc
+#define EIGEN_STACK_ALLOCATION_LIMIT 0
+#define EIGEN_RUNTIME_NO_MALLOC
 #include "main.h"
 #include <Eigen/SVD>
 
@@ -241,6 +244,46 @@ void jacobisvd_inf_nan()
   svd.compute(m, ComputeFullU | ComputeFullV);
 }
 
+void jacobisvd_preallocate()
+{
+  Vector3f v(3.f, 2.f, 1.f);
+  MatrixXf m = v.asDiagonal();
+
+  internal::set_is_malloc_allowed(false);
+  VERIFY_RAISES_ASSERT(VectorXf v(10);)
+  JacobiSVD<MatrixXf> svd;
+  internal::set_is_malloc_allowed(true);
+  svd.compute(m);
+  VERIFY_IS_APPROX(svd.singularValues(), v);
+
+  JacobiSVD<MatrixXf> svd2(3,3);
+  internal::set_is_malloc_allowed(false);
+  svd2.compute(m);
+  internal::set_is_malloc_allowed(true);
+  VERIFY_IS_APPROX(svd2.singularValues(), v);
+  VERIFY_RAISES_ASSERT(svd2.matrixU());
+  VERIFY_RAISES_ASSERT(svd2.matrixV());
+  svd2.compute(m, ComputeFullU | ComputeFullV);
+  VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
+  VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
+  internal::set_is_malloc_allowed(false);
+  svd2.compute(m);
+  internal::set_is_malloc_allowed(true);
+
+  JacobiSVD<MatrixXf> svd3(3,3,ComputeFullU|ComputeFullV);
+  internal::set_is_malloc_allowed(false);
+  svd2.compute(m);
+  internal::set_is_malloc_allowed(true);
+  VERIFY_IS_APPROX(svd2.singularValues(), v);
+  VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity());
+  VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity());
+  internal::set_is_malloc_allowed(false);
+  svd2.compute(m, ComputeFullU|ComputeFullV);
+  internal::set_is_malloc_allowed(true);
+
+  
+}
+
 void test_jacobisvd()
 {
   CALL_SUBTEST_3(( jacobisvd_verify_assert(Matrix3f()) ));
@@ -290,4 +333,7 @@ void test_jacobisvd()
 
   // Test problem size constructors
   CALL_SUBTEST_7( JacobiSVD<MatrixXf>(10,10) );
+
+  // Check that preallocation avoids subsequent mallocs
+  CALL_SUBTEST_9( jacobisvd_preallocate() );
 }

test/map.cpp

@@ -50,7 +50,8 @@ template<typename VectorType> void map_class_vector(const VectorType& m)
   VERIFY_IS_EQUAL(ma1, ma2);
   VERIFY_IS_EQUAL(ma1, ma3);
   #ifdef EIGEN_VECTORIZE
-  VERIFY_RAISES_ASSERT((Map<VectorType,Aligned>(array3unaligned, size)))
+  if(internal::packet_traits<Scalar>::Vectorizable)
+    VERIFY_RAISES_ASSERT((Map<VectorType,Aligned>(array3unaligned, size)))
   #endif
 
   internal::aligned_delete(array1, size);

test/nullary.cpp

@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2010-2011 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -27,17 +27,18 @@
 template<typename MatrixType>
 bool equalsIdentity(const MatrixType& A)
 {
+  typedef typename MatrixType::Index Index;
   typedef typename MatrixType::Scalar Scalar;
   Scalar zero = static_cast<Scalar>(0);
 
   bool offDiagOK = true;
-  for (int i = 0; i < A.rows(); ++i) {
-    for (int j = i+1; j < A.cols(); ++j) {
+  for (Index i = 0; i < A.rows(); ++i) {
+    for (Index j = i+1; j < A.cols(); ++j) {
       offDiagOK = offDiagOK && (A(i,j) == zero);
     }
   }
-  for (int i = 0; i < A.rows(); ++i) {
-    for (int j = 0; j < i; ++j) {
+  for (Index i = 0; i < A.rows(); ++i) {
+    for (Index j = 0; j < std::min(i, A.cols()); ++j) {
       offDiagOK = offDiagOK && (A(i,j) == zero);
     }
   }
@@ -65,22 +66,22 @@ void testVectorType(const VectorType& base)
   for (int i=0; i<size; ++i)
     n(i) = low+i*step;
 
-  VERIFY( (m-n).norm() < std::numeric_limits<Scalar>::epsilon()*10e3 );
+  VERIFY_IS_APPROX(m,n);
 
   // random access version
   m = VectorType::LinSpaced(size,low,high);
-  VERIFY( (m-n).norm() < std::numeric_limits<Scalar>::epsilon()*10e3 );
+  VERIFY_IS_APPROX(m,n);
 
   // Assignment of a RowVectorXd to a MatrixXd (regression test for bug #79).
   VERIFY( (MatrixXd(RowVectorXd::LinSpaced(3, 0, 1)) - RowVector3d(0, 0.5, 1)).norm() < std::numeric_limits<Scalar>::epsilon() );
 
   // These guys sometimes fail! This is not good. Any ideas how to fix them!?
-  //VERIFY( m(m.size()-1) == high );
-  //VERIFY( m(0) == low );
+//   VERIFY( m(m.size()-1) == high );
+//   VERIFY( m(0) == low );
 
   // sequential access version
   m = VectorType::LinSpaced(Sequential,size,low,high);
-  VERIFY( (m-n).norm() < std::numeric_limits<Scalar>::epsilon()*10e3 );
+  VERIFY_IS_APPROX(m,n);
 
   // These guys sometimes fail! This is not good. Any ideas how to fix them!?
   //VERIFY( m(m.size()-1) == high );
@@ -114,12 +115,14 @@ void testMatrixType(const MatrixType& m)
 void test_nullary()
 {
   CALL_SUBTEST_1( testMatrixType(Matrix2d()) );
-  CALL_SUBTEST_2( testMatrixType(MatrixXcf(50,50)) );
-  CALL_SUBTEST_3( testMatrixType(MatrixXf(5,7)) );
-  CALL_SUBTEST_4( testVectorType(VectorXd(51)) );
-  CALL_SUBTEST_5( testVectorType(VectorXd(41)) );
-  CALL_SUBTEST_6( testVectorType(Vector3d()) );
-  CALL_SUBTEST_7( testVectorType(VectorXf(51)) );
-  CALL_SUBTEST_8( testVectorType(VectorXf(41)) );
-  CALL_SUBTEST_9( testVectorType(Vector3f()) );
+  CALL_SUBTEST_2( testMatrixType(MatrixXcf(internal::random<int>(1,300),internal::random<int>(1,300))) );
+  CALL_SUBTEST_3( testMatrixType(MatrixXf(internal::random<int>(1,300),internal::random<int>(1,300))) );
+  
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST_4( testVectorType(VectorXd(internal::random<int>(1,300))) );
+    CALL_SUBTEST_5( testVectorType(Vector4d()) );  // regression test for bug 232
+    CALL_SUBTEST_6( testVectorType(Vector3d()) );
+    CALL_SUBTEST_7( testVectorType(VectorXf(internal::random<int>(1,300))) );
+    CALL_SUBTEST_8( testVectorType(Vector3f()) );
+  }
 }

test/packetmath.cpp

@@ -44,7 +44,7 @@ template<typename Scalar> bool areApproxAbs(const Scalar* a, const Scalar* b, in
   {
     if (!isApproxAbs(a[i],b[i],refvalue))
     {
-      std::cout << "a[" << i << "]: " << a[i] << " != b[" << i << "]: " << b[i] << std::endl;
+      std::cout << "[" << Map<const Matrix<Scalar,1,Dynamic> >(a,size) << "]" << " != " << Map<const Matrix<Scalar,1,Dynamic> >(b,size) << "\n";
       return false;
     }
   }
@@ -57,7 +57,7 @@ template<typename Scalar> bool areApprox(const Scalar* a, const Scalar* b, int s
   {
     if (!internal::isApprox(a[i],b[i]))
     {
-      std::cout << "a[" << i << "]: " << a[i] << " != b[" << i << "]: " << b[i] << std::endl;
+      std::cout << "[" << Map<const Matrix<Scalar,1,Dynamic> >(a,size) << "]" << " != " << Map<const Matrix<Scalar,1,Dynamic> >(b,size) << "\n";
       return false;
     }
   }
@@ -180,7 +180,7 @@ template<typename Scalar> void packetmath()
     internal::pstore(data2, internal::pset1<Packet>(data1[offset]));
     VERIFY(areApprox(ref, data2, PacketSize) && "internal::pset1");
   }
-
+  
   VERIFY(internal::isApprox(data1[0], internal::pfirst(internal::pload<Packet>(data1))) && "internal::pfirst");
   
   if(PacketSize>1)
@@ -275,6 +275,11 @@ template<typename Scalar> void packetmath_real()
   for (int i=0; i<PacketSize; ++i)
     ref[0] = std::max(ref[0],data1[i]);
   VERIFY(internal::isApprox(ref[0], internal::predux_max(internal::pload<Packet>(data1))) && "internal::predux_max");
+  
+  for (int i=0; i<PacketSize; ++i)
+    ref[i] = data1[0]+Scalar(i);
+  internal::pstore(data2, internal::plset(data1[0]));
+  VERIFY(areApprox(ref, data2, PacketSize) && "internal::plset");
 }
 
 template<typename Scalar,bool ConjLhs,bool ConjRhs> void test_conj_helper(Scalar* data1, Scalar* data2, Scalar* ref, Scalar* pval)

test/vectorization_logic.cpp

@@ -123,7 +123,8 @@ template<typename Scalar, bool Enable = internal::packet_traits<Scalar>::Vectori
         (PacketSize==8 ? 4 : PacketSize==4 ? 6 : PacketSize==2 ? ((Matrix11::Flags&RowMajorBit)?2:3) : /*PacketSize==1 ?*/ 1),
         (PacketSize==8 ? 6 : PacketSize==4 ? 2 : PacketSize==2 ? ((Matrix11::Flags&RowMajorBit)?3:2) : /*PacketSize==1 ?*/ 3)
       > Matrix3;
-      
+    
+    #if !EIGEN_GCC_AND_ARCH_DOESNT_WANT_STACK_ALIGNMENT
     VERIFY(test_assign(Vector1(),Vector1(),
       InnerVectorizedTraversal,CompleteUnrolling));
     VERIFY(test_assign(Vector1(),Vector1()+Vector1(),
@@ -173,25 +174,7 @@ template<typename Scalar, bool Enable = internal::packet_traits<Scalar>::Vectori
       VERIFY(test_assign(Matrix11(),Matrix<Scalar,17,17>().template block<PacketSize,PacketSize>(2,3)+Matrix<Scalar,17,17>().template block<PacketSize,PacketSize>(10,4),
       DefaultTraversal,CompleteUnrolling));
     }
-
-    VERIFY(test_assign(MatrixXX(10,10),MatrixXX(20,20).block(10,10,2,3),
-      SliceVectorizedTraversal,NoUnrolling));
-
-    VERIFY((test_assign<
-            Map<Matrix22, Aligned, OuterStride<3*PacketSize> >,
-            Matrix22
-            >(InnerVectorizedTraversal,CompleteUnrolling)));
-
-    VERIFY((test_assign<
-            Map<Matrix22, Aligned, InnerStride<3*PacketSize> >,
-            Matrix22
-            >(DefaultTraversal,CompleteUnrolling)));
-
-    VERIFY((test_assign(Matrix11(), Matrix11()*Matrix11(), InnerVectorizedTraversal, CompleteUnrolling)));
-
-    VERIFY(test_redux(VectorX(10),
-      LinearVectorizedTraversal,NoUnrolling));
-
+    
     VERIFY(test_redux(Matrix3(),
       LinearVectorizedTraversal,CompleteUnrolling));
 
@@ -206,6 +189,27 @@ template<typename Scalar, bool Enable = internal::packet_traits<Scalar>::Vectori
 
     VERIFY(test_redux(Matrix44r().template block<1,2*PacketSize>(2,1),
       LinearVectorizedTraversal,CompleteUnrolling));
+    
+    VERIFY((test_assign<
+            Map<Matrix22, Aligned, OuterStride<3*PacketSize> >,
+            Matrix22
+            >(InnerVectorizedTraversal,CompleteUnrolling)));
+
+    VERIFY((test_assign<
+            Map<Matrix22, Aligned, InnerStride<3*PacketSize> >,
+            Matrix22
+            >(DefaultTraversal,CompleteUnrolling)));
+
+    VERIFY((test_assign(Matrix11(), Matrix11()*Matrix11(), InnerVectorizedTraversal, CompleteUnrolling)));
+    #endif
+
+    VERIFY(test_assign(MatrixXX(10,10),MatrixXX(20,20).block(10,10,2,3),
+      SliceVectorizedTraversal,NoUnrolling));
+
+    VERIFY(test_redux(VectorX(10),
+      LinearVectorizedTraversal,NoUnrolling));
+
+    
   }
 };
 
@@ -219,10 +223,10 @@ void test_vectorization_logic()
 
 #ifdef EIGEN_VECTORIZE
 
-  vectorization_logic<float>::run();
-  vectorization_logic<double>::run();
-  vectorization_logic<std::complex<float> >::run();
-  vectorization_logic<std::complex<double> >::run();
+  CALL_SUBTEST( vectorization_logic<float>::run() );
+  CALL_SUBTEST( vectorization_logic<double>::run() );
+  CALL_SUBTEST( vectorization_logic<std::complex<float> >::run() );
+  CALL_SUBTEST( vectorization_logic<std::complex<double> >::run() );
   
   if(internal::packet_traits<float>::Vectorizable)
   {

unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h

@@ -26,7 +26,7 @@
 #define EIGEN_MATRIX_EXPONENTIAL
 
 #ifdef _MSC_VER
-  template <typename Scalar> Scalar log2(Scalar v) { return std::log(v)/std::log(Scalar(2)); }
+  template <typename Scalar> Scalar log2(Scalar v) { using std::log; return log(v)/log(Scalar(2)); }
 #endif
 
 
@@ -250,14 +250,17 @@ EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade13(const MatrixType
 template <typename MatrixType>
 void MatrixExponential<MatrixType>::computeUV(float)
 {
+  using std::max;
+  using std::pow;
+  using std::ceil;
   if (m_l1norm < 4.258730016922831e-001) {
     pade3(m_M);
   } else if (m_l1norm < 1.880152677804762e+000) {
     pade5(m_M);
   } else {
     const float maxnorm = 3.925724783138660f;
-    m_squarings = std::max(0, (int)ceil(log2(m_l1norm / maxnorm)));
-    MatrixType A = m_M / std::pow(Scalar(2), Scalar(static_cast<RealScalar>(m_squarings)));
+    m_squarings = max(0, (int)ceil(log2(m_l1norm / maxnorm)));
+    MatrixType A = m_M / pow(Scalar(2), Scalar(static_cast<RealScalar>(m_squarings)));
     pade7(A);
   }
 }
@@ -265,6 +268,9 @@ void MatrixExponential<MatrixType>::computeUV(float)
 template <typename MatrixType>
 void MatrixExponential<MatrixType>::computeUV(double)
 {
+  using std::max;
+  using std::pow;
+  using std::ceil;
   if (m_l1norm < 1.495585217958292e-002) {
     pade3(m_M);
   } else if (m_l1norm < 2.539398330063230e-001) {
@@ -275,8 +281,8 @@ void MatrixExponential<MatrixType>::computeUV(double)
     pade9(m_M);
   } else {
     const double maxnorm = 5.371920351148152;
-    m_squarings = std::max(0, (int)ceil(log2(m_l1norm / maxnorm)));
-    MatrixType A = m_M / std::pow(Scalar(2), Scalar(m_squarings));
+    m_squarings = max(0, (int)ceil(log2(m_l1norm / maxnorm)));
+    MatrixType A = m_M / pow(Scalar(2), Scalar(m_squarings));
     pade13(A);
   }
 }

unsupported/Eigen/src/SparseExtra/CholmodSupportLegacy.h

@@ -441,7 +441,7 @@ void SparseLDLT<_MatrixType,Cholmod>::compute(const _MatrixType& a)
   m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
   //m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
   // TODO
-  if (this->m_flags & Base::IncompleteFactorization)
+  if (this->m_flags & IncompleteFactorization)
   {
     m_cholmod.nmethods = 1;
     //m_cholmod.method[0].ordering = CHOLMOD_NATURAL;

unsupported/Eigen/src/SparseExtra/SimplicialCholesky.h

@@ -295,7 +295,7 @@ void SimplicialCholesky<_MatrixType,_UpLo>::analyzePattern(const MatrixType& a)
   m_parent.resize(size);
   m_nonZerosPerCol.resize(size);
   
-  Index* tags = ei_aligned_stack_new(Index, size);
+  ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
   
   // TODO allows to configure the permutation
   {
@@ -338,9 +338,6 @@ void SimplicialCholesky<_MatrixType,_UpLo>::analyzePattern(const MatrixType& a)
     }
   }
   
-  // release workspace
-  ei_aligned_stack_delete(Index, tags, size);
-  
   /* construct Lp index array from m_nonZerosPerCol column counts */
   Index* Lp = m_matrix._outerIndexPtr();
   Lp[0] = 0;
@@ -369,9 +366,9 @@ void SimplicialCholesky<_MatrixType,_UpLo>::factorize(const MatrixType& a)
   Index* Li = m_matrix._innerIndexPtr();
   Scalar* Lx = m_matrix._valuePtr();
 
-  Scalar* y       = ei_aligned_stack_new(Scalar, size);
-  Index* pattern  = ei_aligned_stack_new(Index, size);
-  Index* tags     = ei_aligned_stack_new(Index, size);
+  ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
+  ei_declare_aligned_stack_constructed_variable(Index,  pattern, size, 0);
+  ei_declare_aligned_stack_constructed_variable(Index,  tags, size, 0);
 
   SparseMatrix<Scalar,ColMajor,Index> ap(size,size);
   ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
@@ -443,11 +440,6 @@ void SimplicialCholesky<_MatrixType,_UpLo>::factorize(const MatrixType& a)
     }
   }
 
-  // release workspace
-  ei_aligned_stack_delete(Scalar, y, size);
-  ei_aligned_stack_delete(Index, pattern, size);
-  ei_aligned_stack_delete(Index, tags, size);
-  
   m_info = ok ? Success : NumericalIssue;
   m_factorizationIsOk = true;
 }

unsupported/Eigen/src/SparseExtra/SparseLDLTLegacy.h

@@ -245,7 +245,8 @@ void SparseLDLT<_MatrixType,Backend>::_symbolic(const _MatrixType& a)
   m_matrix.resize(size, size);
   m_parent.resize(size);
   m_nonZerosPerCol.resize(size);
-  Index * tags = ei_aligned_stack_new(Index, size);
+
+  ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0);
 
   const Index* Ap = a._outerIndexPtr();
   const Index* Ai = a._innerIndexPtr();
@@ -298,7 +299,6 @@ void SparseLDLT<_MatrixType,Backend>::_symbolic(const _MatrixType& a)
     Lp[k+1] = Lp[k] + m_nonZerosPerCol[k];
 
   m_matrix.resizeNonZeros(Lp[size]);
-  ei_aligned_stack_delete(Index, tags, size);
 }
 
 template<typename _MatrixType, typename Backend>
@@ -317,9 +317,9 @@ bool SparseLDLT<_MatrixType,Backend>::_numeric(const _MatrixType& a)
   Scalar* Lx = m_matrix._valuePtr();
   m_diag.resize(size);
 
-  Scalar * y = ei_aligned_stack_new(Scalar, size);
-  Index * pattern = ei_aligned_stack_new(Index, size);
-  Index * tags = ei_aligned_stack_new(Index, size);
+  ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0);
+  ei_declare_aligned_stack_constructed_variable(Index,  pattern, size, 0);
+  ei_declare_aligned_stack_constructed_variable(Index,  tags, size, 0);
   
   Index* P = 0;
   Index* Pinv = 0;
@@ -383,10 +383,6 @@ bool SparseLDLT<_MatrixType,Backend>::_numeric(const _MatrixType& a)
     }
   }
 
-  ei_aligned_stack_delete(Scalar, y, size);
-  ei_aligned_stack_delete(Index, pattern, size);
-  ei_aligned_stack_delete(Index, tags, size);
-
   return ok;  /* success, diagonal of D is all nonzero */
 }
 

unsupported/test/openglsupport.cpp

@@ -57,7 +57,7 @@ using namespace Eigen;
   }
   
 #define VERIFY_UNIFORMi(NAME,TYPE) { \
-    TYPE value; value.setRandom(); \
+    TYPE value = TYPE::Random().eval().cast<float>().cast<TYPE::Scalar>(); \
     TYPE data; \
     int loc = glGetUniformLocation(prg_id, #NAME); \
     VERIFY((loc!=-1) && "uniform not found"); \
@@ -320,6 +320,7 @@ void test_openglsupport()
     else
       std::cerr << "Warning: opengl 3.0 was not tested\n";
     
+    #ifdef GLEW_ARB_gpu_shader_fp64
     if(GLEW_ARB_gpu_shader_fp64)
     {
       #ifdef GL_ARB_gpu_shader_fp64
@@ -343,6 +344,9 @@ void test_openglsupport()
     }
     else
       std::cerr << "Warning: GLEW_ARB_gpu_shader_fp64 was not tested\n";
+    #else
+      std::cerr << "Warning: GLEW_ARB_gpu_shader_fp64 was not tested\n";
+    #endif
   }
   
 }