From 745bdbe7e00d32b1e8bb4d8816c7d50522059579 Mon Sep 17 00:00:00 2001
From: Gennady Pospelov <g.pospelov@fz-juelich.de>
Date: Thu, 28 Nov 2013 12:01:41 +0100
Subject: [PATCH] Update to Eigen 3.2

---
 ThirdParty/eigen3/Eigen/Core                  |  46 ++-
 .../Eigen/src/CholmodSupport/CholmodSupport.h |  43 ++-
 ThirdParty/eigen3/Eigen/src/Core/Array.h      |  32 +-
 ThirdParty/eigen3/Eigen/src/Core/ArrayBase.h  |   2 +-
 .../eigen3/Eigen/src/Core/ArrayWrapper.h      |  54 +--
 ThirdParty/eigen3/Eigen/src/Core/Assign.h     |  32 +-
 ThirdParty/eigen3/Eigen/src/Core/Assign_MKL.h |   2 +-
 ThirdParty/eigen3/Eigen/src/Core/Block.h      | 182 ++++++----
 .../eigen3/Eigen/src/Core/BooleanRedux.h      |  28 +-
 .../eigen3/Eigen/src/Core/CommaInitializer.h  |   2 +
 .../src/{SparseCore => Core}/CoreIterators.h  |   0
 .../eigen3/Eigen/src/Core/CwiseBinaryOp.h     |  22 +-
 .../eigen3/Eigen/src/Core/CwiseNullaryOp.h    | 134 +++----
 .../eigen3/Eigen/src/Core/CwiseUnaryOp.h      |   8 +-
 .../eigen3/Eigen/src/Core/CwiseUnaryView.h    |   7 +-
 ThirdParty/eigen3/Eigen/src/Core/DenseBase.h  |  72 ++--
 .../eigen3/Eigen/src/Core/arch/NEON/Complex.h |  10 +-
 .../eigen3/Eigen/src/Core/arch/SSE/Complex.h  |  22 +-
 .../src/Core/products/CoeffBasedProduct.h     |   2 +-
 .../eigen3/Eigen/src/Core/util/BlasUtil.h     |  10 +-
 .../eigen3/Eigen/src/Core/util/Constants.h    |  15 +-
 .../src/Eigen2Support/Geometry/AlignedBox.h   |   4 +-
 .../src/Eigen2Support/Geometry/AngleAxis.h    |   2 +-
 .../src/Eigenvalues/ComplexEigenSolver.h      |  26 +-
 .../Eigen/src/Eigenvalues/ComplexSchur.h      |  98 ++++--
 .../Eigen/src/Eigenvalues/ComplexSchur_MKL.h  |   2 +-
 .../eigen3/Eigen/src/Geometry/AlignedBox.h    |  12 +-
 .../eigen3/Eigen/src/Geometry/AngleAxis.h     |  13 +-
 .../src/IterativeLinearSolvers/BiCGSTAB.h     |  31 +-
 .../ConjugateGradient.h                       |  26 +-
 .../eigen3/Eigen/src/OrderingMethods/Amd.h    |   6 +-
 .../eigen3/Eigen/src/QR/ColPivHouseholderQR.h |  91 ++++-
 .../Eigen/src/QR/ColPivHouseholderQR_MKL.h    |   5 +-
 .../eigen3/Eigen/src/SparseCore/AmbiVector.h  |  10 +-
 .../Eigen/src/SparseCore/CompressedStorage.h  |   8 +-
 .../ConservativeSparseSparseProduct.h         |   6 +-
 .../Eigen/src/plugins/ArrayCwiseBinaryOps.h   |  14 +-
 .../eigen3/Eigen/src/plugins/BlockMethods.h   | 331 +++++++++++++++++-
 38 files changed, 1001 insertions(+), 409 deletions(-)
 rename ThirdParty/eigen3/Eigen/src/{SparseCore => Core}/CoreIterators.h (100%)

diff --git a/ThirdParty/eigen3/Eigen/Core b/ThirdParty/eigen3/Eigen/Core
index 34a6bcef456..9131cc3fc9d 100644
--- a/ThirdParty/eigen3/Eigen/Core
+++ b/ThirdParty/eigen3/Eigen/Core
@@ -19,6 +19,12 @@
 // defined e.g. EIGEN_DONT_ALIGN) so it needs to be done before we do anything with vectorization.
 #include "src/Core/util/Macros.h"
 
+// Disable the ipa-cp-clone optimization flag with MinGW 6.x or newer (enabled by default with -O3)
+// See http://eigen.tuxfamily.org/bz/show_bug.cgi?id=556 for details.
+#if defined(__MINGW32__) && EIGEN_GNUC_AT_LEAST(4,6)
+  #pragma GCC optimize ("-fno-ipa-cp-clone")
+#endif
+
 #include <complex>
 
 // this include file manages BLAS and MKL related macros
@@ -87,19 +93,25 @@
     // so, to avoid compile errors when windows.h is included after Eigen/Core, ensure intrinsics are extern "C" here too.
     // notice that since these are C headers, the extern "C" is theoretically needed anyways.
     extern "C" {
-      #include <emmintrin.h>
-      #include <xmmintrin.h>
-      #ifdef  EIGEN_VECTORIZE_SSE3
-      #include <pmmintrin.h>
-      #endif
-      #ifdef EIGEN_VECTORIZE_SSSE3
-      #include <tmmintrin.h>
-      #endif
-      #ifdef EIGEN_VECTORIZE_SSE4_1
-      #include <smmintrin.h>
-      #endif
-      #ifdef EIGEN_VECTORIZE_SSE4_2
-      #include <nmmintrin.h>
+      // In theory we should only include immintrin.h and not the other *mmintrin.h header files directly.
+      // Doing so triggers some issues with ICC. However old gcc versions seems to not have this file, thus:
+      #ifdef __INTEL_COMPILER
+        #include <immintrin.h>
+      #else
+        #include <emmintrin.h>
+        #include <xmmintrin.h>
+        #ifdef  EIGEN_VECTORIZE_SSE3
+        #include <pmmintrin.h>
+        #endif
+        #ifdef EIGEN_VECTORIZE_SSSE3
+        #include <tmmintrin.h>
+        #endif
+        #ifdef EIGEN_VECTORIZE_SSE4_1
+        #include <smmintrin.h>
+        #endif
+        #ifdef EIGEN_VECTORIZE_SSE4_2
+        #include <nmmintrin.h>
+        #endif
       #endif
     } // end extern "C"
   #elif defined __ALTIVEC__
@@ -236,15 +248,11 @@ using std::ptrdiff_t;
   * \endcode
   */
 
-/** \defgroup Support_modules Support modules [category]
-  * Category of modules which add support for external libraries.
-  */
-
 #include "src/Core/util/Constants.h"
 #include "src/Core/util/ForwardDeclarations.h"
 #include "src/Core/util/Meta.h"
-#include "src/Core/util/XprHelper.h"
 #include "src/Core/util/StaticAssert.h"
+#include "src/Core/util/XprHelper.h"
 #include "src/Core/util/Memory.h"
 
 #include "src/Core/NumTraits.h"
@@ -297,6 +305,7 @@ using std::ptrdiff_t;
 #include "src/Core/Map.h"
 #include "src/Core/Block.h"
 #include "src/Core/VectorBlock.h"
+#include "src/Core/Ref.h"
 #include "src/Core/Transpose.h"
 #include "src/Core/DiagonalMatrix.h"
 #include "src/Core/Diagonal.h"
@@ -330,6 +339,7 @@ using std::ptrdiff_t;
 #include "src/Core/products/TriangularSolverMatrix.h"
 #include "src/Core/products/TriangularSolverVector.h"
 #include "src/Core/BandMatrix.h"
+#include "src/Core/CoreIterators.h"
 
 #include "src/Core/BooleanRedux.h"
 #include "src/Core/Select.h"
diff --git a/ThirdParty/eigen3/Eigen/src/CholmodSupport/CholmodSupport.h b/ThirdParty/eigen3/Eigen/src/CholmodSupport/CholmodSupport.h
index 37f142150ff..783324b0b2d 100644
--- a/ThirdParty/eigen3/Eigen/src/CholmodSupport/CholmodSupport.h
+++ b/ThirdParty/eigen3/Eigen/src/CholmodSupport/CholmodSupport.h
@@ -51,7 +51,6 @@ void cholmod_configure_matrix(CholmodType& mat)
 template<typename _Scalar, int _Options, typename _Index>
 cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat)
 {
-  typedef SparseMatrix<_Scalar,_Options,_Index> MatrixType;
   cholmod_sparse res;
   res.nzmax   = mat.nonZeros();
   res.nrow    = mat.rows();;
@@ -77,9 +76,13 @@ cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat)
   {
     res.itype = CHOLMOD_INT;
   }
+  else if (internal::is_same<_Index,UF_long>::value)
+  {
+    res.itype = CHOLMOD_LONG;
+  }
   else
   {
-    eigen_assert(false && "Index type different than int is not supported yet");
+    eigen_assert(false && "Index type not supported yet");
   }
 
   // setup res.xtype
@@ -123,7 +126,7 @@ cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
   res.ncol   = mat.cols();
   res.nzmax  = res.nrow * res.ncol;
   res.d      = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
-  res.x      = mat.derived().data();
+  res.x      = (void*)(mat.derived().data());
   res.z      = 0;
 
   internal::cholmod_configure_matrix<Scalar>(res);
@@ -137,8 +140,8 @@ template<typename Scalar, int Flags, typename Index>
 MappedSparseMatrix<Scalar,Flags,Index> viewAsEigen(cholmod_sparse& cm)
 {
   return MappedSparseMatrix<Scalar,Flags,Index>
-         (cm.nrow, cm.ncol, reinterpret_cast<Index*>(cm.p)[cm.ncol],
-          reinterpret_cast<Index*>(cm.p), reinterpret_cast<Index*>(cm.i),reinterpret_cast<Scalar*>(cm.x) );
+         (cm.nrow, cm.ncol, static_cast<Index*>(cm.p)[cm.ncol],
+          static_cast<Index*>(cm.p), static_cast<Index*>(cm.i),static_cast<Scalar*>(cm.x) );
 }
 
 enum CholmodMode {
@@ -173,6 +176,7 @@ class CholmodBase : internal::noncopyable
     CholmodBase(const MatrixType& matrix)
       : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
     {
+      m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
       cholmod_start(&m_cholmod);
       compute(matrix);
     }
@@ -269,9 +273,10 @@ class CholmodBase : internal::noncopyable
     {
       eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
       cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
-      cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
+      cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod);
       
-      this->m_info = Success;
+      // If the factorization failed, minor is the column at which it did. On success minor == n.
+      this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue);
       m_factorizationIsOk = true;
     }
     
@@ -286,10 +291,12 @@ class CholmodBase : internal::noncopyable
     {
       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
       const Index size = m_cholmodFactor->n;
+      EIGEN_UNUSED_VARIABLE(size);
       eigen_assert(size==b.rows());
 
       // note: cd stands for Cholmod Dense
-      cholmod_dense b_cd = viewAsCholmod(b.const_cast_derived());
+      Rhs& b_ref(b.const_cast_derived());
+      cholmod_dense b_cd = viewAsCholmod(b_ref);
       cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
       if(!x_cd)
       {
@@ -306,6 +313,7 @@ class CholmodBase : internal::noncopyable
     {
       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
       const Index size = m_cholmodFactor->n;
+      EIGEN_UNUSED_VARIABLE(size);
       eigen_assert(size==b.rows());
 
       // note: cs stands for Cholmod Sparse
@@ -321,13 +329,30 @@ class CholmodBase : internal::noncopyable
     }
     #endif // EIGEN_PARSED_BY_DOXYGEN
     
+    
+    /** Sets the shift parameter that will be used to adjust the diagonal coefficients during the numerical factorization.
+      *
+      * During the numerical factorization, an offset term is added to the diagonal coefficients:\n
+      * \c d_ii = \a offset + \c d_ii
+      *
+      * The default is \a offset=0.
+      *
+      * \returns a reference to \c *this.
+      */
+    Derived& setShift(const RealScalar& offset)
+    {
+      m_shiftOffset[0] = offset;
+      return derived();
+    }
+    
     template<typename Stream>
-    void dumpMemory(Stream& s)
+    void dumpMemory(Stream& /*s*/)
     {}
     
   protected:
     mutable cholmod_common m_cholmod;
     cholmod_factor* m_cholmodFactor;
+    RealScalar m_shiftOffset[2];
     mutable ComputationInfo m_info;
     bool m_isInitialized;
     int m_factorizationIsOk;
diff --git a/ThirdParty/eigen3/Eigen/src/Core/Array.h b/ThirdParty/eigen3/Eigen/src/Core/Array.h
index aaa38997838..497efff66f2 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/Array.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/Array.h
@@ -107,10 +107,10 @@ class Array
       *
       * \sa resize(Index,Index)
       */
-    EIGEN_STRONG_INLINE explicit Array() : Base()
+    EIGEN_STRONG_INLINE Array() : Base()
     {
       Base::_check_template_params();
-      EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
+      EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED
     }
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
@@ -120,7 +120,7 @@ class Array
       : Base(internal::constructor_without_unaligned_array_assert())
     {
       Base::_check_template_params();
-      EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
+      EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED
     }
 #endif
 
@@ -137,15 +137,15 @@ class Array
       EIGEN_STATIC_ASSERT_VECTOR_ONLY(Array)
       eigen_assert(dim >= 0);
       eigen_assert(SizeAtCompileTime == Dynamic || SizeAtCompileTime == dim);
-      EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
+      EIGEN_INITIALIZE_COEFFS_IF_THAT_OPTION_IS_ENABLED
     }
 
     #ifndef EIGEN_PARSED_BY_DOXYGEN
     template<typename T0, typename T1>
-    EIGEN_STRONG_INLINE Array(const T0& x, const T1& y)
+    EIGEN_STRONG_INLINE Array(const T0& val0, const T1& val1)
     {
       Base::_check_template_params();
-      this->template _init2<T0,T1>(x, y);
+      this->template _init2<T0,T1>(val0, val1);
     }
     #else
     /** constructs an uninitialized matrix with \a rows rows and \a cols columns.
@@ -155,27 +155,27 @@ class Array
       * Matrix() instead. */
     Array(Index rows, Index cols);
     /** constructs an initialized 2D vector with given coefficients */
-    Array(const Scalar& x, const Scalar& y);
+    Array(const Scalar& val0, const Scalar& val1);
     #endif
 
     /** constructs an initialized 3D vector with given coefficients */
-    EIGEN_STRONG_INLINE Array(const Scalar& x, const Scalar& y, const Scalar& z)
+    EIGEN_STRONG_INLINE Array(const Scalar& val0, const Scalar& val1, const Scalar& val2)
     {
       Base::_check_template_params();
       EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Array, 3)
-      m_storage.data()[0] = x;
-      m_storage.data()[1] = y;
-      m_storage.data()[2] = z;
+      m_storage.data()[0] = val0;
+      m_storage.data()[1] = val1;
+      m_storage.data()[2] = val2;
     }
     /** constructs an initialized 4D vector with given coefficients */
-    EIGEN_STRONG_INLINE Array(const Scalar& x, const Scalar& y, const Scalar& z, const Scalar& w)
+    EIGEN_STRONG_INLINE Array(const Scalar& val0, const Scalar& val1, const Scalar& val2, const Scalar& val3)
     {
       Base::_check_template_params();
       EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(Array, 4)
-      m_storage.data()[0] = x;
-      m_storage.data()[1] = y;
-      m_storage.data()[2] = z;
-      m_storage.data()[3] = w;
+      m_storage.data()[0] = val0;
+      m_storage.data()[1] = val1;
+      m_storage.data()[2] = val2;
+      m_storage.data()[3] = val3;
     }
 
     explicit Array(const Scalar *data);
diff --git a/ThirdParty/eigen3/Eigen/src/Core/ArrayBase.h b/ThirdParty/eigen3/Eigen/src/Core/ArrayBase.h
index 004b117c933..38852600dc2 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/ArrayBase.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/ArrayBase.h
@@ -143,7 +143,7 @@ template<typename Derived> class ArrayBase
     ArrayBase<Derived>& array() { return *this; }
     const ArrayBase<Derived>& array() const { return *this; }
 
-    /** \returns an \link MatrixBase Matrix \endlink expression of this array
+    /** \returns an \link Eigen::MatrixBase Matrix \endlink expression of this array
       * \sa MatrixBase::array() */
     MatrixWrapper<Derived> matrix() { return derived(); }
     const MatrixWrapper<const Derived> matrix() const { return derived(); }
diff --git a/ThirdParty/eigen3/Eigen/src/Core/ArrayWrapper.h b/ThirdParty/eigen3/Eigen/src/Core/ArrayWrapper.h
index 65ffd64cabd..a791bc3581a 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/ArrayWrapper.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/ArrayWrapper.h
@@ -55,22 +55,22 @@ 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
+    inline CoeffReturnType coeff(Index rowId, Index colId) const
     {
-      return m_expression.coeff(row, col);
+      return m_expression.coeff(rowId, colId);
     }
 
-    inline Scalar& coeffRef(Index row, Index col)
+    inline Scalar& coeffRef(Index rowId, Index colId)
     {
-      return m_expression.const_cast_derived().coeffRef(row, col);
+      return m_expression.const_cast_derived().coeffRef(rowId, colId);
     }
 
-    inline const Scalar& coeffRef(Index row, Index col) const
+    inline const Scalar& coeffRef(Index rowId, Index colId) const
     {
-      return m_expression.const_cast_derived().coeffRef(row, col);
+      return m_expression.const_cast_derived().coeffRef(rowId, colId);
     }
 
     inline CoeffReturnType coeff(Index index) const
@@ -89,15 +89,15 @@ class ArrayWrapper : public ArrayBase<ArrayWrapper<ExpressionType> >
     }
 
     template<int LoadMode>
-    inline const PacketScalar packet(Index row, Index col) const
+    inline const PacketScalar packet(Index rowId, Index colId) const
     {
-      return m_expression.template packet<LoadMode>(row, col);
+      return m_expression.template packet<LoadMode>(rowId, colId);
     }
 
     template<int LoadMode>
