bump to 3.1.3

(transplanted from b5d8299ee7
)

CTestConfig.cmake

@@ -4,10 +4,10 @@
 ## # The following are required to uses Dart and the Cdash dashboard
 ##   ENABLE_TESTING()
 ##   INCLUDE(CTest)
-set(CTEST_PROJECT_NAME "Eigen")
+set(CTEST_PROJECT_NAME "Eigen3.1")
 set(CTEST_NIGHTLY_START_TIME "00:00:00 UTC")
 
 set(CTEST_DROP_METHOD "http")
 set(CTEST_DROP_SITE "manao.inria.fr")
-set(CTEST_DROP_LOCATION "/CDash/submit.php?project=Eigen")
-set(CTEST_DROP_SITE_CDASH TRUE)
+set(CTEST_DROP_LOCATION "/CDash/submit.php?project=Eigen3.1")
+set(CTEST_DROP_SITE_CDASH TRUE)

Eigen/Core

@@ -44,7 +44,7 @@
   #endif
 #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) )
+  #if (defined __SSE2__) && ( (!defined __GNUC__) || (defined __INTEL_COMPILER) || EIGEN_GNUC_AT_LEAST(4,2) )
     #define EIGEN_SSE2_ON_NON_MSVC_BUT_NOT_OLD_GCC
   #endif
 #endif

Eigen/src/Cholesky/LDLT.h

@@ -534,8 +534,7 @@ template<typename Derived>
 bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
 {
   eigen_assert(m_isInitialized && "LDLT is not initialized.");
-  const Index size = m_matrix.rows();
-  eigen_assert(size == bAndX.rows());
+  eigen_assert(m_matrix.rows() == bAndX.rows());
 
   bAndX = this->solve(bAndX);
 

Eigen/src/Core/CwiseUnaryView.h

@@ -44,9 +44,10 @@ struct traits<CwiseUnaryView<ViewOp, MatrixType> >
     // "error: no integral type can represent all of the enumerator values
     InnerStrideAtCompileTime = MatrixTypeInnerStride == Dynamic
                              ? int(Dynamic)
-                             : int(MatrixTypeInnerStride)
-                               * int(sizeof(typename traits<MatrixType>::Scalar) / sizeof(Scalar)),
-    OuterStrideAtCompileTime = outer_stride_at_compile_time<MatrixType>::ret
+                             : int(MatrixTypeInnerStride) * int(sizeof(typename traits<MatrixType>::Scalar) / sizeof(Scalar)),
+    OuterStrideAtCompileTime = outer_stride_at_compile_time<MatrixType>::ret == Dynamic
+                             ? int(Dynamic)
+                             : outer_stride_at_compile_time<MatrixType>::ret * int(sizeof(typename traits<MatrixType>::Scalar) / sizeof(Scalar))
   };
 };
 }
@@ -106,7 +107,7 @@ class CwiseUnaryViewImpl<ViewOp,MatrixType,Dense>
 
     inline Index outerStride() const
     {
-      return derived().nestedExpression().outerStride();
+      return derived().nestedExpression().outerStride() * sizeof(typename internal::traits<MatrixType>::Scalar) / sizeof(Scalar);
     }
 
     EIGEN_STRONG_INLINE CoeffReturnType coeff(Index row, Index col) const

Eigen/src/Core/DenseStorage.h

@@ -39,13 +39,24 @@ struct plain_array
   plain_array(constructor_without_unaligned_array_assert) {}
 };
 
-#ifdef EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT
+#if defined(EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT)
   #define EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(sizemask)
+#elif EIGEN_GNUC_AT_LEAST(4,7) 
+  // GCC 4.7 is too aggressive in its optimizations and remove the alignement test based on the fact the array is declared to be aligned.
+  // See this bug report: http://gcc.gnu.org/bugzilla/show_bug.cgi?id=53900
+  // Hiding the origin of the array pointer behind a function argument seems to do the trick even if the function is inlined:
+  template<typename PtrType>
+  EIGEN_ALWAYS_INLINE PtrType eigen_unaligned_array_assert_workaround_gcc47(PtrType array) { return array; }
+  #define EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(sizemask) \
+    eigen_assert((reinterpret_cast<size_t>(eigen_unaligned_array_assert_workaround_gcc47(array)) & sizemask) == 0 \
+              && "this assertion is explained here: " \
+              "http://eigen.tuxfamily.org/dox-devel/group__TopicUnalignedArrayAssert.html" \
+              " **** READ THIS WEB PAGE !!! ****");
 #else
   #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-devel/TopicUnalignedArrayAssert.html" \
