bump to 3.1.4

(transplanted from 4020d4286f
)

Eigen/SparseCholesky

@@ -2,6 +2,7 @@
 #define EIGEN_SPARSECHOLESKY_MODULE_H
 
 #include "SparseCore"
+#include "OrderingMethods"
 
 #include "src/Core/util/DisableStupidWarnings.h"
 

Eigen/src/Cholesky/LDLT.h

@@ -277,15 +277,13 @@ template<> struct ldlt_inplace<Lower>
         // are compared; if any diagonal is negligible compared
         // to the largest overall, the algorithm bails.
         cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
-
-        if(sign)
-          *sign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
       }
 
       // Finish early if the matrix is not full rank.
       if(biggest_in_corner < cutoff)
       {
         for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
+        if(sign) *sign = 0;
         break;
       }
 
@@ -326,6 +324,16 @@ template<> struct ldlt_inplace<Lower>
       }
       if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
         A21 /= mat.coeffRef(k,k);
+      
+      if(sign)
+      {
+        // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right
+        int newSign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0;
+        if(k == 0)
+          *sign = newSign;
+        else if(*sign != newSign)
+          *sign = 0;
+      }
     }
 
     return true;

Eigen/src/Core/ArrayWrapper.h

@@ -55,7 +55,7 @@ class ArrayWrapper : public ArrayBase<ArrayWrapper<ExpressionType> >
     inline Index outerStride() const { return m_expression.outerStride(); }
     inline Index innerStride() const { return m_expression.innerStride(); }
 
-    inline ScalarWithConstIfNotLvalue* data() { return m_expression.data(); }
+    inline ScalarWithConstIfNotLvalue* data() { return m_expression.const_cast_derived().data(); }
     inline const Scalar* data() const { return m_expression.data(); }
 
     inline CoeffReturnType coeff(Index row, Index col) const
@@ -175,7 +175,7 @@ class MatrixWrapper : public MatrixBase<MatrixWrapper<ExpressionType> >
     inline Index outerStride() const { return m_expression.outerStride(); }
     inline Index innerStride() const { return m_expression.innerStride(); }
 
-    inline ScalarWithConstIfNotLvalue* data() { return m_expression.data(); }
+    inline ScalarWithConstIfNotLvalue* data() { return m_expression.const_cast_derived().data(); }
     inline const Scalar* data() const { return m_expression.data(); }
 
     inline CoeffReturnType coeff(Index row, Index col) const

Eigen/src/Core/Assign_MKL.h

@@ -210,7 +210,7 @@ EIGEN_MKL_VML_DECLARE_UNARY_CALLS_LA(sqrt, Sqrt)
 EIGEN_MKL_VML_DECLARE_UNARY_CALLS_REAL(square, Sqr)
 
 // The vm*powx functions are not avaibale in the windows version of MKL.
-#ifdef _WIN32
+#ifndef _WIN32
 EIGEN_MKL_VML_DECLARE_POW_CALL(pow, vmspowx_, float, float)
 EIGEN_MKL_VML_DECLARE_POW_CALL(pow, vmdpowx_, double, double)
 EIGEN_MKL_VML_DECLARE_POW_CALL(pow, vmcpowx_, scomplex, MKL_Complex8)

Eigen/src/Core/CwiseUnaryView.h

@@ -56,8 +56,7 @@ template<typename ViewOp, typename MatrixType, typename StorageKind>
 class CwiseUnaryViewImpl;
 
 template<typename ViewOp, typename MatrixType>
