bump that too

that eigen2 owes all its performance to nonstandard restrict keyword.
well, this can also improve portability in case some compiler doesn't have __restrict.

CMakeLists.txt

@@ -1,7 +1,7 @@
 project(Eigen)
 cmake_minimum_required(VERSION 2.6.2)
 
-set(EIGEN_VERSION_NUMBER "2.0.5")
+set(EIGEN_VERSION_NUMBER "2.0.6")
 set(EIGEN_VERSION "${EIGEN_VERSION_NUMBER}")
 
 set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)

Eigen/NewStdVector

@@ -0,0 +1,168 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2009 Hauke Heibel <hauke.heibel@googlemail.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/>.
+
+#ifndef EIGEN_STDVECTOR_MODULE_H
+#define EIGEN_STDVECTOR_MODULE_H
+
+#include "Core"
+#include <vector>
+
+namespace Eigen {
+
+// This one is needed to prevent reimplementing the whole std::vector.
+template <class T>
+class aligned_allocator_indirection : public aligned_allocator<T>
+{
+public:
+  typedef size_t    size_type;
+  typedef ptrdiff_t difference_type;
+  typedef T*        pointer;
+  typedef const T*  const_pointer;
+  typedef T&        reference;
+  typedef const T&  const_reference;
+  typedef T         value_type;
+
+  template<class U>
+  struct rebind
+  {
+    typedef aligned_allocator_indirection<U> other;
+  };
+
+  aligned_allocator_indirection() throw() {}
+  aligned_allocator_indirection(const aligned_allocator_indirection& ) throw() : aligned_allocator<T>() {}
+  aligned_allocator_indirection(const aligned_allocator<T>& ) throw() {}
+  template<class U>
+  aligned_allocator_indirection(const aligned_allocator_indirection<U>& ) throw() {}
+  template<class U>
+  aligned_allocator_indirection(const aligned_allocator<U>& ) throw() {}
+  ~aligned_allocator_indirection() throw() {}
+};
+
+#ifdef _MSC_VER
+
+  // sometimes, MSVC detects, at compile time, that the argument x
+  // in std::vector::resize(size_t s,T x) won't be aligned and generate an error
+  // even if this function is never called. Whence this little wrapper.
+  #define EIGEN_WORKAROUND_MSVC_STD_VECTOR(T) Eigen::ei_workaround_msvc_std_vector<T>
+  template<typename T> struct ei_workaround_msvc_std_vector : public T
+  {
+    inline ei_workaround_msvc_std_vector() : T() {}
+    inline ei_workaround_msvc_std_vector(const T& other) : T(other) {}
+    inline operator T& () { return *static_cast<T*>(this); }
+    inline operator const T& () const { return *static_cast<const T*>(this); }
+    template<typename OtherT>
+    inline T& operator=(const OtherT& other)
+    { T::operator=(other); return *this; }
+    inline ei_workaround_msvc_std_vector& operator=(const ei_workaround_msvc_std_vector& other)
+    { T::operator=(other); return *this; }
+  };
+
+#else
+
+  #define EIGEN_WORKAROUND_MSVC_STD_VECTOR(T) T
+
+#endif
+
+}
+
+namespace std {
+
+#define EIGEN_STD_VECTOR_SPECIALIZATION_BODY \
+  public:  \
+    typedef T value_type; \
+    typedef typename vector_base::allocator_type allocator_type; \
+    typedef typename vector_base::size_type size_type;  \
+    typedef typename vector_base::iterator iterator;  \
+    typedef typename vector_base::const_iterator const_iterator;  \
+    explicit vector(const allocator_type& a = allocator_type()) : vector_base(a) {}  \
+    template<typename InputIterator> \
+    vector(InputIterator first, InputIterator last, const allocator_type& a = allocator_type()) \
+    : vector_base(first, last, a) {} \
+    vector(const vector& c) : vector_base(c) {}  \
+    explicit vector(size_type num, const value_type& val = value_type()) : vector_base(num, val) {} \
+    vector(iterator start, iterator end) : vector_base(start, end) {}  \
+    vector& operator=(const vector& x) {  \
+      vector_base::operator=(x);  \
+      return *this;  \
+    }
+
+template<typename T>
+class vector<T,Eigen::aligned_allocator<T> >
+  : public vector<EIGEN_WORKAROUND_MSVC_STD_VECTOR(T),
+                  Eigen::aligned_allocator_indirection<EIGEN_WORKAROUND_MSVC_STD_VECTOR(T)> >
+{
+  typedef vector<EIGEN_WORKAROUND_MSVC_STD_VECTOR(T),
+                 Eigen::aligned_allocator_indirection<EIGEN_WORKAROUND_MSVC_STD_VECTOR(T)> > vector_base;
+  EIGEN_STD_VECTOR_SPECIALIZATION_BODY
+
+  void resize(size_type new_size)
+  { resize(new_size, T()); }
+
+#if defined(_VECTOR_)
+  // workaround MSVC std::vector implementation
+  void resize(size_type new_size, const value_type& x)
+  {
+    if (vector_base::size() < new_size)
+      vector_base::_Insert_n(vector_base::end(), new_size - vector_base::size(), x);
+    else if (new_size < vector_base::size())
+      vector_base::erase(vector_base::begin() + new_size, vector_base::end());
+  }
+  void push_back(const value_type& x)
+  { vector_base::push_back(x); } 
+  using vector_base::insert;  
+  iterator insert(const_iterator position, const value_type& x)
+  { return vector_base::insert(position,x); }
+  void insert(const_iterator position, size_type new_size, const value_type& x)
+  { vector_base::insert(position, new_size, x); }
+#elif defined(_GLIBCXX_VECTOR) && EIGEN_GNUC_AT_LEAST(4,2)
+  // workaround GCC std::vector implementation
+  void resize(size_type new_size, const value_type& x)
+  {
+    if (new_size < vector_base::size())
+      vector_base::_M_erase_at_end(this->_M_impl._M_start + new_size);
+    else
+      vector_base::insert(vector_base::end(), new_size - vector_base::size(), x);
+  }
+#elif defined(_GLIBCXX_VECTOR) && (!EIGEN_GNUC_AT_LEAST(4,1))
+  // Note that before gcc-4.1 we already have: std::vector::resize(size_type,const T&),
+  // no no need to workaround !
+  using vector_base::resize;
+#else
+  // either GCC 4.1 or non-GCC
+  // default implementation which should always work.
+  void resize(size_type new_size, const value_type& x)
+  {
+    if (new_size < vector_base::size())
+      vector_base::erase(vector_base::begin() + new_size, vector_base::end());
+    else if (new_size > vector_base::size())
+      vector_base::insert(vector_base::end(), new_size - vector_base::size(), x);
+  }
+#endif
+
+};
+
+}
+
+#endif // EIGEN_STDVECTOR_MODULE_H

