bump

expand a little the vectorization_logic test and backport EIGEN_DEBUG_ASSIGN.

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.13")
+set(EIGEN_VERSION_NUMBER "2.0.15")
 set(EIGEN_VERSION "${EIGEN_VERSION_NUMBER}")
 
 set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)

Eigen/src/Array/BooleanRedux.h

@@ -139,7 +139,7 @@ inline bool MatrixBase<Derived>::any() const
 template<typename Derived>
 inline int MatrixBase<Derived>::count() const
 {
-  return this->cast<bool>().cast<int>().sum();
+  return this->cast<bool>().template cast<int>().sum();
 }
 
 #endif // EIGEN_ALLANDANY_H

Eigen/src/Cholesky/LLT.h

@@ -132,11 +132,12 @@ void LLT<MatrixType>::compute(const MatrixType& a)
     m_isInitialized = true;
     return;
   }
-  m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
+  if(ei_real(m_matrix.coeff(0,0))>0)
+    m_matrix.col(0).end(size-1) = a.row(0).end(size-1).adjoint() / ei_real(m_matrix.coeff(0,0));
   for (int j = 1; j < size; ++j)
   {
     x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm();
-    if (x < cutoff)
+    if (x <= cutoff)
     {
       m_isPositiveDefinite = false;
       continue;

Eigen/src/Core/Assign.h

@@ -90,6 +90,28 @@ public:
               ? ( int(MayUnrollCompletely) && int(DstIsAligned) ? int(CompleteUnrolling) : int(NoUnrolling) )
               : int(NoUnrolling)
   };
+
+#ifdef EIGEN_DEBUG_ASSIGN
+#define EIGEN_DEBUG_VAR(x) std::cerr << #x << " = " << x << std::endl;
+  static void debug()
+  {
+    EIGEN_DEBUG_VAR(DstIsAligned)
+    EIGEN_DEBUG_VAR(SrcIsAligned)
+    EIGEN_DEBUG_VAR(SrcAlignment)
+    EIGEN_DEBUG_VAR(InnerSize)
+    EIGEN_DEBUG_VAR(InnerMaxSize)
+    EIGEN_DEBUG_VAR(PacketSize)
+    EIGEN_DEBUG_VAR(MightVectorize)
+    EIGEN_DEBUG_VAR(MayInnerVectorize)
+    EIGEN_DEBUG_VAR(MayLinearVectorize)
+    EIGEN_DEBUG_VAR(MaySliceVectorize)
+    EIGEN_DEBUG_VAR(Vectorization)
+    EIGEN_DEBUG_VAR(UnrollingLimit)
+    EIGEN_DEBUG_VAR(MayUnrollCompletely)
+    EIGEN_DEBUG_VAR(MayUnrollInner)
+    EIGEN_DEBUG_VAR(Unrolling)
+  }
+#endif
 };
 
 /***************************************************************************
@@ -400,6 +422,9 @@ template<typename OtherDerived>
 EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>
   ::lazyAssign(const MatrixBase<OtherDerived>& other)
 {
+#ifdef EIGEN_DEBUG_ASSIGN
+  ei_assign_traits<Derived, OtherDerived>::debug();
+#endif
   EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Derived,OtherDerived)
   EIGEN_STATIC_ASSERT((ei_is_same_type<typename Derived::Scalar, typename OtherDerived::Scalar>::ret),
     YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)

Eigen/src/Core/CacheFriendlyProduct.h

@@ -33,7 +33,7 @@ struct ei_L2_block_traits {
 #ifndef EIGEN_EXTERN_INSTANTIATIONS
 
 template<typename Scalar>
-static void ei_cache_friendly_product(
+void ei_cache_friendly_product(
   int _rows, int _cols, int depth,
   bool _lhsRowMajor, const Scalar* _lhs, int _lhsStride,
   bool _rhsRowMajor, const Scalar* _rhs, int _rhsStride,
@@ -352,7 +352,7 @@ static void ei_cache_friendly_product(
  * TODO: since rhs gets evaluated only once, no need to evaluate it
  */
 template<typename Scalar, typename RhsType>
-static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
+EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
   int size,
   const Scalar* lhs, int lhsStride,
   const RhsType& rhs,
@@ -542,7 +542,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
 
 // TODO add peeling to mask unaligned load/stores
 template<typename Scalar, typename ResType>
-static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
+EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
   const Scalar* lhs, int lhsStride,
   const Scalar* rhs, int rhsSize,
   ResType& res)

Eigen/src/Core/CoreInstantiations.cpp

@@ -32,7 +32,7 @@ namespace Eigen
 {
 
 #define EIGEN_INSTANTIATE_PRODUCT(TYPE) \
-template static void ei_cache_friendly_product<TYPE>( \
+template void ei_cache_friendly_product<TYPE>( \
   int _rows, int _cols, int depth, \
   bool _lhsRowMajor, const TYPE* _lhs, int _lhsStride, \
   bool _rhsRowMajor, const TYPE* _rhs, int _rhsStride, \

Eigen/src/Core/MathFunctions.h

@@ -35,16 +35,6 @@ template<typename T> inline T ei_random_amplitude()
   else return static_cast<T>(10);
 }
 
-template<typename T> inline T ei_hypot(T x, T y)
-{
-  T _x = ei_abs(x);
-  T _y = ei_abs(y);
-  T p = std::max(_x, _y);
-  T q = std::min(_x, _y);
-  T qp = q/p;
-  return p * ei_sqrt(T(1) + qp*qp);
-}
-
 /**************
 ***   int   ***
 **************/
@@ -292,4 +282,14 @@ inline bool ei_isApproxOrLessThan(long double a, long double b, long double prec
   return a <= b || ei_isApprox(a, b, prec);
 }
 
+template<typename T> inline T ei_hypot(T x, T y)
+{
+  T _x = ei_abs(x);
+  T _y = ei_abs(y);
+  T p = std::max(_x, _y);
+  T q = std::min(_x, _y);
+  T qp = q/p;
+  return p * ei_sqrt(T(1) + qp*qp);
+}
+
 #endif // EIGEN_MATHFUNCTIONS_H

Eigen/src/Core/Product.h

@@ -614,7 +614,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
     else
     {
       _res = ei_aligned_stack_new(Scalar, res.size());
-      Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1,ColMajor> >(_res, res.size()) = res;
+      Map<Matrix<Scalar,1,DestDerived::SizeAtCompileTime,ColMajor> >(_res, res.size()) = res;
     }
     ei_cache_friendly_product_colmajor_times_vector(res.size(),
       &product.rhs().const_cast_derived().coeffRef(0,0), product.rhs().stride(),
@@ -622,7 +622,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
 
     if (!EvalToRes)
     {
-      res = Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1,ColMajor> >(_res, res.size());
+      res = Map<Matrix<Scalar,1,DestDerived::SizeAtCompileTime,ColMajor> >(_res, res.size());
       ei_aligned_stack_delete(Scalar, _res, res.size());
     }
   }
@@ -677,7 +677,7 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
     else
     {
       _lhs = ei_aligned_stack_new(Scalar, product.lhs().size());
-      Map<Matrix<Scalar,Lhs::SizeAtCompileTime,1,ColMajor> >(_lhs, product.lhs().size()) = product.lhs();
+      Map<Matrix<Scalar,1,Lhs::SizeAtCompileTime,ColMajor> >(_lhs, product.lhs().size()) = product.lhs();
     }
     ei_cache_friendly_product_rowmajor_times_vector(&product.rhs().const_cast_derived().coeffRef(0,0), product.rhs().stride(),
                                                     _lhs, product.lhs().size(), res);

Eigen/src/Core/SolveTriangular.h

@@ -35,7 +35,7 @@ template<typename Lhs, typename Rhs,
                      ? UpperTriangular
                      : -1,
   int StorageOrder = ei_is_part<Lhs>::value ? -1  // this is to solve ambiguous specializations
-                   : int(Lhs::Flags) & (RowMajorBit|SparseBit)
+                   : int(Lhs::Flags) & int(RowMajorBit|SparseBit)
   >
 struct ei_solve_triangular_selector;
 

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

@@ -114,13 +114,13 @@ 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); }
 
-#if defined(_MSC_VER) && (_MSC_VER <= 1500) && defined(_WIN64)
+#if defined(_MSC_VER) && (_MSC_VER <= 1500) && defined(_WIN64) && !defined(__INTEL_COMPILER)
 // 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)
+#elif defined(_MSC_VER) && (_MSC_VER <= 1500) && !defined(__INTEL_COMPILER)
 // 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; }

Eigen/src/Core/util/Macros.h

@@ -30,7 +30,7 @@
 
 #define EIGEN_WORLD_VERSION 2
 #define EIGEN_MAJOR_VERSION 0
-#define EIGEN_MINOR_VERSION 13
+#define EIGEN_MINOR_VERSION 15
 
 #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