-class CwiseUnaryView : internal::no_assignment_operator,
-  public CwiseUnaryViewImpl<ViewOp, MatrixType, typename internal::traits<MatrixType>::StorageKind>
+class CwiseUnaryView : public CwiseUnaryViewImpl<ViewOp, MatrixType, typename internal::traits<MatrixType>::StorageKind>
 {
   public:
 
@@ -99,6 +98,10 @@ class CwiseUnaryViewImpl<ViewOp,MatrixType,Dense>
     typedef typename internal::dense_xpr_base< CwiseUnaryView<ViewOp, MatrixType> >::type Base;
 
     EIGEN_DENSE_PUBLIC_INTERFACE(Derived)
+    EIGEN_INHERIT_ASSIGNMENT_OPERATORS(CwiseUnaryViewImpl)
+    
+    inline Scalar* data() { return &coeffRef(0); }
+    inline const Scalar* data() const { return &coeff(0); }
 
     inline Index innerStride() const
     {

Eigen/src/Core/Transpose.h

@@ -104,6 +104,7 @@ template<typename MatrixType> class TransposeImpl<MatrixType,Dense>
 
     typedef typename internal::TransposeImpl_base<MatrixType>::type Base;
     EIGEN_DENSE_PUBLIC_INTERFACE(Transpose<MatrixType>)
+    EIGEN_INHERIT_ASSIGNMENT_OPERATORS(TransposeImpl)
 
     inline Index innerStride() const { return derived().nestedExpression().innerStride(); }
     inline Index outerStride() const { return derived().nestedExpression().outerStride(); }
@@ -252,7 +253,7 @@ struct inplace_transpose_selector;
 template<typename MatrixType>
 struct inplace_transpose_selector<MatrixType,true> { // square matrix
   static void run(MatrixType& m) {
-    m.template triangularView<StrictlyUpper>().swap(m.transpose());
+    m.matrix().template triangularView<StrictlyUpper>().swap(m.matrix().transpose());
   }
 };
 
@@ -260,7 +261,7 @@ template<typename MatrixType>
 struct inplace_transpose_selector<MatrixType,false> { // non square matrix
   static void run(MatrixType& m) {
     if (m.rows()==m.cols())
-      m.template triangularView<StrictlyUpper>().swap(m.transpose());
+      m.matrix().template triangularView<StrictlyUpper>().swap(m.matrix().transpose());
     else
       m = m.transpose().eval();
   }

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

@@ -374,7 +374,7 @@ Packet4f psqrt<Packet4f>(const Packet4f& _x)
   Packet4f half = pmul(_x, pset1<Packet4f>(.5f));
 
   /* select only the inverse sqrt of non-zero inputs */
-  Packet4f non_zero_mask = _mm_cmpgt_ps(_x, pset1<Packet4f>(std::numeric_limits<float>::epsilon()));
+  Packet4f non_zero_mask = _mm_cmpgt_ps(_x, pset1<Packet4f>((std::numeric_limits<float>::min)()));
   Packet4f x = _mm_and_ps(non_zero_mask, _mm_rsqrt_ps(_x));
 
   x = pmul(x, psub(pset1<Packet4f>(1.5f), pmul(half, pmul(x,x))));

Eigen/src/Core/util/Macros.h

@@ -13,7 +13,7 @@
 
 #define EIGEN_WORLD_VERSION 3
 #define EIGEN_MAJOR_VERSION 1
-#define EIGEN_MINOR_VERSION 3
+#define EIGEN_MINOR_VERSION 4
 
 #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/SVD/JacobiSVD.h

@@ -833,17 +833,13 @@ struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
     // A = U S V^*
     // So A^{-1} = V S^{-1} U^*
 
+    Matrix<Scalar, Dynamic, Rhs::ColsAtCompileTime, 0, _MatrixType::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime> tmp;
     Index diagSize = (std::min)(dec().rows(), dec().cols());
-    typename JacobiSVDType::SingularValuesType invertedSingVals(diagSize);
-
     Index nonzeroSingVals = dec().nonzeroSingularValues();
-    invertedSingVals.head(nonzeroSingVals) = dec().singularValues().head(nonzeroSingVals).array().inverse();
-    invertedSingVals.tail(diagSize - nonzeroSingVals).setZero();
-
-    dst = dec().matrixV().leftCols(diagSize)
-        * invertedSingVals.asDiagonal()
-        * dec().matrixU().leftCols(diagSize).adjoint()
-        * rhs();
+    
+    tmp.noalias() = dec().matrixU().leftCols(nonzeroSingVals).adjoint() * rhs();
+    tmp = dec().singularValues().head(nonzeroSingVals).asDiagonal().inverse() * tmp;
+    dst = dec().matrixV().leftCols(nonzeroSingVals) * tmp;
   }
 };
 } // end namespace internal