Eigen/StdVector

@@ -1,3 +1,7 @@
+#ifdef EIGEN_USE_NEW_STDVECTOR
+#include "NewStdVector"
+#else
+
 #ifndef EIGEN_STDVECTOR_MODULE_H
 #define EIGEN_STDVECTOR_MODULE_H
 
@@ -139,3 +143,5 @@ class vector<T,DummyAlloc,true>
 }
 
 #endif // EIGEN_STDVECTOR_MODULE_H
+
+#endif // EIGEN_USE_NEW_STDVECTOR

Eigen/src/Core/Matrix.h

@@ -137,6 +137,9 @@ class Matrix
     enum { NeedsToAlign = (Options&AutoAlign) == AutoAlign
                           && SizeAtCompileTime!=Dynamic && ((sizeof(Scalar)*SizeAtCompileTime)%16)==0 };
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign)
+    
+    Base& base() { return *static_cast<Base*>(this); }
+    const Base& base() const { return *static_cast<const Base*>(this); }
 
     EIGEN_STRONG_INLINE int rows() const { return m_storage.rows(); }
     EIGEN_STRONG_INLINE int cols() const { return m_storage.cols(); }
@@ -395,13 +398,8 @@ class Matrix
     /** Override MatrixBase::swap() since for dynamic-sized matrices of same type it is enough to swap the
       * data pointers.
       */
-    inline void swap(Matrix& other)
-    {
-      if (Base::SizeAtCompileTime==Dynamic)
-        m_storage.swap(other.m_storage);
-      else
-        this->Base::swap(other);
-    }
+    template<typename OtherDerived>
+    void swap(const MatrixBase<OtherDerived>& other);
 
     /** \name Map
       * These are convenience functions returning Map objects. The Map() static functions return unaligned Map objects,
@@ -540,8 +538,39 @@ class Matrix
                         && (_Options & (AutoAlign|RowMajor)) == _Options),
           INVALID_MATRIX_TEMPLATE_PARAMETERS)
     }
+    
+    template<typename MatrixType, typename OtherDerived, bool IsSameType, bool IsDynamicSize>
+    friend struct ei_matrix_swap_impl;
 };
 
+template<typename MatrixType, typename OtherDerived,
+         bool IsSameType = ei_is_same_type<MatrixType, OtherDerived>::ret,
+         bool IsDynamicSize = MatrixType::SizeAtCompileTime==Dynamic>
+struct ei_matrix_swap_impl
+{
+  static inline void run(MatrixType& matrix, MatrixBase<OtherDerived>& other)
+  {
+    matrix.base().swap(other);
+  }
+};
+
+template<typename MatrixType, typename OtherDerived>
+struct ei_matrix_swap_impl<MatrixType, OtherDerived, true, true>
+{
+  static inline void run(MatrixType& matrix, MatrixBase<OtherDerived>& other)
+  {
+    matrix.m_storage.swap(other.derived().m_storage);
+  }
+};
+
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+template<typename OtherDerived>
+inline void Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::swap(const MatrixBase<OtherDerived>& other)
+{
+  ei_matrix_swap_impl<Matrix, OtherDerived>::run(*this, *const_cast<MatrixBase<OtherDerived>*>(&other));
+}
+
+
 /** \defgroup matrixtypedefs Global matrix typedefs
   *
   * \ingroup Core_Module

Eigen/src/Core/MatrixBase.h

@@ -532,8 +532,11 @@ template<typename Derived> class MatrixBase
     typename ei_traits<Derived>::Scalar minCoeff() const;
     typename ei_traits<Derived>::Scalar maxCoeff() const;
 
-    typename ei_traits<Derived>::Scalar minCoeff(int* row, int* col = 0) const;
-    typename ei_traits<Derived>::Scalar maxCoeff(int* row, int* col = 0) const;
+    typename ei_traits<Derived>::Scalar minCoeff(int* row, int* col) const;
+    typename ei_traits<Derived>::Scalar maxCoeff(int* row, int* col) const;
+
+    typename ei_traits<Derived>::Scalar minCoeff(int* index) const;
+    typename ei_traits<Derived>::Scalar maxCoeff(int* index) const;
 
     template<typename BinaryOp>
     typename ei_result_of<BinaryOp(typename ei_traits<Derived>::Scalar)>::type

Eigen/src/Core/Visitor.h

@@ -164,7 +164,7 @@ struct ei_functor_traits<ei_max_coeff_visitor<Scalar> > {
 /** \returns the minimum of all coefficients of *this
   * and puts in *row and *col its location.
   *
-  * \sa MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff()
+  * \sa MatrixBase::minCoeff(int*), MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff()
   */
 template<typename Derived>
 typename ei_traits<Derived>::Scalar