+              "http://eigen.tuxfamily.org/dox-devel/group__TopicUnalignedArrayAssert.html" \
               " **** READ THIS WEB PAGE !!! ****");
 #endif
 

Eigen/src/Core/Functors.h

@@ -533,8 +533,11 @@ template <typename Scalar, bool RandomAccess> struct linspaced_op_impl;
 // linear access for packet ops:
 // 1) initialization
 //   base = [low, ..., low] + ([step, ..., step] * [-size, ..., 0])
-// 2) each step
+// 2) each step (where size is 1 for coeff access or PacketSize for packet access)
 //   base += [size*step, ..., size*step]
+//
+// TODO: Perhaps it's better to initialize lazily (so not in the constructor but in packetOp)
+//       in order to avoid the padd() in operator() ?
 template <typename Scalar>
 struct linspaced_op_impl<Scalar,false>
 {
@@ -543,10 +546,15 @@ struct linspaced_op_impl<Scalar,false>
   linspaced_op_impl(Scalar low, Scalar step) :
   m_low(low), m_step(step),
   m_packetStep(pset1<Packet>(packet_traits<Scalar>::size*step)),
-  m_base(padd(pset1<Packet>(low),pmul(pset1<Packet>(step),plset<Scalar>(-packet_traits<Scalar>::size)))) {}
+  m_base(padd(pset1<Packet>(low), pmul(pset1<Packet>(step),plset<Scalar>(-packet_traits<Scalar>::size)))) {}
 
   template<typename Index>
-  EIGEN_STRONG_INLINE const Scalar operator() (Index i) const { return m_low+i*m_step; }
+  EIGEN_STRONG_INLINE const Scalar operator() (Index i) const 
+  { 
+    m_base = padd(m_base, pset1<Packet>(m_step));
+    return m_low+i*m_step; 
+  }
+
   template<typename Index>
   EIGEN_STRONG_INLINE const Packet packetOp(Index) const { return m_base = padd(m_base,m_packetStep); }
 

Eigen/src/Core/PermutationMatrix.h

@@ -105,13 +105,13 @@ class PermutationBase : public EigenBase<Derived>
     #endif
 
     /** \returns the number of rows */
-    inline Index rows() const { return indices().size(); }
+    inline Index rows() const { return Index(indices().size()); }
 
     /** \returns the number of columns */
-    inline Index cols() const { return indices().size(); }
+    inline Index cols() const { return Index(indices().size()); }
 
     /** \returns the size of a side of the respective square matrix, i.e., the number of indices */
-    inline Index size() const { return indices().size(); }
+    inline Index size() const { return Index(indices().size()); }
 
     #ifndef EIGEN_PARSED_BY_DOXYGEN
     template<typename DenseDerived>

Eigen/src/Core/PlainObjectBase.h

@@ -551,6 +551,7 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::type
       eigen_assert((this->size()==0 || (IsVectorAtCompileTime ? (this->size() == other.size())
                  : (rows() == other.rows() && cols() == other.cols())))
         && "Size mismatch. Automatic resizing is disabled because EIGEN_NO_AUTOMATIC_RESIZING is defined");
+      EIGEN_ONLY_USED_FOR_DEBUG(other);
       #else
       resizeLike(other);
       #endif

Eigen/src/Core/StableNorm.h

@@ -13,6 +13,7 @@
 namespace Eigen { 
 
 namespace internal {
+
 template<typename ExpressionType, typename Scalar>
 inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale)
 {
@@ -76,21 +77,20 @@ MatrixBase<Derived>::blueNorm() const
   using std::pow;
   using std::min;
   using std::max;
-  static Index nmax = -1;
+  static bool initialized = false;
   static RealScalar b1, b2, s1m, s2m, overfl, rbig, relerr;
-  if(nmax <= 0)
+  if(!initialized)
   {
-    int nbig, ibeta, it, iemin, iemax, iexp;
+    int ibeta, it, iemin, iemax, iexp;
     RealScalar abig, eps;
     // This program calculates the machine-dependent constants
-    // bl, b2, slm, s2m, relerr overfl, nmax
+    // bl, b2, slm, s2m, relerr overfl
     // from the "basic" machine-dependent numbers
-    // nbig, ibeta, it, iemin, iemax, rbig.
+    // ibeta, it, iemin, iemax, rbig.
     // The following define the basic machine-dependent constants.
     // For portability, the PORT subprograms "ilmaeh" and "rlmach"
     // are used. For any specific computer, each of the assignment
     // statements can be replaced
-    nbig  = (std::numeric_limits<Index>::max)();            // largest integer
     ibeta = std::numeric_limits<RealScalar>::radix;         // base for floating-point numbers
     it    = std::numeric_limits<RealScalar>::digits;        // number of base-beta digits in mantissa
     iemin = std::numeric_limits<RealScalar>::min_exponent;  // minimum exponent
@@ -111,8 +111,7 @@ MatrixBase<Derived>::blueNorm() const
     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
-    else                        nmax = nbig;
+    initialized = true;
   }
   Index n = size();
   RealScalar ab2 = b2 / RealScalar(n);

Eigen/src/Core/Transpose.h

@@ -353,7 +353,7 @@ struct check_transpose_aliasing_run_time_selector
 {
   static bool run(const Scalar* dest, const OtherDerived& src)
   {
-    return (bool(blas_traits<OtherDerived>::IsTransposed) != DestIsTransposed) && (dest!=0 && dest==(Scalar*)extract_data(src));
+    return (bool(blas_traits<OtherDerived>::IsTransposed) != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src));
   }
 };
 
@@ -362,8 +362,8 @@ struct check_transpose_aliasing_run_time_selector<Scalar,DestIsTransposed,CwiseB
 {
   static bool run(const Scalar* dest, const CwiseBinaryOp<BinOp,DerivedA,DerivedB>& src)
   {
-    return ((blas_traits<DerivedA>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(Scalar*)extract_data(src.lhs())))
-        || ((blas_traits<DerivedB>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(Scalar*)extract_data(src.rhs())));
+    return ((blas_traits<DerivedA>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src.lhs())))
+        || ((blas_traits<DerivedB>::IsTransposed != DestIsTransposed) && (dest!=0 && dest==(const Scalar*)extract_data(src.rhs())));
   }
 };
 

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

@@ -31,7 +31,8 @@ Packet4f plog<Packet4f>(const Packet4f& _x)
 
   /* the smallest non denormalized float number */
   _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(min_norm_pos,  0x00800000);
-
+  _EIGEN_DECLARE_CONST_Packet4f_FROM_INT(minus_inf,     0xff800000);//-1.f/0.f);
+  
   /* natural logarithm computed for 4 simultaneous float
     return NaN for x <= 0
   */
@@ -51,7 +52,8 @@ Packet4f plog<Packet4f>(const Packet4f& _x)
 
   Packet4i emm0;
 
-  Packet4f invalid_mask = _mm_cmple_ps(x, _mm_setzero_ps());
+  Packet4f invalid_mask = _mm_cmplt_ps(x, _mm_setzero_ps());
+  Packet4f iszero_mask = _mm_cmpeq_ps(x, _mm_setzero_ps());
 
   x = pmax(x, p4f_min_norm_pos);  /* cut off denormalized stuff */
   emm0 = _mm_srli_epi32(_mm_castps_si128(x), 23);
@@ -96,7 +98,9 @@ Packet4f plog<Packet4f>(const Packet4f& _x)
   y2 = pmul(e, p4f_cephes_log_q2);
   x = padd(x, y);
   x = padd(x, y2);
-  return _mm_or_ps(x, invalid_mask); // negative arg will be NAN
+  // negative arg will be NAN, 0 will be -INF
+  return _mm_or_ps(_mm_andnot_ps(iszero_mask, _mm_or_ps(x, invalid_mask)),
+                   _mm_and_ps(iszero_mask, p4f_minus_inf));
 }
 
 template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED

Eigen/src/Core/products/GeneralBlockPanelKernel.h

@@ -69,8 +69,8 @@ inline void manage_caching_sizes(Action action, std::ptrdiff_t* l1=0, std::ptrdi
   * - the number of scalars that fit into a packet (when vectorization is enabled).
   *
   * \sa setCpuCacheSizes */
-template<typename LhsScalar, typename RhsScalar, int KcFactor>
-void computeProductBlockingSizes(std::ptrdiff_t& k, std::ptrdiff_t& m, std::ptrdiff_t& n)
+template<typename LhsScalar, typename RhsScalar, int KcFactor, typename SizeType>
+void computeProductBlockingSizes(SizeType& k, SizeType& m, SizeType& n)
 {
   EIGEN_UNUSED_VARIABLE(n);
   // Explanations:
@@ -91,13 +91,13 @@ void computeProductBlockingSizes(std::ptrdiff_t& k, std::ptrdiff_t& m, std::ptrd
   };
 
   manage_caching_sizes(GetAction, &l1, &l2);
-  k = std::min<std::ptrdiff_t>(k, l1/kdiv);
-  std::ptrdiff_t _m = k>0 ? l2/(4 * sizeof(LhsScalar) * k) : 0;
+  k = std::min<SizeType>(k, l1/kdiv);
+  SizeType _m = k>0 ? l2/(4 * sizeof(LhsScalar) * k) : 0;
   if(_m<m) m = _m & mr_mask;
 }
 
-template<typename LhsScalar, typename RhsScalar>
-inline void computeProductBlockingSizes(std::ptrdiff_t& k, std::ptrdiff_t& m, std::ptrdiff_t& n)
+template<typename LhsScalar, typename RhsScalar, typename SizeType>
+inline void computeProductBlockingSizes(SizeType& k, SizeType& m, SizeType& n)
 {
   computeProductBlockingSizes<LhsScalar,RhsScalar,1>(k, m, n);
 }

Eigen/src/Core/util/Macros.h

@@ -13,7 +13,7 @@
 
 #define EIGEN_WORLD_VERSION 3
 #define EIGEN_MAJOR_VERSION 1
-#define EIGEN_MINOR_VERSION 2
+#define EIGEN_MINOR_VERSION 3
 
 #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/src/Core/util/Memory.h

@@ -88,11 +88,11 @@ inline void throw_std_bad_alloc()
 /** \internal Like malloc, but the returned pointer is guaranteed to be 16-byte aligned.
   * Fast, but wastes 16 additional bytes of memory. Does not throw any exception.
   */
-inline void* handmade_aligned_malloc(size_t size)
+inline void* handmade_aligned_malloc(std::size_t size)
 {
   void *original = std::malloc(size+16);
   if (original == 0) return 0;
-  void *aligned = reinterpret_cast<void*>((reinterpret_cast<size_t>(original) & ~(size_t(15))) + 16);
+  void *aligned = reinterpret_cast<void*>((reinterpret_cast<std::size_t>(original) & ~(std::size_t(15))) + 16);
   *(reinterpret_cast<void**>(aligned) - 1) = original;
   return aligned;
 }
@@ -108,13 +108,18 @@ inline void handmade_aligned_free(void *ptr)
   * Since we know that our handmade version is based on std::realloc
   * we can use std::realloc to implement efficient reallocation.
   */
-inline void* handmade_aligned_realloc(void* ptr, size_t size, size_t = 0)
+inline void* handmade_aligned_realloc(void* ptr, std::size_t size, std::size_t = 0)
 {
   if (ptr == 0) return handmade_aligned_malloc(size);
   void *original = *(reinterpret_cast<void**>(ptr) - 1);
+  std::ptrdiff_t previous_offset = static_cast<char *>(ptr)-static_cast<char *>(original);
   original = std::realloc(original,size+16);
   if (original == 0) return 0;
-  void *aligned = reinterpret_cast<void*>((reinterpret_cast<size_t>(original) & ~(size_t(15))) + 16);
+  void *aligned = reinterpret_cast<void*>((reinterpret_cast<std::size_t>(original) & ~(std::size_t(15))) + 16);
+  void *previous_aligned = static_cast<char *>(original)+previous_offset;
+  if(aligned!=previous_aligned)
+    std::memmove(aligned, previous_aligned, size);
+  
   *(reinterpret_cast<void**>(aligned) - 1) = original;
   return aligned;
 }
@@ -123,7 +128,7 @@ inline void* handmade_aligned_realloc(void* ptr, size_t size, size_t = 0)
 *** Implementation of generic aligned realloc (when no realloc can be used)***
 *****************************************************************************/
 
-void* aligned_malloc(size_t size);
+void* aligned_malloc(std::size_t size);
 void  aligned_free(void *ptr);
 
 /** \internal
@@ -227,7 +232,7 @@ inline void aligned_free(void *ptr)
     std::free(ptr);
   #elif EIGEN_HAS_MM_MALLOC
     _mm_free(ptr);
-  #elif defined(_MSC_VER)
+  #elif defined(_MSC_VER) && (!defined(_WIN32_WCE))
     _aligned_free(ptr);
   #else
     handmade_aligned_free(ptr);

Eigen/src/Geometry/Quaternion.h

@@ -193,7 +193,8 @@ public:
   *
   * \brief The quaternion class used to represent 3D orientations and rotations
   *
-  * \param _Scalar the scalar type, i.e., the type of the coefficients
+  * \tparam _Scalar the scalar type, i.e., the type of the coefficients
+  * \tparam _Options controls the memory alignement of the coeffecients. Can be \# AutoAlign or \# DontAlign. Default is AutoAlign.
   *
   * This class represents a quaternion \f$ w+xi+yj+zk \f$ that is a convenient representation of
   * orientations and rotations of objects in three dimensions. Compared to other representations
@@ -304,41 +305,29 @@ typedef Quaternion<double> Quaterniond;
 
 namespace internal {
   template<typename _Scalar, int _Options>
-  struct traits<Map<Quaternion<_Scalar>, _Options> >:
-  traits<Quaternion<_Scalar, _Options> >
+  struct traits<Map<Quaternion<_Scalar>, _Options> > : traits<Quaternion<_Scalar, (int(_Options)&Aligned)==Aligned ? AutoAlign : DontAlign> >
   {
-    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> >
+  struct traits<Map<const Quaternion<_Scalar>, _Options> > : traits<Quaternion<_Scalar, (int(_Options)&Aligned)==Aligned ? AutoAlign : DontAlign> >
   {
-    typedef _Scalar Scalar;
     typedef Map<const Matrix<_Scalar,4,1>, _Options> Coefficients;
-
-    typedef traits<Quaternion<_Scalar, _Options> > TraitsBase;
+    typedef traits<Quaternion<_Scalar, (int(_Options)&Aligned)==Aligned ? AutoAlign : DontAlign> > TraitsBase;
     enum {
-      IsAligned = TraitsBase::IsAligned,
       Flags = TraitsBase::Flags & ~LvalueBit
     };
   };
 }
 
-/** \brief Quaternion expression mapping a constant memory buffer
+/** \ingroup Geometry_Module
+  * \brief Quaternion expression mapping a constant memory buffer
   *
-  * \param _Scalar the type of the Quaternion coefficients
-  * \param _Options see class Map
+  * \tparam _Scalar the type of the Quaternion coefficients
+  * \tparam _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.
@@ -371,10 +360,11 @@ class Map<const Quaternion<_Scalar>, _Options >
     const Coefficients m_coeffs;
 };
 
-/** \brief Expression of a quaternion from a memory buffer
+/** \ingroup Geometry_Module
+  * \brief Expression of a quaternion from a memory buffer
   *
-  * \param _Scalar the type of the Quaternion coefficients
-  * \param _Options see class Map
+  * \tparam _Scalar the type of the Quaternion coefficients
+  * \tparam _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.

Eigen/src/LU/FullPivLU.h

@@ -577,7 +577,7 @@ struct kernel_retval<FullPivLU<_MatrixType> >
     RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
     Index p = 0;
     for(Index i = 0; i < dec().nonzeroPivots(); ++i)
-      if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
+      if(internal::abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
         pivots.coeffRef(p++) = i;
     eigen_internal_assert(p == rank());
 
@@ -645,7 +645,7 @@ struct image_retval<FullPivLU<_MatrixType> >
     RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
     Index p = 0;
     for(Index i = 0; i < dec().nonzeroPivots(); ++i)
-      if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
+      if(internal::abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
         pivots.coeffRef(p++) = i;
     eigen_internal_assert(p == rank());
 

Eigen/src/QR/ColPivHouseholderQR.h

@@ -56,6 +56,12 @@ template<typename _MatrixType> class ColPivHouseholderQR
     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
     typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
     typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
+    
+  private:
+    
+    typedef typename PermutationType::Index PermIndexType;
+    
+  public:
 
     /**
     * \brief Default Constructor.
@@ -81,7 +87,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
     ColPivHouseholderQR(Index rows, Index cols)
       : m_qr(rows, cols),
         m_hCoeffs((std::min)(rows,cols)),
-        m_colsPermutation(cols),
+        m_colsPermutation(PermIndexType(cols)),
         m_colsTranspositions(cols),
         m_temp(cols),
         m_colSqNorms(cols),
@@ -91,7 +97,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
     ColPivHouseholderQR(const MatrixType& matrix)
       : m_qr(matrix.rows(), matrix.cols()),
         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
-        m_colsPermutation(matrix.cols()),
+        m_colsPermutation(PermIndexType(matrix.cols())),
         m_colsTranspositions(matrix.cols()),
         m_temp(matrix.cols()),
         m_colSqNorms(matrix.cols()),
@@ -436,9 +442,9 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
     m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
   }
 
-  m_colsPermutation.setIdentity(cols);
-  for(Index k = 0; k < m_nonzero_pivots; ++k)
-    m_colsPermutation.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
+  m_colsPermutation.setIdentity(PermIndexType(cols));
+  for(PermIndexType k = 0; k < m_nonzero_pivots; ++k)
+    m_colsPermutation.applyTranspositionOnTheRight(PermIndexType(k), PermIndexType(m_colsTranspositions.coeff(k)));
 
   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
   m_isInitialized = true;

Eigen/src/SparseCore/SparseSelfAdjointView.h

@@ -209,6 +209,7 @@ class SparseSelfAdjointTimeDenseProduct
 
     template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const
     {
+      EIGEN_ONLY_USED_FOR_DEBUG(alpha);
       // TODO use alpha
       eigen_assert(alpha==Scalar(1) && "alpha != 1 is not implemented yet, sorry");
       typedef typename internal::remove_all<Lhs>::type _Lhs;

Eigen/src/SparseCore/SparseVector.h

@@ -202,7 +202,7 @@ class SparseVector
     }
 
     inline SparseVector(const SparseVector& other)
-      : m_size(0)
+      : SparseBase(other), m_size(0)
     {
       *this = other.derived();
     }
@@ -230,7 +230,8 @@ class SparseVector
     template<typename OtherDerived>
     inline SparseVector& operator=(const SparseMatrixBase<OtherDerived>& other)
     {
-      if (int(RowsAtCompileTime)!=int(OtherDerived::RowsAtCompileTime))
+      if ( (bool(OtherDerived::IsVectorAtCompileTime) && int(RowsAtCompileTime)!=int(OtherDerived::RowsAtCompileTime))
+          || ((!bool(OtherDerived::IsVectorAtCompileTime)) && ( bool(IsColVector) ? other.cols()>1 : other.rows()>1 )))
         return assign(other.transpose());
       else
         return assign(other);

cmake/language_support.cmake

@@ -24,6 +24,8 @@ function(workaround_9220 language language_works)
   set(text
     "project(test NONE)
     cmake_minimum_required(VERSION 2.6.0)
+    set (CMAKE_Fortran_FLAGS \"${CMAKE_Fortran_FLAGS}\")
+    set (CMAKE_EXE_LINKER_FLAGS \"${CMAKE_EXE_LINKER_FLAGS}\")
     enable_language(${language} OPTIONAL)
   ")
   file(REMOVE_RECURSE ${CMAKE_BINARY_DIR}/language_tests/${language})

doc/C08_TutorialGeometry.dox

@@ -178,7 +178,7 @@ matNxN = t.linear();
 \endcode</td></tr>
 <tr><td>
 extract the rotation matrix</td><td>\code
-matNxN = t.extractRotation();
+matNxN = t.rotation();
 \endcode</td></tr>
 </table>
 

doc/TopicLinearAlgebraDecompositions.dox

@@ -6,7 +6,6 @@ namespace Eigen {
 \section TopicLinAlgBigTable Catalogue of decompositions offered by Eigen
 
 <table class="manual-vl">
-
     <tr>
         <th class="meta"></th>
         <th class="meta" colspan="5">Generic information, not Eigen-specific</th>

test/array_for_matrix.cpp

@@ -42,10 +42,10 @@ template<typename MatrixType> void array_for_matrix(const MatrixType& m)
   VERIFY_IS_APPROX(m3, (m1.array() - s1).matrix());
 
   // reductions
-  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_MUCH_SMALLER_THAN(m1.colwise().sum().sum() - m1.sum(), m1.squaredNorm());
+  VERIFY_IS_MUCH_SMALLER_THAN(m1.rowwise().sum().sum() - m1.sum(), m1.squaredNorm());
+  VERIFY_IS_MUCH_SMALLER_THAN(m1.colwise().sum() + m2.colwise().sum() - (m1+m2).colwise().sum(), (m1+m2).squaredNorm());
+  VERIFY_IS_MUCH_SMALLER_THAN(m1.rowwise().sum() - m2.rowwise().sum() - (m1-m2).rowwise().sum(), (m1-m2).squaredNorm());
   VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op<Scalar>()));
 
   // vector-wise ops

test/diagonal.cpp

@@ -56,12 +56,12 @@ template<typename MatrixType> void diagonal(const MatrixType& m)
     VERIFY_IS_APPROX(m2.template diagonal<N2>()[0], static_cast<Scalar>(6) * m1.template diagonal<N2>()[0]);
 
     m2.diagonal(N1) = 2 * m1.diagonal(N1);
-    VERIFY_IS_APPROX(m2.diagonal<N1>(), static_cast<Scalar>(2) * m1.diagonal(N1));
+    VERIFY_IS_APPROX(m2.template diagonal<N1>(), static_cast<Scalar>(2) * m1.diagonal(N1));
     m2.diagonal(N1)[0] *= 3;
     VERIFY_IS_APPROX(m2.diagonal(N1)[0], static_cast<Scalar>(6) * m1.diagonal(N1)[0]);
 
     m2.diagonal(N2) = 2 * m1.diagonal(N2);
-    VERIFY_IS_APPROX(m2.diagonal<N2>(), static_cast<Scalar>(2) * m1.diagonal(N2));
+    VERIFY_IS_APPROX(m2.template diagonal<N2>(), static_cast<Scalar>(2) * m1.diagonal(N2));
     m2.diagonal(N2)[0] *= 3;
     VERIFY_IS_APPROX(m2.diagonal(N2)[0], static_cast<Scalar>(6) * m1.diagonal(N2)[0]);
   }

test/geo_quaternion.cpp

@@ -171,23 +171,36 @@ template<typename Scalar, int Options> void quaternion(void)
 
 template<typename Scalar> void mapQuaternion(void){
   typedef Map<Quaternion<Scalar>, Aligned> MQuaternionA;
+  typedef Map<const Quaternion<Scalar>, Aligned> MCQuaternionA;
   typedef Map<Quaternion<Scalar> > MQuaternionUA;
   typedef Map<const Quaternion<Scalar> > MCQuaternionUA;
   typedef Quaternion<Scalar> Quaternionx;
+  typedef Matrix<Scalar,3,1> Vector3;
+  typedef AngleAxis<Scalar> AngleAxisx;
+  
+  Vector3 v0 = Vector3::Random(),
+          v1 = Vector3::Random();
+  Scalar  a = internal::random<Scalar>(-Scalar(M_PI), Scalar(M_PI));
 
   EIGEN_ALIGN16 Scalar array1[4];
   EIGEN_ALIGN16 Scalar array2[4];
   EIGEN_ALIGN16 Scalar array3[4+1];
   Scalar* array3unaligned = array3+1;
+  
+  MQuaternionA    mq1(array1);
+  MCQuaternionA   mcq1(array1);
+  MQuaternionA    mq2(array2);
+  MQuaternionUA   mq3(array3unaligned);
+  MCQuaternionUA  mcq3(array3unaligned);
 
 //  std::cerr << array1 << " " << array2 << " " << array3 << "\n";
-  MQuaternionA(array1).coeffs().setRandom();
-  (MQuaternionA(array2)) = MQuaternionA(array1);
-  (MQuaternionUA(array3unaligned)) = MQuaternionA(array1);
+  mq1 = AngleAxisx(a, v0.normalized());
+  mq2 = mq1;
+  mq3 = mq1;
 
-  Quaternionx q1 = MQuaternionA(array1);
-  Quaternionx q2 = MQuaternionA(array2);
-  Quaternionx q3 = MQuaternionUA(array3unaligned);
+  Quaternionx q1 = mq1;
+  Quaternionx q2 = mq2;
+  Quaternionx q3 = mq3;
   Quaternionx q4 = MCQuaternionUA(array3unaligned);
 
   VERIFY_IS_APPROX(q1.coeffs(), q2.coeffs());
@@ -197,6 +210,23 @@ template<typename Scalar> void mapQuaternion(void){
   if(internal::packet_traits<Scalar>::Vectorizable)
     VERIFY_RAISES_ASSERT((MQuaternionA(array3unaligned)));
   #endif
+    
+  VERIFY_IS_APPROX(mq1 * (mq1.inverse() * v1), v1);
+  VERIFY_IS_APPROX(mq1 * (mq1.conjugate() * v1), v1);
+  
+  VERIFY_IS_APPROX(mcq1 * (mcq1.inverse() * v1), v1);
+  VERIFY_IS_APPROX(mcq1 * (mcq1.conjugate() * v1), v1);
+  
+  VERIFY_IS_APPROX(mq3 * (mq3.inverse() * v1), v1);
+  VERIFY_IS_APPROX(mq3 * (mq3.conjugate() * v1), v1);
+  
+  VERIFY_IS_APPROX(mcq3 * (mcq3.inverse() * v1), v1);
+  VERIFY_IS_APPROX(mcq3 * (mcq3.conjugate() * v1), v1);
+  
+  VERIFY_IS_APPROX(mq1*mq2, q1*q2);
+  VERIFY_IS_APPROX(mq3*mq2, q3*q2);
+  VERIFY_IS_APPROX(mcq1*mq2, q1*q2);
+  VERIFY_IS_APPROX(mcq3*mq2, q3*q2);
 }
 
 template<typename Scalar> void quaternionAlignment(void){

test/nullary.cpp

@@ -91,6 +91,12 @@ void testVectorType(const VectorType& base)
   scalar.setLinSpaced(1,low,high);
   VERIFY_IS_APPROX( scalar, ScalarMatrix::Constant(high) );
   VERIFY_IS_APPROX( ScalarMatrix::LinSpaced(1,low,high), ScalarMatrix::Constant(high) );
+
+  // regression test for bug 526 (linear vectorized transversal)
+  if (size > 1) {
+    m.tail(size-1).setLinSpaced(low, high);
+    VERIFY_IS_APPROX(m(size-1), high);
+  }
 }
 
 template<typename MatrixType>

test/packetmath.cpp

@@ -40,7 +40,7 @@ template<typename Scalar> bool areApprox(const Scalar* a, const Scalar* b, int s
 {
   for (int i=0; i<size; ++i)
   {
-    if (!internal::isApprox(a[i],b[i]))
+    if (a[i]!=b[i] && !internal::isApprox(a[i],b[i]))
     {
       std::cout << "[" << Map<const Matrix<Scalar,1,Dynamic> >(a,size) << "]" << " != " << Map<const Matrix<Scalar,1,Dynamic> >(b,size) << "\n";
       return false;

test/sparse.h

@@ -178,5 +178,30 @@ initSparse(double density,
   }
 }
 
+template<typename Scalar> void
+initSparse(double density,
+           Matrix<Scalar,1,Dynamic>& refVec,
+           SparseVector<Scalar,RowMajor>& sparseVec,
+           std::vector<int>* zeroCoords = 0,
+           std::vector<int>* nonzeroCoords = 0)
+{
+  sparseVec.reserve(int(refVec.size()*density));
+  sparseVec.setZero();
+  for(int i=0; i<refVec.size(); i++)
+  {
+    Scalar v = (internal::random<double>(0,1) < density) ? internal::random<Scalar>() : Scalar(0);
+    if (v!=Scalar(0))
+    {
+      sparseVec.insertBack(i) = v;
+      if (nonzeroCoords)
+        nonzeroCoords->push_back(i);
+    }
+    else if (zeroCoords)
+        zeroCoords->push_back(i);
+    refVec[i] = v;
+  }
+}
+
+
 #include <unsupported/Eigen/SparseExtra>
 #endif // EIGEN_TESTSPARSE_H

test/sparse_product.cpp

@@ -46,8 +46,10 @@ template<typename SparseMatrixType> void sparse_product()
   double density = (std::max)(8./(rows*cols), 0.01);
   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
   typedef Matrix<Scalar,Dynamic,1> DenseVector;
-
-  Scalar s1 = internal::random<Scalar>();
+  typedef Matrix<Scalar,1,Dynamic> RowDenseVector;
+  typedef SparseVector<Scalar,0,Index> ColSpVector;
+  typedef SparseVector<Scalar,RowMajor,Index> RowSpVector;Scalar s1 = internal::random<Scalar>();
+  
   Scalar s2 = internal::random<Scalar>();
 
   // test matrix-matrix product
@@ -117,6 +119,21 @@ template<typename SparseMatrixType> void sparse_product()
     test_outer<SparseMatrixType,DenseMatrix>::run(m2,m4,refMat2,refMat4);
 
     VERIFY_IS_APPROX(m6=m6*m6, refMat6=refMat6*refMat6);
+    
+    // sparse matrix * sparse vector
+    ColSpVector cv0(cols), cv1;
+    DenseVector dcv0(cols), dcv1;
+    initSparse(2*density,dcv0, cv0);
+    
+    RowSpVector rv0(depth), rv1;
+    RowDenseVector drv0(depth), drv1(rv1);
+    initSparse(2*density,drv0, rv0);
+    
+    VERIFY_IS_APPROX(cv1=rv0*m3, dcv1=drv0*refMat3);
+    VERIFY_IS_APPROX(rv1=rv0*m3, drv1=drv0*refMat3);
+    VERIFY_IS_APPROX(cv1=m3*cv0, dcv1=refMat3*dcv0);
+    VERIFY_IS_APPROX(cv1=m3t.adjoint()*cv0, dcv1=refMat3t.adjoint()*dcv0);
+    VERIFY_IS_APPROX(rv1=m3*cv0, drv1=refMat3*dcv0);
   }
 
   // test matrix - diagonal product

test/sparse_vector.cpp

@@ -82,6 +82,12 @@ template<typename Scalar> void sparse_vector(int rows, int cols)
   VERIFY_IS_APPROX((v1 = -v1), (refV1 = -refV1));
   VERIFY_IS_APPROX((v1 = v1.transpose()), (refV1 = refV1.transpose().eval()));
   VERIFY_IS_APPROX((v1 += -v1), (refV1 += -refV1));
+  
+  // sparse matrix to sparse vector
+  SparseMatrixType mv1;
+  VERIFY_IS_APPROX((mv1=v1),v1);
+  VERIFY_IS_APPROX(mv1,(v1=mv1));
+  VERIFY_IS_APPROX(mv1,(v1=mv1.transpose()));
 
 }