bump

conditional jump. This makes Part considered an "expensive" xpr
 that must be evaluated in operations such as Product.
This fixes bug #80.

CMakeLists.txt

@@ -7,7 +7,7 @@ set(INCLUDE_INSTALL_DIR
     "The directory where we install the header files"
     FORCE)
 
-set(EIGEN_VERSION_NUMBER "2.0.9")
+set(EIGEN_VERSION_NUMBER "2.0.11")
 set(EIGEN_VERSION "${EIGEN_VERSION_NUMBER}")
 
 set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
@@ -88,7 +88,7 @@ include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR})
 if(EIGEN_BUILD_PKGCONFIG)
     configure_file(eigen2.pc.in eigen2.pc) # uses INCLUDE_INSTALL_DIR
     install(FILES ${CMAKE_CURRENT_BINARY_DIR}/eigen2.pc
-        DESTINATION lib/pkgconfig
+        DESTINATION share/pkgconfig
         )
 endif(EIGEN_BUILD_PKGCONFIG)
 

Eigen/src/Core/Block.h

@@ -222,7 +222,7 @@ class Block<MatrixType,BlockRows,BlockCols,PacketAccess,HasDirectAccess>
 
     class InnerIterator;
     typedef typename ei_traits<Block>::AlignedDerivedType AlignedDerivedType;
-    friend class Block<MatrixType,BlockRows,BlockCols,PacketAccess==AsRequested?ForceAligned:AsRequested,HasDirectAccess>;
+    friend class Block<MatrixType,BlockRows,BlockCols,PacketAccess==int(AsRequested)?ForceAligned:AsRequested,HasDirectAccess>;
 
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Block)
 

Eigen/src/Core/Matrix.h

@@ -25,6 +25,11 @@
 #ifndef EIGEN_MATRIX_H
 #define EIGEN_MATRIX_H
 
+#ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO
+# define EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED for(int i=0;i<base().size();++i) coeffRef(i)=Scalar(0);
+#else
+# define EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
+#endif
 
 /** \class Matrix
   *
@@ -233,7 +238,14 @@ class Matrix
              && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows)
              && (MaxColsAtCompileTime == Dynamic || MaxColsAtCompileTime >= cols)
              && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols));
-      m_storage.resize(rows * cols, rows, cols);
+      #ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO
+        int size = rows*cols;
+        bool size_changed = size != this->size();
+        m_storage.resize(size, rows, cols);
+        if(size_changed) EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
+      #else
+        m_storage.resize(rows*cols, rows, cols);
+      #endif
     }
 
     /** Resizes \c *this to a vector of length \a size
@@ -243,10 +255,17 @@ class Matrix
     inline void resize(int size)
     {
       EIGEN_STATIC_ASSERT_VECTOR_ONLY(Matrix)
+      ei_assert(SizeAtCompileTime == Dynamic || SizeAtCompileTime == size);
+      #ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO
+        bool size_changed = size != this->size();
+      #endif
       if(RowsAtCompileTime == 1)
         m_storage.resize(size, 1, size);
       else
         m_storage.resize(size, size, 1);
+      #ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO
+        if(size_changed) EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
+      #endif
     }
 
     /** Copies the value of the expression \a other into \c *this with automatic resizing.
@@ -290,13 +309,14 @@ class Matrix
     EIGEN_STRONG_INLINE explicit Matrix() : m_storage()
     {
       _check_template_params();
+      EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
     }
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
     /** \internal */
     Matrix(ei_constructor_without_unaligned_array_assert)
       : m_storage(ei_constructor_without_unaligned_array_assert())
-    {}
+    {EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED}
 #endif
 
     /** Constructs a vector or row-vector with given dimension. \only_for_vectors
@@ -312,6 +332,7 @@ class Matrix
       EIGEN_STATIC_ASSERT_VECTOR_ONLY(Matrix)
       ei_assert(dim > 0);
       ei_assert(SizeAtCompileTime == Dynamic || SizeAtCompileTime == dim);
+      EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
     }
 
     /** This constructor has two very different behaviors, depending on the type of *this.
@@ -337,6 +358,7 @@ class Matrix
       {
         ei_assert(x > 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == x)
                && y > 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == y));
+        EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
       }
     }
     /** constructs an initialized 2D vector with given coefficients */
@@ -402,9 +424,10 @@ class Matrix
     void swap(const MatrixBase<OtherDerived>& other);
 
     /** \name Map
-      * These are convenience functions returning Map objects. The Map() static functions return unaligned Map objects,
-      * while the AlignedMap() functions return aligned Map objects and thus should be called only with 16-byte-aligned
-      * \a data pointers.
+      * These are convenience functions returning Map objects.
+      *
+      * \warning Do not use MapAligned in the Eigen 2.0. Mapping aligned arrays will be fully
+      * supported in Eigen 3.0 (already implemented in the development branch)
       *
       * \see class Map
       */

Eigen/src/Core/MatrixBase.h

@@ -136,12 +136,6 @@ template<typename Derived> class MatrixBase
           */
     };
 
-    /** Default constructor. Just checks at compile-time for self-consistency of the flags. */
-    MatrixBase()
-    {
-      ei_assert(ei_are_flags_consistent<Flags>::ret);
-    }
-
 #ifndef EIGEN_PARSED_BY_DOXYGEN
     /** This is the "real scalar" type; if the \a Scalar type is already real numbers
       * (e.g. int, float or double) then \a RealScalar is just the same as \a Scalar. If
@@ -165,7 +159,7 @@ template<typename Derived> class MatrixBase
     inline int size() const { return rows() * cols(); }
     /** \returns the number of nonzero coefficients which is in practice the number
       * of stored coefficients. */