-    inline void writePacket(Index row, Index col, const PacketScalar& x)
+    inline void writePacket(Index rowId, Index colId, const PacketScalar& val)
     {
-      m_expression.const_cast_derived().template writePacket<LoadMode>(row, col, x);
+      m_expression.const_cast_derived().template writePacket<LoadMode>(rowId, colId, val);
     }
 
     template<int LoadMode>
@@ -107,9 +107,9 @@ class ArrayWrapper : public ArrayBase<ArrayWrapper<ExpressionType> >
     }
 
     template<int LoadMode>
-    inline void writePacket(Index index, const PacketScalar& x)
+    inline void writePacket(Index index, const PacketScalar& val)
     {
-      m_expression.const_cast_derived().template writePacket<LoadMode>(index, x);
+      m_expression.const_cast_derived().template writePacket<LoadMode>(index, val);
     }
 
     template<typename Dest>
@@ -168,29 +168,29 @@ class MatrixWrapper : public MatrixBase<MatrixWrapper<ExpressionType> >
 
     typedef typename internal::nested<ExpressionType>::type NestedExpressionType;
 
-    inline MatrixWrapper(ExpressionType& matrix) : m_expression(matrix) {}
+    inline MatrixWrapper(ExpressionType& a_matrix) : m_expression(a_matrix) {}
 
     inline Index rows() const { return m_expression.rows(); }
     inline Index cols() const { return m_expression.cols(); }
     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
+    inline CoeffReturnType coeff(Index rowId, Index colId) const
     {
-      return m_expression.coeff(row, col);
+      return m_expression.coeff(rowId, colId);
     }
 
-    inline Scalar& coeffRef(Index row, Index col)
+    inline Scalar& coeffRef(Index rowId, Index colId)
     {
-      return m_expression.const_cast_derived().coeffRef(row, col);
+      return m_expression.const_cast_derived().coeffRef(rowId, colId);
     }
 
-    inline const Scalar& coeffRef(Index row, Index col) const
+    inline const Scalar& coeffRef(Index rowId, Index colId) const
     {
-      return m_expression.derived().coeffRef(row, col);
+      return m_expression.derived().coeffRef(rowId, colId);
     }
 
     inline CoeffReturnType coeff(Index index) const
@@ -209,15 +209,15 @@ class MatrixWrapper : public MatrixBase<MatrixWrapper<ExpressionType> >
     }
 
     template<int LoadMode>
-    inline const PacketScalar packet(Index row, Index col) const
+    inline const PacketScalar packet(Index rowId, Index colId) const
     {
-      return m_expression.template packet<LoadMode>(row, col);
+      return m_expression.template packet<LoadMode>(rowId, colId);
     }
 
     template<int LoadMode>
-    inline void writePacket(Index row, Index col, const PacketScalar& x)
+    inline void writePacket(Index rowId, Index colId, const PacketScalar& val)
     {
-      m_expression.const_cast_derived().template writePacket<LoadMode>(row, col, x);
+      m_expression.const_cast_derived().template writePacket<LoadMode>(rowId, colId, val);
     }
 
     template<int LoadMode>
@@ -227,9 +227,9 @@ class MatrixWrapper : public MatrixBase<MatrixWrapper<ExpressionType> >
     }
 
     template<int LoadMode>
-    inline void writePacket(Index index, const PacketScalar& x)
+    inline void writePacket(Index index, const PacketScalar& val)
     {
-      m_expression.const_cast_derived().template writePacket<LoadMode>(index, x);
+      m_expression.const_cast_derived().template writePacket<LoadMode>(index, val);
     }
 
     const typename internal::remove_all<NestedExpressionType>::type& 
diff --git a/ThirdParty/eigen3/Eigen/src/Core/Assign.h b/ThirdParty/eigen3/Eigen/src/Core/Assign.h
index cd29a88f0da..1dccc2f4212 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/Assign.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/Assign.h
@@ -155,7 +155,7 @@ struct assign_DefaultTraversal_CompleteUnrolling<Derived1, Derived2, Stop, Stop>
 template<typename Derived1, typename Derived2, int Index, int Stop>
 struct assign_DefaultTraversal_InnerUnrolling
 {
-  static EIGEN_STRONG_INLINE void run(Derived1 &dst, const Derived2 &src, int outer)
+  static EIGEN_STRONG_INLINE void run(Derived1 &dst, const Derived2 &src, typename Derived1::Index outer)
   {
     dst.copyCoeffByOuterInner(outer, Index, src);
     assign_DefaultTraversal_InnerUnrolling<Derived1, Derived2, Index+1, Stop>::run(dst, src, outer);
@@ -165,7 +165,7 @@ struct assign_DefaultTraversal_InnerUnrolling
 template<typename Derived1, typename Derived2, int Stop>
 struct assign_DefaultTraversal_InnerUnrolling<Derived1, Derived2, Stop, Stop>
 {
-  static EIGEN_STRONG_INLINE void run(Derived1 &, const Derived2 &, int) {}
+  static EIGEN_STRONG_INLINE void run(Derived1 &, const Derived2 &, typename Derived1::Index) {}
 };
 
 /***********************
@@ -218,7 +218,7 @@ struct assign_innervec_CompleteUnrolling<Derived1, Derived2, Stop, Stop>
 template<typename Derived1, typename Derived2, int Index, int Stop>
 struct assign_innervec_InnerUnrolling
 {
-  static EIGEN_STRONG_INLINE void run(Derived1 &dst, const Derived2 &src, int outer)
+  static EIGEN_STRONG_INLINE void run(Derived1 &dst, const Derived2 &src, typename Derived1::Index outer)
   {
     dst.template copyPacketByOuterInner<Derived2, Aligned, Aligned>(outer, Index, src);
     assign_innervec_InnerUnrolling<Derived1, Derived2,
@@ -229,7 +229,7 @@ struct assign_innervec_InnerUnrolling
 template<typename Derived1, typename Derived2, int Stop>
 struct assign_innervec_InnerUnrolling<Derived1, Derived2, Stop, Stop>
 {
-  static EIGEN_STRONG_INLINE void run(Derived1 &, const Derived2 &, int) {}
+  static EIGEN_STRONG_INLINE void run(Derived1 &, const Derived2 &, typename Derived1::Index) {}
 };
 
 /***************************************************************************
@@ -507,19 +507,19 @@ EIGEN_STRONG_INLINE Derived& DenseBase<Derived>
 namespace internal {
 
 template<typename Derived, typename OtherDerived,
-         bool EvalBeforeAssigning = (int(OtherDerived::Flags) & EvalBeforeAssigningBit) != 0,
-         bool NeedToTranspose = Derived::IsVectorAtCompileTime
-                && OtherDerived::IsVectorAtCompileTime
-                && ((int(Derived::RowsAtCompileTime) == 1 && int(OtherDerived::ColsAtCompileTime) == 1)
-                      |  // FIXME | instead of || to please GCC 4.4.0 stupid warning "suggest parentheses around &&".
-                         // revert to || as soon as not needed anymore.
-                    (int(Derived::ColsAtCompileTime) == 1 && int(OtherDerived::RowsAtCompileTime) == 1))
-                && int(Derived::SizeAtCompileTime) != 1>
+         bool EvalBeforeAssigning = (int(internal::traits<OtherDerived>::Flags) & EvalBeforeAssigningBit) != 0,
+         bool NeedToTranspose = ((int(Derived::RowsAtCompileTime) == 1 && int(OtherDerived::ColsAtCompileTime) == 1)
+                              |   // FIXME | instead of || to please GCC 4.4.0 stupid warning "suggest parentheses around &&".
+                                  // revert to || as soon as not needed anymore.
+                                  (int(Derived::ColsAtCompileTime) == 1 && int(OtherDerived::RowsAtCompileTime) == 1))
+                              && int(Derived::SizeAtCompileTime) != 1>
 struct assign_selector;
 
 template<typename Derived, typename OtherDerived>
 struct assign_selector<Derived,OtherDerived,false,false> {
   static EIGEN_STRONG_INLINE Derived& run(Derived& dst, const OtherDerived& other) { return dst.lazyAssign(other.derived()); }
+  template<typename ActualDerived, typename ActualOtherDerived>
+  static EIGEN_STRONG_INLINE Derived& evalTo(ActualDerived& dst, const ActualOtherDerived& other) { other.evalTo(dst); return dst; }
 };
 template<typename Derived, typename OtherDerived>
 struct assign_selector<Derived,OtherDerived,true,false> {
@@ -528,6 +528,8 @@ struct assign_selector<Derived,OtherDerived,true,false> {
 template<typename Derived, typename OtherDerived>
 struct assign_selector<Derived,OtherDerived,false,true> {
   static EIGEN_STRONG_INLINE Derived& run(Derived& dst, const OtherDerived& other) { return dst.lazyAssign(other.transpose()); }
+  template<typename ActualDerived, typename ActualOtherDerived>
+  static EIGEN_STRONG_INLINE Derived& evalTo(ActualDerived& dst, const ActualOtherDerived& other) { Transpose<ActualDerived> dstTrans(dst); other.evalTo(dstTrans); return dst; }
 };
 template<typename Derived, typename OtherDerived>
 struct assign_selector<Derived,OtherDerived,true,true> {
@@ -566,16 +568,14 @@ template<typename Derived>
 template <typename OtherDerived>
 EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const EigenBase<OtherDerived>& other)
 {
-  other.derived().evalTo(derived());
-  return derived();
+  return internal::assign_selector<Derived,OtherDerived,false>::evalTo(derived(), other.derived());
 }
 
 template<typename Derived>
 template<typename OtherDerived>
 EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const ReturnByValue<OtherDerived>& other)
 {
-  other.evalTo(derived());
-  return derived();
+  return internal::assign_selector<Derived,OtherDerived,false>::evalTo(derived(), other.derived());
 }
 
 } // end namespace Eigen
diff --git a/ThirdParty/eigen3/Eigen/src/Core/Assign_MKL.h b/ThirdParty/eigen3/Eigen/src/Core/Assign_MKL.h
index 428c6367b92..7772951b915 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/Assign_MKL.h
+++ b/ThirdParty/eigen3/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)
diff --git a/ThirdParty/eigen3/Eigen/src/Core/Block.h b/ThirdParty/eigen3/Eigen/src/Core/Block.h
index 5f29cb3d1b3..358b3188b38 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/Block.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/Block.h
@@ -21,7 +21,6 @@ namespace Eigen {
   * \param XprType the type of the expression in which we are taking a block
   * \param BlockRows the number of rows of the block we are taking at compile time (optional)
   * \param BlockCols the number of columns of the block we are taking at compile time (optional)
-  * \param _DirectAccessStatus \internal used for partial specialization
   *
   * This class represents an expression of either a fixed-size or dynamic-size block. It is the return
   * type of DenseBase::block(Index,Index,Index,Index) and DenseBase::block<int,int>(Index,Index) and
@@ -47,8 +46,8 @@ namespace Eigen {
   */
 
 namespace internal {
-template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel, bool HasDirectAccess>
-struct traits<Block<XprType, BlockRows, BlockCols, InnerPanel, HasDirectAccess> > : traits<XprType>
+template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel>
+struct traits<Block<XprType, BlockRows, BlockCols, InnerPanel> > : traits<XprType>
 {
   typedef typename traits<XprType>::Scalar Scalar;
   typedef typename traits<XprType>::StorageKind StorageKind;
@@ -92,30 +91,27 @@ struct traits<Block<XprType, BlockRows, BlockCols, InnerPanel, HasDirectAccess>
     Flags = Flags0 | FlagsLinearAccessBit | FlagsLvalueBit | FlagsRowMajorBit
   };
 };
-}
-
-template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel, bool HasDirectAccess> class Block
-  : public internal::dense_xpr_base<Block<XprType, BlockRows, BlockCols, InnerPanel, HasDirectAccess> >::type
-{
-  public:
 
-    typedef typename internal::dense_xpr_base<Block>::type Base;
-    EIGEN_DENSE_PUBLIC_INTERFACE(Block)
+template<typename XprType, int BlockRows=Dynamic, int BlockCols=Dynamic, bool InnerPanel = false,
+         bool HasDirectAccess = internal::has_direct_access<XprType>::ret> class BlockImpl_dense;
+         
+} // end namespace internal
 
-    class InnerIterator;
+template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel, typename StorageKind> class BlockImpl;
 
+template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel> class Block
+  : public BlockImpl<XprType, BlockRows, BlockCols, InnerPanel, typename internal::traits<XprType>::StorageKind>
+{
+    typedef BlockImpl<XprType, BlockRows, BlockCols, InnerPanel, typename internal::traits<XprType>::StorageKind> Impl;
+  public:
+    //typedef typename Impl::Base Base;
+    typedef Impl Base;
+    EIGEN_GENERIC_PUBLIC_INTERFACE(Block)
+    EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Block)
+  
     /** Column or Row constructor
       */
-    inline Block(XprType& xpr, Index i)
-      : m_xpr(xpr),
-        // It is a row if and only if BlockRows==1 and BlockCols==XprType::ColsAtCompileTime,
-        // and it is a column if and only if BlockRows==XprType::RowsAtCompileTime and BlockCols==1,
-        // all other cases are invalid.
-        // The case a 1x1 matrix seems ambiguous, but the result is the same anyway.
-        m_startRow( (BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) ? i : 0),
-        m_startCol( (BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) ? i : 0),
-        m_blockRows(BlockRows==1 ? 1 : xpr.rows()),
-        m_blockCols(BlockCols==1 ? 1 : xpr.cols())
+    inline Block(XprType& xpr, Index i) : Impl(xpr,i)
     {
       eigen_assert( (i>=0) && (
           ((BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) && i<xpr.rows())
@@ -124,50 +120,109 @@ template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel, bool H
 
     /** Fixed-size constructor
       */
-    inline Block(XprType& xpr, Index startRow, Index startCol)
-      : m_xpr(xpr), m_startRow(startRow), m_startCol(startCol),
-        m_blockRows(BlockRows), m_blockCols(BlockCols)
+    inline Block(XprType& xpr, Index a_startRow, Index a_startCol)
+      : Impl(xpr, a_startRow, a_startCol)
     {
       EIGEN_STATIC_ASSERT(RowsAtCompileTime!=Dynamic && ColsAtCompileTime!=Dynamic,THIS_METHOD_IS_ONLY_FOR_FIXED_SIZE)
-      eigen_assert(startRow >= 0 && BlockRows >= 1 && startRow + BlockRows <= xpr.rows()
-             && startCol >= 0 && BlockCols >= 1 && startCol + BlockCols <= xpr.cols());
+      eigen_assert(a_startRow >= 0 && BlockRows >= 1 && a_startRow + BlockRows <= xpr.rows()
+             && a_startCol >= 0 && BlockCols >= 1 && a_startCol + BlockCols <= xpr.cols());
     }
 
     /** Dynamic-size constructor
       */
     inline Block(XprType& xpr,
-          Index startRow, Index startCol,
+          Index a_startRow, Index a_startCol,
           Index blockRows, Index blockCols)
-      : m_xpr(xpr), m_startRow(startRow), m_startCol(startCol),
-                          m_blockRows(blockRows), m_blockCols(blockCols)
+      : Impl(xpr, a_startRow, a_startCol, blockRows, blockCols)
     {
       eigen_assert((RowsAtCompileTime==Dynamic || RowsAtCompileTime==blockRows)
           && (ColsAtCompileTime==Dynamic || ColsAtCompileTime==blockCols));
-      eigen_assert(startRow >= 0 && blockRows >= 0 && startRow + blockRows <= xpr.rows()
-          && startCol >= 0 && blockCols >= 0 && startCol + blockCols <= xpr.cols());
+      eigen_assert(a_startRow >= 0 && blockRows >= 0 && a_startRow  <= xpr.rows() - blockRows
+          && a_startCol >= 0 && blockCols >= 0 && a_startCol <= xpr.cols() - blockCols);
     }
+};
+         
+// The generic default implementation for dense block simplu forward to the internal::BlockImpl_dense
+// that must be specialized for direct and non-direct access...
+template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel>
+class BlockImpl<XprType, BlockRows, BlockCols, InnerPanel, Dense>
+  : public internal::BlockImpl_dense<XprType, BlockRows, BlockCols, InnerPanel>
+{
+    typedef internal::BlockImpl_dense<XprType, BlockRows, BlockCols, InnerPanel> Impl;
+    typedef typename XprType::Index Index;
+  public:
+    typedef Impl Base;
+    EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl)
+    inline BlockImpl(XprType& xpr, Index i) : Impl(xpr,i) {}
+    inline BlockImpl(XprType& xpr, Index a_startRow, Index a_startCol) : Impl(xpr, a_startRow, a_startCol) {}
+    inline BlockImpl(XprType& xpr, Index a_startRow, Index a_startCol, Index blockRows, Index blockCols)
+      : Impl(xpr, a_startRow, a_startCol, blockRows, blockCols) {}
+};
 
-    EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Block)
+namespace internal {
+
+/** \internal Internal implementation of dense Blocks in the general case. */
+template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel, bool HasDirectAccess> class BlockImpl_dense
+  : public internal::dense_xpr_base<Block<XprType, BlockRows, BlockCols, InnerPanel> >::type
+{
+    typedef Block<XprType, BlockRows, BlockCols, InnerPanel> BlockType;
+  public:
+
+    typedef typename internal::dense_xpr_base<BlockType>::type Base;
+    EIGEN_DENSE_PUBLIC_INTERFACE(BlockType)
+    EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl_dense)
+
+    class InnerIterator;
+
+    /** Column or Row constructor
+      */
+    inline BlockImpl_dense(XprType& xpr, Index i)
+      : m_xpr(xpr),
+        // It is a row if and only if BlockRows==1 and BlockCols==XprType::ColsAtCompileTime,
+        // and it is a column if and only if BlockRows==XprType::RowsAtCompileTime and BlockCols==1,
+        // all other cases are invalid.
+        // The case a 1x1 matrix seems ambiguous, but the result is the same anyway.
+        m_startRow( (BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) ? i : 0),
+        m_startCol( (BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) ? i : 0),
+        m_blockRows(BlockRows==1 ? 1 : xpr.rows()),
+        m_blockCols(BlockCols==1 ? 1 : xpr.cols())
+    {}
+
+    /** Fixed-size constructor
+      */
+    inline BlockImpl_dense(XprType& xpr, Index a_startRow, Index a_startCol)
+      : m_xpr(xpr), m_startRow(a_startRow), m_startCol(a_startCol),
+                    m_blockRows(BlockRows), m_blockCols(BlockCols)
+    {}
+
+    /** Dynamic-size constructor
+      */
+    inline BlockImpl_dense(XprType& xpr,
+          Index a_startRow, Index a_startCol,
+          Index blockRows, Index blockCols)
+      : m_xpr(xpr), m_startRow(a_startRow), m_startCol(a_startCol),
+                    m_blockRows(blockRows), m_blockCols(blockCols)
+    {}
 
     inline Index rows() const { return m_blockRows.value(); }
     inline Index cols() const { return m_blockCols.value(); }
 
-    inline Scalar& coeffRef(Index row, Index col)
+    inline Scalar& coeffRef(Index rowId, Index colId)
     {
       EIGEN_STATIC_ASSERT_LVALUE(XprType)
       return m_xpr.const_cast_derived()
-               .coeffRef(row + m_startRow.value(), col + m_startCol.value());
+               .coeffRef(rowId + m_startRow.value(), colId + m_startCol.value());
     }
 
-    inline const Scalar& coeffRef(Index row, Index col) const
+    inline const Scalar& coeffRef(Index rowId, Index colId) const
     {
       return m_xpr.derived()
-               .coeffRef(row + m_startRow.value(), col + m_startCol.value());
+               .coeffRef(rowId + m_startRow.value(), colId + m_startCol.value());
     }
 
-    EIGEN_STRONG_INLINE const CoeffReturnType coeff(Index row, Index col) const
+    EIGEN_STRONG_INLINE const CoeffReturnType coeff(Index rowId, Index colId) const
     {
-      return m_xpr.coeff(row + m_startRow.value(), col + m_startCol.value());
+      return m_xpr.coeff(rowId + m_startRow.value(), colId + m_startCol.value());
     }
 
     inline Scalar& coeffRef(Index index)
@@ -193,17 +248,17 @@ template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel, bool H
     }
 
     template<int LoadMode>
-    inline PacketScalar packet(Index row, Index col) const
+    inline PacketScalar packet(Index rowId, Index colId) const
     {
       return m_xpr.template packet<Unaligned>
-              (row + m_startRow.value(), col + m_startCol.value());
+              (rowId + m_startRow.value(), colId + m_startCol.value());
     }
 
     template<int LoadMode>
-    inline void writePacket(Index row, Index col, const PacketScalar& x)
+    inline void writePacket(Index rowId, Index colId, const PacketScalar& val)
     {
       m_xpr.const_cast_derived().template writePacket<Unaligned>
-              (row + m_startRow.value(), col + m_startCol.value(), x);
+              (rowId + m_startRow.value(), colId + m_startCol.value(), val);
     }
 
     template<int LoadMode>
@@ -215,11 +270,11 @@ template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel, bool H
     }
 
     template<int LoadMode>
-    inline void writePacket(Index index, const PacketScalar& x)
+    inline void writePacket(Index index, const PacketScalar& val)
     {
       m_xpr.const_cast_derived().template writePacket<Unaligned>
          (m_startRow.value() + (RowsAtCompileTime == 1 ? 0 : index),
-          m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0), x);
+          m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0), val);
     }
 
     #ifdef EIGEN_PARSED_BY_DOXYGEN
@@ -253,21 +308,21 @@ template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel, bool H
     const internal::variable_if_dynamic<Index, ColsAtCompileTime> m_blockCols;
 };
 
-/** \internal */
+/** \internal Internal implementation of dense Blocks in the direct access case.*/
 template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel>
-class Block<XprType,BlockRows,BlockCols, InnerPanel,true>
-  : public MapBase<Block<XprType, BlockRows, BlockCols, InnerPanel, true> >
+class BlockImpl_dense<XprType,BlockRows,BlockCols, InnerPanel,true>
+  : public MapBase<Block<XprType, BlockRows, BlockCols, InnerPanel> >
 {
+    typedef Block<XprType, BlockRows, BlockCols, InnerPanel> BlockType;
   public:
 
-    typedef MapBase<Block> Base;
-    EIGEN_DENSE_PUBLIC_INTERFACE(Block)
-
-    EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Block)
+    typedef MapBase<BlockType> Base;
+    EIGEN_DENSE_PUBLIC_INTERFACE(BlockType)
+    EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl_dense)
 
     /** Column or Row constructor
       */
-    inline Block(XprType& xpr, Index i)
+    inline BlockImpl_dense(XprType& xpr, Index i)
       : Base(internal::const_cast_ptr(&xpr.coeffRef(
               (BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) ? i : 0,
               (BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) ? i : 0)),
@@ -275,34 +330,25 @@ class Block<XprType,BlockRows,BlockCols, InnerPanel,true>
              BlockCols==1 ? 1 : xpr.cols()),
         m_xpr(xpr)
     {
-      eigen_assert( (i>=0) && (
-          ((BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) && i<xpr.rows())
-        ||((BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) && i<xpr.cols())));
       init();
     }
 
     /** Fixed-size constructor
       */
-    inline Block(XprType& xpr, Index startRow, Index startCol)
+    inline BlockImpl_dense(XprType& xpr, Index startRow, Index startCol)
       : Base(internal::const_cast_ptr(&xpr.coeffRef(startRow,startCol))), m_xpr(xpr)
     {
-      eigen_assert(startRow >= 0 && BlockRows >= 1 && startRow + BlockRows <= xpr.rows()
-             && startCol >= 0 && BlockCols >= 1 && startCol + BlockCols <= xpr.cols());
       init();
     }
 
     /** Dynamic-size constructor
       */
-    inline Block(XprType& xpr,
+    inline BlockImpl_dense(XprType& xpr,
           Index startRow, Index startCol,
           Index blockRows, Index blockCols)
       : Base(internal::const_cast_ptr(&xpr.coeffRef(startRow,startCol)), blockRows, blockCols),
         m_xpr(xpr)
     {
-      eigen_assert((RowsAtCompileTime==Dynamic || RowsAtCompileTime==blockRows)
-             && (ColsAtCompileTime==Dynamic || ColsAtCompileTime==blockCols));
-      eigen_assert(startRow >= 0 && blockRows >= 0 && startRow + blockRows <= xpr.rows()
-             && startCol >= 0 && blockCols >= 0 && startCol + blockCols <= xpr.cols());
       init();
     }
 
@@ -314,7 +360,7 @@ class Block<XprType,BlockRows,BlockCols, InnerPanel,true>
     /** \sa MapBase::innerStride() */
     inline Index innerStride() const
     {
-      return internal::traits<Block>::HasSameStorageOrderAsXprType
+      return internal::traits<BlockType>::HasSameStorageOrderAsXprType
              ? m_xpr.innerStride()
              : m_xpr.outerStride();
     }
@@ -333,7 +379,7 @@ class Block<XprType,BlockRows,BlockCols, InnerPanel,true>
 
     #ifndef EIGEN_PARSED_BY_DOXYGEN
     /** \internal used by allowAligned() */
-    inline Block(XprType& xpr, const Scalar* data, Index blockRows, Index blockCols)
+    inline BlockImpl_dense(XprType& xpr, const Scalar* data, Index blockRows, Index blockCols)
       : Base(data, blockRows, blockCols), m_xpr(xpr)
     {
       init();
@@ -343,7 +389,7 @@ class Block<XprType,BlockRows,BlockCols, InnerPanel,true>
   protected:
     void init()
     {
-      m_outerStride = internal::traits<Block>::HasSameStorageOrderAsXprType
+      m_outerStride = internal::traits<BlockType>::HasSameStorageOrderAsXprType
                     ? m_xpr.outerStride()
                     : m_xpr.innerStride();
     }
@@ -352,6 +398,8 @@ class Block<XprType,BlockRows,BlockCols, InnerPanel,true>
     Index m_outerStride;
 };
 
+} // end namespace internal
+
 } // end namespace Eigen
 
 #endif // EIGEN_BLOCK_H
diff --git a/ThirdParty/eigen3/Eigen/src/Core/BooleanRedux.h b/ThirdParty/eigen3/Eigen/src/Core/BooleanRedux.h
index 57efd8e6953..6e37e031a8c 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/BooleanRedux.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/BooleanRedux.h
@@ -85,9 +85,7 @@ inline bool DenseBase<Derived>::all() const
           && SizeAtCompileTime * (CoeffReadCost + NumTraits<Scalar>::AddCost) <= EIGEN_UNROLLING_LIMIT
   };
   if(unroll)
-    return internal::all_unroller<Derived,
-                           unroll ? int(SizeAtCompileTime) : Dynamic
-     >::run(derived());
+    return internal::all_unroller<Derived, unroll ? int(SizeAtCompileTime) : Dynamic>::run(derived());
   else
   {
     for(Index j = 0; j < cols(); ++j)
@@ -111,9 +109,7 @@ inline bool DenseBase<Derived>::any() const
           && SizeAtCompileTime * (CoeffReadCost + NumTraits<Scalar>::AddCost) <= EIGEN_UNROLLING_LIMIT
   };
   if(unroll)
-    return internal::any_unroller<Derived,
-                           unroll ? int(SizeAtCompileTime) : Dynamic
-           >::run(derived());
+    return internal::any_unroller<Derived, unroll ? int(SizeAtCompileTime) : Dynamic>::run(derived());
   else
   {
     for(Index j = 0; j < cols(); ++j)
@@ -133,6 +129,26 @@ inline typename DenseBase<Derived>::Index DenseBase<Derived>::count() const
   return derived().template cast<bool>().template cast<Index>().sum();
 }
 
+/** \returns true is \c *this contains at least one Not A Number (NaN).
+  *
+  * \sa allFinite()
+  */
+template<typename Derived>
+inline bool DenseBase<Derived>::hasNaN() const
+{
+  return !((derived().array()==derived().array()).all());
+}
+
+/** \returns true if \c *this contains only finite numbers, i.e., no NaN and no +/-INF values.
+  *
+  * \sa hasNaN()
+  */
+template<typename Derived>
+inline bool DenseBase<Derived>::allFinite() const
+{
+  return !((derived()-derived()).hasNaN());
+}
+    
 } // end namespace Eigen
 
 #endif // EIGEN_ALLANDANY_H
diff --git a/ThirdParty/eigen3/Eigen/src/Core/CommaInitializer.h b/ThirdParty/eigen3/Eigen/src/Core/CommaInitializer.h
index f20c1774c4c..a96867af4d5 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/CommaInitializer.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/CommaInitializer.h
@@ -118,6 +118,8 @@ struct CommaInitializer
   *
   * Example: \include MatrixBase_set.cpp
   * Output: \verbinclude MatrixBase_set.out
+  * 
+  * \note According the c++ standard, the argument expressions of this comma initializer are evaluated in arbitrary order.
   *
   * \sa CommaInitializer::finished(), class CommaInitializer
   */
diff --git a/ThirdParty/eigen3/Eigen/src/SparseCore/CoreIterators.h b/ThirdParty/eigen3/Eigen/src/Core/CoreIterators.h
similarity index 100%
rename from ThirdParty/eigen3/Eigen/src/SparseCore/CoreIterators.h
rename to ThirdParty/eigen3/Eigen/src/Core/CoreIterators.h
diff --git a/ThirdParty/eigen3/Eigen/src/Core/CwiseBinaryOp.h b/ThirdParty/eigen3/Eigen/src/Core/CwiseBinaryOp.h
index 1b93af31b60..586f77aaf32 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/CwiseBinaryOp.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/CwiseBinaryOp.h
@@ -94,8 +94,8 @@ struct traits<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
 // So allowing mixing different types gives very unexpected errors when enabling vectorization, when the user tries to
 // add together a float matrix and a double matrix.
 #define EIGEN_CHECK_BINARY_COMPATIBILIY(BINOP,LHS,RHS) \
-  EIGEN_STATIC_ASSERT((internal::functor_allows_mixing_real_and_complex<BINOP>::ret \
-                        ? int(internal::is_same<typename NumTraits<LHS>::Real, typename NumTraits<RHS>::Real>::value) \
+  EIGEN_STATIC_ASSERT((internal::functor_is_product_like<BINOP>::ret \
+                        ? int(internal::scalar_product_traits<LHS, RHS>::Defined) \
                         : int(internal::is_same<LHS, RHS>::value)), \
     YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
 
@@ -122,13 +122,13 @@ class CwiseBinaryOp : internal::no_assignment_operator,
     typedef typename internal::remove_reference<LhsNested>::type _LhsNested;
     typedef typename internal::remove_reference<RhsNested>::type _RhsNested;
 
-    EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp())
-      : m_lhs(lhs), m_rhs(rhs), m_functor(func)
+    EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& aLhs, const Rhs& aRhs, const BinaryOp& func = BinaryOp())
+      : m_lhs(aLhs), m_rhs(aRhs), m_functor(func)
     {
       EIGEN_CHECK_BINARY_COMPATIBILIY(BinaryOp,typename Lhs::Scalar,typename Rhs::Scalar);
       // require the sizes to match
       EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs)
-      eigen_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols());
+      eigen_assert(aLhs.rows() == aRhs.rows() && aLhs.cols() == aRhs.cols());
     }
 
     EIGEN_STRONG_INLINE Index rows() const {
@@ -169,17 +169,17 @@ class CwiseBinaryOpImpl<BinaryOp, Lhs, Rhs, Dense>
     typedef typename internal::dense_xpr_base<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >::type Base;
     EIGEN_DENSE_PUBLIC_INTERFACE( Derived )
 
-    EIGEN_STRONG_INLINE const Scalar coeff(Index row, Index col) const
+    EIGEN_STRONG_INLINE const Scalar coeff(Index rowId, Index colId) const
     {
-      return derived().functor()(derived().lhs().coeff(row, col),
-                                 derived().rhs().coeff(row, col));
+      return derived().functor()(derived().lhs().coeff(rowId, colId),
+                                 derived().rhs().coeff(rowId, colId));
     }
 
     template<int LoadMode>
-    EIGEN_STRONG_INLINE PacketScalar packet(Index row, Index col) const
+    EIGEN_STRONG_INLINE PacketScalar packet(Index rowId, Index colId) const
     {
-      return derived().functor().packetOp(derived().lhs().template packet<LoadMode>(row, col),
-                                          derived().rhs().template packet<LoadMode>(row, col));
+      return derived().functor().packetOp(derived().lhs().template packet<LoadMode>(rowId, colId),
+                                          derived().rhs().template packet<LoadMode>(rowId, colId));
     }
 
     EIGEN_STRONG_INLINE const Scalar coeff(Index index) const
diff --git a/ThirdParty/eigen3/Eigen/src/Core/CwiseNullaryOp.h b/ThirdParty/eigen3/Eigen/src/Core/CwiseNullaryOp.h
index 2635a62b07b..a93bab2d0f9 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/CwiseNullaryOp.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/CwiseNullaryOp.h
@@ -54,27 +54,27 @@ class CwiseNullaryOp : internal::no_assignment_operator,
     typedef typename internal::dense_xpr_base<CwiseNullaryOp>::type Base;
     EIGEN_DENSE_PUBLIC_INTERFACE(CwiseNullaryOp)
 
-    CwiseNullaryOp(Index rows, Index cols, const NullaryOp& func = NullaryOp())
-      : m_rows(rows), m_cols(cols), m_functor(func)
+    CwiseNullaryOp(Index nbRows, Index nbCols, const NullaryOp& func = NullaryOp())
+      : m_rows(nbRows), m_cols(nbCols), m_functor(func)
     {
-      eigen_assert(rows >= 0
-            && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows)
-            &&  cols >= 0
-            && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols));
+      eigen_assert(nbRows >= 0
+            && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == nbRows)
+            &&  nbCols >= 0
+            && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == nbCols));
     }
 
     EIGEN_STRONG_INLINE Index rows() const { return m_rows.value(); }
     EIGEN_STRONG_INLINE Index cols() const { return m_cols.value(); }
 
-    EIGEN_STRONG_INLINE const Scalar coeff(Index rows, Index cols) const
+    EIGEN_STRONG_INLINE const Scalar coeff(Index rowId, Index colId) const
     {
-      return m_functor(rows, cols);
+      return m_functor(rowId, colId);
     }
 
     template<int LoadMode>
-    EIGEN_STRONG_INLINE PacketScalar packet(Index row, Index col) const
+    EIGEN_STRONG_INLINE PacketScalar packet(Index rowId, Index colId) const
     {
-      return m_functor.packetOp(row, col);
+      return m_functor.packetOp(rowId, colId);
     }
 
     EIGEN_STRONG_INLINE const Scalar coeff(Index index) const
@@ -163,11 +163,11 @@ DenseBase<Derived>::NullaryExpr(const CustomNullaryOp& func)
 
 /** \returns an expression of a constant matrix of value \a value
   *
-  * The parameters \a rows and \a cols are the number of rows and of columns of
+  * The parameters \a nbRows and \a nbCols are the number of rows and of columns of
   * the returned matrix. Must be compatible with this DenseBase type.
   *
   * This variant is meant to be used for dynamic-size matrix types. For fixed-size types,
-  * it is redundant to pass \a rows and \a cols as arguments, so Zero() should be used
+  * it is redundant to pass \a nbRows and \a nbCols as arguments, so Zero() should be used
   * instead.
   *
   * The template parameter \a CustomNullaryOp is the type of the functor.
@@ -176,9 +176,9 @@ DenseBase<Derived>::NullaryExpr(const CustomNullaryOp& func)
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE const typename DenseBase<Derived>::ConstantReturnType
-DenseBase<Derived>::Constant(Index rows, Index cols, const Scalar& value)
+DenseBase<Derived>::Constant(Index nbRows, Index nbCols, const Scalar& value)
 {
-  return DenseBase<Derived>::NullaryExpr(rows, cols, internal::scalar_constant_op<Scalar>(value));
+  return DenseBase<Derived>::NullaryExpr(nbRows, nbCols, internal::scalar_constant_op<Scalar>(value));
 }
 
 /** \returns an expression of a constant matrix of value \a value
@@ -292,14 +292,14 @@ DenseBase<Derived>::LinSpaced(const Scalar& low, const Scalar& high)
   return DenseBase<Derived>::NullaryExpr(Derived::SizeAtCompileTime, internal::linspaced_op<Scalar,true>(low,high,Derived::SizeAtCompileTime));
 }
 
-/** \returns true if all coefficients in this matrix are approximately equal to \a value, to within precision \a prec */
+/** \returns true if all coefficients in this matrix are approximately equal to \a val, to within precision \a prec */
 template<typename Derived>
 bool DenseBase<Derived>::isApproxToConstant