@@ -177,6 +177,22 @@ MatrixBase<Derived>::minCoeff(int* row, int* col) const
   return minVisitor.res;
 }
 
+/** \returns the minimum of all coefficients of *this
+  * and puts in *index its location.
+  *
+  * \sa MatrixBase::minCoeff(int*,int*), MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff()
+  */
+template<typename Derived>
+typename ei_traits<Derived>::Scalar
+MatrixBase<Derived>::minCoeff(int* index) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  ei_min_coeff_visitor<Scalar> minVisitor;
+  this->visit(minVisitor);
+  *index = (RowsAtCompileTime==1) ? minVisitor.col : minVisitor.row;
+  return minVisitor.res;
+}
+
 /** \returns the maximum of all coefficients of *this
   * and puts in *row and *col its location.
   *
@@ -193,5 +209,20 @@ MatrixBase<Derived>::maxCoeff(int* row, int* col) const
   return maxVisitor.res;
 }
 
+/** \returns the maximum of all coefficients of *this
+  * and puts in *index its location.
+  *
+  * \sa MatrixBase::maxCoeff(int*,int*), MatrixBase::minCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::maxCoeff()
+  */
+template<typename Derived>
+typename ei_traits<Derived>::Scalar
+MatrixBase<Derived>::maxCoeff(int* index) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  ei_max_coeff_visitor<Scalar> maxVisitor;
+  this->visit(maxVisitor);
+  *index = (RowsAtCompileTime==1) ? maxVisitor.col : maxVisitor.row;
+  return maxVisitor.res;
+}
 
 #endif // EIGEN_VISITOR_H

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

@@ -315,4 +315,7 @@ struct ei_palign_impl<Offset,__m128d>
 };
 #endif
 
+#define ei_vec4f_swizzle1(v,p,q,r,s) \
+  (_mm_castsi128_ps(_mm_shuffle_epi32( _mm_castps_si128(v), ((s)<<6|(r)<<4|(q)<<2|(p)))))
+
 #endif // EIGEN_PACKET_MATH_SSE_H

Eigen/src/Core/util/Constants.h

@@ -239,4 +239,16 @@ enum {
   HasDirectAccess = DirectAccessBit
 };
 
+const int EiArch_Generic = 0x0;
+const int EiArch_SSE     = 0x1;
+const int EiArch_AltiVec = 0x2;
+
+#if defined EIGEN_VECTORIZE_SSE
+  const int EiArch = EiArch_SSE;
+#elif defined EIGEN_VECTORIZE_ALTIVEC
+  const int EiArch = EiArch_AltiVec;
+#else
+  const int EiArch = EiArch_Generic;
+#endif
+
 #endif // EIGEN_CONSTANTS_H

Eigen/src/Core/util/Macros.h

@@ -30,30 +30,48 @@
 
 #define EIGEN_WORLD_VERSION 2
 #define EIGEN_MAJOR_VERSION 0
-#define EIGEN_MINOR_VERSION 5
+#define EIGEN_MINOR_VERSION 6
 
 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
                                       (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \
                                                                  EIGEN_MINOR_VERSION>=z))))
 
-// if the compiler is GNUC, disable 16 byte alignment on exotic archs that probably don't need it, and on which
-// it may be extra trouble to get aligned memory allocation to work (example: on ARM, overloading new[] is a PITA
-// because extra memory must be allocated for bookkeeping).
-// if the compiler is not GNUC, just cross fingers that the architecture isn't too exotic, because we don't want
-// to keep track of all the different preprocessor symbols for all compilers.
-#if !defined(__GNUC__) || defined(__i386__) || defined(__x86_64__) || defined(__ppc__) || defined(__ia64__)
+// 16 byte alignment is only useful for vectorization. Since it affects the ABI, we need to enable 16 byte alignment on all
+// platforms where vectorization might be enabled. In theory we could always enable alignment, but it can be a cause of problems
+// on some platforms, so we just disable it in certain common platform (compiler+architecture combinations) to avoid these problems.
+#if defined(__GNUC__) && !(defined(__i386__) || defined(__x86_64__) || defined(__powerpc__) || defined(__ia64__))
+#define EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT 1
+#else
+#define EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT 0
+#endif
+
+#if defined(__GNUC__) && (__GNUC__ <= 3)
+#define EIGEN_GCC3_OR_OLDER 1
+#else
+#define EIGEN_GCC3_OR_OLDER 0
+#endif
+
+// FIXME vectorization + alignment is completely disabled with sun studio
+#if !EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT && !EIGEN_GCC3_OR_OLDER && !defined(__SUNPRO_CC)
   #define EIGEN_ARCH_WANTS_ALIGNMENT 1
 #else
-  #ifdef EIGEN_VECTORIZE
-    #error Vectorization enabled, but the architecture is not listed among those for which we require 16 byte alignment. If you added vectorization for another architecture, you also need to edit this list.
-  #endif
   #define EIGEN_ARCH_WANTS_ALIGNMENT 0
+#endif
+
+// EIGEN_ALIGN is the true test whether we want to align or not. It takes into account both the user choice to explicitly disable
+// alignment (EIGEN_DONT_ALIGN) and the architecture config (EIGEN_ARCH_WANTS_ALIGNMENT). Henceforth, only EIGEN_ALIGN should be used.
+#if EIGEN_ARCH_WANTS_ALIGNMENT && !defined(EIGEN_DONT_ALIGN)
+  #define EIGEN_ALIGN 1
+#else
+  #define EIGEN_ALIGN 0
+  #ifdef EIGEN_VECTORIZE
+    #error "Vectorization enabled, but our platform checks say that we don't do 16 byte alignment on this platform. If you added vectorization for another architecture, you also need to edit this platform check."
+  #endif
   #ifndef EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT
     #define EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT
   #endif
 #endif
 
-
 #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR
 #define EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION RowMajor
 #else
@@ -165,7 +183,7 @@ using Eigen::ei_cos;
  * 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_ARCH_WANTS_ALIGNMENT
+#if !EIGEN_ALIGN
 #define EIGEN_ALIGN_128
 #elif (defined __GNUC__)
 #define EIGEN_ALIGN_128 __attribute__((aligned(16)))
@@ -175,7 +193,12 @@ using Eigen::ei_cos;
 #error Please tell me what is the equivalent of __attribute__((aligned(16))) for your compiler
 #endif
 
-#define EIGEN_RESTRICT __restrict
+#ifdef EIGEN_DONT_USE_RESTRICT_KEYWORD
+  #define EIGEN_RESTRICT
+#endif
+#ifndef EIGEN_RESTRICT
+  #define EIGEN_RESTRICT __restrict
+#endif
 
 #ifndef EIGEN_STACK_ALLOCATION_LIMIT
 #define EIGEN_STACK_ALLOCATION_LIMIT 1000000

Eigen/src/Core/util/Memory.h

@@ -73,30 +73,23 @@ inline void* ei_aligned_malloc(size_t size)
     ei_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
   #endif
 
-  void *result;
-  #if EIGEN_HAS_POSIX_MEMALIGN && EIGEN_ARCH_WANTS_ALIGNMENT && !EIGEN_MALLOC_ALREADY_ALIGNED
-    #ifdef EIGEN_EXCEPTIONS
-      const int failed =
-    #endif
-    posix_memalign(&result, 16, size);
+  void *result;  
+  #if !EIGEN_ALIGN
+    result = malloc(size);
+  #elif EIGEN_MALLOC_ALREADY_ALIGNED
+    result = malloc(size);
+  #elif EIGEN_HAS_POSIX_MEMALIGN
+    if(posix_memalign(&result, 16, size)) result = 0;
+  #elif EIGEN_HAS_MM_MALLOC
+    result = _mm_malloc(size, 16);
+  #elif (defined _MSC_VER)
+    result = _aligned_malloc(size, 16);
   #else
-    #if !EIGEN_ARCH_WANTS_ALIGNMENT
-      result = malloc(size);
-    #elif EIGEN_MALLOC_ALREADY_ALIGNED
-      result = malloc(size);
-    #elif EIGEN_HAS_MM_MALLOC
-      result = _mm_malloc(size, 16);
-    #elif (defined _MSC_VER)
-      result = _aligned_malloc(size, 16);
-    #else
-      result = ei_handmade_aligned_malloc(size);
-    #endif
-    #ifdef EIGEN_EXCEPTIONS
-      const int failed = (result == 0);
-    #endif
+    result = ei_handmade_aligned_malloc(size);
   #endif
+    
   #ifdef EIGEN_EXCEPTIONS
-    if(failed)
+    if(result == 0)
       throw std::bad_alloc();
   #endif
   return result;
@@ -143,7 +136,7 @@ template<typename T, bool Align> inline T* ei_conditional_aligned_new(size_t siz
   */
 inline void ei_aligned_free(void *ptr)
 {
-  #if !EIGEN_ARCH_WANTS_ALIGNMENT
+  #if !EIGEN_ALIGN
     free(ptr);
   #elif EIGEN_MALLOC_ALREADY_ALIGNED
     free(ptr);
@@ -237,7 +230,7 @@ inline static int ei_alignmentOffset(const Scalar* ptr, int maxOffset)
                                                    ei_aligned_stack_free(PTR,sizeof(TYPE)*SIZE);} while(0)
 
 
-#if EIGEN_ARCH_WANTS_ALIGNMENT
+#if EIGEN_ALIGN
   #ifdef EIGEN_EXCEPTIONS
     #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \
       void* operator new(size_t size, const std::nothrow_t&) throw() { \

Eigen/src/Geometry/Quaternion.h

@@ -224,17 +224,45 @@ typedef Quaternion<float> Quaternionf;
   * double precision quaternion type */
 typedef Quaternion<double> Quaterniond;
 
+// Generic Quaternion * Quaternion product
+template<int Arch,typename Scalar> inline Quaternion<Scalar>
+ei_quaternion_product(const Quaternion<Scalar>& a, const Quaternion<Scalar>& b)
+{
+  return Quaternion<Scalar>
+  (
+    a.w() * b.w() - a.x() * b.x() - a.y() * b.y() - a.z() * b.z(),
+    a.w() * b.x() + a.x() * b.w() + a.y() * b.z() - a.z() * b.y(),
+    a.w() * b.y() + a.y() * b.w() + a.z() * b.x() - a.x() * b.z(),
+    a.w() * b.z() + a.z() * b.w() + a.x() * b.y() - a.y() * b.x()
+  );
+}
+
+#ifdef EIGEN_VECTORIZE_SSE
+template<> inline Quaternion<float>
+ei_quaternion_product<EiArch_SSE,float>(const Quaternion<float>& _a, const Quaternion<float>& _b)
+{
+  const __m128 mask = _mm_castsi128_ps(_mm_setr_epi32(0,0,0,0x80000000));
+  Quaternion<float> res;
+  __m128 a = _a.coeffs().packet<Aligned>(0);
+  __m128 b = _b.coeffs().packet<Aligned>(0);
+  __m128 flip1 = _mm_xor_ps(_mm_mul_ps(ei_vec4f_swizzle1(a,1,2,0,2),
+                                       ei_vec4f_swizzle1(b,2,0,1,2)),mask);
+  __m128 flip2 = _mm_xor_ps(_mm_mul_ps(ei_vec4f_swizzle1(a,3,3,3,1),
+                                       ei_vec4f_swizzle1(b,0,1,2,1)),mask);
+  ei_pstore(&res.x(),
+            _mm_add_ps(_mm_sub_ps(_mm_mul_ps(a,ei_vec4f_swizzle1(b,3,3,3,3)),
+                                  _mm_mul_ps(ei_vec4f_swizzle1(a,2,0,1,0),
+                                             ei_vec4f_swizzle1(b,1,2,0,0))),
+                       _mm_add_ps(flip1,flip2)));
+  return res;
+}
+#endif
+
 /** \returns the concatenation of two rotations as a quaternion-quaternion product */
 template <typename Scalar>
 inline Quaternion<Scalar> Quaternion<Scalar>::operator* (const Quaternion& other) const
 {
-  return Quaternion
-  (
-    this->w() * other.w() - this->x() * other.x() - this->y() * other.y() - this->z() * other.z(),
-    this->w() * other.x() + this->x() * other.w() + this->y() * other.z() - this->z() * other.y(),
-    this->w() * other.y() + this->y() * other.w() + this->z() * other.x() - this->x() * other.z(),
-    this->w() * other.z() + this->z() * other.w() + this->x() * other.y() - this->y() * other.x()
-  );
+  return ei_quaternion_product<EiArch>(*this,other);
 }
 
 /** \sa operator*(Quaternion) */

Eigen/src/Geometry/Transform.h

@@ -198,6 +198,10 @@ public:
 
   /** \sa MatrixBase::setIdentity() */
   void setIdentity() { m_matrix.setIdentity(); }
+  static const typename MatrixType::IdentityReturnType Identity()
+  {
+    return MatrixType::Identity();
+  }
 
   template<typename OtherDerived>
   inline Transform& scale(const MatrixBase<OtherDerived> &other);

doc/StlContainers.dox

@@ -6,6 +6,8 @@ namespace Eigen {
   - \ref summary
   - \ref allocator
   - \ref vector
+  - \ref newvector
+
 
 \section summary Executive summary
 
@@ -44,6 +46,18 @@ std::vector<Eigen::Vector4f>
 \endcode
 without having to worry about anything.
 
+\section newvector The new and improved workaround for std::vector
+
+Well, except that in Eigen 2.0 the <Eigen/StdVector> header causes some compatibility issues as it reimplements the std::vector<T> container in a not-fully-compatible way. If you want to avoid these issues, starting in Eigen 2.0.6 a new implementation is available, which will become default in the next major version of Eigen. You can enable it by defining EIGEN_USE_NEW_STDVECTOR:
+
+\code
+#define EIGEN_USE_NEW_STDVECTOR
+#include<Eigen/StdVector>
+\endcode
+
+This new implementation <a href="http://eigen.tuxfamily.org/dox-devel/StlContainers.html#vector">is documented here</a>. In particular, note that if you use it, you must specify Eigen::aligned_allocator<T> as the allocator type, otherwise it doesn't make any difference from the standard std::vector. This new std::vector implementation \b only overrides the standard one if used with this allocator, which guarantees that it doesn't break existing non-Eigen code.
+
+
 <span class="note">\b Explanation: The resize() method of std::vector takes a value_type argument (defaulting to value_type()). So with std::vector<Eigen::Vector4f>, some Eigen::Vector4f objects will be passed by value, which discards any alignment modifiers, so a Eigen::Vector4f can be created at an unaligned location. In order to avoid that, the only solution we saw was to specialize std::vector to make it work on a slight modification of, here, Eigen::Vector4f, that is able to deal properly with this situation.
 </span>
 

doc/TutorialGeometry.dox

@@ -152,7 +152,7 @@ OpenGL compatibility \b 3D </td><td>\code
 glLoadMatrixf(t.data());\endcode</td></tr>
 <tr><td>
 OpenGL compatibility \b 2D </td><td>\code
-Transform3f aux(Transform3f::Identity);
+Transform3f aux(Transform3f::Identity());
 aux.linear().corner<2,2>(TopLeft) = t.linear();
 aux.translation().start<2>() = t.translation();
 glLoadMatrixf(aux.data());\endcode</td></tr>

doc/UnalignedArrayAssert.dox

@@ -12,7 +12,7 @@ is explained here: http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html
 **** READ THIS WEB PAGE !!! ****"' failed.
 </pre>
 
-There are 3 known causes for this issue. Please read on to understand them and learn how to fix them.
+There are 4 known causes for this issue. Please read on to understand them and learn how to fix them.
 
 \b Table \b of \b contents
  - \ref where
@@ -21,6 +21,7 @@ There are 3 known causes for this issue. Please read on to understand them and l
  - \ref c3
  - \ref c4
  - \ref explanation
+ - \ref getrid
 
 \section where Where in my own code is the cause of the problem?
 
@@ -101,6 +102,16 @@ Eigen normally takes care of these alignment issues for you, by setting an align
 
 However there are a few corner cases where these alignment settings get overridden: they are the possible causes for this assertion.
 
+\section getrid I don't care about vectorization, how do I get rid of that stuff?
+
+Two possibilities:
+<ul>
+  <li>Define EIGEN_DONT_ALIGN. That disables all 128-bit alignment code, and in particular everything vectorization-related. But do note that this in particular breaks ABI compatibility with vectorized code.</li>
+  <li>Or define both EIGEN_DONT_VECTORIZE and EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT. This keeps the 128-bit alignment code and thus preserves ABI compatibility.</li>
+</ul>
+
+For more information, see <a href="http://eigen.tuxfamily.org/index.php?title=FAQ#I_disabled_vectorization.2C_but_I.27m_still_getting_annoyed_about_alignment_issues.21">this FAQ</a>.
+
 */
 
 }

test/CMakeLists.txt

@@ -211,6 +211,7 @@ ei_add_test(parametrizedline)
 ei_add_test(alignedbox)
 ei_add_test(regression)
 ei_add_test(stdvector)
+ei_add_test(newstdvector)
 if(QT4_FOUND)
   ei_add_test(qtvector " " "${QT_QTCORE_LIBRARY}")
 endif(QT4_FOUND)
@@ -218,6 +219,8 @@ ei_add_test(sparse_vector)
 ei_add_test(sparse_basic)
 ei_add_test(sparse_solvers " " "${SPARSE_LIBS}")
 ei_add_test(sparse_product)
+ei_add_test(swap)
+ei_add_test(visitor)
 
 # print a summary of the different options
 message("************************************************************")

test/geometry.cpp

@@ -151,7 +151,13 @@ template<typename Scalar> void geometry(void)
     a = ei_random<Scalar>(-Scalar(0.4)*Scalar(M_PI), Scalar(0.4)*Scalar(M_PI));
   q1 = AngleAxisx(a, v0.normalized());
   Transform3 t0, t1, t2;
+  // first test setIdentity() and Identity()
   t0.setIdentity();
+  VERIFY_IS_APPROX(t0.matrix(), Transform3::MatrixType::Identity());
+  t0.matrix().setZero();
+  t0 = Transform3::Identity();
+  VERIFY_IS_APPROX(t0.matrix(), Transform3::MatrixType::Identity());
+
   t0.linear() = q1.toRotationMatrix();
   t1.setIdentity();
   t1.linear() = q1.toRotationMatrix();

test/newstdvector.cpp

@@ -0,0 +1,164 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 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/>.
+
+#define EIGEN_USE_NEW_STDVECTOR
+#include "main.h"
+#include <Eigen/StdVector>
+#include <Eigen/Geometry>
+
+template<typename MatrixType>
+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,Eigen::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]);
+  v = w;
+  for(int i = 0; i < 20; i++)
+  {
+    VERIFY_IS_APPROX(w[i], v[i]);
+  }
+
+  v.resize(21);
+  v[20] = x;
+  VERIFY_IS_APPROX(v[20], x);
+  v.resize(22,y);
+  VERIFY_IS_APPROX(v[21], y);
+  v.push_back(x);
+  VERIFY_IS_APPROX(v[22], x);
+  VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(MatrixType));
+
+  // do a lot of push_back such that the vector gets internally resized
+  // (with memory reallocation)
+  MatrixType* ref = &w[0];
+  for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i)
+    v.push_back(w[i%w.size()]);
+  for(unsigned int i=23; i<v.size(); ++i)
+  {
+    VERIFY(v[i]==w[(i-23)%w.size()]);
+  }
+}
+
+template<typename TransformType>
+void check_stdvector_transform(const TransformType&)
+{
+  typedef typename TransformType::MatrixType MatrixType;
+  TransformType x(MatrixType::Random()), y(MatrixType::Random());
+  std::vector<TransformType,Eigen::aligned_allocator<TransformType> > v(10), w(20, y);
+  v[5] = x;
+  w[6] = v[5];
+  VERIFY_IS_APPROX(w[6], v[5]);
+  v = w;
+  for(int i = 0; i < 20; i++)
+  {
+    VERIFY_IS_APPROX(w[i], v[i]);
+  }
+
+  v.resize(21);
+  v[20] = x;
+  VERIFY_IS_APPROX(v[20], x);
+  v.resize(22,y);
+  VERIFY_IS_APPROX(v[21], y);
+  v.push_back(x);
+  VERIFY_IS_APPROX(v[22], x);
+  VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(TransformType));
+
+  // do a lot of push_back such that the vector gets internally resized
+  // (with memory reallocation)
+  TransformType* ref = &w[0];
+  for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i)
+    v.push_back(w[i%w.size()]);
+  for(unsigned int i=23; i<v.size(); ++i)
+  {
+    VERIFY(v[i].matrix()==w[(i-23)%w.size()].matrix());
+  }
+}
+
+template<typename QuaternionType>
+void check_stdvector_quaternion(const QuaternionType&)
+{
+  typedef typename QuaternionType::Coefficients Coefficients;
+  QuaternionType x(Coefficients::Random()), y(Coefficients::Random());
+  std::vector<QuaternionType,Eigen::aligned_allocator<QuaternionType> > v(10), w(20, y);
+  v[5] = x;
+  w[6] = v[5];
+  VERIFY_IS_APPROX(w[6], v[5]);
+  v = w;
+  for(int i = 0; i < 20; i++)
+  {
+    VERIFY_IS_APPROX(w[i], v[i]);
+  }
+
+  v.resize(21);
+  v[20] = x;
+  VERIFY_IS_APPROX(v[20], x);
+  v.resize(22,y);
+  VERIFY_IS_APPROX(v[21], y);
+  v.push_back(x);
+  VERIFY_IS_APPROX(v[22], x);
+  VERIFY((size_t)&(v[22]) == (size_t)&(v[21]) + sizeof(QuaternionType));
+
+  // do a lot of push_back such that the vector gets internally resized
+  // (with memory reallocation)
+  QuaternionType* ref = &w[0];
+  for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i)
+    v.push_back(w[i%w.size()]);
+  for(unsigned int i=23; i<v.size(); ++i)
+  {
+    VERIFY(v[i].coeffs()==w[(i-23)%w.size()].coeffs());
+  }
+}
+
+void test_newstdvector()
+{
+  // some non vectorizable fixed sizes
+  CALL_SUBTEST(check_stdvector_matrix(Vector2f()));
+  CALL_SUBTEST(check_stdvector_matrix(Matrix3f()));
+  CALL_SUBTEST(check_stdvector_matrix(Matrix3d()));
+
+  // some vectorizable fixed sizes
+  CALL_SUBTEST(check_stdvector_matrix(Matrix2f()));
+  CALL_SUBTEST(check_stdvector_matrix(Vector4f()));
+  CALL_SUBTEST(check_stdvector_matrix(Matrix4f()));
+  CALL_SUBTEST(check_stdvector_matrix(Matrix4d()));
+
+  // some dynamic sizes
+  CALL_SUBTEST(check_stdvector_matrix(MatrixXd(1,1)));
+  CALL_SUBTEST(check_stdvector_matrix(VectorXd(20)));
+  CALL_SUBTEST(check_stdvector_matrix(RowVectorXf(20)));
+  CALL_SUBTEST(check_stdvector_matrix(MatrixXcf(10,10)));
+
+  // some Transform
+  CALL_SUBTEST(check_stdvector_transform(Transform2f()));
+  CALL_SUBTEST(check_stdvector_transform(Transform3f()));
+  CALL_SUBTEST(check_stdvector_transform(Transform3d()));
+  //CALL_SUBTEST(check_stdvector_transform(Transform4d()));
+
+  // some Quaternion
+  CALL_SUBTEST(check_stdvector_quaternion(Quaternionf()));
+  CALL_SUBTEST(check_stdvector_quaternion(Quaterniond()));
+}