-    inline int nonZeros() const { return derived.nonZeros(); }
+    inline int nonZeros() const { return size(); }
     /** \returns true if either the number of rows or the number of columns is equal to 1.
       * In other words, this function returns
       * \code rows()==1 || cols()==1 \endcode
@@ -627,6 +621,24 @@ template<typename Derived> class MatrixBase
     #ifdef EIGEN_MATRIXBASE_PLUGIN
     #include EIGEN_MATRIXBASE_PLUGIN
     #endif
+
+  protected:
+    /** Default constructor. Do nothing. */
+    MatrixBase()
+    {
+      /* Just checks for self-consistency of the flags.
+       * Only do it when debugging Eigen, as this borders on paranoiac and could slow compilation down
+       */
+#ifdef EIGEN_INTERNAL_DEBUGGING
+      EIGEN_STATIC_ASSERT(ei_are_flags_consistent<Flags>::ret,
+                          INVALID_MATRIXBASE_TEMPLATE_PARAMETERS)
+#endif
+    }
+
+  private:
+    explicit MatrixBase(int);
+    MatrixBase(int,int);
+    template<typename OtherDerived> explicit MatrixBase(const MatrixBase<OtherDerived>&);
 };
 
 #endif // EIGEN_MATRIXBASE_H

Eigen/src/Core/Part.h

@@ -50,7 +50,7 @@ struct ei_traits<Part<MatrixType, Mode> > : ei_traits<MatrixType>
   typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
   enum {
     Flags = (_MatrixTypeNested::Flags & (HereditaryBits) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit))) | Mode,
-    CoeffReadCost = _MatrixTypeNested::CoeffReadCost
+    CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ConditionalJumpCost
   };
 };
 

Eigen/src/Core/Transpose.h

@@ -69,7 +69,6 @@ template<typename MatrixType> class Transpose
 
     inline int rows() const { return m_matrix.cols(); }
     inline int cols() const { return m_matrix.rows(); }
-    inline int nonZeros() const { return m_matrix.nonZeros(); }
     inline int stride(void) const { return m_matrix.stride(); }
 
     inline Scalar& coeffRef(int row, int col)

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