-(const Scalar& value, RealScalar prec) const
+(const Scalar& val, const RealScalar& prec) const
 {
   for(Index j = 0; j < cols(); ++j)
     for(Index i = 0; i < rows(); ++i)
-      if(!internal::isApprox(this->coeff(i, j), value, prec))
+      if(!internal::isApprox(this->coeff(i, j), val, prec))
         return false;
   return true;
 }
@@ -309,19 +309,19 @@ bool DenseBase<Derived>::isApproxToConstant
   * \returns true if all coefficients in this matrix are approximately equal to \a value, to within precision \a prec */
 template<typename Derived>
 bool DenseBase<Derived>::isConstant
-(const Scalar& value, RealScalar prec) const
+(const Scalar& val, const RealScalar& prec) const
 {
-  return isApproxToConstant(value, prec);
+  return isApproxToConstant(val, prec);
 }
 
-/** Alias for setConstant(): sets all coefficients in this expression to \a value.
+/** Alias for setConstant(): sets all coefficients in this expression to \a val.
   *
   * \sa setConstant(), Constant(), class CwiseNullaryOp
   */
 template<typename Derived>
-EIGEN_STRONG_INLINE void DenseBase<Derived>::fill(const Scalar& value)
+EIGEN_STRONG_INLINE void DenseBase<Derived>::fill(const Scalar& val)
 {
-  setConstant(value);
+  setConstant(val);
 }
 
 /** Sets all coefficients in this expression to \a value.
@@ -329,9 +329,9 @@ EIGEN_STRONG_INLINE void DenseBase<Derived>::fill(const Scalar& value)
   * \sa fill(), setConstant(Index,const Scalar&), setConstant(Index,Index,const Scalar&), setZero(), setOnes(), Constant(), class CwiseNullaryOp, setZero(), setOnes()
   */
 template<typename Derived>
-EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setConstant(const Scalar& value)
+EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setConstant(const Scalar& val)
 {
-  return derived() = Constant(rows(), cols(), value);
+  return derived() = Constant(rows(), cols(), val);
 }
 
 /** Resizes to the given \a size, and sets all coefficients in this expression to the given \a value.
@@ -345,17 +345,17 @@ EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setConstant(const Scalar& value
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE Derived&
-PlainObjectBase<Derived>::setConstant(Index size, const Scalar& value)
+PlainObjectBase<Derived>::setConstant(Index size, const Scalar& val)
 {
   resize(size);
-  return setConstant(value);
+  return setConstant(val);
 }
 
 /** Resizes to the given size, and sets all coefficients in this expression to the given \a value.
   *
-  * \param rows the new number of rows
-  * \param cols the new number of columns
-  * \param value the value to which all coefficients are set
+  * \param nbRows the new number of rows
+  * \param nbCols the new number of columns
+  * \param val the value to which all coefficients are set
   *
   * Example: \include Matrix_setConstant_int_int.cpp
   * Output: \verbinclude Matrix_setConstant_int_int.out
@@ -364,10 +364,10 @@ PlainObjectBase<Derived>::setConstant(Index size, const Scalar& value)
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE Derived&
-PlainObjectBase<Derived>::setConstant(Index rows, Index cols, const Scalar& value)
+PlainObjectBase<Derived>::setConstant(Index nbRows, Index nbCols, const Scalar& val)
 {
-  resize(rows, cols);
-  return setConstant(value);
+  resize(nbRows, nbCols);
+  return setConstant(val);
 }
 
 /**
@@ -384,10 +384,10 @@ PlainObjectBase<Derived>::setConstant(Index rows, Index cols, const Scalar& valu
   * \sa CwiseNullaryOp
   */
 template<typename Derived>
-EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setLinSpaced(Index size, const Scalar& low, const Scalar& high)
+EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setLinSpaced(Index newSize, const Scalar& low, const Scalar& high)
 {
   EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return derived() = Derived::NullaryExpr(size, internal::linspaced_op<Scalar,false>(low,high,size));
+  return derived() = Derived::NullaryExpr(newSize, internal::linspaced_op<Scalar,false>(low,high,newSize));
 }
 
 /**
@@ -425,9 +425,9 @@ EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setLinSpaced(const Scalar& low,
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE const typename DenseBase<Derived>::ConstantReturnType
-DenseBase<Derived>::Zero(Index rows, Index cols)
+DenseBase<Derived>::Zero(Index nbRows, Index nbCols)
 {
-  return Constant(rows, cols, Scalar(0));
+  return Constant(nbRows, nbCols, Scalar(0));
 }
 
 /** \returns an expression of a zero vector.
@@ -479,7 +479,7 @@ DenseBase<Derived>::Zero()
   * \sa class CwiseNullaryOp, Zero()
   */
 template<typename Derived>
-bool DenseBase<Derived>::isZero(RealScalar prec) const
+bool DenseBase<Derived>::isZero(const RealScalar& prec) const
 {
   for(Index j = 0; j < cols(); ++j)
     for(Index i = 0; i < rows(); ++i)
@@ -512,16 +512,16 @@ EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setZero()
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE Derived&
-PlainObjectBase<Derived>::setZero(Index size)
+PlainObjectBase<Derived>::setZero(Index newSize)
 {
-  resize(size);
+  resize(newSize);
   return setConstant(Scalar(0));
 }
 
 /** Resizes to the given size, and sets all coefficients in this expression to zero.
   *
-  * \param rows the new number of rows
-  * \param cols the new number of columns
+  * \param nbRows the new number of rows
+  * \param nbCols the new number of columns
   *
   * Example: \include Matrix_setZero_int_int.cpp
   * Output: \verbinclude Matrix_setZero_int_int.out
@@ -530,9 +530,9 @@ PlainObjectBase<Derived>::setZero(Index size)
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE Derived&
-PlainObjectBase<Derived>::setZero(Index rows, Index cols)
+PlainObjectBase<Derived>::setZero(Index nbRows, Index nbCols)
 {
-  resize(rows, cols);
+  resize(nbRows, nbCols);
   return setConstant(Scalar(0));
 }
 
@@ -540,7 +540,7 @@ PlainObjectBase<Derived>::setZero(Index rows, Index cols)
 
 /** \returns an expression of a matrix where all coefficients equal one.
   *
-  * The parameters \a rows and \a cols are the number of rows and of columns of
+  * The parameters \a nbRows and \a nbCols are the number of rows and of columns of
   * the returned matrix. Must be compatible with this MatrixBase type.
   *
   * This variant is meant to be used for dynamic-size matrix types. For fixed-size types,
@@ -554,14 +554,14 @@ PlainObjectBase<Derived>::setZero(Index rows, Index cols)
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE const typename DenseBase<Derived>::ConstantReturnType
-DenseBase<Derived>::Ones(Index rows, Index cols)
+DenseBase<Derived>::Ones(Index nbRows, Index nbCols)
 {
-  return Constant(rows, cols, Scalar(1));
+  return Constant(nbRows, nbCols, Scalar(1));
 }
 
 /** \returns an expression of a vector where all coefficients equal one.
   *
-  * The parameter \a size is the size of the returned vector.
+  * The parameter \a newSize is the size of the returned vector.
   * Must be compatible with this MatrixBase type.
   *
   * \only_for_vectors
@@ -577,9 +577,9 @@ DenseBase<Derived>::Ones(Index rows, Index cols)
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE const typename DenseBase<Derived>::ConstantReturnType
-DenseBase<Derived>::Ones(Index size)
+DenseBase<Derived>::Ones(Index newSize)
 {
-  return Constant(size, Scalar(1));
+  return Constant(newSize, Scalar(1));
 }
 
 /** \returns an expression of a fixed-size matrix or vector where all coefficients equal one.
@@ -609,7 +609,7 @@ DenseBase<Derived>::Ones()
   */
 template<typename Derived>
 bool DenseBase<Derived>::isOnes
-(RealScalar prec) const
+(const RealScalar& prec) const
 {
   return isApproxToConstant(Scalar(1), prec);
 }
@@ -627,7 +627,7 @@ EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setOnes()
   return setConstant(Scalar(1));
 }
 
-/** Resizes to the given \a size, and sets all coefficients in this expression to one.
+/** Resizes to the given \a newSize, and sets all coefficients in this expression to one.
   *
   * \only_for_vectors
   *
@@ -638,16 +638,16 @@ EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::setOnes()
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE Derived&
-PlainObjectBase<Derived>::setOnes(Index size)
+PlainObjectBase<Derived>::setOnes(Index newSize)
 {
-  resize(size);
+  resize(newSize);
   return setConstant(Scalar(1));
 }
 
 /** Resizes to the given size, and sets all coefficients in this expression to one.
   *
-  * \param rows the new number of rows
-  * \param cols the new number of columns
+  * \param nbRows the new number of rows
+  * \param nbCols the new number of columns
   *
   * Example: \include Matrix_setOnes_int_int.cpp
   * Output: \verbinclude Matrix_setOnes_int_int.out
@@ -656,9 +656,9 @@ PlainObjectBase<Derived>::setOnes(Index size)
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE Derived&
-PlainObjectBase<Derived>::setOnes(Index rows, Index cols)
+PlainObjectBase<Derived>::setOnes(Index nbRows, Index nbCols)
 {
-  resize(rows, cols);
+  resize(nbRows, nbCols);
   return setConstant(Scalar(1));
 }
 
@@ -666,7 +666,7 @@ PlainObjectBase<Derived>::setOnes(Index rows, Index cols)
 
 /** \returns an expression of the identity matrix (not necessarily square).
   *
-  * The parameters \a rows and \a cols are the number of rows and of columns of
+  * The parameters \a nbRows and \a nbCols are the number of rows and of columns of
   * the returned matrix. Must be compatible with this MatrixBase type.
   *
   * This variant is meant to be used for dynamic-size matrix types. For fixed-size types,
@@ -680,9 +680,9 @@ PlainObjectBase<Derived>::setOnes(Index rows, Index cols)
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE const typename MatrixBase<Derived>::IdentityReturnType
-MatrixBase<Derived>::Identity(Index rows, Index cols)
+MatrixBase<Derived>::Identity(Index nbRows, Index nbCols)
 {
-  return DenseBase<Derived>::NullaryExpr(rows, cols, internal::scalar_identity_op<Scalar>());
+  return DenseBase<Derived>::NullaryExpr(nbRows, nbCols, internal::scalar_identity_op<Scalar>());
 }
 
 /** \returns an expression of the identity matrix (not necessarily square).
@@ -714,7 +714,7 @@ MatrixBase<Derived>::Identity()
   */
 template<typename Derived>
 bool MatrixBase<Derived>::isIdentity