@@ -188,7 +188,8 @@ template<> inline void ei_conditional_aligned_free<false>(void *ptr)
 template<typename T> inline void ei_destruct_elements_of_array(T *ptr, std::size_t size)
 {
   // always destruct an array starting from the end.
-  while(size) ptr[--size].~T();
+  if(ptr)
+    while(size) ptr[--size].~T();
 }
 
 /** \internal delete objects constructed with ei_aligned_new

Eigen/src/Core/util/XprHelper.h

@@ -94,10 +94,16 @@ class ei_compute_matrix_flags
 {
     enum {
       row_major_bit = Options&RowMajor ? RowMajorBit : 0,
-      inner_max_size = row_major_bit ? MaxCols : MaxRows,
+      inner_max_size = int(MaxRows==1) ? int(MaxCols)
+                     : int(MaxCols==1) ? int(MaxRows)
+                     : int(row_major_bit) ? int(MaxCols) : int(MaxRows),
       is_big = inner_max_size == Dynamic,
-      is_packet_size_multiple = (Cols*Rows) % ei_packet_traits<Scalar>::size == 0,
-      aligned_bit = ((Options&AutoAlign) && (is_big || is_packet_size_multiple)) ? AlignedBit : 0,
+      storage_has_fixed_size = MaxRows != Dynamic && MaxCols != Dynamic,
+      storage_has_aligned_fixed_size = storage_has_fixed_size
+                     && ( (MaxCols*MaxRows) % ei_packet_traits<Scalar>::size == 0 ),
+      aligned_bit = (    (Options&AutoAlign)
+                      && (is_big || storage_has_aligned_fixed_size)
+                    ) ? AlignedBit : 0,
       packet_access_bit = ei_packet_traits<Scalar>::size > 1 && aligned_bit ? PacketAccessBit : 0
     };
 

Eigen/src/LU/LU.h

@@ -63,10 +63,10 @@ template<typename MatrixType> class LU
 
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
-    typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
-    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
-    typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType;
-    typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType;
+    typedef Matrix<int, 1, MatrixType::ColsAtCompileTime, MatrixType::Options, 1, MatrixType::MaxColsAtCompileTime> IntRowVectorType;
+    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, 1> IntColVectorType;
+    typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime, MatrixType::Options, 1, MatrixType::MaxColsAtCompileTime> RowVectorType;
+    typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, 1> ColVectorType;
 
     enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN(
              MatrixType::MaxColsAtCompileTime,
@@ -515,9 +515,10 @@ bool LU<MatrixType>::solve(
         if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c, m_precision))
           return false;
   }
-  m_lu.corner(TopLeft, m_rank, m_rank)
-      .template marked<UpperTriangular>()
-      .solveTriangularInPlace(c.corner(TopLeft, m_rank, c.cols()));
+  if(m_rank>0)
+    m_lu.corner(TopLeft, m_rank, m_rank)
+        .template marked<UpperTriangular>()
+        .solveTriangularInPlace(c.corner(TopLeft, m_rank, c.cols()));
 
   // Step 4
   result->resize(m_lu.cols(), b.cols());

Eigen/src/QR/QR.h

@@ -270,12 +270,18 @@ bool QR<MatrixType>::solve(
   ei_assert(m_isInitialized && "QR is not initialized.");
   const int rows = m_qr.rows();
   ei_assert(b.rows() == rows);
-  result->resize(rows, b.cols());
+  // enforce the computation of the rank
+  rank();
+  
+  result->resize(m_qr.cols(), b.cols());
   
   // TODO(keir): There is almost certainly a faster way to multiply by
   // Q^T without explicitly forming matrixQ(). Investigate.
   *result = matrixQ().transpose()*b;
-  
+
+  if(m_rank==0)
+    return result->isZero();
+
   if(!isSurjective())
   {
     // is result is in the image of R ?

test/cholesky.cpp

@@ -124,4 +124,8 @@ void test_cholesky()
     CALL_SUBTEST( cholesky(MatrixXf(17,17)) );
     CALL_SUBTEST( cholesky(MatrixXd(33,33)) );
   }
+
+  MatrixXf m = MatrixXf::Zero(10,10);
+  VectorXf b = VectorXf::Zero(10);
+  VERIFY(!m.llt().isPositiveDefinite());
 }

test/lu.cpp

@@ -132,4 +132,10 @@ void test_lu()
     CALL_SUBTEST( lu_invertible<MatrixXcf>() );
     CALL_SUBTEST( lu_invertible<MatrixXcd>() );
   }
+
+  MatrixXf m = MatrixXf::Zero(10,10);
+  VectorXf b = VectorXf::Zero(10);
+  VectorXf x = VectorXf::Random(10);
+  VERIFY(m.lu().solve(b,&x));
+  VERIFY(x.isZero());
 }

test/product_large.cpp

@@ -48,4 +48,11 @@ void test_product_large()
     MatrixXf a = MatrixXf::Random(10,4), b = MatrixXf::Random(4,10), c = a;
     VERIFY_IS_APPROX((a = a * b), (c * b).eval());
   }
+
+  {
+    MatrixXf mat1(10,10); mat1.setRandom();
+    MatrixXf mat2(32,10); mat2.setRandom();
+    MatrixXf result = mat1.row(2)*mat2.transpose();
+    VERIFY_IS_APPROX(result, (mat1.row(2)*mat2.transpose()).eval());
+  }
 }

test/qr.cpp

@@ -75,4 +75,11 @@ void test_qr()
     mat << 1, 1, 1, 2, 2, 2, 1, 2, 3;
     VERIFY(!mat.qr().isFullRank());
   }
+  {
+    MatrixXf m = MatrixXf::Zero(10,10);
+    VectorXf b = VectorXf::Zero(10);
+    VectorXf x = VectorXf::Random(10);
+    VERIFY(m.qr().solve(b,&x));
+    VERIFY(x.isZero());
+  }
 }

test/vectorization_logic.cpp

@@ -85,6 +85,8 @@ void test_vectorization_logic()
   VERIFY(test_assign(MatrixXf(10,10),MatrixXf(20,20).block(10,10,2,3),
     SliceVectorization,NoUnrolling));
 
+  VERIFY(test_assign(VectorXf(10),VectorXf(10)+VectorXf(10),
+    LinearVectorization,NoUnrolling));
 
   VERIFY(test_sum(VectorXf(10),
     LinearVectorization,NoUnrolling));