@@ -114,9 +114,15 @@ template<> EIGEN_STRONG_INLINE void ei_pstoreu<float>(float*  to, const __m128&
 template<> EIGEN_STRONG_INLINE void ei_pstoreu<double>(double* to, const __m128d& from) { _mm_storeu_pd(to, from); }
 template<> EIGEN_STRONG_INLINE void ei_pstoreu<int>(int*    to, const __m128i& from) { _mm_storeu_si128(reinterpret_cast<__m128i*>(to), from); }
 
-#ifdef _MSC_VER
-// this fix internal compilation error
-template<> EIGEN_STRONG_INLINE float  ei_pfirst<__m128>(const __m128&  a) { float x = _mm_cvtss_f32(a); return x; }
+#if defined(_MSC_VER) && (_MSC_VER <= 1500) && defined(_WIN64)
+// The temporary variable fixes an internal compilation error.
+// Direct of the struct members fixed bug #62.
+template<> EIGEN_STRONG_INLINE float  ei_pfirst<__m128>(const __m128& a) { return a.m128_f32[0]; }
+template<> EIGEN_STRONG_INLINE double ei_pfirst<__m128d>(const __m128d& a) { return a.m128d_f64[0]; }
+template<> EIGEN_STRONG_INLINE int    ei_pfirst<__m128i>(const __m128i& a) { int x = _mm_cvtsi128_si32(a); return x; }
+#elif defined(_MSC_VER) && (_MSC_VER <= 1500)
+// The temporary variable fixes an internal compilation error.
+template<> EIGEN_STRONG_INLINE float  ei_pfirst<__m128>(const __m128& a) { float x = _mm_cvtss_f32(a); return x; }
 template<> EIGEN_STRONG_INLINE double ei_pfirst<__m128d>(const __m128d& a) { double x = _mm_cvtsd_f64(a); return x; }
 template<> EIGEN_STRONG_INLINE int    ei_pfirst<__m128i>(const __m128i& a) { int x = _mm_cvtsi128_si32(a); return x; }
 #else

Eigen/src/Core/util/Macros.h

@@ -30,7 +30,7 @@
 
 #define EIGEN_WORLD_VERSION 2
 #define EIGEN_MAJOR_VERSION 0
-#define EIGEN_MINOR_VERSION 9
+#define EIGEN_MINOR_VERSION 11
 
 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
                                       (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \
@@ -39,7 +39,7 @@
 // 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__))
+#if defined(__GNUC__) && !(defined(__i386__) || defined(__x86_64__) || defined(__powerpc__) || defined(__ppc__) || defined(__ia64__))
 #define EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT 1
 #else
 #define EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT 0

Eigen/src/Core/util/Memory.h

@@ -209,18 +209,47 @@ template<typename T, bool Align> inline void ei_conditional_aligned_delete(T *pt
   ei_conditional_aligned_free<Align>(ptr);
 }
 
-/** \internal \returns the number of elements which have to be skipped such that data are 16 bytes aligned */
-template<typename Scalar>
-inline static int ei_alignmentOffset(const Scalar* ptr, int maxOffset)
+/** \internal \returns the index of the first element of the array that is well aligned for vectorization.
+  *
+  * \param array the address of the start of the array
+  * \param size the size of the array
+  *
+  * \note If no element of the array is well aligned, the size of the array is returned. Typically,
+  * for example with SSE, "well aligned" means 16-byte-aligned. If vectorization is disabled or if the
+  * packet size for the given scalar type is 1, then everything is considered well-aligned.
+  *
+  * \note If the scalar type is vectorizable, we rely on the following assumptions: sizeof(Scalar) is a
+  * power of 2, the packet size in bytes is also a power of 2, and is a multiple of sizeof(Scalar). On the
+  * other hand, we do not assume that the array address is a multiple of sizeof(Scalar), as that fails for
+  * example with Scalar=double on certain 32-bit platforms, see bug #79.
+  *
+  * There is also the variant ei_first_aligned(const MatrixBase&, Integer) defined in Coeffs.h.
+  */
+template<typename Scalar, typename Integer>
+inline static Integer ei_alignmentOffset(const Scalar* array, Integer size)
 {
   typedef typename ei_packet_traits<Scalar>::type Packet;
-  const int PacketSize = ei_packet_traits<Scalar>::size;
-  const int PacketAlignedMask = PacketSize-1;
-  const bool Vectorized = PacketSize>1;
-  return Vectorized
-          ? std::min<int>( (PacketSize - (int((size_t(ptr)/sizeof(Scalar))) & PacketAlignedMask))
-                           & PacketAlignedMask, maxOffset)
-          : 0;
+  enum { PacketSize = ei_packet_traits<Scalar>::size,
+         PacketAlignedMask = PacketSize-1
+  };
+  
+  if(PacketSize==1)
+  {
+    // Either there is no vectorization, or a packet consists of exactly 1 scalar so that all elements
+    // of the array have the same aligment.
+    return 0;
+  }
+  else if(size_t(array) & (sizeof(Scalar)-1))
+  {
+    // There is vectorization for this scalar type, but the array is not aligned to the size of a single scalar.
+    // Consequently, no element of the array is well aligned.
+    return size;
+  }
+  else
+  {
+    return std::min<Integer>( (PacketSize - (Integer((size_t(array)/sizeof(Scalar))) & PacketAlignedMask))
+                           & PacketAlignedMask, size);
+  }
 }
 
 /** \internal

Eigen/src/Core/util/StaticAssert.h

@@ -74,6 +74,7 @@
         THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES,
         THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES,
         INVALID_MATRIX_TEMPLATE_PARAMETERS,
+        INVALID_MATRIXBASE_TEMPLATE_PARAMETERS,
         BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER,
         THIS_METHOD_IS_ONLY_FOR_DIAGONAL_MATRIX
       };

Eigen/src/Geometry/Hyperplane.h

@@ -52,9 +52,9 @@ public:
   typedef _Scalar Scalar;
   typedef typename NumTraits<Scalar>::Real RealScalar;
   typedef Matrix<Scalar,AmbientDimAtCompileTime,1> VectorType;
-  typedef Matrix<Scalar,AmbientDimAtCompileTime==Dynamic
+  typedef Matrix<Scalar,int(AmbientDimAtCompileTime)==Dynamic
                         ? Dynamic
-                        : AmbientDimAtCompileTime+1,1> Coefficients;
+                        : int(AmbientDimAtCompileTime)+1,1> Coefficients;
   typedef Block<Coefficients,AmbientDimAtCompileTime,1> NormalReturnType;
 
   /** Default constructor without initialization */

Eigen/src/Geometry/Quaternion.h

@@ -450,22 +450,31 @@ inline Scalar Quaternion<Scalar>::angularDistance(const Quaternion& other) const
 template <typename Scalar>
 Quaternion<Scalar> Quaternion<Scalar>::slerp(Scalar t, const Quaternion& other) const
 {
-  static const Scalar one = Scalar(1) - precision<Scalar>();
+  static const Scalar one = Scalar(1) - machine_epsilon<Scalar>();
   Scalar d = this->dot(other);
   Scalar absD = ei_abs(d);
+
+  Scalar scale0;
+  Scalar scale1;
+
   if (absD>=one)
-    return *this;
+  {
+    scale0 = Scalar(1) - t;
+    scale1 = t;
+  }
+  else
+  {
+    // theta is the angle between the 2 quaternions
+    Scalar theta = std::acos(absD);
+    Scalar sinTheta = ei_sin(theta);
 
-  // theta is the angle between the 2 quaternions
-  Scalar theta = std::acos(absD);
-  Scalar sinTheta = ei_sin(theta);
+    scale0 = ei_sin( ( Scalar(1) - t ) * theta) / sinTheta;
+    scale1 = ei_sin( ( t * theta) ) / sinTheta;
+    if (d<0)
+      scale1 = -scale1;
+  }
 
-  Scalar scale0 = ei_sin( ( Scalar(1) - t ) * theta) / sinTheta;
-  Scalar scale1 = ei_sin( ( t * theta) ) / sinTheta;
-  if (d<0)
-    scale1 = -scale1;
-
-  return Quaternion(scale0 * m_coeffs + scale1 * other.m_coeffs);
+  return Quaternion<Scalar>(scale0 * coeffs() + scale1 * other.coeffs());
 }
 
 // set from a rotation matrix

Eigen/src/LU/Inverse.h

@@ -29,8 +29,8 @@
 *** Part 1 : optimized implementations for fixed-size 2,3,4 cases ***
 ********************************************************************/
 
-template<typename MatrixType>
-void ei_compute_inverse_in_size2_case(const MatrixType& matrix, MatrixType* result)
+template<typename XprType, typename MatrixType>
+void ei_compute_inverse_in_size2_case(const XprType& matrix, MatrixType* result)
 {
   typedef typename MatrixType::Scalar Scalar;
   const Scalar invdet = Scalar(1) / matrix.determinant();
@@ -75,91 +75,155 @@ void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* resu
   result->coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet;
 }
 
-template<typename MatrixType>
-bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixType* result)
+template<typename MatrixType, typename Scalar = typename MatrixType::Scalar>
+struct ei_compute_inverse_in_size4_case
 {
-  /* Let's split M into four 2x2 blocks:
-    * (P Q)
-    * (R S)
-    * If P is invertible, with inverse denoted by P_inverse, and if
-    * (S - R*P_inverse*Q) is also invertible, then the inverse of M is
-    * (P' Q')
-    * (R' S')
-    * where
-    * S' = (S - R*P_inverse*Q)^(-1)
-    * P' = P1 + (P1*Q) * S' *(R*P_inverse)
-    * Q' = -(P_inverse*Q) * S'
-    * R' = -S' * (R*P_inverse)
-    */
-  typedef Block<MatrixType,2,2> XprBlock22;
-  typedef typename MatrixBase<XprBlock22>::PlainMatrixType Block22;
-  Block22 P_inverse;
-  if(ei_compute_inverse_in_size2_case_with_check(matrix.template block<2,2>(0,0), &P_inverse))
+  static void run(const MatrixType& matrix, MatrixType& result)
   {
-    const Block22 Q = matrix.template block<2,2>(0,2);
-    const Block22 P_inverse_times_Q = P_inverse * Q;
-    const XprBlock22 R = matrix.template block<2,2>(2,0);
-    const Block22 R_times_P_inverse = R * P_inverse;
-    const Block22 R_times_P_inverse_times_Q = R_times_P_inverse * Q;
-    const XprBlock22 S = matrix.template block<2,2>(2,2);
-    const Block22 X = S - R_times_P_inverse_times_Q;
-    Block22 Y;
-    ei_compute_inverse_in_size2_case(X, &Y);
-    result->template block<2,2>(2,2) = Y;
-    result->template block<2,2>(2,0) = - Y * R_times_P_inverse;
-    const Block22 Z = P_inverse_times_Q * Y;
-    result->template block<2,2>(0,2) = - Z;
-    result->template block<2,2>(0,0) = P_inverse + Z * R_times_P_inverse;
-    return true;
+    result.coeffRef(0,0) = matrix.minor(0,0).determinant();
+    result.coeffRef(1,0) = -matrix.minor(0,1).determinant();
+    result.coeffRef(2,0) = matrix.minor(0,2).determinant();
+    result.coeffRef(3,0) = -matrix.minor(0,3).determinant();
+    result.coeffRef(0,2) = matrix.minor(2,0).determinant();
+    result.coeffRef(1,2) = -matrix.minor(2,1).determinant();
+    result.coeffRef(2,2) = matrix.minor(2,2).determinant();
+    result.coeffRef(3,2) = -matrix.minor(2,3).determinant();
+    result.coeffRef(0,1) = -matrix.minor(1,0).determinant();
+    result.coeffRef(1,1) = matrix.minor(1,1).determinant();
+    result.coeffRef(2,1) = -matrix.minor(1,2).determinant();
+    result.coeffRef(3,1) = matrix.minor(1,3).determinant();
+    result.coeffRef(0,3) = -matrix.minor(3,0).determinant();
+    result.coeffRef(1,3) = matrix.minor(3,1).determinant();
+    result.coeffRef(2,3) = -matrix.minor(3,2).determinant();
+    result.coeffRef(3,3) = matrix.minor(3,3).determinant();
+    result /= (matrix.col(0).cwise()*result.row(0).transpose()).sum();
   }
-  else
-  {
-    return false;
-  }
-}
+};
 