test/swap.cpp

@@ -0,0 +1,98 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 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/>.
+
+#define EIGEN_NO_STATIC_ASSERT
+#include "main.h"
+
+template<typename T>
+struct other_matrix_type
+{
+  typedef int type;
+};
+
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+struct other_matrix_type<Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> >
+{
+  typedef Matrix<_Scalar, _Rows, _Cols, _Options^RowMajor, _MaxRows, _MaxCols> type;
+};
+
+template<typename MatrixType> void swap(const MatrixType& m)
+{
+  typedef typename other_matrix_type<MatrixType>::type OtherMatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+
+  ei_assert((!ei_is_same_type<MatrixType,OtherMatrixType>::ret));
+  int rows = m.rows();
+  int cols = m.cols();
+  
+  // construct 3 matrix guaranteed to be distinct
+  MatrixType m1 = MatrixType::Random(rows,cols);
+  MatrixType m2 = MatrixType::Random(rows,cols) + Scalar(100) * MatrixType::Identity(rows,cols);
+  OtherMatrixType m3 = OtherMatrixType::Random(rows,cols) + Scalar(200) * OtherMatrixType::Identity(rows,cols);
+  
+  MatrixType m1_copy = m1;
+  MatrixType m2_copy = m2;
+  OtherMatrixType m3_copy = m3;
+  
+  // test swapping 2 matrices of same type
+  m1.swap(m2);
+  VERIFY_IS_APPROX(m1,m2_copy);
+  VERIFY_IS_APPROX(m2,m1_copy);
+  m1 = m1_copy;
+  m2 = m2_copy;
+  
+  // test swapping 2 matrices of different types
+  m1.swap(m3);
+  VERIFY_IS_APPROX(m1,m3_copy);
+  VERIFY_IS_APPROX(m3,m1_copy);
+  m1 = m1_copy;
+  m3 = m3_copy;
+  
+  // test swapping matrix with expression
+  m1.swap(m2.block(0,0,rows,cols));
+  VERIFY_IS_APPROX(m1,m2_copy);
+  VERIFY_IS_APPROX(m2,m1_copy);
+  m1 = m1_copy;
+  m2 = m2_copy;
+
+  // test swapping two expressions of different types
+  m1.transpose().swap(m3.transpose());
+  VERIFY_IS_APPROX(m1,m3_copy);
+  VERIFY_IS_APPROX(m3,m1_copy);
+  m1 = m1_copy;
+  m3 = m3_copy;
+  
+  // test assertion on mismatching size -- matrix case
+  VERIFY_RAISES_ASSERT(m1.swap(m1.row(0)));
+  // test assertion on mismatching size -- xpr case
+  VERIFY_RAISES_ASSERT(m1.row(0).swap(m1));
+}
+
+void test_swap()
+{
+  CALL_SUBTEST( swap(Matrix3f()) ); // fixed size, no vectorization 
+  CALL_SUBTEST( swap(Matrix4d()) ); // fixed size, possible vectorization 
+  CALL_SUBTEST( swap(MatrixXd(3,3)) ); // dyn size, no vectorization 
+  CALL_SUBTEST( swap(MatrixXf(30,30)) ); // dyn size, possible vectorization 
+}

test/visitor.cpp

@@ -0,0 +1,131 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+
+template<typename MatrixType> void matrixVisitor(const MatrixType& p)
+{
+  typedef typename MatrixType::Scalar Scalar;
+
+  int rows = p.rows();
+  int cols = p.cols();
+
+  // construct a random matrix where all coefficients are different
+  MatrixType m;
+  m = MatrixType::Random(rows, cols);
+  for(int i = 0; i < m.size(); i++)
+    for(int i2 = 0; i2 < i; i2++)
+      while(m(i) == m(i2)) // yes, ==
+        m(i) = ei_random<Scalar>();
+  
+  Scalar minc = Scalar(1000), maxc = Scalar(-1000);
+  int minrow,mincol,maxrow,maxcol;
+  for(int j = 0; j < cols; j++)
+  for(int i = 0; i < rows; i++)
+  {
+    if(m(i,j) < minc)
+    {
+      minc = m(i,j);
+      minrow = i;
+      mincol = j;
+    }
+    if(m(i,j) > maxc)
+    {
+      maxc = m(i,j);
+      maxrow = i;
+      maxcol = j;
+    }
+  }
+  int eigen_minrow, eigen_mincol, eigen_maxrow, eigen_maxcol;
+  Scalar eigen_minc, eigen_maxc;
+  eigen_minc = m.minCoeff(&eigen_minrow,&eigen_mincol);
+  eigen_maxc = m.maxCoeff(&eigen_maxrow,&eigen_maxcol);
+  VERIFY(minrow == eigen_minrow);
+  VERIFY(maxrow == eigen_maxrow);
+  VERIFY(mincol == eigen_mincol);
+  VERIFY(maxcol == eigen_maxcol);
+  VERIFY_IS_APPROX(minc, eigen_minc);
+  VERIFY_IS_APPROX(maxc, eigen_maxc);
+  VERIFY_IS_APPROX(minc, m.minCoeff());
+  VERIFY_IS_APPROX(maxc, m.maxCoeff());
+}
+
+template<typename VectorType> void vectorVisitor(const VectorType& w)
+{
+  typedef typename VectorType::Scalar Scalar;
+
+  int size = w.size();
+
+  // construct a random vector where all coefficients are different
+  VectorType v;
+  v = VectorType::Random(size);
+  for(int i = 0; i < size; i++)
+    for(int i2 = 0; i2 < i; i2++)
+      while(v(i) == v(i2)) // yes, ==
+        v(i) = ei_random<Scalar>();
+  
+  Scalar minc = Scalar(1000), maxc = Scalar(-1000);
+  int minidx,maxidx;
+  for(int i = 0; i < size; i++)
+  {
+    if(v(i) < minc)
+    {
+      minc = v(i);
+      minidx = i;
+    }
+    if(v(i) > maxc)
+    {
+      maxc = v(i);
+      maxidx = i;
+    }
+  }
+  int eigen_minidx, eigen_maxidx;
+  Scalar eigen_minc, eigen_maxc;
+  eigen_minc = v.minCoeff(&eigen_minidx);
+  eigen_maxc = v.maxCoeff(&eigen_maxidx);
+  VERIFY(minidx == eigen_minidx);
+  VERIFY(maxidx == eigen_maxidx);
+  VERIFY_IS_APPROX(minc, eigen_minc);
+  VERIFY_IS_APPROX(maxc, eigen_maxc);
+  VERIFY_IS_APPROX(minc, v.minCoeff());
+  VERIFY_IS_APPROX(maxc, v.maxCoeff());
+}
+
+void test_visitor()
+{
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST( matrixVisitor(Matrix<float, 1, 1>()) );
+    CALL_SUBTEST( matrixVisitor(Matrix2f()) );
+    CALL_SUBTEST( matrixVisitor(Matrix4d()) );
+    CALL_SUBTEST( matrixVisitor(MatrixXd(8, 12)) );
+    CALL_SUBTEST( matrixVisitor(Matrix<double,Dynamic,Dynamic,RowMajor>(20, 20)) );
+    CALL_SUBTEST( matrixVisitor(MatrixXi(8, 12)) );
+  }
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST( vectorVisitor(Vector4f()) );
+    CALL_SUBTEST( vectorVisitor(VectorXd(10)) );
+    CALL_SUBTEST( vectorVisitor(RowVectorXd(10)) );
+    CALL_SUBTEST( vectorVisitor(VectorXf(33)) );
+  }
+}