-(RealScalar prec) const
+(const RealScalar& prec) const
 {
   for(Index j = 0; j < cols(); ++j)
   {
@@ -776,8 +776,8 @@ EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::setIdentity()
 
 /** \brief Resizes to the given size, and writes the identity expression (not necessarily square) into *this.
   *
-  * \param rows the new number of rows
-  * \param cols the new number of columns
+  * \param nbRows the new number of rows
+  * \param nbCols the new number of columns
   *
   * Example: \include Matrix_setIdentity_int_int.cpp
   * Output: \verbinclude Matrix_setIdentity_int_int.out
@@ -785,9 +785,9 @@ EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::setIdentity()
   * \sa MatrixBase::setIdentity(), class CwiseNullaryOp, MatrixBase::Identity()
   */
 template<typename Derived>
-EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::setIdentity(Index rows, Index cols)
+EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::setIdentity(Index nbRows, Index nbCols)
 {
-  derived().resize(rows, cols);
+  derived().resize(nbRows, nbCols);
   return setIdentity();
 }
 
@@ -798,10 +798,10 @@ EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::setIdentity(Index rows, Index
   * \sa MatrixBase::Unit(Index), MatrixBase::UnitX(), MatrixBase::UnitY(), MatrixBase::UnitZ(), MatrixBase::UnitW()
   */
 template<typename Derived>
-EIGEN_STRONG_INLINE const typename MatrixBase<Derived>::BasisReturnType MatrixBase<Derived>::Unit(Index size, Index i)
+EIGEN_STRONG_INLINE const typename MatrixBase<Derived>::BasisReturnType MatrixBase<Derived>::Unit(Index newSize, Index i)
 {
   EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
-  return BasisReturnType(SquareMatrixType::Identity(size,size), i);
+  return BasisReturnType(SquareMatrixType::Identity(newSize,newSize), i);
 }
 
 /** \returns an expression of the i-th unit (basis) vector.
diff --git a/ThirdParty/eigen3/Eigen/src/Core/CwiseUnaryOp.h b/ThirdParty/eigen3/Eigen/src/Core/CwiseUnaryOp.h
index 063355ae521..f2de749f92b 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/CwiseUnaryOp.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/CwiseUnaryOp.h
@@ -98,15 +98,15 @@ class CwiseUnaryOpImpl<UnaryOp,XprType,Dense>
     typedef typename internal::dense_xpr_base<CwiseUnaryOp<UnaryOp, XprType> >::type Base;
     EIGEN_DENSE_PUBLIC_INTERFACE(Derived)
 
-    EIGEN_STRONG_INLINE const Scalar coeff(Index row, Index col) const
+    EIGEN_STRONG_INLINE const Scalar coeff(Index rowId, Index colId) const
     {
-      return derived().functor()(derived().nestedExpression().coeff(row, col));
+      return derived().functor()(derived().nestedExpression().coeff(rowId, colId));
     }
 
     template<int LoadMode>
-    EIGEN_STRONG_INLINE PacketScalar packet(Index row, Index col) const
+    EIGEN_STRONG_INLINE PacketScalar packet(Index rowId, Index colId) const
     {
-      return derived().functor().packetOp(derived().nestedExpression().template packet<LoadMode>(row, col));
+      return derived().functor().packetOp(derived().nestedExpression().template packet<LoadMode>(rowId, colId));
     }
 
     EIGEN_STRONG_INLINE const Scalar coeff(Index index) const
diff --git a/ThirdParty/eigen3/Eigen/src/Core/CwiseUnaryView.h b/ThirdParty/eigen3/Eigen/src/Core/CwiseUnaryView.h
index 57fbd7247b7..b2638d3265b 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/CwiseUnaryView.h
+++ b/ThirdParty/eigen3/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
     {
diff --git a/ThirdParty/eigen3/Eigen/src/Core/DenseBase.h b/ThirdParty/eigen3/Eigen/src/Core/DenseBase.h
index 1cc0314ef0b..c5800f6c8c8 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/DenseBase.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/DenseBase.h
@@ -13,6 +13,16 @@
 
 namespace Eigen {
 
+namespace internal {
+  
+// The index type defined by EIGEN_DEFAULT_DENSE_INDEX_TYPE must be a signed type.
+// This dummy function simply aims at checking that at compile time.
+static inline void check_DenseIndex_is_signed() {
+  EIGEN_STATIC_ASSERT(NumTraits<DenseIndex>::IsSigned,THE_INDEX_TYPE_MUST_BE_A_SIGNED_TYPE); 
+}
+
+} // end namespace internal
+  
 /** \class DenseBase
   * \ingroup Core_Module
   *
@@ -204,21 +214,21 @@ template<typename Derived> class DenseBase
       * Matrix::resize() and Array::resize(). The present method only asserts that the new size equals the old size, and does
       * nothing else.
       */
-    void resize(Index size)
+    void resize(Index newSize)
     {
-      EIGEN_ONLY_USED_FOR_DEBUG(size);
-      eigen_assert(size == this->size()
+      EIGEN_ONLY_USED_FOR_DEBUG(newSize);
+      eigen_assert(newSize == this->size()
                 && "DenseBase::resize() does not actually allow to resize.");
     }
     /** Only plain matrices/arrays, not expressions, may be resized; therefore the only useful resize methods are
       * Matrix::resize() and Array::resize(). The present method only asserts that the new size equals the old size, and does
       * nothing else.
       */
-    void resize(Index rows, Index cols)
+    void resize(Index nbRows, Index nbCols)
     {
-      EIGEN_ONLY_USED_FOR_DEBUG(rows);
-      EIGEN_ONLY_USED_FOR_DEBUG(cols);
-      eigen_assert(rows == this->rows() && cols == this->cols()
+      EIGEN_ONLY_USED_FOR_DEBUG(nbRows);
+      EIGEN_ONLY_USED_FOR_DEBUG(nbCols);
+      eigen_assert(nbRows == this->rows() && nbCols == this->cols()
                 && "DenseBase::resize() does not actually allow to resize.");
     }
 
@@ -271,7 +281,7 @@ template<typename Derived> class DenseBase
     CommaInitializer<Derived> operator<< (const DenseBase<OtherDerived>& other);
 
     Eigen::Transpose<Derived> transpose();
-    typedef const Transpose<const Derived> ConstTransposeReturnType;
+	typedef typename internal::add_const<Transpose<const Derived> >::type ConstTransposeReturnType;
     ConstTransposeReturnType transpose() const;
     void transposeInPlace();
 #ifndef EIGEN_NO_DEBUG
@@ -281,29 +291,6 @@ template<typename Derived> class DenseBase
   public:
 #endif
 
-    typedef VectorBlock<Derived> SegmentReturnType;
-    typedef const VectorBlock<const Derived> ConstSegmentReturnType;
-    template<int Size> struct FixedSegmentReturnType { typedef VectorBlock<Derived, Size> Type; };
-    template<int Size> struct ConstFixedSegmentReturnType { typedef const VectorBlock<const Derived, Size> Type; };
-    
-    // Note: The "DenseBase::" prefixes are added to help MSVC9 to match these declarations with the later implementations.
-    SegmentReturnType segment(Index start, Index size);
-    typename DenseBase::ConstSegmentReturnType segment(Index start, Index size) const;
-
-    SegmentReturnType head(Index size);
-    typename DenseBase::ConstSegmentReturnType head(Index size) const;
-
-    SegmentReturnType tail(Index size);
-    typename DenseBase::ConstSegmentReturnType tail(Index size) const;
-
-    template<int Size> typename FixedSegmentReturnType<Size>::Type head();
-    template<int Size> typename ConstFixedSegmentReturnType<Size>::Type head() const;
-
-    template<int Size> typename FixedSegmentReturnType<Size>::Type tail();
-    template<int Size> typename ConstFixedSegmentReturnType<Size>::Type tail() const;
-
-    template<int Size> typename FixedSegmentReturnType<Size>::Type segment(Index start);
-    template<int Size> typename ConstFixedSegmentReturnType<Size>::Type segment(Index start) const;
 
     static const ConstantReturnType
     Constant(Index rows, Index cols, const Scalar& value);
@@ -348,17 +335,20 @@ template<typename Derived> class DenseBase
 
     template<typename OtherDerived>
     bool isApprox(const DenseBase<OtherDerived>& other,
-                  RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+                  const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
     bool isMuchSmallerThan(const RealScalar& other,
-                           RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+                           const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
     template<typename OtherDerived>
     bool isMuchSmallerThan(const DenseBase<OtherDerived>& other,
-                           RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+                           const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
 
-    bool isApproxToConstant(const Scalar& value, RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
-    bool isConstant(const Scalar& value, RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
-    bool isZero(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
-    bool isOnes(RealScalar prec = NumTraits<Scalar>::dummy_precision()) const;
+    bool isApproxToConstant(const Scalar& value, const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
+    bool isConstant(const Scalar& value, const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
+    bool isZero(const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
+    bool isOnes(const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const;
+    
+    inline bool hasNaN() const;
+    inline bool allFinite() const;
 
     inline Derived& operator*=(const Scalar& other);
     inline Derived& operator/=(const Scalar& other);
@@ -438,8 +428,6 @@ template<typename Derived> class DenseBase
       return derived().coeff(0,0);
     }
 
-/////////// Array module ///////////
-
     bool all(void) const;
     bool any(void) const;
     Index count() const;
@@ -465,11 +453,11 @@ template<typename Derived> class DenseBase
 
     template<typename ThenDerived>
     inline const Select<Derived,ThenDerived, typename ThenDerived::ConstantReturnType>
-    select(const DenseBase<ThenDerived>& thenMatrix, typename ThenDerived::Scalar elseScalar) const;
+    select(const DenseBase<ThenDerived>& thenMatrix, const typename ThenDerived::Scalar& elseScalar) const;
 
     template<typename ElseDerived>
     inline const Select<Derived, typename ElseDerived::ConstantReturnType, ElseDerived >
-    select(typename ElseDerived::Scalar thenScalar, const DenseBase<ElseDerived>& elseMatrix) const;
+    select(const typename ElseDerived::Scalar& thenScalar, const DenseBase<ElseDerived>& elseMatrix) const;
 
     template<int p> RealScalar lpNorm() const;
 
diff --git a/ThirdParty/eigen3/Eigen/src/Core/arch/NEON/Complex.h b/ThirdParty/eigen3/Eigen/src/Core/arch/NEON/Complex.h
index 795b4be7303..f183d31de2a 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/arch/NEON/Complex.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/arch/NEON/Complex.h
@@ -68,7 +68,6 @@ template<> EIGEN_STRONG_INLINE Packet2cf pconj(const Packet2cf& a)
 template<> EIGEN_STRONG_INLINE Packet2cf pmul<Packet2cf>(const Packet2cf& a, const Packet2cf& b)
 {
   Packet4f v1, v2;
-  float32x2_t a_lo, a_hi;
 
   // Get the real values of a | a1_re | a1_re | a2_re | a2_re |
   v1 = vcombine_f32(vdup_lane_f32(vget_low_f32(a.v), 0), vdup_lane_f32(vget_high_f32(a.v), 0));
@@ -81,9 +80,7 @@ template<> EIGEN_STRONG_INLINE Packet2cf pmul<Packet2cf>(const Packet2cf& a, con
   // Conjugate v2 
   v2 = vreinterpretq_f32_u32(veorq_u32(vreinterpretq_u32_f32(v2), p4ui_CONJ_XOR));
   // Swap real/imag elements in v2.
-  a_lo = vrev64_f32(vget_low_f32(v2));
-  a_hi = vrev64_f32(vget_high_f32(v2));
-  v2 = vcombine_f32(a_lo, a_hi);
+  v2 = vrev64q_f32(v2);
   // Add and return the result
   return Packet2cf(vaddq_f32(v1, v2));
 }
@@ -241,13 +238,10 @@ template<> EIGEN_STRONG_INLINE Packet2cf pdiv<Packet2cf>(const Packet2cf& a, con
   // TODO optimize it for AltiVec
   Packet2cf res = conj_helper<Packet2cf,Packet2cf,false,true>().pmul(a,b);
   Packet4f s, rev_s;
-  float32x2_t a_lo, a_hi;
 
   // this computes the norm
   s = vmulq_f32(b.v, b.v);
-  a_lo = vrev64_f32(vget_low_f32(s));
-  a_hi = vrev64_f32(vget_high_f32(s));
-  rev_s = vcombine_f32(a_lo, a_hi);
+  rev_s = vrev64q_f32(s);
 
   return Packet2cf(pdiv(res.v, vaddq_f32(s,rev_s)));
 }
diff --git a/ThirdParty/eigen3/Eigen/src/Core/arch/SSE/Complex.h b/ThirdParty/eigen3/Eigen/src/Core/arch/SSE/Complex.h
index 12df987754c..91bba5e38fe 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/arch/SSE/Complex.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/arch/SSE/Complex.h
@@ -81,25 +81,31 @@ template<> EIGEN_STRONG_INLINE Packet2cf por    <Packet2cf>(const Packet2cf& a,
 template<> EIGEN_STRONG_INLINE Packet2cf pxor   <Packet2cf>(const Packet2cf& a, const Packet2cf& b) { return Packet2cf(_mm_xor_ps(a.v,b.v)); }
 template<> EIGEN_STRONG_INLINE Packet2cf pandnot<Packet2cf>(const Packet2cf& a, const Packet2cf& b) { return Packet2cf(_mm_andnot_ps(a.v,b.v)); }
 
-template<> EIGEN_STRONG_INLINE Packet2cf pload <Packet2cf>(const std::complex<float>* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet2cf(pload<Packet4f>(&real_ref(*from))); }
-template<> EIGEN_STRONG_INLINE Packet2cf ploadu<Packet2cf>(const std::complex<float>* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet2cf(ploadu<Packet4f>(&real_ref(*from))); }
+template<> EIGEN_STRONG_INLINE Packet2cf pload <Packet2cf>(const std::complex<float>* from) { EIGEN_DEBUG_ALIGNED_LOAD return Packet2cf(pload<Packet4f>(&numext::real_ref(*from))); }
+template<> EIGEN_STRONG_INLINE Packet2cf ploadu<Packet2cf>(const std::complex<float>* from) { EIGEN_DEBUG_UNALIGNED_LOAD return Packet2cf(ploadu<Packet4f>(&numext::real_ref(*from))); }
 
 template<> EIGEN_STRONG_INLINE Packet2cf pset1<Packet2cf>(const std::complex<float>&  from)
 {
   Packet2cf res;
-  #if EIGEN_GNUC_AT_MOST(4,2)
-  // workaround annoying "may be used uninitialized in this function" warning with gcc 4.2
+#if EIGEN_GNUC_AT_MOST(4,2)
+  // Workaround annoying "may be used uninitialized in this function" warning with gcc 4.2
   res.v = _mm_loadl_pi(_mm_set1_ps(0.0f), reinterpret_cast<const __m64*>(&from));
-  #else
+#elif EIGEN_GNUC_AT_LEAST(4,6)
+  // Suppress annoying "may be used uninitialized in this function" warning with gcc >= 4.6
+  #pragma GCC diagnostic push
+  #pragma GCC diagnostic ignored "-Wuninitialized"
   res.v = _mm_loadl_pi(res.v, (const __m64*)&from);
-  #endif
+  #pragma GCC diagnostic pop
+#else
+  res.v = _mm_loadl_pi(res.v, (const __m64*)&from);
+#endif
   return Packet2cf(_mm_movelh_ps(res.v,res.v));
 }
 
 template<> EIGEN_STRONG_INLINE Packet2cf ploaddup<Packet2cf>(const std::complex<float>* from) { return pset1<Packet2cf>(*from); }
 
-template<> EIGEN_STRONG_INLINE void pstore <std::complex<float> >(std::complex<float> *   to, const Packet2cf& from) { EIGEN_DEBUG_ALIGNED_STORE pstore(&real_ref(*to), from.v); }
-template<> EIGEN_STRONG_INLINE void pstoreu<std::complex<float> >(std::complex<float> *   to, const Packet2cf& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu(&real_ref(*to), from.v); }
+template<> EIGEN_STRONG_INLINE void pstore <std::complex<float> >(std::complex<float> *   to, const Packet2cf& from) { EIGEN_DEBUG_ALIGNED_STORE pstore(&numext::real_ref(*to), from.v); }
+template<> EIGEN_STRONG_INLINE void pstoreu<std::complex<float> >(std::complex<float> *   to, const Packet2cf& from) { EIGEN_DEBUG_UNALIGNED_STORE pstoreu(&numext::real_ref(*to), from.v); }
 
 template<> EIGEN_STRONG_INLINE void prefetch<std::complex<float> >(const std::complex<float> *   addr) { _mm_prefetch((const char*)(addr), _MM_HINT_T0); }
 
diff --git a/ThirdParty/eigen3/Eigen/src/Core/products/CoeffBasedProduct.h b/ThirdParty/eigen3/Eigen/src/Core/products/CoeffBasedProduct.h
index 403d25fa9eb..c06a0df1c21 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/products/CoeffBasedProduct.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/products/CoeffBasedProduct.h
@@ -150,7 +150,7 @@ class CoeffBasedProduct
     {
       // we don't allow taking products of matrices of different real types, as that wouldn't be vectorizable.
       // We still allow to mix T and complex<T>.
-      EIGEN_STATIC_ASSERT((internal::is_same<typename Lhs::RealScalar, typename Rhs::RealScalar>::value),
+      EIGEN_STATIC_ASSERT((internal::scalar_product_traits<typename Lhs::RealScalar, typename Rhs::RealScalar>::Defined),
         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
       eigen_assert(lhs.cols() == rhs.rows()
         && "invalid matrix product"
diff --git a/ThirdParty/eigen3/Eigen/src/Core/util/BlasUtil.h b/ThirdParty/eigen3/Eigen/src/Core/util/BlasUtil.h
index 91496651c82..a28f16fa04b 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/util/BlasUtil.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/util/BlasUtil.h
@@ -42,7 +42,7 @@ template<bool Conjugate> struct conj_if;
 
 template<> struct conj_if<true> {
   template<typename T>
-  inline T operator()(const T& x) { return conj(x); }
+  inline T operator()(const T& x) { return numext::conj(x); }
   template<typename T>
   inline T pconj(const T& x) { return internal::pconj(x); }
 };
@@ -67,7 +67,7 @@ template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::
   { return c + pmul(x,y); }
 
   EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
-  { return Scalar(real(x)*real(y) + imag(x)*imag(y), imag(x)*real(y) - real(x)*imag(y)); }
+  { return Scalar(numext::real(x)*numext::real(y) + numext::imag(x)*numext::imag(y), numext::imag(x)*numext::real(y) - numext::real(x)*numext::imag(y)); }
 };
 
 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, true,false>
@@ -77,7 +77,7 @@ template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::
   { return c + pmul(x,y); }
 
   EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
-  { return Scalar(real(x)*real(y) + imag(x)*imag(y), real(x)*imag(y) - imag(x)*real(y)); }
+  { return Scalar(numext::real(x)*numext::real(y) + numext::imag(x)*numext::imag(y), numext::real(x)*numext::imag(y) - numext::imag(x)*numext::real(y)); }
 };
 
 template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::complex<RealScalar>, true,true>
@@ -87,7 +87,7 @@ template<typename RealScalar> struct conj_helper<std::complex<RealScalar>, std::
   { return c + pmul(x,y); }
 
   EIGEN_STRONG_INLINE Scalar pmul(const Scalar& x, const Scalar& y) const
-  { return Scalar(real(x)*real(y) - imag(x)*imag(y), - real(x)*imag(y) - imag(x)*real(y)); }
+  { return Scalar(numext::real(x)*numext::real(y) - numext::imag(x)*numext::imag(y), - numext::real(x)*numext::imag(y) - numext::imag(x)*numext::real(y)); }
 };
 
 template<typename RealScalar,bool Conj> struct conj_helper<std::complex<RealScalar>, RealScalar, Conj,false>
@@ -113,7 +113,7 @@ template<typename From,typename To> struct get_factor {
 };
 
 template<typename Scalar> struct get_factor<Scalar,typename NumTraits<Scalar>::Real> {
-  static EIGEN_STRONG_INLINE typename NumTraits<Scalar>::Real run(const Scalar& x) { return real(x); }
+  static EIGEN_STRONG_INLINE typename NumTraits<Scalar>::Real run(const Scalar& x) { return numext::real(x); }
 };
 
 // Lightweight helper class to access matrix coefficients.
diff --git a/ThirdParty/eigen3/Eigen/src/Core/util/Constants.h b/ThirdParty/eigen3/Eigen/src/Core/util/Constants.h
index 3fd45e84f8e..14b9624e1d9 100644
--- a/ThirdParty/eigen3/Eigen/src/Core/util/Constants.h
+++ b/ThirdParty/eigen3/Eigen/src/Core/util/Constants.h
@@ -13,13 +13,18 @@
 
 namespace Eigen {
 
-/** This value means that a quantity is not known at compile-time, and that instead the value is
+/** This value means that a positive quantity (e.g., a size) is not known at compile-time, and that instead the value is
   * stored in some runtime variable.
   *
   * Changing the value of Dynamic breaks the ABI, as Dynamic is often used as a template parameter for Matrix.
   */
 const int Dynamic = -1;
 
+/** This value means that a signed quantity (e.g., a signed index) is not known at compile-time, and that instead its value
+  * has to be specified at runtime.
+  */
+const int DynamicIndex = 0xffffff;
+
 /** This value means +Infinity; it is currently used only as the p parameter to MatrixBase::lpNorm<int>().
   * The value Infinity there means the L-infinity norm.
   */
@@ -227,7 +232,9 @@ enum {
     * scalar loops to handle the unaligned boundaries */
   SliceVectorizedTraversal,
   /** \internal Special case to properly handle incompatible scalar types or other defecting cases*/
-  InvalidTraversal
+  InvalidTraversal,
+  /** \internal Evaluate all entries at once */
+  AllAtOnceTraversal
 };
 
 /** \internal \ingroup enums
@@ -257,9 +264,9 @@ enum {
   ColMajor = 0,
   /** Storage order is row major (see \ref TopicStorageOrders). */
   RowMajor = 0x1,  // it is only a coincidence that this is equal to RowMajorBit -- don't rely on that
-  /** \internal Align the matrix itself if it is vectorizable fixed-size */
+  /** Align the matrix itself if it is vectorizable fixed-size */
   AutoAlign = 0,
-  /** \internal Don't require alignment for the matrix itself (the array of coefficients, if dynamically allocated, may still be requested to be aligned) */ // FIXME --- clarify the situation
+  /** Don't require alignment for the matrix itself (the array of coefficients, if dynamically allocated, may still be requested to be aligned) */ // FIXME --- clarify the situation
   DontAlign = 0x2
 };
 
diff --git a/ThirdParty/eigen3/Eigen/src/Eigen2Support/Geometry/AlignedBox.h b/ThirdParty/eigen3/Eigen/src/Eigen2Support/Geometry/AlignedBox.h
index 5c928e8fc2d..2e4309dd945 100644
--- a/ThirdParty/eigen3/Eigen/src/Eigen2Support/Geometry/AlignedBox.h
+++ b/ThirdParty/eigen3/Eigen/src/Eigen2Support/Geometry/AlignedBox.h
@@ -1,5 +1,5 @@
 // This file is part of Eigen, a lightweight C++ template library
-// for linear algebra. Eigen itself is part of the KDE project.
+// for linear algebra.
 //
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
 //
@@ -34,7 +34,7 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim==
   typedef Matrix<Scalar,AmbientDimAtCompileTime,1> VectorType;
 
   /** Default constructor initializing a null box. */
-  inline explicit AlignedBox()
+  inline AlignedBox()
   { if (AmbientDimAtCompileTime!=Dynamic) setNull(); }
 
   /** Constructs a null box with \a _dim the dimension of the ambient space. */
diff --git a/ThirdParty/eigen3/Eigen/src/Eigen2Support/Geometry/AngleAxis.h b/ThirdParty/eigen3/Eigen/src/Eigen2Support/Geometry/AngleAxis.h
index 20f1fceeb19..af598a40313 100644
--- a/ThirdParty/eigen3/Eigen/src/Eigen2Support/Geometry/AngleAxis.h
+++ b/ThirdParty/eigen3/Eigen/src/Eigen2Support/Geometry/AngleAxis.h
@@ -1,5 +1,5 @@
 // This file is part of Eigen, a lightweight C++ template library
-// for linear algebra. Eigen itself is part of the KDE project.
+// for linear algebra.
 //
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
 //
diff --git a/ThirdParty/eigen3/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/ThirdParty/eigen3/Eigen/src/Eigenvalues/ComplexEigenSolver.h
index c4b8a308cee..af434bc9bd6 100644
--- a/ThirdParty/eigen3/Eigen/src/Eigenvalues/ComplexEigenSolver.h
+++ b/ThirdParty/eigen3/Eigen/src/Eigenvalues/ComplexEigenSolver.h
@@ -3,7 +3,7 @@
 //
 // Copyright (C) 2009 Claire Maurice
 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -220,6 +220,19 @@ template<typename _MatrixType> class ComplexEigenSolver
       return m_schur.info();
     }
 
+    /** \brief Sets the maximum number of iterations allowed. */
+    ComplexEigenSolver& setMaxIterations(Index maxIters)
+    {
+      m_schur.setMaxIterations(maxIters);
+      return *this;
+    }
+
+    /** \brief Returns the maximum number of iterations. */
+    Index getMaxIterations()
+    {
+      return m_schur.getMaxIterations();
+    }
+
   protected:
     EigenvectorType m_eivec;
     EigenvalueType m_eivalues;
@@ -229,16 +242,17 @@ template<typename _MatrixType> class ComplexEigenSolver
     EigenvectorType m_matX;
 
   private:
-    void doComputeEigenvectors(RealScalar matrixnorm);
+    void doComputeEigenvectors(const RealScalar& matrixnorm);
     void sortEigenvalues(bool computeEigenvectors);
 };
 
 
 template<typename MatrixType>
-ComplexEigenSolver<MatrixType>& ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
+ComplexEigenSolver<MatrixType>& 
+ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
 {
   // this code is inspired from Jampack
-  assert(matrix.cols() == matrix.rows());
+  eigen_assert(matrix.cols() == matrix.rows());
 
   // Do a complex Schur decomposition, A = U T U^*
   // The eigenvalues are on the diagonal of T.
@@ -259,7 +273,7 @@ ComplexEigenSolver<MatrixType>& ComplexEigenSolver<MatrixType>::compute(const Ma
 
 
 template<typename MatrixType>
-void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm)
+void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(const RealScalar& matrixnorm)
 {
   const Index n = m_eivalues.size();
 
@@ -280,7 +294,7 @@ void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm
       {
         // If the i-th and k-th eigenvalue are equal, then z equals 0.
         // Use a small value instead, to prevent division by zero.
-        internal::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
+        numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
       }
       m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
     }
diff --git a/ThirdParty/eigen3/Eigen/src/Eigenvalues/ComplexSchur.h b/ThirdParty/eigen3/Eigen/src/Eigenvalues/ComplexSchur.h
index 55aeedb90e8..89e6cade334 100644
--- a/ThirdParty/eigen3/Eigen/src/Eigenvalues/ComplexSchur.h
+++ b/ThirdParty/eigen3/Eigen/src/Eigenvalues/ComplexSchur.h
@@ -3,7 +3,7 @@
 //
 // Copyright (C) 2009 Claire Maurice
 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
-// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -96,7 +96,8 @@ template<typename _MatrixType> class ComplexSchur
         m_matU(size,size),
         m_hess(size),
         m_isInitialized(false),
-        m_matUisUptodate(false)
+        m_matUisUptodate(false),
+        m_maxIters(-1)
     {}
 
     /** \brief Constructor; computes Schur decomposition of given matrix. 
@@ -109,11 +110,12 @@ template<typename _MatrixType> class ComplexSchur
       * \sa matrixT() and matrixU() for examples.
       */
     ComplexSchur(const MatrixType& matrix, bool computeU = true)
-            : m_matT(matrix.rows(),matrix.cols()),
-              m_matU(matrix.rows(),matrix.cols()),
-              m_hess(matrix.rows()),
-              m_isInitialized(false),
-              m_matUisUptodate(false)
+      : m_matT(matrix.rows(),matrix.cols()),
+        m_matU(matrix.rows(),matrix.cols()),
+        m_hess(matrix.rows()),
+        m_isInitialized(false),
+        m_matUisUptodate(false),
+        m_maxIters(-1)
     {
       compute(matrix, computeU);
     }
@@ -166,6 +168,7 @@ template<typename _MatrixType> class ComplexSchur
       * 
       * \param[in]  matrix  Square matrix whose Schur decomposition is to be computed.
       * \param[in]  computeU  If true, both T and U are computed; if false, only T is computed.
+
       * \returns    Reference to \c *this
       *
       * The Schur decomposition is computed by first reducing the
@@ -180,8 +183,30 @@ template<typename _MatrixType> class ComplexSchur
       *
       * Example: \include ComplexSchur_compute.cpp
       * Output: \verbinclude ComplexSchur_compute.out
+      *
+      * \sa compute(const MatrixType&, bool, Index)
       */
     ComplexSchur& compute(const MatrixType& matrix, bool computeU = true);
+    
+    /** \brief Compute Schur decomposition from a given Hessenberg matrix
+     *  \param[in] matrixH Matrix in Hessenberg form H
+     *  \param[in] matrixQ orthogonal matrix Q that transform a matrix A to H : A = Q H Q^T
+     *  \param computeU Computes the matriX U of the Schur vectors
+     * \return Reference to \c *this
+     * 
+     *  This routine assumes that the matrix is already reduced in Hessenberg form matrixH
+     *  using either the class HessenbergDecomposition or another mean. 
+     *  It computes the upper quasi-triangular matrix T of the Schur decomposition of H
+     *  When computeU is true, this routine computes the matrix U such that 
+     *  A = U T U^T =  (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix
+     * 
+     * NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix
+     * is not available, the user should give an identity matrix (Q.setIdentity())
+     * 
+     * \sa compute(const MatrixType&, bool)
+     */
+    template<typename HessMatrixType, typename OrthMatrixType>
+    ComplexSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ,  bool computeU=true);
 
     /** \brief Reports whether previous computation was successful.
       *
@@ -189,15 +214,33 @@ template<typename _MatrixType> class ComplexSchur
       */
     ComputationInfo info() const
     {
-      eigen_assert(m_isInitialized && "RealSchur is not initialized.");
+      eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
       return m_info;
     }
 
-    /** \brief Maximum number of iterations.
+    /** \brief Sets the maximum number of iterations allowed. 
       *
-      * Maximum number of iterations allowed for an eigenvalue to converge. 
+      * If not specified by the user, the maximum number of iterations is m_maxIterationsPerRow times the size
+      * of the matrix.
       */
-    static const int m_maxIterations = 30;
+    ComplexSchur& setMaxIterations(Index maxIters)
+    {
+      m_maxIters = maxIters;
+      return *this;
+    }
+
+    /** \brief Returns the maximum number of iterations. */
+    Index getMaxIterations()
+    {
+      return m_maxIters;
+    }
+
+    /** \brief Maximum number of iterations per row.
+      *
+      * If not otherwise specified, the maximum number of iterations is this number times the size of the
+      * matrix. It is currently set to 30.
+      */
+    static const int m_maxIterationsPerRow = 30;
 
   protected:
     ComplexMatrixType m_matT, m_matU;
@@ -205,6 +248,7 @@ template<typename _MatrixType> class ComplexSchur
     ComputationInfo m_info;
     bool m_isInitialized;
     bool m_matUisUptodate;
+    Index m_maxIters;
 
   private:  
     bool subdiagonalEntryIsNeglegible(Index i);
@@ -219,8 +263,8 @@ template<typename _MatrixType> class ComplexSchur
 template<typename MatrixType>
 inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
 {
-  RealScalar d = internal::norm1(m_matT.coeff(i,i)) + internal::norm1(m_matT.coeff(i+1,i+1));
-  RealScalar sd = internal::norm1(m_matT.coeff(i+1,i));
+  RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1));
+  RealScalar sd = numext::norm1(m_matT.coeff(i+1,i));
   if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
   {
     m_matT.coeffRef(i+1,i) = ComplexScalar(0);
@@ -234,10 +278,11 @@ inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
 template<typename MatrixType>
 typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
 {
+  using std::abs;
   if (iter == 10 || iter == 20) 
   {
     // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
-    return internal::abs(internal::real(m_matT.coeff(iu,iu-1))) + internal::abs(internal::real(m_matT.coeff(iu-1,iu-2)));
+    return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2)));
   }
 
   // compute the shift as one of the eigenvalues of t, the 2x2
@@ -254,13 +299,13 @@ typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::compu
   ComplexScalar eival1 = (trace + disc) / RealScalar(2);
   ComplexScalar eival2 = (trace - disc) / RealScalar(2);
 
-  if(internal::norm1(eival1) > internal::norm1(eival2))
+  if(numext::norm1(eival1) > numext::norm1(eival2))
     eival2 = det / eival1;
   else
     eival1 = det / eival2;
 
   // choose the eigenvalue closest to the bottom entry of the diagonal
-  if(internal::norm1(eival1-t.coeff(1,1)) < internal::norm1(eival2-t.coeff(1,1)))
+  if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1)))
     return normt * eival1;
   else
     return normt * eival2;
@@ -284,10 +329,20 @@ ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& ma
   }
 
   internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU);
-  reduceToTriangularForm(computeU);
+  computeFromHessenberg(m_matT, m_matU, computeU);
   return *this;
 }
 
+template<typename MatrixType>
+template<typename HessMatrixType, typename OrthMatrixType>
+ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
+{
+  m_matT = matrixH;
+  if(computeU)
+    m_matU = matrixQ;
+  reduceToTriangularForm(computeU);
+  return *this;
+}
 namespace internal {
 
 /* Reduce given matrix to Hessenberg form */
@@ -309,7 +364,6 @@ struct complex_schur_reduce_to_hessenberg<MatrixType, false>
   static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
   {
     typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
-    typedef typename ComplexSchur<MatrixType>::ComplexMatrixType ComplexMatrixType;
 
     // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
     _this.m_hess.compute(matrix);
@@ -329,6 +383,10 @@ struct complex_schur_reduce_to_hessenberg<MatrixType, false>
 template<typename MatrixType>
 void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
 {  
+  Index maxIters = m_maxIters;
+  if (maxIters == -1)
+    maxIters = m_maxIterationsPerRow * m_matT.rows();
+
   // The matrix m_matT is divided in three parts. 
   // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 
   // Rows il,...,iu is the part we are working on (the active submatrix).
@@ -354,7 +412,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
     // if we spent too many iterations, we give up
     iter++;
     totalIter++;
-    if(totalIter > m_maxIterations * m_matT.cols()) break;
+    if(totalIter > maxIters) break;
 
     // find il, the top row of the active submatrix
     il = iu-1;
@@ -384,7 +442,7 @@ void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
     }
   }
 
-  if(totalIter <= m_maxIterations * m_matT.cols())
+  if(totalIter <= maxIters)
     m_info = Success;
   else
     m_info = NoConvergence;
diff --git a/ThirdParty/eigen3/Eigen/src/Eigenvalues/ComplexSchur_MKL.h b/ThirdParty/eigen3/Eigen/src/Eigenvalues/ComplexSchur_MKL.h
index ada7a24e34f..91496ae5bdb 100644
--- a/ThirdParty/eigen3/Eigen/src/Eigenvalues/ComplexSchur_MKL.h
+++ b/ThirdParty/eigen3/Eigen/src/Eigenvalues/ComplexSchur_MKL.h
@@ -49,7 +49,7 @@ ComplexSchur<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW> >::compute(const Matri
   typedef MatrixType::RealScalar RealScalar; \
   typedef std::complex<RealScalar> ComplexScalar; \
 \
-  assert(matrix.cols() == matrix.rows()); \
+  eigen_assert(matrix.cols() == matrix.rows()); \
 \
   m_matUisUptodate = false; \
   if(matrix.cols() == 1) \
diff --git a/ThirdParty/eigen3/Eigen/src/Geometry/AlignedBox.h b/ThirdParty/eigen3/Eigen/src/Geometry/AlignedBox.h
index c0f97300c43..8e186d57a34 100644
--- a/ThirdParty/eigen3/Eigen/src/Geometry/AlignedBox.h
+++ b/ThirdParty/eigen3/Eigen/src/Geometry/AlignedBox.h
@@ -56,7 +56,7 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim)
 
 
   /** Default constructor initializing a null box. */
-  inline explicit AlignedBox()
+  inline AlignedBox()
   { if (AmbientDimAtCompileTime!=Dynamic) setEmpty(); }
 
   /** Constructs a null box with \a _dim the dimension of the ambient space. */
@@ -71,7 +71,7 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim)
   template<typename Derived>
   inline explicit AlignedBox(const MatrixBase<Derived>& a_p)
   {
-    const typename internal::nested<Derived,2>::type p(a_p.derived());
+    typename internal::nested<Derived,2>::type p(a_p.derived());
     m_min = p;
     m_max = p;
   }
@@ -248,14 +248,14 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim)
     */
   template<typename Derived>
   inline NonInteger exteriorDistance(const MatrixBase<Derived>& p) const
-  { return internal::sqrt(NonInteger(squaredExteriorDistance(p))); }
+  { using std::sqrt; return sqrt(NonInteger(squaredExteriorDistance(p))); }
 
   /** \returns the distance between the boxes \a b and \c *this,
     * and zero if the boxes intersect.
     * \sa squaredExteriorDistance()
     */
   inline NonInteger exteriorDistance(const AlignedBox& b) const
-  { return internal::sqrt(NonInteger(squaredExteriorDistance(b))); }
+  { using std::sqrt; return sqrt(NonInteger(squaredExteriorDistance(b))); }
 
   /** \returns \c *this with scalar type casted to \a NewScalarType
     *
@@ -282,7 +282,7 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim)
     * determined by \a prec.
     *
     * \sa MatrixBase::isApprox() */
-  bool isApprox(const AlignedBox& other, RealScalar prec = ScalarTraits::dummy_precision()) const
+  bool isApprox(const AlignedBox& other, const RealScalar& prec = ScalarTraits::dummy_precision()) const
   { return m_min.isApprox(other.m_min, prec) && m_max.isApprox(other.m_max, prec); }
 
 protected:
@@ -296,7 +296,7 @@ template<typename Scalar,int AmbientDim>
 template<typename Derived>
 inline Scalar AlignedBox<Scalar,AmbientDim>::squaredExteriorDistance(const MatrixBase<Derived>& a_p) const
 {
-  const typename internal::nested<Derived,2*AmbientDim>::type p(a_p.derived());
+  typename internal::nested<Derived,2*AmbientDim>::type p(a_p.derived());
   Scalar dist2(0);
   Scalar aux;
   for (Index k=0; k<dim(); ++k)
diff --git a/ThirdParty/eigen3/Eigen/src/Geometry/AngleAxis.h b/ThirdParty/eigen3/Eigen/src/Geometry/AngleAxis.h
index 67197ac78c3..553d38c7449 100644
--- a/ThirdParty/eigen3/Eigen/src/Geometry/AngleAxis.h
+++ b/ThirdParty/eigen3/Eigen/src/Geometry/AngleAxis.h
@@ -76,7 +76,7 @@ public:
     * \warning If the \a axis vector is not normalized, then the angle-axis object
     *          represents an invalid rotation. */
   template<typename Derived>
-  inline AngleAxis(Scalar angle, const MatrixBase<Derived>& axis) : m_axis(axis), m_angle(angle) {}
+  inline AngleAxis(const Scalar& angle, const MatrixBase<Derived>& axis) : m_axis(axis), m_angle(angle) {}
   /** Constructs and initialize the angle-axis rotation from a quaternion \a q. */
   template<typename QuatDerived> inline explicit AngleAxis(const QuaternionBase<QuatDerived>& q) { *this = q; }
   /** Constructs and initialize the angle-axis rotation from a 3x3 rotation matrix. */
@@ -137,7 +137,7 @@ public:
     * determined by \a prec.
     *
     * \sa MatrixBase::isApprox() */
-  bool isApprox(const AngleAxis& other, typename NumTraits<Scalar>::Real prec = NumTraits<Scalar>::dummy_precision()) const
+  bool isApprox(const AngleAxis& other, const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::dummy_precision()) const
   { return m_axis.isApprox(other.m_axis, prec) && internal::isApprox(m_angle,other.m_angle, prec); }
 };
 
@@ -161,6 +161,7 @@ AngleAxis<Scalar>& AngleAxis<Scalar>::operator=(const QuaternionBase<QuatDerived
   using std::acos;
   using std::min;
   using std::max;
+  using std::sqrt;
   Scalar n2 = q.vec().squaredNorm();
   if (n2 < NumTraits<Scalar>::dummy_precision()*NumTraits<Scalar>::dummy_precision())
   {
@@ -170,7 +171,7 @@ AngleAxis<Scalar>& AngleAxis<Scalar>::operator=(const QuaternionBase<QuatDerived
   else
   {
     m_angle = Scalar(2)*acos((min)((max)(Scalar(-1),q.w()),Scalar(1)));
-    m_axis = q.vec() / internal::sqrt(n2);
+    m_axis = q.vec() / sqrt(n2);
   }
   return *this;
 }
@@ -202,9 +203,11 @@ template<typename Scalar>
 typename AngleAxis<Scalar>::Matrix3
 AngleAxis<Scalar>::toRotationMatrix(void) const
 {
+  using std::sin;
+  using std::cos;
   Matrix3 res;
-  Vector3 sin_axis  = internal::sin(m_angle) * m_axis;
-  Scalar c = internal::cos(m_angle);
+  Vector3 sin_axis  = sin(m_angle) * m_axis;
+  Scalar c = cos(m_angle);
   Vector3 cos1_axis = (Scalar(1)-c) * m_axis;
 
   Scalar tmp;
diff --git a/ThirdParty/eigen3/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/ThirdParty/eigen3/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
index 126341be8d8..6fc6ab85225 100644
--- a/ThirdParty/eigen3/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
+++ b/ThirdParty/eigen3/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
@@ -39,10 +39,17 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
   int maxIters = iters;
 
   int n = mat.cols();
+  x = precond.solve(x);
   VectorType r  = rhs - mat * x;
   VectorType r0 = r;
   
   RealScalar r0_sqnorm = r0.squaredNorm();
+  RealScalar rhs_sqnorm = rhs.squaredNorm();
+  if(rhs_sqnorm == 0)
+  {
+    x.setZero();
+    return true;
+  }
   Scalar rho    = 1;
   Scalar alpha  = 1;
   Scalar w      = 1;
@@ -55,13 +62,22 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
 
   RealScalar tol2 = tol*tol;
   int i = 0;
+  int restarts = 0;
 
-  while ( r.squaredNorm()/r0_sqnorm > tol2 && i<maxIters )
+  while ( r.squaredNorm()/rhs_sqnorm > tol2 && i<maxIters )
   {
     Scalar rho_old = rho;
 
     rho = r0.dot(r);
-    if (rho == Scalar(0)) return false; /* New search directions cannot be found */
+    if (internal::isMuchSmallerThan(rho,r0_sqnorm))
+    {
+      // The new residual vector became too orthogonal to the arbitrarily choosen direction r0
+      // Let's restart with a new r0:
+      r0 = r;
+      rho = r0_sqnorm = r.squaredNorm();
+      if(restarts++ == 0)
+        i = 0;
+    }
     Scalar beta = (rho/rho_old) * (alpha / w);
     p = r + beta * (p - w * v);
     
@@ -75,12 +91,16 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
     z = precond.solve(s);
     t.noalias() = mat * z;
 
-    w = t.dot(s) / t.squaredNorm();
+    RealScalar tmp = t.squaredNorm();
+    if(tmp>RealScalar(0))
+      w = t.dot(s) / tmp;
+    else
+      w = Scalar(0);
     x += alpha * y + w * z;
     r = s - w * t;
     ++i;
   }
-  tol_error = sqrt(r.squaredNorm()/r0_sqnorm);
+  tol_error = sqrt(r.squaredNorm()/rhs_sqnorm);
   iters = i;
   return true; 
 }
@@ -223,7 +243,8 @@ public:
   template<typename Rhs,typename Dest>
   void _solve(const Rhs& b, Dest& x) const
   {
-    x.setZero();
+//     x.setZero();
+  x = b;
     _solveWithGuess(b,x);
   }
 
diff --git a/ThirdParty/eigen3/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h b/ThirdParty/eigen3/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
index f64f2534d28..a74a8155e68 100644
--- a/ThirdParty/eigen3/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
+++ b/ThirdParty/eigen3/Eigen/src/IterativeLinearSolvers/ConjugateGradient.h
@@ -41,15 +41,29 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
   int n = mat.cols();
 
   VectorType residual = rhs - mat * x; //initial residual
-  VectorType p(n);
 
+  RealScalar rhsNorm2 = rhs.squaredNorm();
+  if(rhsNorm2 == 0) 
+  {
+    x.setZero();
+    iters = 0;
+    tol_error = 0;
+    return;
+  }
+  RealScalar threshold = tol*tol*rhsNorm2;
+  RealScalar residualNorm2 = residual.squaredNorm();
+  if (residualNorm2 < threshold)
+  {
+    iters = 0;
+    tol_error = sqrt(residualNorm2 / rhsNorm2);
+    return;
+  }
+  
+  VectorType p(n);
   p = precond.solve(residual);      //initial search direction
 
   VectorType z(n), tmp(n);
-  RealScalar absNew = internal::real(residual.dot(p));  // the square of the absolute value of r scaled by invM
-  RealScalar rhsNorm2 = rhs.squaredNorm();
-  RealScalar residualNorm2 = 0;
-  RealScalar threshold = tol*tol*rhsNorm2;
+  RealScalar absNew = numext::real(residual.dot(p));  // the square of the absolute value of r scaled by invM
   int i = 0;
   while(i < maxIters)
   {
@@ -66,7 +80,7 @@ void conjugate_gradient(const MatrixType& mat, const Rhs& rhs, Dest& x,
     z = precond.solve(residual);          // approximately solve for "A z = residual"
 
     RealScalar absOld = absNew;
-    absNew = internal::real(residual.dot(z));     // update the absolute value of r
+    absNew = numext::real(residual.dot(z));     // update the absolute value of r
     RealScalar beta = absNew / absOld;            // calculate the Gram-Schmidt value used to create the new search direction
     p = z + beta * p;                             // update search direction
     i++;
diff --git a/ThirdParty/eigen3/Eigen/src/OrderingMethods/Amd.h b/ThirdParty/eigen3/Eigen/src/OrderingMethods/Amd.h
index ce04852b872..41b4fd7e392 100644
--- a/ThirdParty/eigen3/Eigen/src/OrderingMethods/Amd.h
+++ b/ThirdParty/eigen3/Eigen/src/OrderingMethods/Amd.h
@@ -2,10 +2,6 @@
 // for linear algebra.
 //
 // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr>
-//
-// This Source Code Form is subject to the terms of the Mozilla
-// Public License v. 2.0. If a copy of the MPL was not distributed
-// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 /*
 
@@ -86,6 +82,7 @@ Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Ind
 
 
 /** \internal
+  * \ingroup OrderingMethods_Module 
   * Approximate minimum degree ordering algorithm.
   * \returns the permutation P reducing the fill-in of the input matrix \a C
   * The input matrix \a C must be a selfadjoint compressed column major SparseMatrix object. Both the upper and lower parts have to be stored, but the diagonal entries are optional.
@@ -94,7 +91,6 @@ template<typename Scalar, typename Index>
 void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, PermutationMatrix<Dynamic,Dynamic,Index>& perm)
 {
   using std::sqrt;
-  typedef SparseMatrix<Scalar,ColMajor,Index> CCS;
   
   int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
       k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
diff --git a/ThirdParty/eigen3/Eigen/src/QR/ColPivHouseholderQR.h b/ThirdParty/eigen3/Eigen/src/QR/ColPivHouseholderQR.h
index 726b9fa99b7..8b01f8179e4 100644
--- a/ThirdParty/eigen3/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/ThirdParty/eigen3/Eigen/src/QR/ColPivHouseholderQR.h
@@ -55,7 +55,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
     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;
+    typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
     
   private:
     
@@ -94,6 +94,18 @@ template<typename _MatrixType> class ColPivHouseholderQR
         m_isInitialized(false),
         m_usePrescribedThreshold(false) {}
 
+    /** \brief Constructs a QR factorization from a given matrix
+      *
+      * This constructor computes the QR factorization of the matrix \a matrix by calling
+      * the method compute(). It is a short cut for:
+      * 
+      * \code
+      * ColPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols());
+      * qr.compute(matrix);
+      * \endcode
+      * 
+      * \sa compute()
+      */
     ColPivHouseholderQR(const MatrixType& matrix)
       : m_qr(matrix.rows(), matrix.cols()),
         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
@@ -133,6 +145,10 @@ template<typename _MatrixType> class ColPivHouseholderQR
     }
 
     HouseholderSequenceType householderQ(void) const;
+    HouseholderSequenceType matrixQ(void) const
+    {
+      return householderQ(); 
+    }
 
     /** \returns a reference to the matrix where the Householder QR decomposition is stored
       */
@@ -141,9 +157,25 @@ template<typename _MatrixType> class ColPivHouseholderQR
       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
       return m_qr;
     }
-
+    
+    /** \returns a reference to the matrix where the result Householder QR is stored 
+     * \warning The strict lower part of this matrix contains internal values. 
+     * Only the upper triangular part should be referenced. To get it, use
+     * \code matrixR().template triangularView<Upper>() \endcode
+     * For rank-deficient matrices, use 
+     * \code 
+     * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>() 
+     * \endcode
+     */
+    const MatrixType& matrixR() const
+    {
+      eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
+      return m_qr;
+    }
+    
     ColPivHouseholderQR& compute(const MatrixType& matrix);
 
+    /** \returns a const reference to the column permutation matrix */
     const PermutationType& colsPermutation() const
     {
       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
@@ -187,11 +219,12 @@ template<typename _MatrixType> class ColPivHouseholderQR
       */
     inline Index rank() const
     {
+      using std::abs;
       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
-      RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
+      RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
       Index result = 0;
       for(Index i = 0; i < m_nonzero_pivots; ++i)
-        result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);
+        result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
       return result;
     }
 
@@ -261,6 +294,11 @@ template<typename _MatrixType> class ColPivHouseholderQR
 
     inline Index rows() const { return m_qr.rows(); }
     inline Index cols() const { return m_qr.cols(); }
+    
+    /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q.
+      * 
+      * For advanced uses only.
+      */
     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
 
     /** Allows to prescribe a threshold to be used by certain methods, such as rank(),
@@ -331,6 +369,18 @@ template<typename _MatrixType> class ColPivHouseholderQR
       *          diagonal coefficient of R.
       */
     RealScalar maxPivot() const { return m_maxpivot; }
+    
+    /** \brief Reports whether the QR factorization was succesful.
+      *
+      * \note This function always returns \c Success. It is provided for compatibility 
+      * with other factorization routines.
+      * \returns \c Success 
+      */
+    ComputationInfo info() const
+    {
+      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
+      return Success;
+    }
 
   protected:
     MatrixType m_qr;
@@ -348,9 +398,10 @@ template<typename _MatrixType> class ColPivHouseholderQR
 template<typename MatrixType>
 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
 {
+  using std::abs;
   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
-  return internal::abs(m_qr.diagonal().prod());
+  return abs(m_qr.diagonal().prod());
 }
 
 template<typename MatrixType>
@@ -361,12 +412,22 @@ typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDetermina
   return m_qr.diagonal().cwiseAbs().array().log().sum();
 }
 
+/** Performs the QR factorization of the given matrix \a matrix. The result of
+  * the factorization is stored into \c *this, and a reference to \c *this
+  * is returned.
+  *
+  * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&)
+  */
 template<typename MatrixType>
 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
 {
+  using std::abs;
   Index rows = matrix.rows();
   Index cols = matrix.cols();
   Index size = matrix.diagonalSize();
+  
+  // the column permutation is stored as int indices, so just to be sure:
+  eigen_assert(cols<=NumTraits<int>::highest());
 
   m_qr = matrix;
   m_hCoeffs.resize(size);
@@ -380,7 +441,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
   for(Index k = 0; k < cols; ++k)
     m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
 
-  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * internal::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
+  RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
 
   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
   m_maxpivot = RealScalar(0);
@@ -432,7 +493,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
     m_qr.coeffRef(k,k) = beta;
 
     // remember the maximum absolute value of diagonal coefficients
-    if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta);
+    if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
 
     // apply the householder transformation
     m_qr.bottomRightCorner(rows-k, cols-k-1)
@@ -444,7 +505,7 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
 
   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_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
 
   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
   m_isInitialized = true;
@@ -464,8 +525,8 @@ struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
   {
     eigen_assert(rhs().rows() == dec().rows());
 
-    const int cols = dec().cols(),
-    nonzero_pivots = dec().nonzeroPivots();
+    const Index cols = dec().cols(),
+				nonzero_pivots = dec().nonzeroPivots();
 
     if(nonzero_pivots == 0)
     {
@@ -481,19 +542,11 @@ struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
 		     .transpose()
       );
 
-    dec().matrixQR()
+    dec().matrixR()
        .topLeftCorner(nonzero_pivots, nonzero_pivots)
        .template triangularView<Upper>()
        .solveInPlace(c.topRows(nonzero_pivots));
 
-
-    typename Rhs::PlainObject d(c);
-    d.topRows(nonzero_pivots)
-      = dec().matrixQR()
-       .topLeftCorner(nonzero_pivots, nonzero_pivots)
-       .template triangularView<Upper>()
-       * c.topRows(nonzero_pivots);
-
     for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
     for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
   }
diff --git a/ThirdParty/eigen3/Eigen/src/QR/ColPivHouseholderQR_MKL.h b/ThirdParty/eigen3/Eigen/src/QR/ColPivHouseholderQR_MKL.h
index ebcafe7dae5..b5b19832652 100644
--- a/ThirdParty/eigen3/Eigen/src/QR/ColPivHouseholderQR_MKL.h
+++ b/ThirdParty/eigen3/Eigen/src/QR/ColPivHouseholderQR_MKL.h
@@ -47,6 +47,7 @@ ColPivHouseholderQR<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynami
               const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>& matrix) \
 \
 { \
+  using std::abs; \
   typedef Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> MatrixType; \
   typedef MatrixType::Scalar Scalar; \
   typedef MatrixType::RealScalar RealScalar; \
@@ -71,10 +72,10 @@ ColPivHouseholderQR<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynami
   m_isInitialized = true; \
   m_maxpivot=m_qr.diagonal().cwiseAbs().maxCoeff(); \
   m_hCoeffs.adjointInPlace(); \
-  RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold(); \
+  RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); \
   lapack_int *perm = m_colsPermutation.indices().data(); \
   for(i=0;i<size;i++) { \
-    m_nonzero_pivots += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);\
+    m_nonzero_pivots += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);\
   } \
   for(i=0;i<cols;i++) perm[i]--;\
 \
diff --git a/ThirdParty/eigen3/Eigen/src/SparseCore/AmbiVector.h b/ThirdParty/eigen3/Eigen/src/SparseCore/AmbiVector.h
index 6cfaadbaa9a..17fff96a78b 100644
--- a/ThirdParty/eigen3/Eigen/src/SparseCore/AmbiVector.h
+++ b/ThirdParty/eigen3/Eigen/src/SparseCore/AmbiVector.h
@@ -288,9 +288,10 @@ class AmbiVector<_Scalar,_Index>::Iterator
       * In practice, all coefficients having a magnitude smaller than \a epsilon
       * are skipped.
       */
-    Iterator(const AmbiVector& vec, RealScalar epsilon = 0)
+    Iterator(const AmbiVector& vec, const RealScalar& epsilon = 0)
       : m_vector(vec)
     {
+      using std::abs;
       m_epsilon = epsilon;
       m_isDense = m_vector.m_mode==IsDense;
       if (m_isDense)
@@ -304,7 +305,7 @@ class AmbiVector<_Scalar,_Index>::Iterator
       {
         ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_vector.m_buffer);
         m_currentEl = m_vector.m_llStart;
-        while (m_currentEl>=0 && internal::abs(llElements[m_currentEl].value)<=m_epsilon)
+        while (m_currentEl>=0 && abs(llElements[m_currentEl].value)<=m_epsilon)
           m_currentEl = llElements[m_currentEl].next;
         if (m_currentEl<0)
         {
@@ -326,11 +327,12 @@ class AmbiVector<_Scalar,_Index>::Iterator
 
     Iterator& operator++()
     {
+      using std::abs;
       if (m_isDense)
       {
         do {
           ++m_cachedIndex;
-        } while (m_cachedIndex<m_vector.m_end && internal::abs(m_vector.m_buffer[m_cachedIndex])<m_epsilon);
+        } while (m_cachedIndex<m_vector.m_end && abs(m_vector.m_buffer[m_cachedIndex])<m_epsilon);
         if (m_cachedIndex<m_vector.m_end)
           m_cachedValue = m_vector.m_buffer[m_cachedIndex];
         else
@@ -341,7 +343,7 @@ class AmbiVector<_Scalar,_Index>::Iterator
         ListEl* EIGEN_RESTRICT llElements = reinterpret_cast<ListEl*>(m_vector.m_buffer);
         do {
           m_currentEl = llElements[m_currentEl].next;
-        } while (m_currentEl>=0 && internal::abs(llElements[m_currentEl].value)<m_epsilon);
+        } while (m_currentEl>=0 && abs(llElements[m_currentEl].value)<m_epsilon);
         if (m_currentEl<0)
         {
           m_cachedIndex = -1;
diff --git a/ThirdParty/eigen3/Eigen/src/SparseCore/CompressedStorage.h b/ThirdParty/eigen3/Eigen/src/SparseCore/CompressedStorage.h
index 85a998aff10..3321fab4a8a 100644
--- a/ThirdParty/eigen3/Eigen/src/SparseCore/CompressedStorage.h
+++ b/ThirdParty/eigen3/Eigen/src/SparseCore/CompressedStorage.h
@@ -139,7 +139,7 @@ class CompressedStorage
 
     /** \returns the stored value at index \a key
       * If the value does not exist, then the value \a defaultValue is returned without any insertion. */
-    inline Scalar at(Index key, Scalar defaultValue = Scalar(0)) const
+    inline Scalar at(Index key, const Scalar& defaultValue = Scalar(0)) const
     {
       if (m_size==0)
         return defaultValue;
@@ -152,7 +152,7 @@ class CompressedStorage
     }
 
     /** Like at(), but the search is performed in the range [start,end) */
-    inline Scalar atInRange(size_t start, size_t end, Index key, Scalar defaultValue = Scalar(0)) const
+    inline Scalar atInRange(size_t start, size_t end, Index key, const Scalar& defaultValue = Scalar(0)) const
     {
       if (start>=end)
         return Scalar(0);
@@ -167,7 +167,7 @@ class CompressedStorage
     /** \returns a reference to the value at index \a key
       * If the value does not exist, then the value \a defaultValue is inserted
       * such that the keys are sorted. */
-    inline Scalar& atWithInsertion(Index key, Scalar defaultValue = Scalar(0))
+    inline Scalar& atWithInsertion(Index key, const Scalar& defaultValue = Scalar(0))
     {
       size_t id = searchLowerIndex(0,m_size,key);
       if (id>=m_size || m_indices[id]!=key)
@@ -184,7 +184,7 @@ class CompressedStorage
       return m_values[id];
     }
 
-    void prune(Scalar reference, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision())
+    void prune(const Scalar& reference, const RealScalar& epsilon = NumTraits<RealScalar>::dummy_precision())
     {
       size_t k = 0;
       size_t n = size();
diff --git a/ThirdParty/eigen3/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h b/ThirdParty/eigen3/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h
index 16b5e1dba6c..4b13f08d4e4 100644
--- a/ThirdParty/eigen3/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h
+++ b/ThirdParty/eigen3/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h
@@ -121,9 +121,9 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r
 namespace internal {
 
 template<typename Lhs, typename Rhs, typename ResultType,
-  int LhsStorageOrder = traits<Lhs>::Flags&RowMajorBit,
-  int RhsStorageOrder = traits<Rhs>::Flags&RowMajorBit,
-  int ResStorageOrder = traits<ResultType>::Flags&RowMajorBit>
+  int LhsStorageOrder = (traits<Lhs>::Flags&RowMajorBit) ? RowMajor : ColMajor,
+  int RhsStorageOrder = (traits<Rhs>::Flags&RowMajorBit) ? RowMajor : ColMajor,
+  int ResStorageOrder = (traits<ResultType>::Flags&RowMajorBit) ? RowMajor : ColMajor>
 struct conservative_sparse_sparse_product_selector;
 
 template<typename Lhs, typename Rhs, typename ResultType>
diff --git a/ThirdParty/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.h b/ThirdParty/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.h
index 1e751ad629d..5c8c476eecf 100644
--- a/ThirdParty/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.h
+++ b/ThirdParty/eigen3/Eigen/src/plugins/ArrayCwiseBinaryOps.h
@@ -35,7 +35,12 @@ EIGEN_MAKE_CWISE_BINARY_OP(min,internal::scalar_min_op)
   */
 EIGEN_STRONG_INLINE const CwiseBinaryOp<internal::scalar_min_op<Scalar>, const Derived,
                                         const CwiseNullaryOp<internal::scalar_constant_op<Scalar>, PlainObject> >
-(min)(const Scalar &other) const
+#ifdef EIGEN_PARSED_BY_DOXYGEN
+min
+#else
+(min)
+#endif
+(const Scalar &other) const
 {
   return (min)(Derived::PlainObject::Constant(rows(), cols(), other));
 }
@@ -55,7 +60,12 @@ EIGEN_MAKE_CWISE_BINARY_OP(max,internal::scalar_max_op)
   */
 EIGEN_STRONG_INLINE const CwiseBinaryOp<internal::scalar_max_op<Scalar>, const Derived,
                                         const CwiseNullaryOp<internal::scalar_constant_op<Scalar>, PlainObject> >
-(max)(const Scalar &other) const
+#ifdef EIGEN_PARSED_BY_DOXYGEN
+max
+#else
+(max)
+#endif
+(const Scalar &other) const
 {
   return (max)(Derived::PlainObject::Constant(rows(), cols(), other));
 }
diff --git a/ThirdParty/eigen3/Eigen/src/plugins/BlockMethods.h b/ThirdParty/eigen3/Eigen/src/plugins/BlockMethods.h
index ef224001a54..6911bedef04 100644
--- a/ThirdParty/eigen3/Eigen/src/plugins/BlockMethods.h
+++ b/ThirdParty/eigen3/Eigen/src/plugins/BlockMethods.h
@@ -8,8 +8,6 @@
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-#ifndef EIGEN_BLOCKMETHODS_H
-#define EIGEN_BLOCKMETHODS_H
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
 
@@ -32,6 +30,10 @@ template<int N> struct ConstNColsBlockXpr { typedef const Block<const Derived, i
 template<int N> struct NRowsBlockXpr { typedef Block<Derived, N, internal::traits<Derived>::ColsAtCompileTime, IsRowMajor> Type; };
 template<int N> struct ConstNRowsBlockXpr { typedef const Block<const Derived, N, internal::traits<Derived>::ColsAtCompileTime, IsRowMajor> Type; };
 
+typedef VectorBlock<Derived> SegmentReturnType;
+typedef const VectorBlock<const Derived> ConstSegmentReturnType;
+template<int Size> struct FixedSegmentReturnType { typedef VectorBlock<Derived, Size> Type; };
+template<int Size> struct ConstFixedSegmentReturnType { typedef const VectorBlock<const Derived, Size> Type; };
 
 #endif // not EIGEN_PARSED_BY_DOXYGEN
 
@@ -88,12 +90,13 @@ inline const Block<const Derived> topRightCorner(Index cRows, Index cCols) const
 
 /** \returns an expression of a fixed-size top-right corner of *this.
   *
-  * The template parameters CRows and CCols are the number of rows and columns in the corner.
+  * \tparam CRows the number of rows in the corner
+  * \tparam CCols the number of columns in the corner
   *
   * Example: \include MatrixBase_template_int_int_topRightCorner.cpp
   * Output: \verbinclude MatrixBase_template_int_int_topRightCorner.out
   *
-  * \sa class Block, block(Index,Index,Index,Index)
+  * \sa class Block, block<int,int>(Index,Index)
   */
 template<int CRows, int CCols>
 inline Block<Derived, CRows, CCols> topRightCorner()
@@ -108,6 +111,35 @@ inline const Block<const Derived, CRows, CCols> topRightCorner() const
   return Block<const Derived, CRows, CCols>(derived(), 0, cols() - CCols);
 }
 
+/** \returns an expression of a top-right corner of *this.
+  *
+  * \tparam CRows number of rows in corner as specified at compile time
+  * \tparam CCols number of columns in corner as specified at compile time
+  * \param  cRows number of rows in corner as specified at run time
+  * \param  cCols number of columns in corner as specified at run time
+  *
+  * This function is mainly useful for corners where the number of rows is specified at compile time
+  * and the number of columns is specified at run time, or vice versa. The compile-time and run-time
+  * information should not contradict. In other words, \a cRows should equal \a CRows unless
+  * \a CRows is \a Dynamic, and the same for the number of columns.
+  *
+  * Example: \include MatrixBase_template_int_int_topRightCorner_int_int.cpp
+  * Output: \verbinclude MatrixBase_template_int_int_topRightCorner_int_int.out
+  *
+  * \sa class Block
+  */
+template<int CRows, int CCols>
+inline Block<Derived, CRows, CCols> topRightCorner(Index cRows, Index cCols)
+{
+  return Block<Derived, CRows, CCols>(derived(), 0, cols() - cCols, cRows, cCols);
+}
+
+/** This is the const version of topRightCorner<int, int>(Index, Index).*/
+template<int CRows, int CCols>
+inline const Block<const Derived, CRows, CCols> topRightCorner(Index cRows, Index cCols) const
+{
+  return Block<const Derived, CRows, CCols>(derived(), 0, cols() - cCols, cRows, cCols);
+}
 
 
 
@@ -154,6 +186,36 @@ inline const Block<const Derived, CRows, CCols> topLeftCorner() const
   return Block<const Derived, CRows, CCols>(derived(), 0, 0);
 }
 
+/** \returns an expression of a top-left corner of *this.
+  *
+  * \tparam CRows number of rows in corner as specified at compile time
+  * \tparam CCols number of columns in corner as specified at compile time
+  * \param  cRows number of rows in corner as specified at run time
+  * \param  cCols number of columns in corner as specified at run time
+  *
+  * This function is mainly useful for corners where the number of rows is specified at compile time
+  * and the number of columns is specified at run time, or vice versa. The compile-time and run-time
+  * information should not contradict. In other words, \a cRows should equal \a CRows unless
+  * \a CRows is \a Dynamic, and the same for the number of columns.
+  *
+  * Example: \include MatrixBase_template_int_int_topLeftCorner_int_int.cpp
+  * Output: \verbinclude MatrixBase_template_int_int_topLeftCorner_int_int.out
+  *
+  * \sa class Block
+  */
+template<int CRows, int CCols>
+inline Block<Derived, CRows, CCols> topLeftCorner(Index cRows, Index cCols)
+{
+  return Block<Derived, CRows, CCols>(derived(), 0, 0, cRows, cCols);
+}
+
+/** This is the const version of topLeftCorner<int, int>(Index, Index).*/
+template<int CRows, int CCols>
+inline const Block<const Derived, CRows, CCols> topLeftCorner(Index cRows, Index cCols) const
+{
+  return Block<const Derived, CRows, CCols>(derived(), 0, 0, cRows, cCols);
+}
+
 
 
 /** \returns a dynamic-size expression of a bottom-right corner of *this.
@@ -199,6 +261,36 @@ inline const Block<const Derived, CRows, CCols> bottomRightCorner() const
   return Block<const Derived, CRows, CCols>(derived(), rows() - CRows, cols() - CCols);
 }
 
+/** \returns an expression of a bottom-right corner of *this.
+  *
+  * \tparam CRows number of rows in corner as specified at compile time
+  * \tparam CCols number of columns in corner as specified at compile time
+  * \param  cRows number of rows in corner as specified at run time
+  * \param  cCols number of columns in corner as specified at run time
+  *
+  * This function is mainly useful for corners where the number of rows is specified at compile time
+  * and the number of columns is specified at run time, or vice versa. The compile-time and run-time
+  * information should not contradict. In other words, \a cRows should equal \a CRows unless
+  * \a CRows is \a Dynamic, and the same for the number of columns.
+  *
+  * Example: \include MatrixBase_template_int_int_bottomRightCorner_int_int.cpp
+  * Output: \verbinclude MatrixBase_template_int_int_bottomRightCorner_int_int.out
+  *
+  * \sa class Block
+  */
+template<int CRows, int CCols>
+inline Block<Derived, CRows, CCols> bottomRightCorner(Index cRows, Index cCols)
+{
+  return Block<Derived, CRows, CCols>(derived(), rows() - cRows, cols() - cCols, cRows, cCols);
+}
+
+/** This is the const version of bottomRightCorner<int, int>(Index, Index).*/
+template<int CRows, int CCols>
+inline const Block<const Derived, CRows, CCols> bottomRightCorner(Index cRows, Index cCols) const
+{
+  return Block<const Derived, CRows, CCols>(derived(), rows() - cRows, cols() - cCols, cRows, cCols);
+}
+
 
 
 /** \returns a dynamic-size expression of a bottom-left corner of *this.
@@ -244,6 +336,36 @@ inline const Block<const Derived, CRows, CCols> bottomLeftCorner() const
   return Block<const Derived, CRows, CCols>(derived(), rows() - CRows, 0);
 }
 
+/** \returns an expression of a bottom-left corner of *this.
+  *
+  * \tparam CRows number of rows in corner as specified at compile time
+  * \tparam CCols number of columns in corner as specified at compile time
+  * \param  cRows number of rows in corner as specified at run time
+  * \param  cCols number of columns in corner as specified at run time
+  *
+  * This function is mainly useful for corners where the number of rows is specified at compile time
+  * and the number of columns is specified at run time, or vice versa. The compile-time and run-time
+  * information should not contradict. In other words, \a cRows should equal \a CRows unless
+  * \a CRows is \a Dynamic, and the same for the number of columns.
+  *
+  * Example: \include MatrixBase_template_int_int_bottomLeftCorner_int_int.cpp
+  * Output: \verbinclude MatrixBase_template_int_int_bottomLeftCorner_int_int.out
+  *
+  * \sa class Block
+  */
+template<int CRows, int CCols>
+inline Block<Derived, CRows, CCols> bottomLeftCorner(Index cRows, Index cCols)
+{
+  return Block<Derived, CRows, CCols>(derived(), rows() - cRows, 0, cRows, cCols);
+}
+
+/** This is the const version of bottomLeftCorner<int, int>(Index, Index).*/
+template<int CRows, int CCols>
+inline const Block<const Derived, CRows, CCols> bottomLeftCorner(Index cRows, Index cCols) const
+{
+  return Block<const Derived, CRows, CCols>(derived(), rows() - cRows, 0, cRows, cCols);
+}
+
 
 
 /** \returns a block consisting of the top rows of *this.
@@ -543,6 +665,40 @@ inline const Block<const Derived, BlockRows, BlockCols> block(Index startRow, In
   return Block<const Derived, BlockRows, BlockCols>(derived(), startRow, startCol);
 }
 
+/** \returns an expression of a block in *this.
+  *
+  * \tparam BlockRows number of rows in block as specified at compile time
+  * \tparam BlockCols number of columns in block as specified at compile time
+  * \param  startRow  the first row in the block
+  * \param  startCol  the first column in the block
+  * \param  blockRows number of rows in block as specified at run time
+  * \param  blockCols number of columns in block as specified at run time
+  *
+  * This function is mainly useful for blocks where the number of rows is specified at compile time
+  * and the number of columns is specified at run time, or vice versa. The compile-time and run-time
+  * information should not contradict. In other words, \a blockRows should equal \a BlockRows unless
+  * \a BlockRows is \a Dynamic, and the same for the number of columns.
+  *
+  * Example: \include MatrixBase_template_int_int_block_int_int_int_int.cpp
+  * Output: \verbinclude MatrixBase_template_int_int_block_int_int_int_int.cpp
+  *
+  * \sa class Block, block(Index,Index,Index,Index)
+  */
+template<int BlockRows, int BlockCols>
+inline Block<Derived, BlockRows, BlockCols> block(Index startRow, Index startCol, 
+                                                  Index blockRows, Index blockCols)
+{
+  return Block<Derived, BlockRows, BlockCols>(derived(), startRow, startCol, blockRows, blockCols);
+}
+
+/** This is the const version of block<>(Index, Index, Index, Index). */
+template<int BlockRows, int BlockCols>
+inline const Block<const Derived, BlockRows, BlockCols> block(Index startRow, Index startCol,
+                                                              Index blockRows, Index blockCols) const
+{
+  return Block<const Derived, BlockRows, BlockCols>(derived(), startRow, startCol, blockRows, blockCols);
+}
+
 /** \returns an expression of the \a i-th column of *this. Note that the numbering starts at 0.
   *
   * Example: \include MatrixBase_col.cpp
@@ -577,4 +733,169 @@ inline ConstRowXpr row(Index i) const
   return ConstRowXpr(derived(), i);
 }
 
-#endif // EIGEN_BLOCKMETHODS_H
+/** \returns a dynamic-size expression of a segment (i.e. a vector block) in *this.
+  *
+  * \only_for_vectors
+  *
+  * \param start the first coefficient in the segment
+  * \param vecSize the number of coefficients in the segment
+  *
+  * Example: \include MatrixBase_segment_int_int.cpp
+  * Output: \verbinclude MatrixBase_segment_int_int.out
+  *
+  * \note Even though the returned expression has dynamic size, in the case
+  * when it is applied to a fixed-size vector, it inherits a fixed maximal size,
+  * which means that evaluating it does not cause a dynamic memory allocation.
+  *
+  * \sa class Block, segment(Index)
+  */
+inline SegmentReturnType segment(Index start, Index vecSize)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return SegmentReturnType(derived(), start, vecSize);
+}
+
+
+/** This is the const version of segment(Index,Index).*/
+inline ConstSegmentReturnType segment(Index start, Index vecSize) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return ConstSegmentReturnType(derived(), start, vecSize);
+}
+
+/** \returns a dynamic-size expression of the first coefficients of *this.
+  *
+  * \only_for_vectors
+  *
+  * \param vecSize the number of coefficients in the block
+  *
+  * Example: \include MatrixBase_start_int.cpp
+  * Output: \verbinclude MatrixBase_start_int.out
+  *
+  * \note Even though the returned expression has dynamic size, in the case
+  * when it is applied to a fixed-size vector, it inherits a fixed maximal size,
+  * which means that evaluating it does not cause a dynamic memory allocation.
+  *
+  * \sa class Block, block(Index,Index)
+  */
+inline SegmentReturnType head(Index vecSize)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return SegmentReturnType(derived(), 0, vecSize);
+}
+
+/** This is the const version of head(Index).*/
+inline ConstSegmentReturnType
+  head(Index vecSize) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return ConstSegmentReturnType(derived(), 0, vecSize);
+}
+
+/** \returns a dynamic-size expression of the last coefficients of *this.
+  *
+  * \only_for_vectors
+  *
+  * \param vecSize the number of coefficients in the block
+  *
+  * Example: \include MatrixBase_end_int.cpp
+  * Output: \verbinclude MatrixBase_end_int.out
+  *
+  * \note Even though the returned expression has dynamic size, in the case
+  * when it is applied to a fixed-size vector, it inherits a fixed maximal size,
+  * which means that evaluating it does not cause a dynamic memory allocation.
+  *
+  * \sa class Block, block(Index,Index)
+  */
+inline SegmentReturnType tail(Index vecSize)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return SegmentReturnType(derived(), this->size() - vecSize, vecSize);
+}
+
+/** This is the const version of tail(Index).*/
+inline ConstSegmentReturnType tail(Index vecSize) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return ConstSegmentReturnType(derived(), this->size() - vecSize, vecSize);
+}
+
+/** \returns a fixed-size expression of a segment (i.e. a vector block) in \c *this
+  *
+  * \only_for_vectors
+  *
+  * The template parameter \a Size is the number of coefficients in the block
+  *
+  * \param start the index of the first element of the sub-vector
+  *
+  * Example: \include MatrixBase_template_int_segment.cpp
+  * Output: \verbinclude MatrixBase_template_int_segment.out
+  *
+  * \sa class Block
+  */
+template<int Size>
+inline typename FixedSegmentReturnType<Size>::Type segment(Index start)
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return typename FixedSegmentReturnType<Size>::Type(derived(), start);
+}
+
+/** This is the const version of segment<int>(Index).*/
+template<int Size>
+inline typename ConstFixedSegmentReturnType<Size>::Type segment(Index start) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return typename ConstFixedSegmentReturnType<Size>::Type(derived(), start);
+}
+
+/** \returns a fixed-size expression of the first coefficients of *this.
+  *
+  * \only_for_vectors
+  *
+  * The template parameter \a Size is the number of coefficients in the block
+  *
+  * Example: \include MatrixBase_template_int_start.cpp
+  * Output: \verbinclude MatrixBase_template_int_start.out
+  *
+  * \sa class Block
+  */
+template<int Size>
+inline typename FixedSegmentReturnType<Size>::Type head()
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return typename FixedSegmentReturnType<Size>::Type(derived(), 0);
+}
+
+/** This is the const version of head<int>().*/
+template<int Size>
+inline typename ConstFixedSegmentReturnType<Size>::Type head() const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return typename ConstFixedSegmentReturnType<Size>::Type(derived(), 0);
+}
+
+/** \returns a fixed-size expression of the last coefficients of *this.
+  *
+  * \only_for_vectors
+  *
+  * The template parameter \a Size is the number of coefficients in the block
+  *
+  * Example: \include MatrixBase_template_int_end.cpp
+  * Output: \verbinclude MatrixBase_template_int_end.out
+  *
+  * \sa class Block
+  */
+template<int Size>
+inline typename FixedSegmentReturnType<Size>::Type tail()
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return typename FixedSegmentReturnType<Size>::Type(derived(), size() - Size);
+}
+
+/** This is the const version of tail<int>.*/
+template<int Size>
+inline typename ConstFixedSegmentReturnType<Size>::Type tail() const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  return typename ConstFixedSegmentReturnType<Size>::Type(derived(), size() - Size);
+}
-- 
GitLab