+#ifdef EIGEN_VECTORIZE_SSE
+// The SSE code for the 4x4 float matrix inverse in this file comes from the file
+//   ftp://download.intel.com/design/PentiumIII/sml/24504301.pdf
+// its copyright information is:
+//   Copyright (C) 1999 Intel Corporation
+// See page ii of that document for legal stuff. Not being lawyers, we just assume
+// here that if Intel makes this document publically available, with source code
+// and detailed explanations, it's because they want their CPUs to be fed with
+// good code, and therefore they presumably don't mind us using it in Eigen.
 template<typename MatrixType>
-void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* result)
+struct ei_compute_inverse_in_size4_case<MatrixType, float>
 {
-  if(ei_compute_inverse_in_size4_case_helper(matrix, result))
+  static void run(const MatrixType& matrix, MatrixType& result)
   {
-    // good ! The topleft 2x2 block was invertible, so the 2x2 blocks approach is successful.
-    return;
+    // Variables (Streaming SIMD Extensions registers) which will contain cofactors and, later, the
+    // lines of the inverted matrix.
+    __m128 minor0, minor1, minor2, minor3;
+
+    // Variables which will contain the lines of the reference matrix and, later (after the transposition),
+    // the columns of the original matrix.
+    __m128 row0, row1, row2, row3;
+
+    // Temporary variables and the variable that will contain the matrix determinant.
+    __m128 det, tmp1;
+
+    // Matrix transposition
+    const float *src = matrix.data();
+    tmp1  = _mm_loadh_pi(_mm_castpd_ps(_mm_load_sd((double*)src)), (__m64*)(src+ 4));
+    row1  = _mm_loadh_pi(_mm_castpd_ps(_mm_load_sd((double*)(src+8))), (__m64*)(src+12));
+    row0  = _mm_shuffle_ps(tmp1, row1, 0x88);
+    row1  = _mm_shuffle_ps(row1, tmp1, 0xDD);
+    tmp1  = _mm_loadh_pi(_mm_castpd_ps(_mm_load_sd((double*)(src+ 2))), (__m64*)(src+ 6));
+    row3  = _mm_loadh_pi(_mm_castpd_ps(_mm_load_sd((double*)(src+10))), (__m64*)(src+14));
+    row2  = _mm_shuffle_ps(tmp1, row3, 0x88);
+    row3  = _mm_shuffle_ps(row3, tmp1, 0xDD);
+
+
+    // Cofactors calculation. Because in the process of cofactor computation some pairs in three-
+    // element products are repeated, it is not reasonable to load these pairs anew every time. The
+    // values in the registers with these pairs are formed using shuffle instruction. Cofactors are
+    // calculated row by row (4 elements are placed in 1 SP FP SIMD floating point register).
+
+    tmp1   = _mm_mul_ps(row2, row3);
+    tmp1   = _mm_shuffle_ps(tmp1, tmp1, 0xB1);
+    minor0  = _mm_mul_ps(row1, tmp1);
+    minor1  = _mm_mul_ps(row0, tmp1);
+    tmp1   = _mm_shuffle_ps(tmp1, tmp1, 0x4E);
+    minor0  = _mm_sub_ps(_mm_mul_ps(row1, tmp1), minor0);
+    minor1  = _mm_sub_ps(_mm_mul_ps(row0, tmp1), minor1);
+    minor1  = _mm_shuffle_ps(minor1, minor1, 0x4E);
+    //    -----------------------------------------------
+    tmp1   = _mm_mul_ps(row1, row2);
+    tmp1   = _mm_shuffle_ps(tmp1, tmp1, 0xB1);
+    minor0  = _mm_add_ps(_mm_mul_ps(row3, tmp1), minor0);
+    minor3  = _mm_mul_ps(row0, tmp1);
+    tmp1   = _mm_shuffle_ps(tmp1, tmp1, 0x4E);
+    minor0  = _mm_sub_ps(minor0, _mm_mul_ps(row3, tmp1));
+    minor3  = _mm_sub_ps(_mm_mul_ps(row0, tmp1), minor3);
+    minor3  = _mm_shuffle_ps(minor3, minor3, 0x4E);
+    //    -----------------------------------------------
+    tmp1   = _mm_mul_ps(_mm_shuffle_ps(row1, row1, 0x4E), row3);
+    tmp1   = _mm_shuffle_ps(tmp1, tmp1, 0xB1);
+    row2   = _mm_shuffle_ps(row2, row2, 0x4E);
+    minor0  = _mm_add_ps(_mm_mul_ps(row2, tmp1), minor0);
+    minor2  = _mm_mul_ps(row0, tmp1);
+    tmp1   = _mm_shuffle_ps(tmp1, tmp1, 0x4E);
+    minor0  = _mm_sub_ps(minor0, _mm_mul_ps(row2, tmp1));
+    minor2  = _mm_sub_ps(_mm_mul_ps(row0, tmp1), minor2);
+    minor2  = _mm_shuffle_ps(minor2, minor2, 0x4E);
+    //    -----------------------------------------------
+    tmp1   = _mm_mul_ps(row0, row1);
+    tmp1   = _mm_shuffle_ps(tmp1, tmp1, 0xB1);
+    minor2 = _mm_add_ps(_mm_mul_ps(row3, tmp1), minor2);
+    minor3 = _mm_sub_ps(_mm_mul_ps(row2, tmp1), minor3);
+    tmp1   = _mm_shuffle_ps(tmp1, tmp1, 0x4E);
+    minor2 = _mm_sub_ps(_mm_mul_ps(row3, tmp1), minor2);
+    minor3 = _mm_sub_ps(minor3, _mm_mul_ps(row2, tmp1));
+    //           -----------------------------------------------
+    tmp1   = _mm_mul_ps(row0, row3);
+    tmp1   = _mm_shuffle_ps(tmp1, tmp1, 0xB1);
+    minor1 = _mm_sub_ps(minor1, _mm_mul_ps(row2, tmp1));
+    minor2 = _mm_add_ps(_mm_mul_ps(row1, tmp1), minor2);
+    tmp1   = _mm_shuffle_ps(tmp1, tmp1, 0x4E);
+    minor1 = _mm_add_ps(_mm_mul_ps(row2, tmp1), minor1);
+    minor2 = _mm_sub_ps(minor2, _mm_mul_ps(row1, tmp1));
+    //           -----------------------------------------------
+    tmp1   = _mm_mul_ps(row0, row2);
+    tmp1   = _mm_shuffle_ps(tmp1, tmp1, 0xB1);
+    minor1 = _mm_add_ps(_mm_mul_ps(row3, tmp1), minor1);
+    minor3 = _mm_sub_ps(minor3, _mm_mul_ps(row1, tmp1));
+    tmp1   = _mm_shuffle_ps(tmp1, tmp1, 0x4E);
+    minor1 = _mm_sub_ps(minor1, _mm_mul_ps(row3, tmp1));
+    minor3 = _mm_add_ps(_mm_mul_ps(row1, tmp1), minor3);
+
+    // Evaluation of determinant and its reciprocal value. In the original Intel document,
+    // 1/det was evaluated using a fast rcpps command with subsequent approximation using
+    // the Newton-Raphson algorithm. Here, we go for a IEEE-compliant division instead,
+    // so as to not compromise precision at all.
+    det    = _mm_mul_ps(row0, minor0);
+    det    = _mm_add_ps(_mm_shuffle_ps(det, det, 0x4E), det);
+    det    = _mm_add_ss(_mm_shuffle_ps(det, det, 0xB1), det);
+//     tmp1= _mm_rcp_ss(det);
+//     det= _mm_sub_ss(_mm_add_ss(tmp1, tmp1), _mm_mul_ss(det, _mm_mul_ss(tmp1, tmp1)));
+    det    = _mm_div_ss(_mm_set_ss(1.0f), det); // <--- yay, one original line not copied from Intel
+    det    = _mm_shuffle_ps(det, det, 0x00);
+    // warning, Intel's variable naming is very confusing: now 'det' is 1/det !
+
+    // Multiplication of cofactors by 1/det. Storing the inverse matrix to the address in pointer src.
+    minor0 = _mm_mul_ps(det, minor0);
+    float *dst = result.data();
+    _mm_storel_pi((__m64*)(dst), minor0);
+    _mm_storeh_pi((__m64*)(dst+2), minor0);
+    minor1 = _mm_mul_ps(det, minor1);
+    _mm_storel_pi((__m64*)(dst+4), minor1);
+    _mm_storeh_pi((__m64*)(dst+6), minor1);
+    minor2 = _mm_mul_ps(det, minor2);
+    _mm_storel_pi((__m64*)(dst+ 8), minor2);
+    _mm_storeh_pi((__m64*)(dst+10), minor2);
+    minor3 = _mm_mul_ps(det, minor3);
+    _mm_storel_pi((__m64*)(dst+12), minor3);
+    _mm_storeh_pi((__m64*)(dst+14), minor3);
   }
-  else
-  {
-    // rare case: the topleft 2x2 block is not invertible (but the matrix itself is assumed to be).
-    // since this is a rare case, we don't need to optimize it. We just want to handle it with little
-    // additional code.
-    MatrixType m(matrix);
-    m.row(0).swap(m.row(2));
-    m.row(1).swap(m.row(3));
-    if(ei_compute_inverse_in_size4_case_helper(m, result))
-    {
-      // good, the topleft 2x2 block of m is invertible. Since m is different from matrix in that some
-      // rows were permuted, the actual inverse of matrix is derived from the inverse of m by permuting
-      // the corresponding columns.
-      result->col(0).swap(result->col(2));
-      result->col(1).swap(result->col(3));
-    }
-    else
-    {
-      // last possible case. Since matrix is assumed to be invertible, this last case has to work.
-      // first, undo the swaps previously made
-      m.row(0).swap(m.row(2));
-      m.row(1).swap(m.row(3));
-      // swap row 0 with the the row among 0 and 1 that has the biggest 2 first coeffs
-      int swap0with = ei_abs(m.coeff(0,0))+ei_abs(m.coeff(0,1))>ei_abs(m.coeff(1,0))+ei_abs(m.coeff(1,1)) ? 0 : 1;
-      m.row(0).swap(m.row(swap0with));
-      // swap row 1 with the the row among 2 and 3 that has the biggest 2 first coeffs
-      int swap1with = ei_abs(m.coeff(2,0))+ei_abs(m.coeff(2,1))>ei_abs(m.coeff(3,0))+ei_abs(m.coeff(3,1)) ? 2 : 3;
-      m.row(1).swap(m.row(swap1with));
-      ei_compute_inverse_in_size4_case_helper(m, result);
-      result->col(1).swap(result->col(swap1with));
-      result->col(0).swap(result->col(swap0with));
-    }
-  }
-}
+};
+#endif
 
 /***********************************************
 *** Part 2 : selector and MatrixBase methods ***
@@ -208,7 +272,7 @@ struct ei_compute_inverse<MatrixType, 4>
 {
   static inline void run(const MatrixType& matrix, MatrixType* result)
   {
-    ei_compute_inverse_in_size4_case(matrix, result);
+    ei_compute_inverse_in_size4_case<MatrixType>::run(matrix, *result);
   }
 };
 

Eigen/src/Sparse/TriangularSolver.h

@@ -43,8 +43,11 @@ struct ei_solve_triangular_selector<Lhs,Rhs,LowerTriangular,RowMajor|IsSparse>
         {
           lastVal = it.value();
           lastIndex = it.index();
+          if(lastIndex == i)
+            break;
           tmp -= lastVal * other.coeff(lastIndex,col);
         }
+
         if (Lhs::Flags & UnitDiagBit)
           other.coeffRef(i,col) = tmp;
         else

doc/QuickStartGuide.dox

@@ -229,17 +229,24 @@ Of course, fixed-size matrices can't be resized.
 
 
 \subsection TutorialMap Map
-Any memory buffer can be mapped as an Eigen expression:
-<table class="tutorial_code"><tr><td>
+Any memory buffer can be mapped as an Eigen expression using the Map() static method:
 \code
 std::vector<float> stlarray(10);
-Map<VectorXf>(&stlarray[0], stlarray.size()).setOnes();
-int data[4] = 1, 2, 3, 4;
-Matrix2i mat2x2(data);
-MatrixXi mat2x2 = Map<Matrix2i>(data);
-MatrixXi mat2x2 = Map<MatrixXi>(data,2,2);
+VectorXf::Map(&stlarray[0], stlarray.size()).squaredNorm();
+\endcode
+Here VectorXf::Map returns an object of class Map<VectorXf>, which behaves like a VectorXf except that it uses the existing array. You can write to this object, that will write to the existing array. You can also construct a named obtect to reuse it:
+\code
+float array[rows*cols];
+Map<MatrixXf> m(array,rows,cols);
+m = othermatrix1 * othermatrix2;
+m.eigenvalues();
+\endcode
+In the fixed-size case, no need to pass sizes:
+\code
+float array[9];
+Map<Matrix3d> m(array);
+Matrix3d::Map(array).setIdentity();
 \endcode
-</td></tr></table>
 
 
 

doc/StlContainers.dox

@@ -11,7 +11,7 @@ namespace Eigen {
 
 \section summary Executive summary
 
-Using STL containers on \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types" requires taking the following two steps:
+Using STL containers on \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types", or classes having members of such types, requires taking the following two steps:
 
 \li A 16-byte-aligned allocator must be used. Eigen does provide one ready for use: aligned_allocator.
 \li If you want to use the std::vector container, you need to \#include <Eigen/StdVector>.

doc/UnalignedArrayAssert.dox

@@ -55,16 +55,17 @@ Note that here, Eigen::Vector2d is only used as an example, more generally the i
 
 \section c2 Cause 2: STL Containers
 
-If you use STL Containers such as std::vector, std::map, ..., with Eigen objects, like this,
+If you use STL Containers such as std::vector, std::map, ..., with Eigen objects, or with classes containing Eigen objects, like this,
 
 \code
 std::vector<Eigen::Matrix2f> my_vector;
-std::map<int, Eigen::Matrix2f> my_map;
+struct my_class { ... Eigen::Matrix2f m; ... };
+std::map<int, my_class> my_map;
 \endcode
 
 then you need to read this separate page: \ref StlContainers "Using STL Containers with Eigen".
 
-Note that here, Eigen::Matrix2f is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types".
+Note that here, Eigen::Matrix2f is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types" and \ref StructHavingEigenMembers "structures having such Eigen objects as member".
 
 \section c3 Cause 3: Passing Eigen objects by value
 

scripts/eigen_gen_docs

@@ -11,7 +11,7 @@ USER=${USER:-'orzel'}
 # step 1 : build
 # todo if 'build is not there, create one:
 #mkdir build
-(cd build && cmake .. && make -j3 doc) || (echo "make failed"; exit 1)
+(cd build && cmake .. && make -j3 doc) || echo "make failed"; exit 1
 #todo: n+1 where n = number of cpus
 
 #step 2 : upload

test/CMakeLists.txt

@@ -180,6 +180,7 @@ ei_add_test(meta)
 ei_add_test(sizeof)
 ei_add_test(dynalloc)
 ei_add_test(nomalloc)
+ei_add_test(first_aligned)
 ei_add_test(mixingtypes)
 ei_add_test(packetmath)
 ei_add_test(unalignedassert)
@@ -222,6 +223,8 @@ ei_add_test(sparse_product)
 ei_add_test(swap)
 ei_add_test(visitor)
 
+ei_add_test(prec_inverse_4x4 ${EI_OFLAG})
+
 # print a summary of the different options
 message("************************************************************")
 message("***    Eigen's unit tests configuration summary          ***")

test/first_aligned.cpp

@@ -0,0 +1,64 @@
+// 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/>.
+
+#include "main.h"
+
+template<typename Scalar>
+void test_first_aligned_helper(Scalar *array, int size)
+{
+  const int packet_size = sizeof(Scalar) * ei_packet_traits<Scalar>::size;
+  VERIFY(((size_t(array) + sizeof(Scalar) * ei_alignmentOffset(array, size)) % packet_size) == 0);
+}
+
+template<typename Scalar>
+void test_none_aligned_helper(Scalar *array, int size)
+{
+  VERIFY(ei_packet_traits<Scalar>::size == 1 || ei_alignmentOffset(array, size) == size);
+}
+
+struct some_non_vectorizable_type { float x; };
+
+void test_first_aligned()
+{
+  EIGEN_ALIGN_128 float array_float[100];
+  test_first_aligned_helper(array_float, 50);
+  test_first_aligned_helper(array_float+1, 50);
+  test_first_aligned_helper(array_float+2, 50);
+  test_first_aligned_helper(array_float+3, 50);
+  test_first_aligned_helper(array_float+4, 50);
+  test_first_aligned_helper(array_float+5, 50);
+  
+  EIGEN_ALIGN_128 double array_double[100];
+  test_first_aligned_helper(array_float, 50);
+  test_first_aligned_helper(array_float+1, 50);
+  test_first_aligned_helper(array_float+2, 50);
+  
+  double *array_double_plus_4_bytes = (double*)(size_t(array_double)+4);
+  test_none_aligned_helper(array_double_plus_4_bytes, 50);
+  test_none_aligned_helper(array_double_plus_4_bytes+1, 50);
+  
+  some_non_vectorizable_type array_nonvec[100];
+  test_first_aligned_helper(array_nonvec, 100);
+  test_none_aligned_helper(array_nonvec, 100);
+}

test/prec_inverse_4x4.cpp

@@ -0,0 +1,99 @@
+// 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/>.
+
+#include "main.h"
+#include <Eigen/LU>
+#include <algorithm>
+
+template<typename T> std::string type_name() { return "other"; }
+template<> std::string type_name<float>() { return "float"; }
+template<> std::string type_name<double>() { return "double"; }
+template<> std::string type_name<int>() { return "int"; }
+template<> std::string type_name<std::complex<float> >() { return "complex<float>"; }
+template<> std::string type_name<std::complex<double> >() { return "complex<double>"; }
+template<> std::string type_name<std::complex<int> >() { return "complex<int>"; }
+
+#define EIGEN_DEBUG_VAR(x) std::cerr << #x << " = " << x << std::endl;
+
+template<typename T> inline typename NumTraits<T>::Real epsilon()
+{
+ return std::numeric_limits<typename NumTraits<T>::Real>::epsilon();
+}
+
+template<typename MatrixType> void inverse_permutation_4x4()
+{
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  Vector4i indices(0,1,2,3);
+  for(int i = 0; i < 24; ++i)
+  {
+    MatrixType m = MatrixType::Zero();
+    m(indices(0),0) = 1;
+    m(indices(1),1) = 1;
+    m(indices(2),2) = 1;
+    m(indices(3),3) = 1;
+    MatrixType inv = m.inverse();
+    double error = double( (m*inv-MatrixType::Identity()).norm() / epsilon<Scalar>() );
+    VERIFY(error == 0.0);
+    std::next_permutation(indices.data(),indices.data()+4);
+  }
+}
+
+template<typename MatrixType> void inverse_general_4x4(int repeat)
+{
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  double error_sum = 0., error_max = 0.;
+  for(int i = 0; i < repeat; ++i)
+  {
+    MatrixType m;
+    RealScalar absdet;
+    do {
+      m = MatrixType::Random();
+      absdet = ei_abs(m.determinant());
+    } while(absdet < 10 * epsilon<Scalar>());
+    MatrixType inv = m.inverse();
+    double error = double( (m*inv-MatrixType::Identity()).norm() * absdet / epsilon<Scalar>() );
+    error_sum += error;
+    error_max = std::max(error_max, error);
+  }
+  std::cerr << "inverse_general_4x4, Scalar = " << type_name<Scalar>() << std::endl;
+  double error_avg = error_sum / repeat;
+  EIGEN_DEBUG_VAR(error_avg);
+  EIGEN_DEBUG_VAR(error_max);
+  VERIFY(error_avg < (NumTraits<Scalar>::IsComplex ? 8.0 : 1.0));
+  VERIFY(error_max < (NumTraits<Scalar>::IsComplex ? 64.0 : 20.0));
+}
+
+void test_prec_inverse_4x4()
+{
+  CALL_SUBTEST((inverse_permutation_4x4<Matrix4f>()));
+  CALL_SUBTEST(( inverse_general_4x4<Matrix4f>(200000 * g_repeat) ));
+
+  CALL_SUBTEST((inverse_permutation_4x4<Matrix<double,4,4,RowMajor> >()));
+  CALL_SUBTEST(( inverse_general_4x4<Matrix<double,4,4,RowMajor> >(200000 * g_repeat) ));
+
+  CALL_SUBTEST((inverse_permutation_4x4<Matrix4cf>()));
+  CALL_SUBTEST((inverse_general_4x4<Matrix4cf>(50000 * g_repeat)));
+}

test/sparse_solvers.cpp

@@ -68,12 +68,20 @@ template<typename Scalar> void sparse_solvers(int rows, int cols)
     VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().solveTriangular(vec2),
                      m2.template marked<LowerTriangular>().solveTriangular(vec3));
 
+    // lower - transpose
+    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
+    VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().transpose().solveTriangular(vec2),
+                     m2.template marked<LowerTriangular>().transpose().solveTriangular(vec3));
+
     // upper
     initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
     VERIFY_IS_APPROX(refMat2.template marked<UpperTriangular>().solveTriangular(vec2),
                      m2.template marked<UpperTriangular>().solveTriangular(vec3));
 
-    // TODO test row major
+    // upper - transpose
+    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
+    VERIFY_IS_APPROX(refMat2.template marked<UpperTriangular>().transpose().solveTriangular(vec2),
+                     m2.template marked<UpperTriangular>().transpose().solveTriangular(vec3));
   }
 
   // test LLT

test/testsuite.cmake

@@ -39,7 +39,7 @@
 # VERSION=opensuse-11.1
 # WORK_DIR=/home/gael/Coding/eigen2/cdash
 # # get the last version of the script
-# wget http://bitbucket.org/eigen/eigen2/raw/tip/test/testsuite.cmake -o $WORK_DIR/testsuite.cmake
+# wget http://bitbucket.org/eigen/eigen/raw/tip/test/testsuite.cmake -o $WORK_DIR/testsuite.cmake
 # COMMON="ctest -S $WORK_DIR/testsuite.cmake,EIGEN_WORK_DIR=$WORK_DIR,EIGEN_SITE=$SITE,EIGEN_MODE=$1,EIGEN_BUILD_STRING=$OS_VERSION-$ARCH"
 # $COMMON-gcc-3.4.6,EIGEN_CXX=g++-3.4
 # $COMMON-gcc-4.0.1,EIGEN_CXX=g++-4.0.1
@@ -133,7 +133,7 @@ endif(NOT EIGEN_MODE)
 ## mandatory variables (the default should be ok in most cases):
 
 SET (CTEST_CVS_COMMAND "hg")
-SET (CTEST_CVS_CHECKOUT "${CTEST_CVS_COMMAND} clone -r 2.0 http://bitbucket.org/eigen/eigen2 \"${CTEST_SOURCE_DIRECTORY}\"")
+SET (CTEST_CVS_CHECKOUT "${CTEST_CVS_COMMAND} clone -r 2.0 http://bitbucket.org/eigen/eigen \"${CTEST_SOURCE_DIRECTORY}\"")
 
 # which ctest command to use for running the dashboard
 SET (CTEST_COMMAND "${EIGEN_CMAKE_DIR}ctest -D ${EIGEN_MODE}")