doc/C01_TutorialMatrixClass.dox

@@ -98,8 +98,8 @@ Matrix3f a;
 MatrixXf b;
 \endcode
 Here,
-\li \c a is a 3x3 matrix, with a static float[9] array of uninitialized coefficients,
-\li \c b is a dynamic-size matrix whose size is currently 0x0, and whose array of
+\li \c a is a 3-by-3 matrix, with a plain float[9] array of uninitialized coefficients,
+\li \c b is a dynamic-size matrix whose size is currently 0-by-0, and whose array of
 coefficients hasn't yet been allocated at all.
 
 Constructors taking sizes are also available. For matrices, the number of rows is always passed first.
@@ -216,7 +216,7 @@ The simple answer is: use fixed
 sizes for very small sizes where you can, and use dynamic sizes for larger sizes or where you have to. For small sizes,
 especially for sizes smaller than (roughly) 16, using fixed sizes is hugely beneficial
 to performance, as it allows Eigen to avoid dynamic memory allocation and to unroll
-loops. Internally, a fixed-size Eigen matrix is just a plain static array, i.e. doing
+loops. Internally, a fixed-size Eigen matrix is just a plain array, i.e. doing
 \code Matrix4f mymatrix; \endcode
 really amounts to just doing
 \code float mymatrix[16]; \endcode
@@ -231,8 +231,9 @@ member variables.
 The limitation of using fixed sizes, of course, is that this is only possible
 when you know the sizes at compile time. Also, for large enough sizes, say for sizes
 greater than (roughly) 32, the performance benefit of using fixed sizes becomes negligible.
-Worse, trying to create a very large matrix using fixed sizes could result in a stack overflow,
-since Eigen will try to allocate the array as a static array, which by default goes on the stack.
+Worse, trying to create a very large matrix using fixed sizes inside a function could result in a
+stack overflow, since Eigen will try to allocate the array automatically as a local variable, and
+this is normally done on the stack.
 Finally, depending on circumstances, Eigen can also be more aggressive trying to vectorize
 (use SIMD instructions) when dynamic sizes are used, see \ref TopicVectorization "Vectorization".
 

doc/C09_TutorialSparse.dox

@@ -130,7 +130,7 @@ Describing the \a buildProblem and \a save functions is out of the scope of this
 
 The SparseMatrix and SparseVector classes take three template arguments:
  * the scalar type (e.g., double)
- * the storage order (ColMajor or RowMajor, the default is RowMajor)
+ * the storage order (ColMajor or RowMajor, the default is ColMajor)
  * the inner index type (default is \c int).
 
 As for dense Matrix objects, constructors takes the size of the object.

doc/QuickReference.dox

@@ -55,7 +55,7 @@ All combinations are allowed: you can have a matrix with a fixed number of rows
 Matrix<double, 6, Dynamic>                  // Dynamic number of columns (heap allocation)
 Matrix<double, Dynamic, 2>                  // Dynamic number of rows (heap allocation)
 Matrix<double, Dynamic, Dynamic, RowMajor>  // Fully dynamic, row major (heap allocation)
-Matrix<double, 13, 3>                       // Fully fixed (static allocation)
+Matrix<double, 13, 3>                       // Fully fixed (usually allocated on stack)
 \endcode
 
 In most cases, you can simply use one of the convenience typedefs for \ref matrixtypedefs "matrices" and \ref arraytypedefs "arrays". Some examples: