bump to 2.0.16

(matrix-vector product did not use the nesting rules....)

.hgignore

@@ -0,0 +1,25 @@
+syntax: glob
+qrc_*cxx
+*.orig
+*.pyc
+*.diff
+diff
+*.save
+save
+*.old
+*.gmo
+*.qm
+core
+core.*
+*.bak
+*~
+build*
+*.moc.*
+*.moc
+ui_*
+CMakeCache.txt
+tags
+.*.swp
+activity.png
+*.out
+*.php*

CMakeLists.txt

@@ -1,22 +1,15 @@
 project(Eigen)
-set(EIGEN_VERSION_NUMBER "2.0.3")
-
-#if the svnversion program is absent, this will leave the SVN_REVISION string empty,
-#but won't stop CMake.
-execute_process(COMMAND svnversion -n ${CMAKE_SOURCE_DIR}
-                OUTPUT_VARIABLE EIGEN_SVNVERSION_OUTPUT)
-
-#we only want EIGEN_SVN_REVISION if it is an actual revision number, not a string like "exported"
-string(REGEX MATCH "^[0-9]+.*" EIGEN_SVN_REVISION "${EIGEN_SVNVERSION_OUTPUT}")
-
-if(EIGEN_SVN_REVISION)
-  set(EIGEN_VERSION "${EIGEN_VERSION_NUMBER} (SVN revision ${EIGEN_SVN_REVISION})")
-else(EIGEN_SVN_REVISION)
-  set(EIGEN_VERSION "${EIGEN_VERSION_NUMBER}")
-endif(EIGEN_SVN_REVISION)
-
 cmake_minimum_required(VERSION 2.6.2)
 
+set(INCLUDE_INSTALL_DIR
+    "${CMAKE_INSTALL_PREFIX}/include/eigen2"
+    CACHE PATH
+    "The directory where we install the header files"
+    FORCE)
+
+set(EIGEN_VERSION_NUMBER "2.0.16")
+set(EIGEN_VERSION "${EIGEN_VERSION_NUMBER}")
+
 set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
 
 option(EIGEN_BUILD_TESTS "Build tests" OFF)
@@ -25,16 +18,56 @@ if(NOT WIN32)
   option(EIGEN_BUILD_LIB "Build the binary shared library" OFF)
 endif(NOT WIN32)
 option(EIGEN_BUILD_BTL "Build benchmark suite" OFF)
+if(NOT WIN32)
+  option(EIGEN_BUILD_PKGCONFIG "Build pkg-config .pc file for Eigen" ON)
+endif(NOT WIN32)
 
 if(EIGEN_BUILD_LIB)
   option(EIGEN_TEST_LIB "Build the unit tests using the library (disable -pedantic)" OFF)
 endif(EIGEN_BUILD_LIB)
 
+
+#############################################################################
+# find how to link to the standard libraries                                #
+#############################################################################
+
+find_package(StandardMathLibrary)
+
+set(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO "")
+
+if(NOT STANDARD_MATH_LIBRARY_FOUND)
+
+  message(FATAL_ERROR
+    "Can't link to the standard math library. Please report to the Eigen developers, telling them about your platform.")
+
+else()
+
+  if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
+    set(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO "${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO} ${STANDARD_MATH_LIBRARY}")
+  else()
+    set(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO "${STANDARD_MATH_LIBRARY}")
+  endif()
+
+endif()
+
+if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
+  message(STATUS "Standard libraries to link to explicitly: ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO}")
+else()
+  message(STATUS "Standard libraries to link to explicitly: none")
+endif()
+
+
+
 set(CMAKE_INCLUDE_CURRENT_DIR ON)
 
 if(CMAKE_COMPILER_IS_GNUCXX)
   if(CMAKE_SYSTEM_NAME MATCHES Linux)
-    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wnon-virtual-dtor -Wno-long-long -ansi -Wundef -Wcast-align -Wchar-subscripts -Wall -W -Wpointer-arith -Wwrite-strings -Wformat-security -Wextra -fno-exceptions -fno-check-new -fno-common -fstrict-aliasing")
+    include(CheckCXXCompilerFlag)
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wnon-virtual-dtor -Wno-long-long -ansi -Wundef -Wcast-align -Wchar-subscripts -Wall -W -Wpointer-arith -Wwrite-strings -Wformat-security -fno-exceptions -fno-check-new -fno-common -fstrict-aliasing")
+    check_cxx_compiler_flag("-Wextra" has_wextra)
+    if(has_wextra)
+      set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wextra")
+    endif()
     if(NOT EIGEN_TEST_LIB)
       set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pedantic")
     endif(NOT EIGEN_TEST_LIB)
@@ -84,6 +117,13 @@ endif(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
 
 include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR})
 
+if(EIGEN_BUILD_PKGCONFIG)
+    configure_file(eigen2.pc.in eigen2.pc) # uses INCLUDE_INSTALL_DIR
+    install(FILES ${CMAKE_CURRENT_BINARY_DIR}/eigen2.pc
+        DESTINATION share/pkgconfig
+        )
+endif(EIGEN_BUILD_PKGCONFIG)
+
 add_subdirectory(Eigen)
 add_subdirectory(unsupported)
 

CTestConfig.cmake

@@ -3,11 +3,11 @@
 ## project to incorporate the testing dashboard.
 ## # The following are required to uses Dart and the Cdash dashboard
 ##   ENABLE_TESTING()
-##   INCLUDE(Dart)
-set(CTEST_PROJECT_NAME "Eigen")
-set(CTEST_NIGHTLY_START_TIME "05:00:00 UTC")
+##   INCLUDE(CTest)
+set(CTEST_PROJECT_NAME "Eigen 2.0")
+set(CTEST_NIGHTLY_START_TIME "06:00:00 UTC")
 
 set(CTEST_DROP_METHOD "http")
-set(CTEST_DROP_SITE "www.cdash.org")
-set(CTEST_DROP_LOCATION "/CDashPublic/submit.php?project=Eigen")
+set(CTEST_DROP_SITE "eigen.tuxfamily.org")
+set(CTEST_DROP_LOCATION "/CDash/submit.php?project=Eigen+2.0")
 set(CTEST_DROP_SITE_CDASH TRUE)

Eigen/CMakeLists.txt

@@ -1,4 +1,7 @@
-set(Eigen_HEADERS Core LU Cholesky QR Geometry Sparse Array SVD LeastSquares QtAlignedMalloc StdVector)
+set(Eigen_HEADERS Core LU Cholesky QR Geometry
+                  Sparse Array SVD LeastSquares
+                  QtAlignedMalloc StdVector NewStdVector
+                  Eigen Dense)
 
 if(EIGEN_BUILD_LIB)
     set(Eigen_SRCS
@@ -20,12 +23,6 @@ if(CMAKE_COMPILER_IS_GNUCXX)
     set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} -g1 -O2")
 endif(CMAKE_COMPILER_IS_GNUCXX)
 
-set(INCLUDE_INSTALL_DIR
-    "${CMAKE_INSTALL_PREFIX}/include/eigen2"
-    CACHE PATH
-    "The directory where we install the header files"
-    FORCE)
-
 install(FILES
   ${Eigen_HEADERS}
   DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen

Eigen/Cholesky

@@ -38,8 +38,8 @@ namespace Eigen {
 } // namespace Eigen
 
 #define EIGEN_CHOLESKY_MODULE_INSTANTIATE_TYPE(MATRIXTYPE,PREFIX) \
-  PREFIX template class Cholesky<MATRIXTYPE>; \
-  PREFIX template class CholeskyWithoutSquareRoot<MATRIXTYPE>
+  PREFIX template class LLT<MATRIXTYPE>; \
+  PREFIX template class LDLT<MATRIXTYPE>
 
 #define EIGEN_CHOLESKY_MODULE_INSTANTIATE(PREFIX) \
   EIGEN_CHOLESKY_MODULE_INSTANTIATE_TYPE(Matrix2f,PREFIX); \

Eigen/Core

@@ -15,7 +15,9 @@
   #endif
 #endif
 
-#ifdef __GNUC__
+// FIXME: this check should not be against __QNXNTO__, which is also defined
+// while compiling with GCC for QNX target. Better solution is welcome!
+#if defined(__GNUC__) && !defined(__QNXNTO__)
   #define EIGEN_GNUC_AT_LEAST(x,y) ((__GNUC__>=x && __GNUC_MINOR__>=y) || __GNUC__>x)
 #else
   #define EIGEN_GNUC_AT_LEAST(x,y) 0
@@ -26,7 +28,7 @@
   #define EIGEN_SSE2_BUT_NOT_OLD_GCC
 #endif
 
-#ifndef EIGEN_DONT_VECTORIZE
+#if !defined(EIGEN_DONT_VECTORIZE) && !defined(EIGEN_DONT_ALIGN)
   #if defined (EIGEN_SSE2_BUT_NOT_OLD_GCC) || defined(EIGEN_SSE2_ON_MSVC_2008_OR_LATER)
     #define EIGEN_VECTORIZE
     #define EIGEN_VECTORIZE_SSE
@@ -50,6 +52,7 @@
   #endif
 #endif
 
+#include <cstddef>
 #include <cstdlib>
 #include <cmath>
 #include <complex>
@@ -77,6 +80,10 @@
 
 namespace Eigen {
 
+// we use size_t frequently and we'll never remember to prepend it with std:: everytime just to
+// ensure QNX/QCC support
+using std::size_t;
+
 /** \defgroup Core_Module Core module
   * This is the main module of Eigen providing dense matrix and vector support
   * (both fixed and dynamic size) with all the features corresponding to a BLAS library

Eigen/NewStdVector

@@ -0,0 +1,168 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2009 Hauke Heibel <hauke.heibel@googlemail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_STDVECTOR_MODULE_H
+#define EIGEN_STDVECTOR_MODULE_H
+
+#include "Core"
+#include <vector>
+
+namespace Eigen {
+
+// This one is needed to prevent reimplementing the whole std::vector.
+template <class T>
+class aligned_allocator_indirection : public aligned_allocator<T>
+{
+public:
+  typedef std::size_t    size_type;
+  typedef std::ptrdiff_t difference_type;
+  typedef T*        pointer;
+  typedef const T*  const_pointer;
+  typedef T&        reference;
+  typedef const T&  const_reference;
+  typedef T         value_type;
+
+  template<class U>
+  struct rebind
+  {
+    typedef aligned_allocator_indirection<U> other;
+  };
+
+  aligned_allocator_indirection() throw() {}
+  aligned_allocator_indirection(const aligned_allocator_indirection& ) throw() : aligned_allocator<T>() {}
+  aligned_allocator_indirection(const aligned_allocator<T>& ) throw() {}
+  template<class U>
+  aligned_allocator_indirection(const aligned_allocator_indirection<U>& ) throw() {}
+  template<class U>
+  aligned_allocator_indirection(const aligned_allocator<U>& ) throw() {}
+  ~aligned_allocator_indirection() throw() {}
+};
+
+#ifdef _MSC_VER
+
+  // sometimes, MSVC detects, at compile time, that the argument x
+  // in std::vector::resize(size_t s,T x) won't be aligned and generate an error
+  // even if this function is never called. Whence this little wrapper.
+  #define EIGEN_WORKAROUND_MSVC_STD_VECTOR(T) Eigen::ei_workaround_msvc_std_vector<T>
+  template<typename T> struct ei_workaround_msvc_std_vector : public T
+  {
+    inline ei_workaround_msvc_std_vector() : T() {}
+    inline ei_workaround_msvc_std_vector(const T& other) : T(other) {}
+    inline operator T& () { return *static_cast<T*>(this); }
+    inline operator const T& () const { return *static_cast<const T*>(this); }
+    template<typename OtherT>
+    inline T& operator=(const OtherT& other)
+    { T::operator=(other); return *this; }
+    inline ei_workaround_msvc_std_vector& operator=(const ei_workaround_msvc_std_vector& other)
+    { T::operator=(other); return *this; }
+  };
+
+#else
+
+  #define EIGEN_WORKAROUND_MSVC_STD_VECTOR(T) T
+
+#endif
+
+}
+
+namespace std {
+
+#define EIGEN_STD_VECTOR_SPECIALIZATION_BODY \
+  public:  \
+    typedef T value_type; \
+    typedef typename vector_base::allocator_type allocator_type; \
+    typedef typename vector_base::size_type size_type;  \
+    typedef typename vector_base::iterator iterator;  \
+    typedef typename vector_base::const_iterator const_iterator;  \
+    explicit vector(const allocator_type& a = allocator_type()) : vector_base(a) {}  \
+    template<typename InputIterator> \
+    vector(InputIterator first, InputIterator last, const allocator_type& a = allocator_type()) \
+    : vector_base(first, last, a) {} \
+    vector(const vector& c) : vector_base(c) {}  \
+    explicit vector(size_type num, const value_type& val = value_type()) : vector_base(num, val) {} \
+    vector(iterator start, iterator end) : vector_base(start, end) {}  \
+    vector& operator=(const vector& x) {  \
+      vector_base::operator=(x);  \
+      return *this;  \
+    }
+
+template<typename T>
+class vector<T,Eigen::aligned_allocator<T> >
+  : public vector<EIGEN_WORKAROUND_MSVC_STD_VECTOR(T),
+                  Eigen::aligned_allocator_indirection<EIGEN_WORKAROUND_MSVC_STD_VECTOR(T)> >
+{
+  typedef vector<EIGEN_WORKAROUND_MSVC_STD_VECTOR(T),
+                 Eigen::aligned_allocator_indirection<EIGEN_WORKAROUND_MSVC_STD_VECTOR(T)> > vector_base;
+  EIGEN_STD_VECTOR_SPECIALIZATION_BODY
+
+  void resize(size_type new_size)
+  { resize(new_size, T()); }
+
+#if defined(_VECTOR_)
+  // workaround MSVC std::vector implementation
+  void resize(size_type new_size, const value_type& x)
+  {
+    if (vector_base::size() < new_size)
+      vector_base::_Insert_n(vector_base::end(), new_size - vector_base::size(), x);
+    else if (new_size < vector_base::size())
+      vector_base::erase(vector_base::begin() + new_size, vector_base::end());
+  }
+  void push_back(const value_type& x)
+  { vector_base::push_back(x); } 
+  using vector_base::insert;  
+  iterator insert(const_iterator position, const value_type& x)
+  { return vector_base::insert(position,x); }
+  void insert(const_iterator position, size_type new_size, const value_type& x)
+  { vector_base::insert(position, new_size, x); }
+#elif defined(_GLIBCXX_VECTOR) && EIGEN_GNUC_AT_LEAST(4,2)
+  // workaround GCC std::vector implementation
+  void resize(size_type new_size, const value_type& x)
+  {
+    if (new_size < vector_base::size())
+      vector_base::_M_erase_at_end(this->_M_impl._M_start + new_size);
+    else
+      vector_base::insert(vector_base::end(), new_size - vector_base::size(), x);
+  }
+#elif defined(_GLIBCXX_VECTOR) && (!EIGEN_GNUC_AT_LEAST(4,1))
+  // Note that before gcc-4.1 we already have: std::vector::resize(size_type,const T&),
+  // no no need to workaround !
+  using vector_base::resize;
+#else
+  // either GCC 4.1 or non-GCC
+  // default implementation which should always work.
+  void resize(size_type new_size, const value_type& x)
+  {
+    if (new_size < vector_base::size())
+      vector_base::erase(vector_base::begin() + new_size, vector_base::end());
+    else if (new_size > vector_base::size())
+      vector_base::insert(vector_base::end(), new_size - vector_base::size(), x);
+  }
+#endif
+
+};
+
+}
+
+#endif // EIGEN_STDVECTOR_MODULE_H

Eigen/QtAlignedMalloc

@@ -1,10 +1,23 @@
-
 #ifndef EIGEN_QTMALLOC_MODULE_H
 #define EIGEN_QTMALLOC_MODULE_H
 
+#if (!EIGEN_MALLOC_ALREADY_ALIGNED)
+
+#ifdef QVECTOR_H
+#error You must include <Eigen/QtAlignedMalloc> before <QtCore/QVector>.
+#endif                                           
+
+#ifdef Q_DECL_IMPORT
+  #define Q_DECL_IMPORT_ORIG Q_DECL_IMPORT
+  #undef Q_DECL_IMPORT
+  #define Q_DECL_IMPORT
+#else
+  #define Q_DECL_IMPORT
+#endif
+
 #include "Core"
 
-#if (!EIGEN_MALLOC_ALREADY_ALIGNED)
+#include <QtCore/QVector>
 
 inline void *qMalloc(size_t size)
 {
@@ -26,4 +39,11 @@ inline void *qRealloc(void *ptr, size_t size)
 
 #endif
 
+#ifdef Q_DECL_IMPORT_ORIG
+  #define Q_DECL_IMPORT Q_DECL_IMPORT_ORIG
+  #undef Q_DECL_IMPORT_ORIG
+#else
+  #undef Q_DECL_IMPORT
+#endif
+
 #endif // EIGEN_QTMALLOC_MODULE_H

Eigen/StdVector

@@ -1,8 +1,12 @@
+#ifdef EIGEN_USE_NEW_STDVECTOR
+#include "NewStdVector"
+#else
+
 #ifndef EIGEN_STDVECTOR_MODULE_H
 #define EIGEN_STDVECTOR_MODULE_H
 
 #if defined(_GLIBCXX_VECTOR) || defined(_VECTOR_)
-#error you must include Eigen/StdVector before std::vector
+#error you must include <Eigen/StdVector> before <vector>. Also note that <Eigen/Sparse> includes <vector>, so it must be included after <Eigen/StdVector> too.
 #endif                                                    
 
 #ifndef EIGEN_GNUC_AT_LEAST
@@ -112,10 +116,8 @@ class vector<T,DummyAlloc,true>
     else if (__new_size < vector_base::size())
       vector_base::erase(vector_base::begin() + __new_size, vector_base::end());
   }
-  #elif defined(_GLIBCXX_VECTOR) && EIGEN_GNUC_AT_LEAST(4,1)
+  #elif defined(_GLIBCXX_VECTOR) && EIGEN_GNUC_AT_LEAST(4,2)
   // workaround GCC std::vector implementation
-  // Note that before gcc-4.1 we already have: std::vector::resize(size_type,const T&),
-  // no no need to workaround !
   void resize(size_type __new_size, const value_type& __x)
   {                                              
     if (__new_size < vector_base::size())               
@@ -123,7 +125,17 @@ class vector<T,DummyAlloc,true>
     else                       
       vector_base::insert(vector_base::end(), __new_size - vector_base::size(), __x); 
   }                                                              
+  #elif defined(_GLIBCXX_VECTOR) && EIGEN_GNUC_AT_LEAST(4,1)
+  void resize(size_type __new_size, const value_type& __x)
+  {
+    if (__new_size < vector_base::size())
+      vector_base::erase(vector_base::begin() + __new_size, vector_base::end());
+    else
+      vector_base::insert(vector_base::end(), __new_size - vector_base::size(), __x);
+  }
   #else
+  // Before gcc-4.1 we already have: std::vector::resize(size_type,const T&),
+  // so no need for a workaround !
   using vector_base::resize;
   #endif  
 };
@@ -131,3 +143,5 @@ class vector<T,DummyAlloc,true>
 }
 
 #endif // EIGEN_STDVECTOR_MODULE_H
+
+#endif // EIGEN_USE_NEW_STDVECTOR

Eigen/src/Array/BooleanRedux.h

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

Eigen/src/Array/Functors.h

@@ -43,6 +43,8 @@ struct ei_scalar_add_op {
   inline const PacketScalar packetOp(const PacketScalar& a) const
   { return ei_padd(a, ei_pset1(m_other)); }
   const Scalar m_other;
+private:
+  ei_scalar_add_op& operator=(const ei_scalar_add_op&);
 };
 template<typename Scalar>
 struct ei_functor_traits<ei_scalar_add_op<Scalar> >
@@ -138,6 +140,8 @@ struct ei_scalar_pow_op {
   inline ei_scalar_pow_op(const Scalar& exponent) : m_exponent(exponent) {}
   inline Scalar operator() (const Scalar& a) const { return ei_pow(a, m_exponent); }
   const Scalar m_exponent;
+private:
+  ei_scalar_pow_op& operator=(const ei_scalar_pow_op&);
 };
 template<typename Scalar>
 struct ei_functor_traits<ei_scalar_pow_op<Scalar> >

Eigen/src/Array/PartialRedux.h

@@ -133,6 +133,8 @@ struct ei_member_redux {
   inline result_type operator()(const MatrixBase<Derived>& mat) const
   { return mat.redux(m_functor); }
   const BinaryOp m_functor;
+private:
+  ei_member_redux& operator=(const ei_member_redux&);
 };
 
 /** \array_module \ingroup Array
@@ -158,13 +160,15 @@ template<typename ExpressionType, int Direction> class PartialRedux
   public:
 
     typedef typename ei_traits<ExpressionType>::Scalar Scalar;
+    typedef typename NumTraits<Scalar>::Real RealScalar;
     typedef typename ei_meta_if<ei_must_nest_by_value<ExpressionType>::ret,
         ExpressionType, const ExpressionType&>::ret ExpressionTypeNested;
 
-    template<template<typename _Scalar> class Functor> struct ReturnType
+    template<template<typename _Scalar> class Functor,
+             typename Scalar = typename ei_traits<ExpressionType>::Scalar> struct ReturnType
     {
       typedef PartialReduxExpr<ExpressionType,
-                               Functor<typename ei_traits<ExpressionType>::Scalar>,
+                               Functor<Scalar>,
                                Direction
                               > Type;
     };
@@ -215,7 +219,7 @@ template<typename ExpressionType, int Direction> class PartialRedux
       * Output: \verbinclude PartialRedux_squaredNorm.out
       *
       * \sa MatrixBase::squaredNorm() */
-    const typename ReturnType<ei_member_squaredNorm>::Type squaredNorm() const
+    const typename ReturnType<ei_member_squaredNorm,RealScalar>::Type squaredNorm() const
     { return _expression(); }
 
     /** \returns a row (or column) vector expression of the norm
@@ -225,7 +229,7 @@ template<typename ExpressionType, int Direction> class PartialRedux
       * Output: \verbinclude PartialRedux_norm.out
       *
       * \sa MatrixBase::norm() */
-    const typename ReturnType<ei_member_norm>::Type norm() const
+    const typename ReturnType<ei_member_norm,RealScalar>::Type norm() const
     { return _expression(); }
 
     /** \returns a row (or column) vector expression of the sum
@@ -290,6 +294,9 @@ template<typename ExpressionType, int Direction> class PartialRedux
 
   protected:
     ExpressionTypeNested m_matrix;
+
+  private:
+    PartialRedux& operator=(const PartialRedux&);
 };
 
 /** \array_module

Eigen/src/Cholesky/LDLT.h

@@ -96,8 +96,7 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
   assert(a.rows()==a.cols());
   const int size = a.rows();
   m_matrix.resize(size, size);
-  m_isPositiveDefinite = true;
-  const RealScalar eps = ei_sqrt(precision<Scalar>());
+  m_isPositiveDefinite = true; // always true. This decomposition is not rank-revealing anyway.
 
   if (size<=1)
   {
@@ -121,12 +120,6 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
     RealScalar tmp = ei_real(a.coeff(j,j) - (m_matrix.row(j).start(j) * m_matrix.col(j).start(j).conjugate()).coeff(0,0));
     m_matrix.coeffRef(j,j) = tmp;
 
-    if (tmp < eps)
-    {
-      m_isPositiveDefinite = false;
-      return;
-    }
-
     int endSize = size-j-1;
     if (endSize>0)
     {
@@ -136,7 +129,8 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
       m_matrix.row(j).end(endSize) = a.row(j).end(endSize).conjugate()
                                    - _temporary.end(endSize).transpose();
 
-      m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp;
+      if(tmp != RealScalar(0))
+        m_matrix.col(j).end(endSize) = m_matrix.row(j).end(endSize) / tmp;
     }
   }
 }
@@ -192,7 +186,7 @@ template<typename Derived>
 inline const LDLT<typename MatrixBase<Derived>::PlainMatrixType>
 MatrixBase<Derived>::ldlt() const
 {
-  return derived();
+  return LDLT<PlainMatrixType>(derived());
 }
 
 #endif // EIGEN_LDLT_H

Eigen/src/Cholesky/LLT.h

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

Eigen/src/Core/Assign.h

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

Eigen/src/Core/Block.h

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

Eigen/src/Core/CacheFriendlyProduct.h

@@ -33,7 +33,7 @@ struct ei_L2_block_traits {
 #ifndef EIGEN_EXTERN_INSTANTIATIONS
 
 template<typename Scalar>
-static void ei_cache_friendly_product(
+void ei_cache_friendly_product(
   int _rows, int _cols, int depth,
   bool _lhsRowMajor, const Scalar* _lhs, int _lhsStride,
   bool _rhsRowMajor, const Scalar* _rhs, int _rhsStride,
@@ -84,7 +84,7 @@ static void ei_cache_friendly_product(
     MaxL2BlockSize = ei_L2_block_traits<EIGEN_TUNE_FOR_CPU_CACHE_SIZE,Scalar>::width
   };
 
-  const bool resIsAligned = (PacketSize==1) || (((resStride%PacketSize) == 0) && (size_t(res)%16==0));
+  const bool resIsAligned = (PacketSize==1) || (((resStride%PacketSize) == 0) && (std::size_t(res)%16==0));
 
   const int remainingSize = depth % PacketSize;
   const int size = depth - remainingSize; // third dimension of the product clamped to packet boundaries
@@ -92,7 +92,7 @@ static void ei_cache_friendly_product(
   const int l2BlockCols = MaxL2BlockSize > cols ? cols : MaxL2BlockSize;
   const int l2BlockSize = MaxL2BlockSize > size ? size : MaxL2BlockSize;
   const int l2BlockSizeAligned = (1 + std::max(l2BlockSize,l2BlockCols)/PacketSize)*PacketSize;
-  const bool needRhsCopy = (PacketSize>1) && ((rhsStride%PacketSize!=0) || (size_t(rhs)%16!=0));
+  const bool needRhsCopy = (PacketSize>1) && ((rhsStride%PacketSize!=0) || (std::size_t(rhs)%16!=0));
   Scalar* EIGEN_RESTRICT block = 0;
   const int allocBlockSize = l2BlockRows*size;
   block = ei_aligned_stack_new(Scalar, allocBlockSize);
@@ -172,7 +172,7 @@ static void ei_cache_friendly_product(
           for(int l1j=l2j; l1j<l2blockColEnd; l1j+=1)
           {
             ei_internal_assert(l2BlockSizeAligned*(l1j-l2j)+(l2blockSizeEnd-l2k) < l2BlockSizeAligned*l2BlockSizeAligned);
-            memcpy(rhsCopy+l2BlockSizeAligned*(l1j-l2j),&(rhs[l1j*rhsStride+l2k]),(l2blockSizeEnd-l2k)*sizeof(Scalar));
+            std::memcpy(rhsCopy+l2BlockSizeAligned*(l1j-l2j),&(rhs[l1j*rhsStride+l2k]),(l2blockSizeEnd-l2k)*sizeof(Scalar));
           }
 
         // for each bw x 1 result's block
@@ -180,7 +180,7 @@ static void ei_cache_friendly_product(
         {
           int offsetblock = l2k * (l2blockRowEnd-l2i) + (l1i-l2i)*(l2blockSizeEnd-l2k) - l2k*MaxBlockRows;
           const Scalar* EIGEN_RESTRICT localB = &block[offsetblock];
-          
+
           for(int l1j=l2j; l1j<l2blockColEnd; l1j+=1)
           {
             const Scalar* EIGEN_RESTRICT rhsColumn;
@@ -352,7 +352,7 @@ static void ei_cache_friendly_product(
  * TODO: since rhs gets evaluated only once, no need to evaluate it
  */
 template<typename Scalar, typename RhsType>
-static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
+EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
   int size,
   const Scalar* lhs, int lhsStride,
   const RhsType& rhs,
@@ -397,7 +397,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
   int skipColumns = 0;
   if (PacketSize>1)
   {
-    ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);
+    ei_internal_assert(std::size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);
 
     while (skipColumns<PacketSize &&
            alignedStart != ((lhsAlignmentOffset + alignmentStep*skipColumns)%PacketSize))
@@ -414,7 +414,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
       // note that the skiped columns are processed later.
     }
 
-    ei_internal_assert((alignmentPattern==NoneAligned) || (size_t(lhs+alignedStart+lhsStride*skipColumns)%sizeof(Packet))==0);
+    ei_internal_assert((alignmentPattern==NoneAligned) || (std::size_t(lhs+alignedStart+lhsStride*skipColumns)%sizeof(Packet))==0);
   }
 
   int offset1 = (FirstAligned && alignmentStep==1?3:1);
@@ -516,7 +516,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
           res[j] += ei_pfirst(ptmp0) * lhs0[j];
 
         // process aligned result's coeffs
-        if ((size_t(lhs0+alignedStart)%sizeof(Packet))==0)
+        if ((std::size_t(lhs0+alignedStart)%sizeof(Packet))==0)
           for (int j = alignedStart;j<alignedSize;j+=PacketSize)
             ei_pstore(&res[j], ei_pmadd(ptmp0,ei_pload(&lhs0[j]),ei_pload(&res[j])));
         else
@@ -542,7 +542,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
 
 // TODO add peeling to mask unaligned load/stores
 template<typename Scalar, typename ResType>
-static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
+EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
   const Scalar* lhs, int lhsStride,
   const Scalar* rhs, int rhsSize,
   ResType& res)
@@ -586,7 +586,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
   int skipRows = 0;
   if (PacketSize>1)
   {
-    ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0  || size<PacketSize);
+    ei_internal_assert(std::size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0  || size<PacketSize);
 
     while (skipRows<PacketSize &&
            alignedStart != ((lhsAlignmentOffset + alignmentStep*skipRows)%PacketSize))
@@ -603,7 +603,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
       // note that the skiped columns are processed later.
     }
     ei_internal_assert((alignmentPattern==NoneAligned) || PacketSize==1
-      || (size_t(lhs+alignedStart+lhsStride*skipRows)%sizeof(Packet))==0);
+      || (std::size_t(lhs+alignedStart+lhsStride*skipRows)%sizeof(Packet))==0);
   }
 
   int offset1 = (FirstAligned && alignmentStep==1?3:1);
@@ -722,7 +722,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
       if (alignedSize>alignedStart)
       {
         // process aligned rhs coeffs
-        if ((size_t(lhs0+alignedStart)%sizeof(Packet))==0)
+        if ((std::size_t(lhs0+alignedStart)%sizeof(Packet))==0)
           for (int j = alignedStart;j<alignedSize;j+=PacketSize)
             ptmp0 = ei_pmadd(ei_pload(&rhs[j]), ei_pload(&lhs0[j]), ptmp0);
         else

Eigen/src/Core/CommaInitializer.h

@@ -116,6 +116,9 @@ struct CommaInitializer
   int m_row;              // current row id
   int m_col;              // current col id
   int m_currentBlockRows; // current block height
+
+private:
+  CommaInitializer& operator=(const CommaInitializer&);
 };
 
 /** \anchor MatrixBaseCommaInitRef

Eigen/src/Core/CoreInstantiations.cpp

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

Eigen/src/Core/Cwise.h

@@ -178,6 +178,9 @@ template<typename ExpressionType> class Cwise
 
   protected:
     ExpressionTypeNested m_matrix;
+
+  private:
+    Cwise& operator=(const Cwise&);
 };
 
 /** \returns a Cwise wrapper of *this providing additional coefficient-wise operations

Eigen/src/Core/DiagonalCoeffs.h

@@ -47,11 +47,11 @@ struct ei_traits<DiagonalCoeffs<MatrixType> >
   typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
   enum {
     RowsAtCompileTime = int(MatrixType::SizeAtCompileTime) == Dynamic ? Dynamic
-                      : EIGEN_ENUM_MIN(MatrixType::RowsAtCompileTime,
+                      : EIGEN_SIZE_MIN(MatrixType::RowsAtCompileTime,
                                        MatrixType::ColsAtCompileTime),
     ColsAtCompileTime = 1,
     MaxRowsAtCompileTime = int(MatrixType::MaxSizeAtCompileTime) == Dynamic ? Dynamic
-                            : EIGEN_ENUM_MIN(MatrixType::MaxRowsAtCompileTime,
+                            : EIGEN_SIZE_MIN(MatrixType::MaxRowsAtCompileTime,
                                              MatrixType::MaxColsAtCompileTime),
     MaxColsAtCompileTime = 1,
     Flags = (unsigned int)_MatrixTypeNested::Flags & (HereditaryBits | LinearAccessBit),

Eigen/src/Core/Flagged.h

@@ -109,6 +109,9 @@ template<typename ExpressionType, unsigned int Added, unsigned int Removed> clas
 
   protected:
     ExpressionTypeNested m_matrix;
+
+private:
+  Flagged& operator=(const Flagged&);
 };
 
 /** \returns an expression of *this with added flags

Eigen/src/Core/Functors.h

@@ -279,6 +279,8 @@ struct ei_scalar_multiple_op {
   EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a) const
   { return ei_pmul(a, ei_pset1(m_other)); }
   const Scalar m_other;
+private:
+  ei_scalar_multiple_op& operator=(const ei_scalar_multiple_op&);
 };
 template<typename Scalar>
 struct ei_functor_traits<ei_scalar_multiple_op<Scalar> >
@@ -294,6 +296,8 @@ struct ei_scalar_quotient1_impl {
   EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a) const
   { return ei_pmul(a, ei_pset1(m_other)); }
   const Scalar m_other;
+private:
+  ei_scalar_quotient1_impl& operator=(const ei_scalar_quotient1_impl&);
 };
 template<typename Scalar>
 struct ei_functor_traits<ei_scalar_quotient1_impl<Scalar,true> >
@@ -306,6 +310,8 @@ struct ei_scalar_quotient1_impl<Scalar,false> {
   EIGEN_STRONG_INLINE ei_scalar_quotient1_impl(const Scalar& other) : m_other(other) {}
   EIGEN_STRONG_INLINE Scalar operator() (const Scalar& a) const { return a / m_other; }
   const Scalar m_other;
+private:
+  ei_scalar_quotient1_impl& operator=(const ei_scalar_quotient1_impl&);
 };
 template<typename Scalar>
 struct ei_functor_traits<ei_scalar_quotient1_impl<Scalar,false> >
@@ -323,6 +329,8 @@ template<typename Scalar>
 struct ei_scalar_quotient1_op : ei_scalar_quotient1_impl<Scalar, NumTraits<Scalar>::HasFloatingPoint > {
   EIGEN_STRONG_INLINE ei_scalar_quotient1_op(const Scalar& other)
     : ei_scalar_quotient1_impl<Scalar, NumTraits<Scalar>::HasFloatingPoint >(other) {}
+private:
+  ei_scalar_quotient1_op& operator=(const ei_scalar_quotient1_op&);
 };
 
 // nullary functors
@@ -335,6 +343,8 @@ struct ei_scalar_constant_op {
   EIGEN_STRONG_INLINE const Scalar operator() (int, int = 0) const { return m_other; }
   EIGEN_STRONG_INLINE const PacketScalar packetOp() const { return ei_pset1(m_other); }
   const Scalar m_other;
+private:
+  ei_scalar_constant_op& operator=(const ei_scalar_constant_op&);
 };
 template<typename Scalar>
 struct ei_functor_traits<ei_scalar_constant_op<Scalar> >

Eigen/src/Core/MapBase.h

@@ -66,11 +66,11 @@ template<typename Derived> class MapBase
     inline const Scalar* data() const { return m_data; }
 
     template<bool IsForceAligned,typename Dummy> struct force_aligned_impl {
-      AlignedDerivedType static run(MapBase& a) { return a.derived(); }
+      static AlignedDerivedType run(MapBase& a) { return a.derived(); }
     };
 
     template<typename Dummy> struct force_aligned_impl<false,Dummy> {
-      AlignedDerivedType static run(MapBase& a) { return a.derived()._convertToForceAligned(); }
+      static AlignedDerivedType run(MapBase& a) { return a.derived()._convertToForceAligned(); }
     };
 
     /** \returns an expression equivalent to \c *this but having the \c PacketAccess constant
@@ -99,7 +99,7 @@ template<typename Derived> class MapBase
     inline const Scalar coeff(int index) const
     {
       ei_assert(Derived::IsVectorAtCompileTime || (ei_traits<Derived>::Flags & LinearAccessBit));
-      if ( ((RowsAtCompileTime == 1) == IsRowMajor) )
+      if ( ((RowsAtCompileTime == 1) == IsRowMajor) || !int(Derived::IsVectorAtCompileTime) )
         return m_data[index];
       else
         return m_data[index*stride()];
@@ -108,7 +108,7 @@ template<typename Derived> class MapBase
     inline Scalar& coeffRef(int index)
     {
       ei_assert(Derived::IsVectorAtCompileTime || (ei_traits<Derived>::Flags & LinearAccessBit));
-      if ( ((RowsAtCompileTime == 1) == IsRowMajor) )
+      if ( ((RowsAtCompileTime == 1) == IsRowMajor)  || !int(Derived::IsVectorAtCompileTime) )
         return const_cast<Scalar*>(m_data)[index];
       else
         return const_cast<Scalar*>(m_data)[index*stride()];

Eigen/src/Core/MathFunctions.h

@@ -35,16 +35,6 @@ template<typename T> inline T ei_random_amplitude()
   else return static_cast<T>(10);
 }
 
-template<typename T> inline T ei_hypot(T x, T y)
-{
-  T _x = ei_abs(x);
-  T _y = ei_abs(y);
-  T p = std::max(_x, _y);
-  T q = std::min(_x, _y);
-  T qp = q/p;
-  return p * ei_sqrt(T(1) + qp*qp);
-}
-
 /**************
 ***   int   ***
 **************/
@@ -54,24 +44,20 @@ template<> inline int machine_epsilon<int>() { return 0; }
 inline int ei_real(int x)  { return x; }
 inline int ei_imag(int)    { return 0; }
 inline int ei_conj(int x)  { return x; }
-inline int ei_abs(int x)   { return abs(x); }
+inline int ei_abs(int x)   { return std::abs(x); }
 inline int ei_abs2(int x)  { return x*x; }
 inline int ei_sqrt(int)  { ei_assert(false); return 0; }
 inline int ei_exp(int)  { ei_assert(false); return 0; }
 inline int ei_log(int)  { ei_assert(false); return 0; }
 inline int ei_sin(int)  { ei_assert(false); return 0; }
 inline int ei_cos(int)  { ei_assert(false); return 0; }
-
-#if EIGEN_GNUC_AT_LEAST(4,3)
-inline int ei_pow(int x, int y) { return int(std::pow(x, y)); }
-#else
+inline int ei_atan2(int, int)  { ei_assert(false); return 0; }
 inline int ei_pow(int x, int y) { return int(std::pow(double(x), y)); }
-#endif
 
 template<> inline int ei_random(int a, int b)
 {
   // We can't just do rand()%n as only the high-order bits are really random
-  return a + static_cast<int>((b-a+1) * (rand() / (RAND_MAX + 1.0)));
+  return a + static_cast<int>((b-a+1) * (std::rand() / (RAND_MAX + 1.0)));
 }
 template<> inline int ei_random()
 {
@@ -106,6 +92,7 @@ inline float ei_exp(float x)   { return std::exp(x); }
 inline float ei_log(float x)   { return std::log(x); }
 inline float ei_sin(float x)   { return std::sin(x); }
 inline float ei_cos(float x)   { return std::cos(x); }
+inline float ei_atan2(float y, float x) { return std::atan2(y,x); }
 inline float ei_pow(float x, float y)  { return std::pow(x, y); }
 
 template<> inline float ei_random(float a, float b)
@@ -153,6 +140,7 @@ inline double ei_exp(double x)   { return std::exp(x); }
 inline double ei_log(double x)   { return std::log(x); }
 inline double ei_sin(double x)   { return std::sin(x); }
 inline double ei_cos(double x)   { return std::cos(x); }
+inline double ei_atan2(double y, double x) { return std::atan2(y,x); }
 inline double ei_pow(double x, double y) { return std::pow(x, y); }
 
 template<> inline double ei_random(double a, double b)
@@ -197,6 +185,7 @@ inline float ei_abs2(const std::complex<float>& x) { return std::norm(x); }
 inline std::complex<float> ei_exp(std::complex<float> x)  { return std::exp(x); }
 inline std::complex<float> ei_sin(std::complex<float> x)  { return std::sin(x); }
 inline std::complex<float> ei_cos(std::complex<float> x)  { return std::cos(x); }
+inline std::complex<float> ei_atan2(std::complex<float>, std::complex<float> )  { ei_assert(false); return 0; }
 
 template<> inline std::complex<float> ei_random()
 {
@@ -231,6 +220,7 @@ inline double ei_abs2(const std::complex<double>& x) { return std::norm(x); }
 inline std::complex<double> ei_exp(std::complex<double> x)  { return std::exp(x); }
 inline std::complex<double> ei_sin(std::complex<double> x)  { return std::sin(x); }
 inline std::complex<double> ei_cos(std::complex<double> x)  { return std::cos(x); }
+inline std::complex<double> ei_atan2(std::complex<double>, std::complex<double>)  { ei_assert(false); return 0; }
 
 template<> inline std::complex<double> ei_random()
 {
@@ -268,6 +258,7 @@ inline long double ei_exp(long double x)   { return std::exp(x); }
 inline long double ei_log(long double x)   { return std::log(x); }
 inline long double ei_sin(long double x)   { return std::sin(x); }
 inline long double ei_cos(long double x)   { return std::cos(x); }
+inline long double ei_atan2(long double y, long double x) { return std::atan2(y,x); }
 inline long double ei_pow(long double x, long double y)  { return std::pow(x, y); }
 
 template<> inline long double ei_random(long double a, long double b)
@@ -291,4 +282,14 @@ inline bool ei_isApproxOrLessThan(long double a, long double b, long double prec
   return a <= b || ei_isApprox(a, b, prec);
 }
 
+template<typename T> inline T ei_hypot(T x, T y)
+{
+  T _x = ei_abs(x);
+  T _y = ei_abs(y);
+  T p = std::max(_x, _y);
+  T q = std::min(_x, _y);
+  T qp = q/p;
+  return p * ei_sqrt(T(1) + qp*qp);
+}
+
 #endif // EIGEN_MATHFUNCTIONS_H

Eigen/src/Core/Matrix.h

@@ -25,6 +25,11 @@
 #ifndef EIGEN_MATRIX_H
 #define EIGEN_MATRIX_H
 
+#ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO
+# define EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED for(int i=0;i<base().size();++i) coeffRef(i)=Scalar(0);
+#else
+# define EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
+#endif
 
 /** \class Matrix
   *
@@ -137,6 +142,9 @@ class Matrix
     enum { NeedsToAlign = (Options&AutoAlign) == AutoAlign
                           && SizeAtCompileTime!=Dynamic && ((sizeof(Scalar)*SizeAtCompileTime)%16)==0 };
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign)
+    
+    Base& base() { return *static_cast<Base*>(this); }
+    const Base& base() const { return *static_cast<const Base*>(this); }
 
     EIGEN_STRONG_INLINE int rows() const { return m_storage.rows(); }
     EIGEN_STRONG_INLINE int cols() const { return m_storage.cols(); }
@@ -230,7 +238,14 @@ class Matrix
              && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows)
              && (MaxColsAtCompileTime == Dynamic || MaxColsAtCompileTime >= cols)
              && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols));
-      m_storage.resize(rows * cols, rows, cols);
+      #ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO
+        int size = rows*cols;
+        bool size_changed = size != this->size();
+        m_storage.resize(size, rows, cols);
+        if(size_changed) EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
+      #else
+        m_storage.resize(rows*cols, rows, cols);
+      #endif
     }
 
     /** Resizes \c *this to a vector of length \a size
@@ -240,10 +255,17 @@ class Matrix
     inline void resize(int size)
     {
       EIGEN_STATIC_ASSERT_VECTOR_ONLY(Matrix)
+      ei_assert(SizeAtCompileTime == Dynamic || SizeAtCompileTime == size);
+      #ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO
+        bool size_changed = size != this->size();
+      #endif
       if(RowsAtCompileTime == 1)
         m_storage.resize(size, 1, size);
       else
         m_storage.resize(size, size, 1);
+      #ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO
+        if(size_changed) EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
+      #endif
     }
 
     /** Copies the value of the expression \a other into \c *this with automatic resizing.
@@ -287,13 +309,14 @@ class Matrix
     EIGEN_STRONG_INLINE explicit Matrix() : m_storage()
     {
       _check_template_params();
+      EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
     }
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN
     /** \internal */
     Matrix(ei_constructor_without_unaligned_array_assert)
       : m_storage(ei_constructor_without_unaligned_array_assert())
-    {}
+    {EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED}
 #endif
 
     /** Constructs a vector or row-vector with given dimension. \only_for_vectors
@@ -309,6 +332,7 @@ class Matrix
       EIGEN_STATIC_ASSERT_VECTOR_ONLY(Matrix)
       ei_assert(dim > 0);
       ei_assert(SizeAtCompileTime == Dynamic || SizeAtCompileTime == dim);
+      EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
     }
 
     /** This constructor has two very different behaviors, depending on the type of *this.
@@ -334,6 +358,7 @@ class Matrix
       {
         ei_assert(x > 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == x)
                && y > 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == y));
+        EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED
       }
     }
     /** constructs an initialized 2D vector with given coefficients */
@@ -395,18 +420,14 @@ class Matrix
     /** Override MatrixBase::swap() since for dynamic-sized matrices of same type it is enough to swap the
       * data pointers.
       */
-    inline void swap(Matrix& other)
-    {
-      if (Base::SizeAtCompileTime==Dynamic)
-        m_storage.swap(other.m_storage);
-      else
-        this->Base::swap(other);
-    }
+    template<typename OtherDerived>
+    void swap(const MatrixBase<OtherDerived>& other);
 
     /** \name Map
-      * These are convenience functions returning Map objects. The Map() static functions return unaligned Map objects,
-      * while the AlignedMap() functions return aligned Map objects and thus should be called only with 16-byte-aligned
-      * \a data pointers.
+      * These are convenience functions returning Map objects.
+      *
+      * \warning Do not use MapAligned in the Eigen 2.0. Mapping aligned arrays will be fully
+      * supported in Eigen 3.0 (already implemented in the development branch)
       *
       * \see class Map
       */
@@ -507,10 +528,18 @@ class Matrix
     template<typename OtherDerived>
     EIGEN_STRONG_INLINE Matrix& _set(const MatrixBase<OtherDerived>& other)
     {
-      _resize_to_match(other);
-      return Base::operator=(other);
+      // this enum introduced to fix compilation with gcc 3.3
+      enum { cond = int(OtherDerived::Flags) & EvalBeforeAssigningBit };
+      _set_selector(other.derived(), typename ei_meta_if<bool(cond), ei_meta_true, ei_meta_false>::ret());
+      return *this;
     }
 
+    template<typename OtherDerived>
+    EIGEN_STRONG_INLINE void _set_selector(const OtherDerived& other, const ei_meta_true&) { _set_noalias(other.eval()); }
+
+    template<typename OtherDerived>
+    EIGEN_STRONG_INLINE void _set_selector(const OtherDerived& other, const ei_meta_false&) { _set_noalias(other); }
+
     /** \internal Like _set() but additionally makes the assumption that no aliasing effect can happen (which
       * is the case when creating a new matrix) so one can enforce lazy evaluation.
       *
@@ -534,8 +563,40 @@ class Matrix
                         && (_Options & (AutoAlign|RowMajor)) == _Options),
           INVALID_MATRIX_TEMPLATE_PARAMETERS)
     }
+    
+    template<typename MatrixType, typename OtherDerived, bool IsSameType, bool IsDynamicSize>
+    friend struct ei_matrix_swap_impl;
 };
 
+template<typename MatrixType, typename OtherDerived,
+         bool IsSameType = ei_is_same_type<MatrixType, OtherDerived>::ret,
+         bool IsDynamicSize = MatrixType::SizeAtCompileTime==Dynamic>
+struct ei_matrix_swap_impl
+{
+  static inline void run(MatrixType& matrix, MatrixBase<OtherDerived>& other)
+  {
+    matrix.base().swap(other);
+  }
+};
+
+template<typename MatrixType, typename OtherDerived>
+struct ei_matrix_swap_impl<MatrixType, OtherDerived, true, true>
+{
+  static inline void run(MatrixType& matrix, MatrixBase<OtherDerived>& other)
+  {
+    matrix.m_storage.swap(other.derived().m_storage);
+  }
+};
+
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+template<typename OtherDerived>
+inline void Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::swap(const MatrixBase<OtherDerived>& other)
+{
+  // the Eigen:: here is to work around a stupid ICC 11.1 bug.
+  Eigen::ei_matrix_swap_impl<Matrix, OtherDerived>::run(*this, *const_cast<MatrixBase<OtherDerived>*>(&other));
+}
+
+
 /** \defgroup matrixtypedefs Global matrix typedefs
   *
   * \ingroup Core_Module

Eigen/src/Core/MatrixBase.h

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

Eigen/src/Core/NestByValue.h

@@ -100,6 +100,9 @@ template<typename ExpressionType> class NestByValue
 
   protected:
     const ExpressionType m_expression;
+
+  private:
+    NestByValue& operator=(const NestByValue&);
 };
 
 /** \returns an expression of the temporary version of *this.

Eigen/src/Core/NumTraits.h

@@ -94,7 +94,7 @@ template<typename _Real> struct NumTraits<std::complex<_Real> >
   enum {
     IsComplex = 1,
     HasFloatingPoint = NumTraits<Real>::HasFloatingPoint,
-    ReadCost = 2,
+    ReadCost = 2 * NumTraits<_Real>::ReadCost,
     AddCost = 2 * NumTraits<Real>::AddCost,
     MulCost = 4 * NumTraits<Real>::MulCost + 2 * NumTraits<Real>::AddCost
   };

Eigen/src/Core/Part.h

@@ -50,7 +50,7 @@ struct ei_traits<Part<MatrixType, Mode> > : ei_traits<MatrixType>
   typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
   enum {
     Flags = (_MatrixTypeNested::Flags & (HereditaryBits) & (~(PacketAccessBit | DirectAccessBit | LinearAccessBit))) | Mode,
-    CoeffReadCost = _MatrixTypeNested::CoeffReadCost
+    CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ConditionalJumpCost
   };
 };
 
@@ -124,8 +124,10 @@ template<typename MatrixType, unsigned int Mode> class Part
     }
 
   protected:
-
     const typename MatrixType::Nested m_matrix;
+
+  private:
+    Part& operator=(const Part&);
 };
 
 /** \nonstableyet

Eigen/src/Core/Product.h

@@ -66,11 +66,8 @@ struct ProductReturnType
 template<typename Lhs, typename Rhs>
 struct ProductReturnType<Lhs,Rhs,CacheFriendlyProduct>
 {
-  typedef typename ei_nested<Lhs,Rhs::ColsAtCompileTime>::type LhsNested;
-
-  typedef typename ei_nested<Rhs,Lhs::RowsAtCompileTime,
-                             typename ei_plain_matrix_type_column_major<Rhs>::type
-                   >::type RhsNested;
+  typedef const Lhs& LhsNested;
+  typedef const Rhs& RhsNested;
 
   typedef Product<LhsNested, RhsNested, CacheFriendlyProduct> Type;
 };
@@ -128,7 +125,7 @@ struct ei_traits<Product<LhsNested, RhsNested, ProductMode> >
 
     RowsAtCompileTime = _LhsNested::RowsAtCompileTime,
     ColsAtCompileTime = _RhsNested::ColsAtCompileTime,
-    InnerSize = EIGEN_ENUM_MIN(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
+    InnerSize = EIGEN_SIZE_MIN(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
 
     MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
     MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,
@@ -144,7 +141,7 @@ struct ei_traits<Product<LhsNested, RhsNested, ProductMode> >
 
     EvalToRowMajor = RhsRowMajor && (ProductMode==(int)CacheFriendlyProduct ? LhsRowMajor : (!CanVectorizeLhs)),
 
-    RemovedBits = ~(EvalToRowMajor ? 0 : RowMajorBit),
+    RemovedBits = ~((EvalToRowMajor ? 0 : RowMajorBit)|DirectAccessBit),
 
     Flags = ((unsigned int)(LhsFlags | RhsFlags) & HereditaryBits & RemovedBits)
           | EvalBeforeAssigningBit
@@ -190,6 +187,17 @@ template<typename LhsNested, typename RhsNested, int ProductMode> class Product
                                   _LhsNested, _RhsNested, Scalar> ScalarCoeffImpl;
 
   public:
+    typedef typename ei_nested<_LhsNested,_RhsNested::ColsAtCompileTime>::type LhsNestedX;
+    typedef typename ei_nested<_RhsNested,_LhsNested::RowsAtCompileTime>::type RhsNestedX;
+    typedef typename ei_cleantype<LhsNestedX>::type _LhsNestedX;
+    typedef typename ei_cleantype<RhsNestedX>::type _RhsNestedX;
+  
+    
+    enum {
+      LhsNestedFlags = _LhsNestedX::Flags,
+      RhsNestedFlags = _RhsNestedX::Flags
+    };
+    
 
     template<typename Lhs, typename Rhs>
     inline Product(const Lhs& lhs, const Rhs& rhs)
@@ -299,7 +307,7 @@ template<typename OtherDerived>
 inline Derived &
 MatrixBase<Derived>::operator*=(const MatrixBase<OtherDerived> &other)
 {
-  return *this = *this * other;
+  return derived() = derived() * other.derived();
 }
 
 /***************************************************************************
@@ -525,11 +533,11 @@ static void ei_cache_friendly_product_rowmajor_times_vector(
 
 template<typename ProductType,
   int LhsRows  = ei_traits<ProductType>::RowsAtCompileTime,
-  int LhsOrder = int(ei_traits<ProductType>::LhsFlags)&RowMajorBit ? RowMajor : ColMajor,
-  int LhsHasDirectAccess = int(ei_traits<ProductType>::LhsFlags)&DirectAccessBit? HasDirectAccess : NoDirectAccess,
-  int RhsCols  = ei_traits<ProductType>::ColsAtCompileTime,
-  int RhsOrder = int(ei_traits<ProductType>::RhsFlags)&RowMajorBit ? RowMajor : ColMajor,
-  int RhsHasDirectAccess = int(ei_traits<ProductType>::RhsFlags)&DirectAccessBit? HasDirectAccess : NoDirectAccess>
+  int LhsOrder = int(ProductType::LhsNestedFlags)&RowMajorBit ? RowMajor : ColMajor,
+  int LhsHasDirectAccess = int(ProductType::LhsNestedFlags)&DirectAccessBit? HasDirectAccess : NoDirectAccess,
+  int RhsCols  = ProductType::ColsAtCompileTime,
+  int RhsOrder = int(ProductType::RhsNestedFlags)&RowMajorBit ? RowMajor : ColMajor,
+  int RhsHasDirectAccess = int(ProductType::RhsNestedFlags)&DirectAccessBit? HasDirectAccess : NoDirectAccess>
 struct ei_cache_friendly_product_selector
 {
   template<typename DestDerived>
@@ -546,9 +554,12 @@ struct ei_cache_friendly_product_selector<ProductType,LhsRows,ColMajor,NoDirectA
   template<typename DestDerived>
   inline static void run(DestDerived& res, const ProductType& product)
   {
-    const int size = product.rhs().rows();
+    typename ProductType::LhsNestedX lhs(product.lhs());
+    typename ProductType::RhsNestedX rhs(product.rhs());
+    
+    const int size = rhs.rows();
     for (int k=0; k<size; ++k)
-        res += product.rhs().coeff(k) * product.lhs().col(k);
+        res += rhs.coeff(k) * lhs.col(k);
   }
 };
 
@@ -562,6 +573,9 @@ struct ei_cache_friendly_product_selector<ProductType,LhsRows,ColMajor,HasDirect
   template<typename DestDerived>
   inline static void run(DestDerived& res, const ProductType& product)
   {
+    typename ProductType::LhsNestedX lhs(product.lhs());
+    typename ProductType::RhsNestedX rhs(product.rhs());
+    
     enum {
       EvalToRes = (ei_packet_traits<Scalar>::size==1)
                 ||((DestDerived::Flags&ActualPacketAccessBit) && (!(DestDerived::Flags & RowMajorBit))) };
@@ -571,15 +585,15 @@ struct ei_cache_friendly_product_selector<ProductType,LhsRows,ColMajor,HasDirect
     else
     {
       _res = ei_aligned_stack_new(Scalar,res.size());
-      Map<Matrix<Scalar,DestDerived::RowsAtCompileTime,1> >(_res, res.size()) = res;
+      Map<Matrix<Scalar,DestDerived::RowsAtCompileTime,1,ColMajor> >(_res, res.size()) = res;
     }
     ei_cache_friendly_product_colmajor_times_vector(res.size(),
-      &product.lhs().const_cast_derived().coeffRef(0,0), product.lhs().stride(),
-      product.rhs(), _res);
+      &lhs.const_cast_derived().coeffRef(0,0), lhs.stride(),
+      rhs, _res);
 
     if (!EvalToRes)
     {
-      res = Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1> >(_res, res.size());
+      res = Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1,ColMajor> >(_res, res.size());
       ei_aligned_stack_delete(Scalar, _res, res.size());
     }
   }
@@ -592,9 +606,12 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
   template<typename DestDerived>
   inline static void run(DestDerived& res, const ProductType& product)
   {
-    const int cols = product.lhs().cols();
+    typename ProductType::LhsNestedX lhs(product.lhs());
+    typename ProductType::RhsNestedX rhs(product.rhs());
+    
+    const int cols = lhs.cols();
     for (int j=0; j<cols; ++j)
-      res += product.lhs().coeff(j) * product.rhs().row(j);
+      res += lhs.coeff(j) * rhs.row(j);
   }
 };
 
@@ -608,6 +625,9 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
   template<typename DestDerived>
   inline static void run(DestDerived& res, const ProductType& product)
   {
+    typename ProductType::LhsNestedX lhs(product.lhs());
+    typename ProductType::RhsNestedX rhs(product.rhs());
+    
     enum {
       EvalToRes = (ei_packet_traits<Scalar>::size==1)
                 ||((DestDerived::Flags & ActualPacketAccessBit) && (DestDerived::Flags & RowMajorBit)) };
@@ -617,15 +637,15 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
     else
     {
       _res = ei_aligned_stack_new(Scalar, res.size());
-      Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1> >(_res, res.size()) = res;
+      Map<Matrix<Scalar,1,DestDerived::SizeAtCompileTime,ColMajor> >(_res, res.size()) = res;
     }
     ei_cache_friendly_product_colmajor_times_vector(res.size(),
-      &product.rhs().const_cast_derived().coeffRef(0,0), product.rhs().stride(),
-      product.lhs().transpose(), _res);
+      &rhs.const_cast_derived().coeffRef(0,0), rhs.stride(),
+      lhs.transpose(), _res);
 
     if (!EvalToRes)
     {
-      res = Map<Matrix<Scalar,DestDerived::SizeAtCompileTime,1> >(_res, res.size());
+      res = Map<Matrix<Scalar,1,DestDerived::SizeAtCompileTime,ColMajor> >(_res, res.size());
       ei_aligned_stack_delete(Scalar, _res, res.size());
     }
   }
@@ -638,24 +658,28 @@ struct ei_cache_friendly_product_selector<ProductType,LhsRows,RowMajor,HasDirect
   typedef typename ProductType::Scalar Scalar;
   typedef typename ei_traits<ProductType>::_RhsNested Rhs;
   enum {
-      UseRhsDirectly = ((ei_packet_traits<Scalar>::size==1) || (Rhs::Flags&ActualPacketAccessBit))
-                     && (!(Rhs::Flags & RowMajorBit)) };
+      UseRhsDirectly = ((ei_packet_traits<Scalar>::size==1) || (ProductType::RhsNestedFlags&ActualPacketAccessBit))
+                     && (ProductType::RhsNestedFlags&DirectAccessBit)
+                     && (!(ProductType::RhsNestedFlags & RowMajorBit)) };
 
   template<typename DestDerived>
   inline static void run(DestDerived& res, const ProductType& product)
   {
+    typename ProductType::LhsNestedX lhs(product.lhs());
+    typename ProductType::RhsNestedX rhs(product.rhs());
+  
     Scalar* EIGEN_RESTRICT _rhs;
     if (UseRhsDirectly)
-       _rhs = &product.rhs().const_cast_derived().coeffRef(0);
+       _rhs = &rhs.const_cast_derived().coeffRef(0);
     else
     {
-      _rhs = ei_aligned_stack_new(Scalar, product.rhs().size());
-      Map<Matrix<Scalar,Rhs::SizeAtCompileTime,1> >(_rhs, product.rhs().size()) = product.rhs();
+      _rhs = ei_aligned_stack_new(Scalar, rhs.size());
+      Map<Matrix<Scalar,Rhs::SizeAtCompileTime,1,ColMajor> >(_rhs, rhs.size()) = rhs;
     }
-    ei_cache_friendly_product_rowmajor_times_vector(&product.lhs().const_cast_derived().coeffRef(0,0), product.lhs().stride(),
-                                                    _rhs, product.rhs().size(), res);
+    ei_cache_friendly_product_rowmajor_times_vector(&lhs.const_cast_derived().coeffRef(0,0), lhs.stride(),
+                                                    _rhs, rhs.size(), res);
 
-    if (!UseRhsDirectly) ei_aligned_stack_delete(Scalar, _rhs, product.rhs().size());
+    if (!UseRhsDirectly) ei_aligned_stack_delete(Scalar, _rhs, rhs.size());
   }
 };
 
@@ -666,24 +690,28 @@ struct ei_cache_friendly_product_selector<ProductType,1,LhsOrder,LhsAccess,RhsCo
   typedef typename ProductType::Scalar Scalar;
   typedef typename ei_traits<ProductType>::_LhsNested Lhs;
   enum {
-      UseLhsDirectly = ((ei_packet_traits<Scalar>::size==1) || (Lhs::Flags&ActualPacketAccessBit))
-                     && (Lhs::Flags & RowMajorBit) };
+      UseLhsDirectly = ((ei_packet_traits<Scalar>::size==1) || (ProductType::LhsNestedFlags&ActualPacketAccessBit))
+                     && (ProductType::LhsNestedFlags&DirectAccessBit)
+                     && (ProductType::LhsNestedFlags & RowMajorBit) };
 
   template<typename DestDerived>
   inline static void run(DestDerived& res, const ProductType& product)
   {
+    typename ProductType::LhsNestedX lhs(product.lhs());
+    typename ProductType::RhsNestedX rhs(product.rhs());
+  
     Scalar* EIGEN_RESTRICT _lhs;
     if (UseLhsDirectly)
-       _lhs = &product.lhs().const_cast_derived().coeffRef(0);
+       _lhs = &lhs.const_cast_derived().coeffRef(0);
     else
     {
-      _lhs = ei_aligned_stack_new(Scalar, product.lhs().size());
-      Map<Matrix<Scalar,Lhs::SizeAtCompileTime,1> >(_lhs, product.lhs().size()) = product.lhs();
+      _lhs = ei_aligned_stack_new(Scalar, lhs.size());
+      Map<Matrix<Scalar,1,Lhs::SizeAtCompileTime,ColMajor> >(_lhs, lhs.size()) = lhs;
     }
-    ei_cache_friendly_product_rowmajor_times_vector(&product.rhs().const_cast_derived().coeffRef(0,0), product.rhs().stride(),
-                                                    _lhs, product.lhs().size(), res);
+    ei_cache_friendly_product_rowmajor_times_vector(&rhs.const_cast_derived().coeffRef(0,0), rhs.stride(),
+                                                    _lhs, lhs.size(), res);
 
-    if(!UseLhsDirectly) ei_aligned_stack_delete(Scalar, _lhs, product.lhs().size());
+    if(!UseLhsDirectly) ei_aligned_stack_delete(Scalar, _lhs, lhs.size());
   }
 };
 
@@ -709,7 +737,17 @@ MatrixBase<Derived>::operator+=(const Flagged<Product<Lhs,Rhs,CacheFriendlyProdu
   if (other._expression()._useCacheFriendlyProduct())
     ei_cache_friendly_product_selector<Product<Lhs,Rhs,CacheFriendlyProduct> >::run(const_cast_derived(), other._expression());
   else
-    lazyAssign(derived() + other._expression());
+  {
+    typedef typename ei_cleantype<Lhs>::type _Lhs;
+    typedef typename ei_cleantype<Rhs>::type _Rhs;
+  
+    typedef typename ei_nested<_Lhs,_Rhs::ColsAtCompileTime>::type LhsNested;
+    typedef typename ei_nested<_Rhs,_Lhs::RowsAtCompileTime>::type RhsNested;
+
+    Product<LhsNested,RhsNested,NormalProduct> prod(other._expression().lhs(),other._expression().rhs());
+    
+    lazyAssign(derived() + prod);
+  }
   return derived();
 }
 
@@ -724,12 +762,21 @@ inline Derived& MatrixBase<Derived>::lazyAssign(const Product<Lhs,Rhs,CacheFrien
   }
   else
   {
-    lazyAssign<Product<Lhs,Rhs,CacheFriendlyProduct> >(product);
+    typedef typename ei_cleantype<Lhs>::type _Lhs;
+    typedef typename ei_cleantype<Rhs>::type _Rhs;
+    
+    typedef typename ei_nested<_Lhs,_Rhs::ColsAtCompileTime>::type LhsNested;
+    typedef typename ei_nested<_Rhs,_Lhs::RowsAtCompileTime>::type RhsNested;
+
+    typedef Product<LhsNested,RhsNested,NormalProduct> NormalProduct;
+    NormalProduct normal_prod(product.lhs(),product.rhs());
+
+    lazyAssign<NormalProduct>(normal_prod);
   }
   return derived();
 }
 
-template<typename T> struct ei_product_copy_rhs
+template<typename T,int StorageOrder> struct ei_product_copy_rhs
 {
   typedef typename ei_meta_if<
          (ei_traits<T>::Flags & RowMajorBit)
@@ -739,11 +786,30 @@ template<typename T> struct ei_product_copy_rhs
     >::ret type;
 };
 
-template<typename T> struct ei_product_copy_lhs
+template<typename T> struct ei_product_copy_rhs<T,RowMajorBit>
+{
+  typedef typename ei_meta_if<
+      (!(ei_traits<T>::Flags & DirectAccessBit)),
+      typename ei_plain_matrix_type<T>::type,
+      const T&
+    >::ret type;
+};
+
+template<typename T,int StorageOrder> struct ei_product_copy_lhs
 {
   typedef typename ei_meta_if<
       (!(int(ei_traits<T>::Flags) & DirectAccessBit)),
-      typename ei_plain_matrix_type<T>::type,
+      typename ei_plain_matrix_type_row_major<T>::type,
+      const T&
+    >::ret type;
+};
+
+template<typename T> struct ei_product_copy_lhs<T,RowMajorBit>
+{
+  typedef typename ei_meta_if<
+         ((ei_traits<T>::Flags & RowMajorBit)==0)
+      || (!(int(ei_traits<T>::Flags) & DirectAccessBit)),
+      typename ei_plain_matrix_type_row_major<T>::type,
       const T&
     >::ret type;
 };
@@ -752,9 +818,9 @@ template<typename Lhs, typename Rhs, int ProductMode>
 template<typename DestDerived>
 inline void Product<Lhs,Rhs,ProductMode>::_cacheFriendlyEvalAndAdd(DestDerived& res) const
 {
-  typedef typename ei_product_copy_lhs<_LhsNested>::type LhsCopy;
+  typedef typename ei_product_copy_lhs<_LhsNested,DestDerived::Flags&RowMajorBit>::type LhsCopy;
   typedef typename ei_unref<LhsCopy>::type _LhsCopy;
-  typedef typename ei_product_copy_rhs<_RhsNested>::type RhsCopy;
+  typedef typename ei_product_copy_rhs<_RhsNested,DestDerived::Flags&RowMajorBit>::type RhsCopy;
   typedef typename ei_unref<RhsCopy>::type _RhsCopy;
   LhsCopy lhs(m_lhs);
   RhsCopy rhs(m_rhs);
@@ -762,8 +828,9 @@ inline void Product<Lhs,Rhs,ProductMode>::_cacheFriendlyEvalAndAdd(DestDerived&
     rows(), cols(), lhs.cols(),
     _LhsCopy::Flags&RowMajorBit, (const Scalar*)&(lhs.const_cast_derived().coeffRef(0,0)), lhs.stride(),
     _RhsCopy::Flags&RowMajorBit, (const Scalar*)&(rhs.const_cast_derived().coeffRef(0,0)), rhs.stride(),
-    Flags&RowMajorBit, (Scalar*)&(res.coeffRef(0,0)), res.stride()
+    DestDerived::Flags&RowMajorBit, (Scalar*)&(res.coeffRef(0,0)), res.stride()
   );
+
 }
 
 #endif // EIGEN_PRODUCT_H

Eigen/src/Core/SolveTriangular.h

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

Eigen/src/Core/Swap.h

@@ -117,6 +117,9 @@ template<typename ExpressionType> class SwapWrapper
 
   protected:
     ExpressionType& m_expression;
+
+  private:
+    SwapWrapper& operator=(const SwapWrapper&);
 };
 
 /** swaps *this with the expression \a other.

Eigen/src/Core/Transpose.h

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

Eigen/src/Core/Visitor.h

@@ -164,7 +164,7 @@ struct ei_functor_traits<ei_max_coeff_visitor<Scalar> > {
 /** \returns the minimum of all coefficients of *this
   * and puts in *row and *col its location.
   *
-  * \sa MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff()
+  * \sa MatrixBase::minCoeff(int*), MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff()
   */
 template<typename Derived>
 typename ei_traits<Derived>::Scalar
@@ -177,6 +177,22 @@ MatrixBase<Derived>::minCoeff(int* row, int* col) const
   return minVisitor.res;
 }
 
+/** \returns the minimum of all coefficients of *this
+  * and puts in *index its location.
+  *
+  * \sa MatrixBase::minCoeff(int*,int*), MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff()
+  */
+template<typename Derived>
+typename ei_traits<Derived>::Scalar
+MatrixBase<Derived>::minCoeff(int* index) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  ei_min_coeff_visitor<Scalar> minVisitor;
+  this->visit(minVisitor);
+  *index = (RowsAtCompileTime==1) ? minVisitor.col : minVisitor.row;
+  return minVisitor.res;
+}
+
 /** \returns the maximum of all coefficients of *this
   * and puts in *row and *col its location.
   *
@@ -193,5 +209,20 @@ MatrixBase<Derived>::maxCoeff(int* row, int* col) const
   return maxVisitor.res;
 }
 
+/** \returns the maximum of all coefficients of *this
+  * and puts in *index its location.
+  *
+  * \sa MatrixBase::maxCoeff(int*,int*), MatrixBase::minCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::maxCoeff()
+  */
+template<typename Derived>
+typename ei_traits<Derived>::Scalar
+MatrixBase<Derived>::maxCoeff(int* index) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  ei_max_coeff_visitor<Scalar> maxVisitor;
+  this->visit(maxVisitor);
+  *index = (RowsAtCompileTime==1) ? maxVisitor.col : maxVisitor.row;
+  return maxVisitor.res;
+}
 
 #endif // EIGEN_VISITOR_H

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

@@ -114,9 +114,15 @@ template<> EIGEN_STRONG_INLINE void ei_pstoreu<float>(float*  to, const __m128&
 template<> EIGEN_STRONG_INLINE void ei_pstoreu<double>(double* to, const __m128d& from) { _mm_storeu_pd(to, from); }
 template<> EIGEN_STRONG_INLINE void ei_pstoreu<int>(int*    to, const __m128i& from) { _mm_storeu_si128(reinterpret_cast<__m128i*>(to), from); }
 
-#ifdef _MSC_VER
-// this fix internal compilation error
-template<> EIGEN_STRONG_INLINE float  ei_pfirst<__m128>(const __m128&  a) { float x = _mm_cvtss_f32(a); return x; }
+#if defined(_MSC_VER) && (_MSC_VER <= 1500) && defined(_WIN64) && !defined(__INTEL_COMPILER)
+// The temporary variable fixes an internal compilation error.
+// Direct of the struct members fixed bug #62.
+template<> EIGEN_STRONG_INLINE float  ei_pfirst<__m128>(const __m128& a) { return a.m128_f32[0]; }
+template<> EIGEN_STRONG_INLINE double ei_pfirst<__m128d>(const __m128d& a) { return a.m128d_f64[0]; }
+template<> EIGEN_STRONG_INLINE int    ei_pfirst<__m128i>(const __m128i& a) { int x = _mm_cvtsi128_si32(a); return x; }
+#elif defined(_MSC_VER) && (_MSC_VER <= 1500) && !defined(__INTEL_COMPILER)
+// The temporary variable fixes an internal compilation error.
+template<> EIGEN_STRONG_INLINE float  ei_pfirst<__m128>(const __m128& a) { float x = _mm_cvtss_f32(a); return x; }
 template<> EIGEN_STRONG_INLINE double ei_pfirst<__m128d>(const __m128d& a) { double x = _mm_cvtsd_f64(a); return x; }
 template<> EIGEN_STRONG_INLINE int    ei_pfirst<__m128i>(const __m128i& a) { int x = _mm_cvtsi128_si32(a); return x; }
 #else
@@ -315,4 +321,7 @@ struct ei_palign_impl<Offset,__m128d>
 };
 #endif
 
+#define ei_vec4f_swizzle1(v,p,q,r,s) \
+  (_mm_castsi128_ps(_mm_shuffle_epi32( _mm_castps_si128(v), ((s)<<6|(r)<<4|(q)<<2|(p)))))
+
 #endif // EIGEN_PACKET_MATH_SSE_H

Eigen/src/Core/util/Constants.h

@@ -239,4 +239,16 @@ enum {
   HasDirectAccess = DirectAccessBit
 };
 
+const int EiArch_Generic = 0x0;
+const int EiArch_SSE     = 0x1;
+const int EiArch_AltiVec = 0x2;
+
+#if defined EIGEN_VECTORIZE_SSE
+  const int EiArch = EiArch_SSE;
+#elif defined EIGEN_VECTORIZE_ALTIVEC
+  const int EiArch = EiArch_AltiVec;
+#else
+  const int EiArch = EiArch_Generic;
+#endif
+
 #endif // EIGEN_CONSTANTS_H

Eigen/src/Core/util/Macros.h

@@ -30,30 +30,48 @@
 
 #define EIGEN_WORLD_VERSION 2
 #define EIGEN_MAJOR_VERSION 0
-#define EIGEN_MINOR_VERSION 3
+#define EIGEN_MINOR_VERSION 16
 
 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
                                       (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \
                                                                  EIGEN_MINOR_VERSION>=z))))
 
-// if the compiler is GNUC, disable 16 byte alignment on exotic archs that probably don't need it, and on which
-// it may be extra trouble to get aligned memory allocation to work (example: on ARM, overloading new[] is a PITA
-// because extra memory must be allocated for bookkeeping).
-// if the compiler is not GNUC, just cross fingers that the architecture isn't too exotic, because we don't want
-// to keep track of all the different preprocessor symbols for all compilers.
-#if !defined(__GNUC__) || defined(__i386__) || defined(__x86_64__) || defined(__ppc__) || defined(__ia64__)
+// 16 byte alignment is only useful for vectorization. Since it affects the ABI, we need to enable 16 byte alignment on all
+// platforms where vectorization might be enabled. In theory we could always enable alignment, but it can be a cause of problems
+// on some platforms, so we just disable it in certain common platform (compiler+architecture combinations) to avoid these problems.
+#if defined(__GNUC__) && !(defined(__i386__) || defined(__x86_64__) || defined(__powerpc__) || defined(__ppc__) || defined(__ia64__))
+#define EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT 1
+#else
+#define EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT 0
+#endif
+
+#if defined(__GNUC__) && (__GNUC__ <= 3)
+#define EIGEN_GCC3_OR_OLDER 1
+#else
+#define EIGEN_GCC3_OR_OLDER 0
+#endif
+
+// FIXME vectorization + alignment is completely disabled with sun studio
+#if !EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT && !EIGEN_GCC3_OR_OLDER && !defined(__SUNPRO_CC) && !defined(__QNXNTO__)
   #define EIGEN_ARCH_WANTS_ALIGNMENT 1
 #else
-  #ifdef EIGEN_VECTORIZE
-    #error Vectorization enabled, but the architecture is not listed among those for which we require 16 byte alignment. If you added vectorization for another architecture, you also need to edit this list.
-  #endif
   #define EIGEN_ARCH_WANTS_ALIGNMENT 0
+#endif
+
+// EIGEN_ALIGN is the true test whether we want to align or not. It takes into account both the user choice to explicitly disable
+// alignment (EIGEN_DONT_ALIGN) and the architecture config (EIGEN_ARCH_WANTS_ALIGNMENT). Henceforth, only EIGEN_ALIGN should be used.
+#if EIGEN_ARCH_WANTS_ALIGNMENT && !defined(EIGEN_DONT_ALIGN)
+  #define EIGEN_ALIGN 1
+#else
+  #define EIGEN_ALIGN 0
+  #ifdef EIGEN_VECTORIZE
+    #error "Vectorization enabled, but our platform checks say that we don't do 16 byte alignment on this platform. If you added vectorization for another architecture, you also need to edit this platform check."
+  #endif
   #ifndef EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT
     #define EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT
   #endif
 #endif
 
-
 #ifdef EIGEN_DEFAULT_TO_ROW_MAJOR
 #define EIGEN_DEFAULT_MATRIX_STORAGE_ORDER_OPTION RowMajor
 #else
@@ -165,7 +183,7 @@ using Eigen::ei_cos;
  * If we made alignment depend on whether or not EIGEN_VECTORIZE is defined, it would be impossible to link
  * vectorized and non-vectorized code.
  */
-#if !EIGEN_ARCH_WANTS_ALIGNMENT
+#if !EIGEN_ALIGN
 #define EIGEN_ALIGN_128
 #elif (defined __GNUC__)
 #define EIGEN_ALIGN_128 __attribute__((aligned(16)))
@@ -175,10 +193,15 @@ using Eigen::ei_cos;
 #error Please tell me what is the equivalent of __attribute__((aligned(16))) for your compiler
 #endif
 
-#define EIGEN_RESTRICT __restrict
+#ifdef EIGEN_DONT_USE_RESTRICT_KEYWORD
+  #define EIGEN_RESTRICT
+#endif
+#ifndef EIGEN_RESTRICT
+  #define EIGEN_RESTRICT __restrict
+#endif
 
 #ifndef EIGEN_STACK_ALLOCATION_LIMIT
-#define EIGEN_STACK_ALLOCATION_LIMIT 16000000
+#define EIGEN_STACK_ALLOCATION_LIMIT 1000000
 #endif
 
 #ifndef EIGEN_DEFAULT_IO_FORMAT
@@ -234,6 +257,9 @@ enum { RowsAtCompileTime = Eigen::ei_traits<Derived>::RowsAtCompileTime, \
 _EIGEN_GENERIC_PUBLIC_INTERFACE(Derived, Eigen::MatrixBase<Derived>)
 
 #define EIGEN_ENUM_MIN(a,b) (((int)a <= (int)b) ? (int)a : (int)b)
+#define EIGEN_SIZE_MIN(a,b) (((int)a == 1 || (int)b == 1) ? 1 \
+                           : ((int)a == Dynamic || (int)b == Dynamic) ? Dynamic \
+                           : ((int)a <= (int)b) ? (int)a : (int)b)
 #define EIGEN_ENUM_MAX(a,b) (((int)a >= (int)b) ? (int)a : (int)b)
 
 // just an empty macro !

Eigen/src/Core/util/Memory.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>
 // Copyright (C) 2008-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
@@ -27,13 +27,23 @@
 #ifndef EIGEN_MEMORY_H
 #define EIGEN_MEMORY_H
 
-#if defined(__APPLE__) || defined(_WIN64)
+// FreeBSD 6 seems to have 16-byte aligned malloc
+// See http://svn.freebsd.org/viewvc/base/stable/6/lib/libc/stdlib/malloc.c?view=markup
+// FreeBSD 7 seems to have 16-byte aligned malloc except on ARM and MIPS architectures
+// See http://svn.freebsd.org/viewvc/base/stable/7/lib/libc/stdlib/malloc.c?view=markup
+#if defined(__FreeBSD__) && !defined(__arm__) && !defined(__mips__)
+#define EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED 1
+#else
+#define EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED 0
+#endif
+
+#if defined(__APPLE__) || defined(_WIN64) || EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED
   #define EIGEN_MALLOC_ALREADY_ALIGNED 1
 #else
   #define EIGEN_MALLOC_ALREADY_ALIGNED 0
 #endif
 
-#if ((defined _GNU_SOURCE) || ((defined _XOPEN_SOURCE) && (_XOPEN_SOURCE >= 600))) && (defined _POSIX_ADVISORY_INFO) && (_POSIX_ADVISORY_INFO > 0)
+#if (defined __QNXNTO__) || (((defined _GNU_SOURCE) || ((defined _XOPEN_SOURCE) && (_XOPEN_SOURCE >= 600))) && (defined _POSIX_ADVISORY_INFO) && (_POSIX_ADVISORY_INFO > 0))
   #define EIGEN_HAS_POSIX_MEMALIGN 1
 #else
   #define EIGEN_HAS_POSIX_MEMALIGN 0
@@ -49,10 +59,10 @@
   * Fast, but wastes 16 additional bytes of memory.
   * Does not throw any exception.
   */
-inline void* ei_handmade_aligned_malloc(size_t size)
+inline void* ei_handmade_aligned_malloc(std::size_t size)
 {
-  void *original = malloc(size+16);
-  void *aligned = reinterpret_cast<void*>((reinterpret_cast<size_t>(original) & ~(size_t(15))) + 16);
+  void *original = std::malloc(size+16);
+  void *aligned = reinterpret_cast<void*>((reinterpret_cast<std::size_t>(original) & ~(std::size_t(15))) + 16);
   *(reinterpret_cast<void**>(aligned) - 1) = original;
   return aligned;
 }
@@ -61,94 +71,96 @@ inline void* ei_handmade_aligned_malloc(size_t size)
 inline void ei_handmade_aligned_free(void *ptr)
 {
   if(ptr)
-    free(*(reinterpret_cast<void**>(ptr) - 1));
+    std::free(*(reinterpret_cast<void**>(ptr) - 1));
 }
 
 /** \internal allocates \a size bytes. The returned pointer is guaranteed to have 16 bytes alignment.
-  * On allocation error, the returned pointer is undefined, but if exceptions are enabled then a std::bad_alloc is thrown.
+  * On allocation error, the returned pointer is null, and if exceptions are enabled then a std::bad_alloc is thrown.
   */
-inline void* ei_aligned_malloc(size_t size)
+inline void* ei_aligned_malloc(std::size_t size)
 {
   #ifdef EIGEN_NO_MALLOC
     ei_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
   #endif
 
-  void *result;
-  #if EIGEN_HAS_POSIX_MEMALIGN && EIGEN_ARCH_WANTS_ALIGNMENT && !EIGEN_MALLOC_ALREADY_ALIGNED
-    #ifdef EIGEN_EXCEPTIONS
-      const int failed =
-    #endif
-    posix_memalign(&result, 16, size);
+  void *result;  
+  #if !EIGEN_ALIGN
+    result = std::malloc(size);
+  #elif EIGEN_MALLOC_ALREADY_ALIGNED
+    result = std::malloc(size);
+  #elif EIGEN_HAS_POSIX_MEMALIGN
+    if(posix_memalign(&result, 16, size)) result = 0;
+  #elif EIGEN_HAS_MM_MALLOC
+    result = _mm_malloc(size, 16);
+  #elif (defined _MSC_VER)
+    result = _aligned_malloc(size, 16);
   #else
-    #if !EIGEN_ARCH_WANTS_ALIGNMENT
-      result = malloc(size);
-    #elif EIGEN_MALLOC_ALREADY_ALIGNED
-      result = malloc(size);
-    #elif EIGEN_HAS_MM_MALLOC
-      result = _mm_malloc(size, 16);
-    #elif (defined _MSC_VER)
-      result = _aligned_malloc(size, 16);
-    #else
-      result = ei_handmade_aligned_malloc(size);
-    #endif
-    #ifdef EIGEN_EXCEPTIONS
-      const int failed = (result == 0);
-    #endif
+    result = ei_handmade_aligned_malloc(size);
   #endif
+    
   #ifdef EIGEN_EXCEPTIONS
-    if(failed)
+    if(result == 0)
       throw std::bad_alloc();
   #endif
   return result;
 }
 
 /** allocates \a size bytes. If Align is true, then the returned ptr is 16-byte-aligned.
-  * On allocation error, the returned pointer is undefined, but if exceptions are enabled then a std::bad_alloc is thrown.
+  * On allocation error, the returned pointer is null, and if exceptions are enabled then a std::bad_alloc is thrown.
   */
-template<bool Align> inline void* ei_conditional_aligned_malloc(size_t size)
+template<bool Align> inline void* ei_conditional_aligned_malloc(std::size_t size)
 {
   return ei_aligned_malloc(size);
 }
 
-template<> inline void* ei_conditional_aligned_malloc<false>(size_t size)
+template<> inline void* ei_conditional_aligned_malloc<false>(std::size_t size)
 {
   #ifdef EIGEN_NO_MALLOC
     ei_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
   #endif
 
-  void *result = malloc(size);
+  void *result = std::malloc(size);
   #ifdef EIGEN_EXCEPTIONS
     if(!result) throw std::bad_alloc();
   #endif
   return result;
 }
 
+/** \internal construct the elements of an array.
+  * The \a size parameter tells on how many objects to call the constructor of T.
+  */
+template<typename T> inline T* ei_construct_elements_of_array(T *ptr, std::size_t size)
+{
+  for (std::size_t i=0; i < size; ++i) ::new (ptr + i) T;
+  return ptr;
+}
+
 /** allocates \a size objects of type T. The returned pointer is guaranteed to have 16 bytes alignment.
   * On allocation error, the returned pointer is undefined, but if exceptions are enabled then a std::bad_alloc is thrown.
   * The default constructor of T is called.
   */
-template<typename T> inline T* ei_aligned_new(size_t size)
+template<typename T> inline T* ei_aligned_new(std::size_t size)
 {
-  void *void_result = ei_aligned_malloc(sizeof(T)*size);
-  return ::new(void_result) T[size];
+  T *result = reinterpret_cast<T*>(ei_aligned_malloc(sizeof(T)*size));
+  return ei_construct_elements_of_array(result, size);
 }
 
-template<typename T, bool Align> inline T* ei_conditional_aligned_new(size_t size)
+template<typename T, bool Align> inline T* ei_conditional_aligned_new(std::size_t size)
 {
-  void *void_result = ei_conditional_aligned_malloc<Align>(sizeof(T)*size);
-  return ::new(void_result) T[size];
+  T *result = reinterpret_cast<T*>(ei_conditional_aligned_malloc<Align>(sizeof(T)*size));
+  return ei_construct_elements_of_array(result, size);
 }
 
 /** \internal free memory allocated with ei_aligned_malloc
   */
 inline void ei_aligned_free(void *ptr)
 {
-  #if !EIGEN_ARCH_WANTS_ALIGNMENT
-    free(ptr);
+  #if !EIGEN_ALIGN
+    std::free(ptr);
   #elif EIGEN_MALLOC_ALREADY_ALIGNED
-    free(ptr);
+    std::free(ptr);
   #elif EIGEN_HAS_POSIX_MEMALIGN
-    free(ptr);
+    std::free(ptr);
   #elif EIGEN_HAS_MM_MALLOC
     _mm_free(ptr);
   #elif defined(_MSC_VER)
@@ -167,48 +179,78 @@ template<bool Align> inline void ei_conditional_aligned_free(void *ptr)
 
 template<> inline void ei_conditional_aligned_free<false>(void *ptr)
 {
-  free(ptr);
+  std::free(ptr);
 }
 
-/** \internal delete the elements of an array.
+/** \internal destruct the elements of an array.
   * The \a size parameters tells on how many objects to call the destructor of T.
   */
-template<typename T> inline void ei_delete_elements_of_array(T *ptr, size_t size)
+template<typename T> inline void ei_destruct_elements_of_array(T *ptr, std::size_t size)
 {
   // always destruct an array starting from the end.
-  while(size) ptr[--size].~T();
+  if(ptr)
+    while(size) ptr[--size].~T();
 }
 
 /** \internal delete objects constructed with ei_aligned_new
   * The \a size parameters tells on how many objects to call the destructor of T.
   */
-template<typename T> inline void ei_aligned_delete(T *ptr, size_t size)
+template<typename T> inline void ei_aligned_delete(T *ptr, std::size_t size)
 {
-  ei_delete_elements_of_array<T>(ptr, size);
+  ei_destruct_elements_of_array<T>(ptr, size);
   ei_aligned_free(ptr);
 }
 
 /** \internal delete objects constructed with ei_conditional_aligned_new
   * The \a size parameters tells on how many objects to call the destructor of T.
   */
-template<typename T, bool Align> inline void ei_conditional_aligned_delete(T *ptr, size_t size)
+template<typename T, bool Align> inline void ei_conditional_aligned_delete(T *ptr, std::size_t size)
 {
-  ei_delete_elements_of_array<T>(ptr, size);
+  ei_destruct_elements_of_array<T>(ptr, size);
   ei_conditional_aligned_free<Align>(ptr);
 }
 
-/** \internal \returns the number of elements which have to be skipped such that data are 16 bytes aligned */
-template<typename Scalar>
-inline static int ei_alignmentOffset(const Scalar* ptr, int maxOffset)
+/** \internal \returns the index of the first element of the array that is well aligned for vectorization.
+  *
+  * \param array the address of the start of the array
+  * \param size the size of the array
+  *
+  * \note If no element of the array is well aligned, the size of the array is returned. Typically,
+  * for example with SSE, "well aligned" means 16-byte-aligned. If vectorization is disabled or if the
+  * packet size for the given scalar type is 1, then everything is considered well-aligned.
+  *
+  * \note If the scalar type is vectorizable, we rely on the following assumptions: sizeof(Scalar) is a
+  * power of 2, the packet size in bytes is also a power of 2, and is a multiple of sizeof(Scalar). On the
+  * other hand, we do not assume that the array address is a multiple of sizeof(Scalar), as that fails for
+  * example with Scalar=double on certain 32-bit platforms, see bug #79.
+  *
+  * There is also the variant ei_first_aligned(const MatrixBase&, Integer) defined in Coeffs.h.
+  */
+template<typename Scalar, typename Integer>
+inline static Integer ei_alignmentOffset(const Scalar* array, Integer size)
 {
   typedef typename ei_packet_traits<Scalar>::type Packet;
-  const int PacketSize = ei_packet_traits<Scalar>::size;
-  const int PacketAlignedMask = PacketSize-1;
-  const bool Vectorized = PacketSize>1;
-  return Vectorized
-          ? std::min<int>( (PacketSize - (int((size_t(ptr)/sizeof(Scalar))) & PacketAlignedMask))
-                           & PacketAlignedMask, maxOffset)
-          : 0;
+  enum { PacketSize = ei_packet_traits<Scalar>::size,
+         PacketAlignedMask = PacketSize-1
+  };
+  
+  if(PacketSize==1)
+  {
+    // Either there is no vectorization, or a packet consists of exactly 1 scalar so that all elements
+    // of the array have the same aligment.
+    return 0;
+  }
+  else if(std::size_t(array) & (sizeof(Scalar)-1))
+  {
+    // There is vectorization for this scalar type, but the array is not aligned to the size of a single scalar.
+    // Consequently, no element of the array is well aligned.
+    return size;
+  }
+  else
+  {
+    return std::min<Integer>( (PacketSize - (Integer((std::size_t(array)/sizeof(Scalar))) & PacketAlignedMask))
+                           & PacketAlignedMask, size);
+  }
 }
 
 /** \internal
@@ -232,22 +274,45 @@ inline static int ei_alignmentOffset(const Scalar* ptr, int maxOffset)
   #define ei_aligned_stack_free(PTR,SIZE) ei_aligned_free(PTR)
 #endif
 
-#define ei_aligned_stack_new(TYPE,SIZE) ::new(ei_aligned_stack_alloc(sizeof(TYPE)*SIZE)) TYPE[SIZE]
-#define ei_aligned_stack_delete(TYPE,PTR,SIZE) do {ei_delete_elements_of_array<TYPE>(PTR, SIZE); \
+#define ei_aligned_stack_new(TYPE,SIZE) ei_construct_elements_of_array(reinterpret_cast<TYPE*>(ei_aligned_stack_alloc(sizeof(TYPE)*SIZE)), SIZE)
+#define ei_aligned_stack_delete(TYPE,PTR,SIZE) do {ei_destruct_elements_of_array<TYPE>(PTR, SIZE); \
                                                    ei_aligned_stack_free(PTR,sizeof(TYPE)*SIZE);} while(0)
 
 
-#if EIGEN_ARCH_WANTS_ALIGNMENT
+#if EIGEN_ALIGN
+  #ifdef EIGEN_EXCEPTIONS
+    #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \
+      void* operator new(std::size_t size, const std::nothrow_t&) throw() { \
+        try { return Eigen::ei_conditional_aligned_malloc<NeedsToAlign>(size); } \
+        catch (...) { return 0; } \
+        return 0; \
+      }
+  #else
+    #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \
+      void* operator new(std::size_t size, const std::nothrow_t&) throw() { \
+        return Eigen::ei_conditional_aligned_malloc<NeedsToAlign>(size); \
+      }
+  #endif
+
   #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign) \
-      void *operator new(size_t size) throw() { \
+      void *operator new(std::size_t size) { \
         return Eigen::ei_conditional_aligned_malloc<NeedsToAlign>(size); \
       } \
-      void *operator new[](size_t size) throw() { \
+      void *operator new[](std::size_t size) { \
         return Eigen::ei_conditional_aligned_malloc<NeedsToAlign>(size); \
       } \
-      void operator delete(void * ptr) { Eigen::ei_conditional_aligned_free<NeedsToAlign>(ptr); } \
-      void operator delete[](void * ptr) { Eigen::ei_conditional_aligned_free<NeedsToAlign>(ptr); } \
-      void *operator new(size_t, void *ptr) throw() { return ptr; } \
+      void operator delete(void * ptr) throw() { Eigen::ei_conditional_aligned_free<NeedsToAlign>(ptr); } \
+      void operator delete[](void * ptr) throw() { Eigen::ei_conditional_aligned_free<NeedsToAlign>(ptr); } \
+      /* in-place new and delete. since (at least afaik) there is no actual   */ \
+      /* memory allocated we can safely let the default implementation handle */ \
+      /* this particular case. */ \
+      static void *operator new(std::size_t size, void *ptr) { return ::operator new(size,ptr); } \
+      void operator delete(void * memory, void *ptr) throw() { return ::operator delete(memory,ptr); } \
+      /* nothrow-new (returns zero instead of std::bad_alloc) */ \
+      EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \
+      void operator delete(void *ptr, const std::nothrow_t&) throw() { \
+        Eigen::ei_conditional_aligned_free<NeedsToAlign>(ptr); \
+      } \
       typedef void ei_operator_new_marker_type;
 #else
   #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign)
@@ -275,8 +340,8 @@ template<class T>
 class aligned_allocator
 {
 public:
-    typedef size_t    size_type;
-    typedef ptrdiff_t difference_type;
+    typedef std::size_t    size_type;
+    typedef std::ptrdiff_t difference_type;
     typedef T*        pointer;
     typedef const T*  const_pointer;
     typedef T&        reference;

Eigen/src/Core/util/StaticAssert.h

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

Eigen/src/Core/util/XprHelper.h

@@ -37,6 +37,10 @@
 //classes inheriting ei_no_assignment_operator don't generate a default operator=.
 class ei_no_assignment_operator
 {
+#if EIGEN_GCC3_OR_OLDER
+  protected:
+    void nevermind_this_is_just_to_work_around_a_stupid_gcc3_warning();
+#endif
   private:
     ei_no_assignment_operator& operator=(const ei_no_assignment_operator&);
 };
@@ -90,10 +94,16 @@ class ei_compute_matrix_flags
 {
     enum {
       row_major_bit = Options&RowMajor ? RowMajorBit : 0,
-      inner_max_size = row_major_bit ? MaxCols : MaxRows,
+      inner_max_size = int(MaxRows==1) ? int(MaxCols)
+                     : int(MaxCols==1) ? int(MaxRows)
+                     : int(row_major_bit) ? int(MaxCols) : int(MaxRows),
       is_big = inner_max_size == Dynamic,
-      is_packet_size_multiple = (Cols*Rows) % ei_packet_traits<Scalar>::size == 0,
-      aligned_bit = ((Options&AutoAlign) && (is_big || is_packet_size_multiple)) ? AlignedBit : 0,
+      storage_has_fixed_size = MaxRows != Dynamic && MaxCols != Dynamic,
+      storage_has_aligned_fixed_size = storage_has_fixed_size
+                     && ( (MaxCols*MaxRows) % ei_packet_traits<Scalar>::size == 0 ),
+      aligned_bit = (    (Options&AutoAlign)
+                      && (is_big || storage_has_aligned_fixed_size)
+                    ) ? AlignedBit : 0,
       packet_access_bit = ei_packet_traits<Scalar>::size > 1 && aligned_bit ? PacketAccessBit : 0
     };
 
@@ -110,7 +120,7 @@ template<int _Rows, int _Cols> struct ei_size_at_compile_time
  * in order to avoid a useless copy
  */
 
-template<typename T, int Sparseness = ei_traits<T>::Flags&SparseBit> class ei_eval;
+template<typename T, int Sparseness = ei_traits<T>::Flags&SparseBit> struct ei_eval;
 
 template<typename T> struct ei_eval<T,IsDense>
 {
@@ -157,6 +167,19 @@ template<typename T> struct ei_plain_matrix_type_column_major
           > type;
 };
 
+/* ei_plain_matrix_type_row_major : same as ei_plain_matrix_type but guaranteed to be row-major
+ */
+template<typename T> struct ei_plain_matrix_type_row_major
+{
+  typedef Matrix<typename ei_traits<T>::Scalar,
+                ei_traits<T>::RowsAtCompileTime,
+                ei_traits<T>::ColsAtCompileTime,
+                AutoAlign | RowMajor,
+                ei_traits<T>::MaxRowsAtCompileTime,
+                ei_traits<T>::MaxColsAtCompileTime
+          > type;
+};
+
 template<typename T> struct ei_must_nest_by_value { enum { ret = false }; };
 template<typename T> struct ei_must_nest_by_value<NestByValue<T> > { enum { ret = true }; };
 

Eigen/src/Geometry/EulerAngles.h

@@ -60,31 +60,31 @@ MatrixBase<Derived>::eulerAngles(int a0, int a1, int a2) const
   if (a0==a2)
   {
     Scalar s = Vector2(coeff(j,i) , coeff(k,i)).norm();
-    res[1] = std::atan2(s, coeff(i,i));
+    res[1] = ei_atan2(s, coeff(i,i));
     if (s > epsilon)
     {
-      res[0] = std::atan2(coeff(j,i), coeff(k,i));
-      res[2] = std::atan2(coeff(i,j),-coeff(i,k));
+      res[0] = ei_atan2(coeff(j,i), coeff(k,i));
+      res[2] = ei_atan2(coeff(i,j),-coeff(i,k));
     }
     else
     {
       res[0] = Scalar(0);
-      res[2] = (coeff(i,i)>0?1:-1)*std::atan2(-coeff(k,j), coeff(j,j));
+      res[2] = (coeff(i,i)>0?1:-1)*ei_atan2(-coeff(k,j), coeff(j,j));
     }
   }
   else
   {
     Scalar c = Vector2(coeff(i,i) , coeff(i,j)).norm();
-    res[1] = std::atan2(-coeff(i,k), c);
+    res[1] = ei_atan2(-coeff(i,k), c);
     if (c > epsilon)
     {
-      res[0] = std::atan2(coeff(j,k), coeff(k,k));
-      res[2] = std::atan2(coeff(i,j), coeff(i,i));
+      res[0] = ei_atan2(coeff(j,k), coeff(k,k));
+      res[2] = ei_atan2(coeff(i,j), coeff(i,i));
     }
     else
     {
       res[0] = Scalar(0);
-      res[2] = (coeff(i,k)>0?1:-1)*std::atan2(-coeff(k,j), coeff(j,j));
+      res[2] = (coeff(i,k)>0?1:-1)*ei_atan2(-coeff(k,j), coeff(j,j));
     }
   }
   if (!odd)

Eigen/src/Geometry/Hyperplane.h

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

Eigen/src/Geometry/Quaternion.h

@@ -224,17 +224,45 @@ typedef Quaternion<float> Quaternionf;
   * double precision quaternion type */
 typedef Quaternion<double> Quaterniond;
 
+// Generic Quaternion * Quaternion product
+template<int Arch,typename Scalar> inline Quaternion<Scalar>
+ei_quaternion_product(const Quaternion<Scalar>& a, const Quaternion<Scalar>& b)
+{
+  return Quaternion<Scalar>
+  (
+    a.w() * b.w() - a.x() * b.x() - a.y() * b.y() - a.z() * b.z(),
+    a.w() * b.x() + a.x() * b.w() + a.y() * b.z() - a.z() * b.y(),
+    a.w() * b.y() + a.y() * b.w() + a.z() * b.x() - a.x() * b.z(),
+    a.w() * b.z() + a.z() * b.w() + a.x() * b.y() - a.y() * b.x()
+  );
+}
+
+#ifdef EIGEN_VECTORIZE_SSE
+template<> inline Quaternion<float>
+ei_quaternion_product<EiArch_SSE,float>(const Quaternion<float>& _a, const Quaternion<float>& _b)
+{
+  const __m128 mask = _mm_castsi128_ps(_mm_setr_epi32(0,0,0,0x80000000));
+  Quaternion<float> res;
+  __m128 a = _a.coeffs().packet<Aligned>(0);
+  __m128 b = _b.coeffs().packet<Aligned>(0);
+  __m128 flip1 = _mm_xor_ps(_mm_mul_ps(ei_vec4f_swizzle1(a,1,2,0,2),
+                                       ei_vec4f_swizzle1(b,2,0,1,2)),mask);
+  __m128 flip2 = _mm_xor_ps(_mm_mul_ps(ei_vec4f_swizzle1(a,3,3,3,1),
+                                       ei_vec4f_swizzle1(b,0,1,2,1)),mask);
+  ei_pstore(&res.x(),
+            _mm_add_ps(_mm_sub_ps(_mm_mul_ps(a,ei_vec4f_swizzle1(b,3,3,3,3)),
+                                  _mm_mul_ps(ei_vec4f_swizzle1(a,2,0,1,0),
+                                             ei_vec4f_swizzle1(b,1,2,0,0))),
+                       _mm_add_ps(flip1,flip2)));
+  return res;
+}
+#endif
+
 /** \returns the concatenation of two rotations as a quaternion-quaternion product */
 template <typename Scalar>
 inline Quaternion<Scalar> Quaternion<Scalar>::operator* (const Quaternion& other) const
 {
-  return Quaternion
-  (
-    this->w() * other.w() - this->x() * other.x() - this->y() * other.y() - this->z() * other.z(),
-    this->w() * other.x() + this->x() * other.w() + this->y() * other.z() - this->z() * other.y(),
-    this->w() * other.y() + this->y() * other.w() + this->z() * other.x() - this->x() * other.z(),
-    this->w() * other.z() + this->z() * other.w() + this->x() * other.y() - this->y() * other.x()
-  );
+  return ei_quaternion_product<EiArch>(*this,other);
 }
 
 /** \sa operator*(Quaternion) */
@@ -346,7 +374,6 @@ inline Quaternion<Scalar>& Quaternion<Scalar>::setFromTwoVectors(const MatrixBas
 {
   Vector3 v0 = a.normalized();
   Vector3 v1 = b.normalized();
-  Vector3 axis = v0.cross(v1);
   Scalar c = v0.dot(v1);
 
   // if dot == 1, vectors are the same
@@ -354,7 +381,17 @@ inline Quaternion<Scalar>& Quaternion<Scalar>::setFromTwoVectors(const MatrixBas
   {
     // set to identity
     this->w() = 1; this->vec().setZero();
+    return *this;
   }
+  // if dot == -1, vectors are opposites
+  if (ei_isApprox(c,Scalar(-1)))
+  {
+    this->vec() = v0.unitOrthogonal();
+    this->w() = 0;
+    return *this;
+  }
+
+  Vector3 axis = v0.cross(v1);
   Scalar s = ei_sqrt((Scalar(1)+c)*Scalar(2));
   Scalar invs = Scalar(1)/s;
   this->vec() = axis * invs;
@@ -413,22 +450,31 @@ inline Scalar Quaternion<Scalar>::angularDistance(const Quaternion& other) const
 template <typename Scalar>
 Quaternion<Scalar> Quaternion<Scalar>::slerp(Scalar t, const Quaternion& other) const
 {
-  static const Scalar one = Scalar(1) - precision<Scalar>();
+  static const Scalar one = Scalar(1) - machine_epsilon<Scalar>();
   Scalar d = this->dot(other);
   Scalar absD = ei_abs(d);
+
+  Scalar scale0;
+  Scalar scale1;
+
   if (absD>=one)
-    return *this;
+  {
+    scale0 = Scalar(1) - t;
+    scale1 = t;
+  }
+  else
+  {
+    // theta is the angle between the 2 quaternions
+    Scalar theta = std::acos(absD);
+    Scalar sinTheta = ei_sin(theta);
 
-  // theta is the angle between the 2 quaternions
-  Scalar theta = std::acos(absD);
-  Scalar sinTheta = ei_sin(theta);
+    scale0 = ei_sin( ( Scalar(1) - t ) * theta) / sinTheta;
+    scale1 = ei_sin( ( t * theta) ) / sinTheta;
+    if (d<0)
+      scale1 = -scale1;
+  }
 
-  Scalar scale0 = ei_sin( ( Scalar(1) - t ) * theta) / sinTheta;
-  Scalar scale1 = ei_sin( ( t * theta) ) / sinTheta;
-  if (d<0)
-    scale1 = -scale1;
-
-  return Quaternion(scale0 * m_coeffs + scale1 * other.m_coeffs);
+  return Quaternion<Scalar>(scale0 * coeffs() + scale1 * other.coeffs());
 }
 
 // set from a rotation matrix

Eigen/src/Geometry/Rotation2D.h

@@ -85,7 +85,7 @@ public:
 
   /** Concatenates two rotations */
   inline Rotation2D& operator*=(const Rotation2D& other)
-  { return m_angle += other.m_angle; }
+  { return m_angle += other.m_angle; return *this; }
 
   /** Applies the rotation to a 2D vector */
   Vector2 operator* (const Vector2& vec) const

Eigen/src/Geometry/Transform.h

@@ -198,6 +198,10 @@ public:
 
   /** \sa MatrixBase::setIdentity() */
   void setIdentity() { m_matrix.setIdentity(); }
+  static const typename MatrixType::IdentityReturnType Identity()
+  {
+    return MatrixType::Identity();
+  }
 
   template<typename OtherDerived>
   inline Transform& scale(const MatrixBase<OtherDerived> &other);
@@ -283,6 +287,10 @@ public:
   bool isApprox(const Transform& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const
   { return m_matrix.isApprox(other.m_matrix, prec); }
 
+  #ifdef EIGEN_TRANSFORM_PLUGIN
+  #include EIGEN_TRANSFORM_PLUGIN
+  #endif
+
 protected:
 
 };
@@ -335,9 +343,9 @@ template<typename Scalar, int Dim>
 QMatrix Transform<Scalar,Dim>::toQMatrix(void) const
 {
   EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE)
-  return QMatrix(other.coeffRef(0,0), other.coeffRef(1,0),
-                 other.coeffRef(0,1), other.coeffRef(1,1),
-                 other.coeffRef(0,2), other.coeffRef(1,2));
+  return QMatrix(m_matrix.coeff(0,0), m_matrix.coeff(1,0),
+                 m_matrix.coeff(0,1), m_matrix.coeff(1,1),
+                 m_matrix.coeff(0,2), m_matrix.coeff(1,2));
 }
 
 /** Initialises \c *this from a QTransform assuming the dimension is 2.
@@ -369,12 +377,12 @@ Transform<Scalar,Dim>& Transform<Scalar,Dim>::operator=(const QTransform& other)
   * This function is available only if the token EIGEN_QT_SUPPORT is defined.
   */
 template<typename Scalar, int Dim>
-QMatrix Transform<Scalar,Dim>::toQTransform(void) const
+QTransform Transform<Scalar,Dim>::toQTransform(void) const
 {
   EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE)
-  return QTransform(other.coeffRef(0,0), other.coeffRef(1,0), other.coeffRef(2,0)
-                    other.coeffRef(0,1), other.coeffRef(1,1), other.coeffRef(2,1)
-                    other.coeffRef(0,2), other.coeffRef(1,2), other.coeffRef(2,2);
+  return QTransform(m_matrix.coeff(0,0), m_matrix.coeff(1,0), m_matrix.coeff(2,0),
+                    m_matrix.coeff(0,1), m_matrix.coeff(1,1), m_matrix.coeff(2,1),
+                    m_matrix.coeff(0,2), m_matrix.coeff(1,2), m_matrix.coeff(2,2));
 }
 #endif
 

Eigen/src/LU/Inverse.h

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

Eigen/src/LU/LU.h

@@ -63,12 +63,12 @@ template<typename MatrixType> class LU
 
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
-    typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
-    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
-    typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType;
-    typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType;
+    typedef Matrix<int, 1, MatrixType::ColsAtCompileTime, MatrixType::Options, 1, MatrixType::MaxColsAtCompileTime> IntRowVectorType;
+    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, 1> IntColVectorType;
+    typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime, MatrixType::Options, 1, MatrixType::MaxColsAtCompileTime> RowVectorType;
+    typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1, MatrixType::Options, MatrixType::MaxRowsAtCompileTime, 1> ColVectorType;
 
-    enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN(
+    enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN(
              MatrixType::MaxColsAtCompileTime,
              MatrixType::MaxRowsAtCompileTime)
     };
@@ -297,7 +297,8 @@ template<typename MatrixType> class LU
       *
       * \sa MatrixBase::computeInverse(), inverse()
       */
-    inline void computeInverse(MatrixType *result) const
+    template<typename ResultType>
+    inline void computeInverse(ResultType *result) const
     {
       solve(MatrixType::Identity(m_lu.rows(), m_lu.cols()), result);
     }
@@ -508,15 +509,16 @@ bool LU<MatrixType>::solve(
   if(!isSurjective())
   {
     // is c is in the image of U ?
-    RealScalar biggest_in_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff();
+    RealScalar biggest_in_c = m_rank>0 ? c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff() : RealScalar(0);
     for(int col = 0; col < c.cols(); ++col)
       for(int row = m_rank; row < c.rows(); ++row)
         if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c, m_precision))
           return false;
   }
-  m_lu.corner(TopLeft, m_rank, m_rank)
-      .template marked<UpperTriangular>()
-      .solveTriangularInPlace(c.corner(TopLeft, m_rank, c.cols()));
+  if(m_rank>0)
+    m_lu.corner(TopLeft, m_rank, m_rank)
+        .template marked<UpperTriangular>()
+        .solveTriangularInPlace(c.corner(TopLeft, m_rank, c.cols()));
 
   // Step 4
   result->resize(m_lu.cols(), b.cols());

Eigen/src/LeastSquares/LeastSquares.h

@@ -54,10 +54,13 @@
   * constants \f$a,b,c\f$ so that the plane of equation \f$y=ax+bz+c\f$ fits
   * best the five above points. To do that, call this function as follows:
   * @code
+    // create a vector of pointers to the points
+    std::vector<Vector3d> points_ptrs(5);
+    for(int k=0; k<5; ++k) points_ptrs[k] = &points[k];
     Vector3d coeffs; // will store the coefficients a, b, c
     linearRegression(
       5,
-      &points,
+      &(points_ptrs[0]),
       &coeffs,
       1 // the coord to express as a function of
         // the other ones. 0 means x, 1 means y, 2 means z.

Eigen/src/QR/QR.h

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

Eigen/src/QR/Tridiagonalization.h

@@ -293,7 +293,7 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
       {
         int starti = i+1;
         int alignedEnd = starti;
-        if (PacketSize>1)
+        if (PacketSize>1 && (int(MatrixType::Flags)&RowMajor) == 0)
         {
           int alignedStart = (starti) + ei_alignmentOffset(&matA.coeffRef(starti,j1), n-starti);
           alignedEnd = alignedStart + ((n-alignedStart)/PacketSize)*PacketSize;
@@ -391,7 +391,7 @@ void Tridiagonalization<MatrixType>::_decomposeInPlace3x3(MatrixType& mat, Diago
 {
   diag[0] = ei_real(mat(0,0));
   RealScalar v1norm2 = ei_abs2(mat(0,2));
-  if (ei_isMuchSmallerThan(v1norm2, RealScalar(1)))
+  if (v1norm2==RealScalar(0))
   {
     diag[1] = ei_real(mat(1,1));
     diag[2] = ei_real(mat(2,2));

Eigen/src/SVD/SVD.h

@@ -49,7 +49,7 @@ template<typename MatrixType> class SVD
     enum {
       PacketSize = ei_packet_traits<Scalar>::size,
       AlignmentMask = int(PacketSize)-1,
-      MinSize = EIGEN_ENUM_MIN(MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime)
+      MinSize = EIGEN_SIZE_MIN(MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime)
     };
 
     typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVector;
@@ -61,6 +61,8 @@ template<typename MatrixType> class SVD
 
   public:
 
+    SVD() {} // a user who relied on compiler-generated default compiler reported problems with MSVC in 2.0.7
+    
     SVD(const MatrixType& matrix)
       : m_matU(matrix.rows(), std::min(matrix.rows(), matrix.cols())),
         m_matV(matrix.cols(),matrix.cols()),
@@ -107,6 +109,8 @@ void SVD<MatrixType>::compute(const MatrixType& matrix)
   const int m = matrix.rows();
   const int n = matrix.cols();
   const int nu = std::min(m,n);
+  ei_assert(m>=n && "In Eigen 2.0, SVD only works for MxN matrices with M>=N. Sorry!");
+  ei_assert(m>1 && "In Eigen 2.0, SVD doesn't work on 1x1 matrices");
 
   m_matU.resize(m, nu);
   m_matU.setZero();

Eigen/src/Sparse/AmbiVector.h

@@ -98,7 +98,9 @@ template<typename _Scalar> class AmbiVector
       int allocSize = m_allocatedElements * sizeof(ListEl);
       allocSize = allocSize/sizeof(Scalar) + (allocSize%sizeof(Scalar)>0?1:0);
       Scalar* newBuffer = new Scalar[allocSize];
-      memcpy(newBuffer,  m_buffer,  copyElements * sizeof(ListEl));
+	  std::memcpy(newBuffer,  m_buffer,  copyElements * sizeof(ListEl));
+      delete[] m_buffer;
+      m_buffer = newBuffer;
     }
 
   protected:
@@ -238,8 +240,11 @@ Scalar& AmbiVector<Scalar>::coeffRef(int i)
       else
       {
         if (m_llSize>=m_allocatedElements)
+        {
           reallocateSparse();
-        ei_internal_assert(m_llSize<m_size && "internal error: overflow in sparse mode");
+          llElements = reinterpret_cast<ListEl*>(m_buffer);
+        }
+        ei_internal_assert(m_llSize<m_allocatedElements && "internal error: overflow in sparse mode");
         // let's insert a new coefficient
         ListEl& el = llElements[m_llSize];
         el.value = Scalar(0);
@@ -365,6 +370,9 @@ class AmbiVector<_Scalar>::Iterator
     int m_cachedIndex;          // current coordinate
     Scalar m_cachedValue;       // current value
     bool m_isDense;             // mode of the vector
+
+  private:
+    Iterator& operator=(const Iterator&);
 };
 
 

Eigen/src/Sparse/CompressedStorage.h

@@ -37,7 +37,7 @@ class CompressedStorage
       : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0)
     {}
 
-    CompressedStorage(size_t size)
+    CompressedStorage(std::size_t size)
       : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0)
     {
       resize(size);
@@ -52,8 +52,8 @@ class CompressedStorage
     CompressedStorage& operator=(const CompressedStorage& other)
     {
       resize(other.size());
-      memcpy(m_values, other.m_values, m_size * sizeof(Scalar));
-      memcpy(m_indices, other.m_indices, m_size * sizeof(int));
+	  std::memcpy(m_values, other.m_values, m_size * sizeof(Scalar));
+	  std::memcpy(m_indices, other.m_indices, m_size * sizeof(int));
       return *this;
     }
 
@@ -71,9 +71,9 @@ class CompressedStorage
       delete[] m_indices;
     }
 
-    void reserve(size_t size)
+    void reserve(std::size_t size)
     {
-      size_t newAllocatedSize = m_size + size;
+      std::size_t newAllocatedSize = m_size + size;
       if (newAllocatedSize > m_allocatedSize)
         reallocate(newAllocatedSize);
     }
@@ -84,10 +84,10 @@ class CompressedStorage
         reallocate(m_size);
     }
 
-    void resize(size_t size, float reserveSizeFactor = 0)
+    void resize(std::size_t size, float reserveSizeFactor = 0)
     {
       if (m_allocatedSize<size)
-        reallocate(size + size_t(reserveSizeFactor*size));
+        reallocate(size + std::size_t(reserveSizeFactor*size));
       m_size = size;
     }
 
@@ -99,17 +99,17 @@ class CompressedStorage
       m_indices[id] = i;
     }
 
-    inline size_t size() const { return m_size; }
-    inline size_t allocatedSize() const { return m_allocatedSize; }
+    inline std::size_t size() const { return m_size; }
+    inline std::size_t allocatedSize() const { return m_allocatedSize; }
     inline void clear() { m_size = 0; }
 
-    inline Scalar& value(size_t i) { return m_values[i]; }
-    inline const Scalar& value(size_t i) const { return m_values[i]; }
+    inline Scalar& value(std::size_t i) { return m_values[i]; }
+    inline const Scalar& value(std::size_t i) const { return m_values[i]; }
 
-    inline int& index(size_t i) { return m_indices[i]; }
-    inline const int& index(size_t i) const { return m_indices[i]; }
+    inline int& index(std::size_t i) { return m_indices[i]; }
+    inline const int& index(std::size_t i) const { return m_indices[i]; }
 
-    static CompressedStorage Map(int* indices, Scalar* values, size_t size)
+    static CompressedStorage Map(int* indices, Scalar* values, std::size_t size)
     {
       CompressedStorage res;
       res.m_indices = indices;
@@ -125,11 +125,11 @@ class CompressedStorage
     }
     
     /** \returns the largest \c k in [start,end) such that for all \c j in [start,k) index[\c j]\<\a key */
-    inline int searchLowerIndex(size_t start, size_t end, int key) const
+    inline int searchLowerIndex(std::size_t start, std::size_t end, int key) const
     {
       while(end>start)
       {
-        size_t mid = (end+start)>>1;
+        std::size_t mid = (end+start)>>1;
         if (m_indices[mid]<key)
           start = mid+1;
         else
@@ -148,12 +148,12 @@ class CompressedStorage
         return m_values[m_size-1];
       // ^^  optimization: let's first check if it is the last coefficient
       // (very common in high level algorithms)
-      const size_t id = searchLowerIndex(0,m_size-1,key);
+      const std::size_t id = searchLowerIndex(0,m_size-1,key);
       return ((id<m_size) && (m_indices[id]==key)) ? m_values[id] : defaultValue;
     }
     
     /** Like at(), but the search is performed in the range [start,end) */
-    inline Scalar atInRange(size_t start, size_t end, int key, Scalar defaultValue = Scalar(0)) const
+    inline Scalar atInRange(std::size_t start, std::size_t end, int key, Scalar defaultValue = Scalar(0)) const
     {
       if (start==end)
         return Scalar(0);
@@ -161,7 +161,7 @@ class CompressedStorage
         return m_values[end-1];
       // ^^  optimization: let's first check if it is the last coefficient
       // (very common in high level algorithms)
-      const size_t id = searchLowerIndex(start,end-1,key);
+      const std::size_t id = searchLowerIndex(start,end-1,key);
       return ((id<end) && (m_indices[id]==key)) ? m_values[id] : defaultValue;
     }
     
@@ -170,11 +170,11 @@ class CompressedStorage
       * such that the keys are sorted. */
     inline Scalar& atWithInsertion(int key, Scalar defaultValue = Scalar(0))
     {
-      size_t id = searchLowerIndex(0,m_size,key);
+      std::size_t id = searchLowerIndex(0,m_size,key);
       if (id>=m_size || m_indices[id]!=key)
       {
         resize(m_size+1,1);
-        for (size_t j=m_size-1; j>id; --j)
+        for (std::size_t j=m_size-1; j>id; --j)
         {
           m_indices[j] = m_indices[j-1];
           m_values[j] = m_values[j-1];
@@ -187,9 +187,9 @@ class CompressedStorage
     
     void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>())
     {
-      size_t k = 0;
-      size_t n = size();
-      for (size_t i=0; i<n; ++i)
+      std::size_t k = 0;
+      std::size_t n = size();
+      for (std::size_t i=0; i<n; ++i)
       {
         if (!ei_isMuchSmallerThan(value(i), reference, epsilon))
         {
@@ -203,14 +203,14 @@ class CompressedStorage
 
   protected:
 
-    inline void reallocate(size_t size)
+    inline void reallocate(std::size_t size)
     {
       Scalar* newValues  = new Scalar[size];
       int* newIndices = new int[size];
-      size_t copySize = std::min(size, m_size);
+      std::size_t copySize = std::min(size, m_size);
       // copy
-      memcpy(newValues,  m_values,  copySize * sizeof(Scalar));
-      memcpy(newIndices, m_indices, copySize * sizeof(int));
+	  std::memcpy(newValues,  m_values,  copySize * sizeof(Scalar));
+	  std::memcpy(newIndices, m_indices, copySize * sizeof(int));
       // delete old stuff
       delete[] m_values;
       delete[] m_indices;
@@ -222,8 +222,8 @@ class CompressedStorage
   protected:
     Scalar* m_values;
     int* m_indices;
-    size_t m_size;
-    size_t m_allocatedSize;
+    std::size_t m_size;
+    std::size_t m_allocatedSize;
 
 };
 

Eigen/src/Sparse/DynamicSparseMatrix.h

@@ -212,7 +212,7 @@ class DynamicSparseMatrix
         // remove all coefficients with innerCoord>=innerSize
         // TODO
         std::cerr << "not implemented yet\n";
-        exit(2);
+        std::exit(2);
       }
       if (m_data.size() != outerSize)
       {
@@ -289,9 +289,11 @@ class DynamicSparseMatrix<Scalar,_Flags>::InnerIterator : public SparseVector<Sc
     inline int row() const { return IsRowMajor ? m_outer : Base::index(); }
     inline int col() const { return IsRowMajor ? Base::index() : m_outer; }
 
-
   protected:
     const int m_outer;
+
+  private:
+    InnerIterator& operator=(const InnerIterator&);
 };
 
 #endif // EIGEN_DYNAMIC_SPARSEMATRIX_H

Eigen/src/Sparse/SparseBlock.h

@@ -53,6 +53,9 @@ class SparseInnerVectorSet : ei_no_assignment_operator,
         inline InnerIterator(const SparseInnerVectorSet& xpr, int outer)
           : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer)
         {}
+
+      private:
+        InnerIterator& operator=(const InnerIterator&);
     };
 
     inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize)
@@ -110,6 +113,8 @@ class SparseInnerVectorSet<DynamicSparseMatrix<_Scalar, _Options>, Size>
         inline InnerIterator(const SparseInnerVectorSet& xpr, int outer)
           : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer)
         {}
+      private:
+        InnerIterator& operator=(const InnerIterator&);
     };
 
     inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize)

Eigen/src/Sparse/SparseCwise.h

@@ -156,6 +156,9 @@ template<typename ExpressionType> class SparseCwise
 
   protected:
     ExpressionTypeNested m_matrix;
+
+  private:
+    SparseCwise& operator=(const SparseCwise&);
 };
 
 template<typename Derived>

Eigen/src/Sparse/SparseCwiseBinaryOp.h

@@ -126,6 +126,8 @@ class SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator
     EIGEN_STRONG_INLINE InnerIterator(const SparseCwiseBinaryOp& binOp, int outer)
       : Base(binOp,outer)
     {}
+  private:
+    InnerIterator& operator=(const InnerIterator&);
 };
 
 /***************************************************************************
@@ -197,6 +199,9 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp, Lhs, Rhs, Deri
     const BinaryOp& m_functor;
     Scalar m_value;
     int m_id;
+
+  private:
+    ei_sparse_cwise_binary_op_inner_iterator_selector& operator=(const ei_sparse_cwise_binary_op_inner_iterator_selector&);
 };
 
 // sparse - sparse  (product)
@@ -250,6 +255,9 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>,
     LhsIterator m_lhsIter;
     RhsIterator m_rhsIter;
     const BinaryFunc& m_functor;
+
+  private:
+    ei_sparse_cwise_binary_op_inner_iterator_selector& operator=(const ei_sparse_cwise_binary_op_inner_iterator_selector&);
 };
 
 // sparse - dense  (product)
@@ -290,6 +298,9 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>,
     LhsIterator m_lhsIter;
     const BinaryFunc m_functor;
     const int m_outer;
+
+  private:
+    ei_sparse_cwise_binary_op_inner_iterator_selector& operator=(const ei_sparse_cwise_binary_op_inner_iterator_selector&);
 };
 
 // sparse - dense  (product)

Eigen/src/Sparse/SparseCwiseUnaryOp.h

@@ -90,6 +90,9 @@ class SparseCwiseUnaryOp<UnaryOp,MatrixType>::InnerIterator
   protected:
     MatrixTypeIterator m_iter;
     const UnaryOp m_functor;
+
+  private:
+    InnerIterator& operator=(const InnerIterator&);
 };
 
 template<typename Derived>

Eigen/src/Sparse/SparseDiagonalProduct.h

@@ -120,6 +120,8 @@ class ei_sparse_diagonal_product_inner_iterator_selector
               const SparseDiagonalProductType& expr, int outer)
       : Base(expr.rhs().innerVector(outer) .cwise()* expr.lhs().diagonal(), 0)
     {}
+  private:
+    ei_sparse_diagonal_product_inner_iterator_selector& operator=(const ei_sparse_diagonal_product_inner_iterator_selector&);
 };
 
 template<typename Lhs, typename Rhs, typename SparseDiagonalProductType>

Eigen/src/Sparse/SparseFlagged.h

@@ -64,16 +64,21 @@ template<typename ExpressionType, unsigned int Added, unsigned int Removed> clas
 
   protected:
     ExpressionTypeNested m_matrix;
+
+  private:
+    SparseFlagged& operator=(const SparseFlagged&);
 };
 
 template<typename ExpressionType, unsigned int Added, unsigned int Removed>
   class SparseFlagged<ExpressionType,Added,Removed>::InnerIterator : public ExpressionType::InnerIterator
 {
   public:
-
     EIGEN_STRONG_INLINE InnerIterator(const SparseFlagged& xpr, int outer)
       : ExpressionType::InnerIterator(xpr.m_matrix, outer)
     {}
+
+  private:
+    InnerIterator& operator=(const InnerIterator&);
 };
 
 template<typename ExpressionType, unsigned int Added, unsigned int Removed>

Eigen/src/Sparse/SparseLU.h

@@ -96,8 +96,8 @@ class SparseLU
 
     void setOrderingMethod(int m)
     {
-      ei_assert(m&~OrderingMask == 0 && m!=0 && "invalid ordering method");
-      m_flags = m_flags&~OrderingMask | m&OrderingMask;
+      ei_assert((m&~OrderingMask) == 0 && m!=0 && "invalid ordering method");
+      m_flags = (m_flags&~OrderingMask) | (m&OrderingMask);
     }
 
     int orderingMethod() const

Eigen/src/Sparse/SparseMatrix.h

@@ -122,7 +122,7 @@ class SparseMatrix
     {
       m_data.clear();
       //if (m_outerSize)
-      memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(int));
+	  std::memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(int));
 //       for (int i=0; i<m_outerSize; ++i)
 //         m_outerIndex[i] = 0;
 //       if (m_outerSize)
@@ -164,7 +164,7 @@ class SparseMatrix
       {
         ei_assert(m_data.index(m_data.size()-1)<inner && "wrong sorted insertion");
       }
-      assert(size_t(m_outerIndex[outer+1]) == m_data.size());
+      assert(std::size_t(m_outerIndex[outer+1]) == m_data.size());
       int id = m_outerIndex[outer+1];
       ++m_outerIndex[outer+1];
 
@@ -190,10 +190,10 @@ class SparseMatrix
         }
         m_outerIndex[outer+1] = m_outerIndex[outer];
       }
-      assert(size_t(m_outerIndex[outer+1]) == m_data.size() && "invalid outer index");
-      size_t startId = m_outerIndex[outer];
-      // FIXME let's make sure sizeof(long int) == sizeof(size_t)
-      size_t id = m_outerIndex[outer+1];
+      assert(std::size_t(m_outerIndex[outer+1]) == m_data.size() && "invalid outer index");
+      std::size_t startId = m_outerIndex[outer];
+      // FIXME let's make sure sizeof(long int) == sizeof(std::size_t)
+      std::size_t id = m_outerIndex[outer+1];
       ++m_outerIndex[outer+1];
 
       float reallocRatio = 1;
@@ -259,19 +259,21 @@ class SparseMatrix
       m_data.resize(k,0);
     }
 
+    /** Resizes the matrix to a \a rows x \a cols matrix and initializes it to zero
+      * \sa resizeNonZeros(int), reserve(), setZero()
+      */
     void resize(int rows, int cols)
     {
-//       std::cerr << this << " resize " << rows << "x" << cols << "\n";
       const int outerSize = IsRowMajor ? rows : cols;
       m_innerSize = IsRowMajor ? cols : rows;
       m_data.clear();
-      if (m_outerSize != outerSize)
+      if (m_outerSize != outerSize || m_outerSize==0)
       {
         delete[] m_outerIndex;
         m_outerIndex = new int [outerSize+1];
         m_outerSize = outerSize;
-        memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(int));
       }
+	  std::memset(m_outerIndex, 0, (m_outerSize+1)*sizeof(int));
     }
     void resizeNonZeros(int size)
     {
@@ -322,7 +324,7 @@ class SparseMatrix
       else
       {
         resize(other.rows(), other.cols());
-        memcpy(m_outerIndex, other.m_outerIndex, (m_outerSize+1)*sizeof(int));
+		std::memcpy(m_outerIndex, other.m_outerIndex, (m_outerSize+1)*sizeof(int));
         m_data = other.m_data;
       }
       return *this;
@@ -442,6 +444,9 @@ class SparseMatrix<Scalar,_Flags>::InnerIterator
     int m_id;
     const int m_start;
     const int m_end;
+
+  private:
+    InnerIterator& operator=(const InnerIterator&);
 };
 
 #endif // EIGEN_SPARSEMATRIX_H

Eigen/src/Sparse/SparseProduct.h

@@ -97,7 +97,7 @@ struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >
 
     RowsAtCompileTime = _LhsNested::RowsAtCompileTime,
     ColsAtCompileTime = _RhsNested::ColsAtCompileTime,
-    InnerSize = EIGEN_ENUM_MIN(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
+    InnerSize = EIGEN_SIZE_MIN(_LhsNested::ColsAtCompileTime, _RhsNested::RowsAtCompileTime),
 
     MaxRowsAtCompileTime = _LhsNested::MaxRowsAtCompileTime,
     MaxColsAtCompileTime = _RhsNested::MaxColsAtCompileTime,

Eigen/src/Sparse/SparseTranspose.h

@@ -62,15 +62,20 @@ template<typename MatrixType> class SparseTranspose
 
   protected:
     const typename MatrixType::Nested m_matrix;
+
+  private:
+    SparseTranspose& operator=(const SparseTranspose&);
 };
 
 template<typename MatrixType> class SparseTranspose<MatrixType>::InnerIterator : public MatrixType::InnerIterator
 {
   public:
-
     EIGEN_STRONG_INLINE InnerIterator(const SparseTranspose& trans, int outer)
       : MatrixType::InnerIterator(trans.m_matrix, outer)
     {}
+
+  private:
+    InnerIterator& operator=(const InnerIterator&);
 };
 
 template<typename MatrixType> class SparseTranspose<MatrixType>::ReverseInnerIterator : public MatrixType::ReverseInnerIterator

Eigen/src/Sparse/SparseVector.h

@@ -360,6 +360,9 @@ class SparseVector<Scalar,_Flags>::InnerIterator
     const CompressedStorage<Scalar>& m_data;
     int m_id;
     const int m_end;
+
+  private:
+    InnerIterator& operator=(const InnerIterator&);
 };
 
 #endif // EIGEN_SPARSEVECTOR_H

Eigen/src/Sparse/TriangularSolver.h

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

bench/BenchTimer.h

@@ -1,8 +1,8 @@
 // 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>
-// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -26,8 +26,15 @@
 #ifndef EIGEN_BENCH_TIMER_H
 #define EIGEN_BENCH_TIMER_H
 
-#include <sys/time.h>
+#if defined(_WIN32) || defined(__CYGWIN__)
+#define NOMINMAX
+#define WIN32_LEAN_AND_MEAN
+#include <windows.h>
+#else
+#include <time.h>
 #include <unistd.h>
+#endif
+
 #include <cstdlib>
 #include <numeric>
 
@@ -35,12 +42,25 @@ namespace Eigen
 {
 
 /** Elapsed time timer keeping the best try.
+  *
+  * On POSIX platforms we use clock_gettime with CLOCK_PROCESS_CPUTIME_ID.
+  * On Windows we use QueryPerformanceCounter
+  *
+  * Important: on linux, you must link with -lrt
   */
 class BenchTimer
 {
 public:
 
-  BenchTimer() { reset(); }
+  BenchTimer() 
+  { 
+#if defined(_WIN32) || defined(__CYGWIN__)
+    LARGE_INTEGER freq;
+    QueryPerformanceFrequency(&freq);
+    m_frequency = (double)freq.QuadPart;
+#endif
+    reset(); 
+  }
 
   ~BenchTimer() {}
 
@@ -51,23 +71,34 @@ public:
     m_best = std::min(m_best, getTime() - m_start);
   }
 
-  /** Return the best elapsed time.
+  /** Return the best elapsed time in seconds.
     */
   inline double value(void)
   {
-      return m_best;
+    return m_best;
   }
 
+#if defined(_WIN32) || defined(__CYGWIN__)
+  inline double getTime(void)
+#else
   static inline double getTime(void)
+#endif
   {
-      struct timeval tv;
-      struct timezone tz;
-      gettimeofday(&tv, &tz);
-      return (double)tv.tv_sec + 1.e-6 * (double)tv.tv_usec;
+#ifdef WIN32
+    LARGE_INTEGER query_ticks;
+    QueryPerformanceCounter(&query_ticks);
+    return query_ticks.QuadPart/m_frequency;
+#else
+    timespec ts;
+    clock_gettime(CLOCK_PROCESS_CPUTIME_ID, &ts);
+    return double(ts.tv_sec) + 1e-9 * double(ts.tv_nsec);
+#endif
   }
 
 protected:
-
+#if defined(_WIN32) || defined(__CYGWIN__)
+  double m_frequency;
+#endif
   double m_best, m_start;
 
 };

bench/btl/libs/eigen2/eigen2_interface.hh

@@ -142,7 +142,7 @@ public :
   }
 
   static inline void cholesky(const gene_matrix & X, gene_matrix & C, int N){
-    C = X.cholesky().matrixL();
+    C = X.llt().matrixL();
 //     C = X;
 //     Cholesky<gene_matrix>::computeInPlace(C);
 //     Cholesky<gene_matrix>::computeInPlaceBlock(C);

cmake/FindStandardMathLibrary.cmake

@@ -0,0 +1,64 @@
+# - Try to find how to link to the standard math library, if anything at all is needed to do.
+# On most platforms this is automatic, but for example it's not automatic on QNX.
+#
+# Once done this will define
+#
+#  STANDARD_MATH_LIBRARY_FOUND - we found how to successfully link to the standard math library
+#  STANDARD_MATH_LIBRARY - the name of the standard library that one has to link to.
+#                            -- this will be left empty if it's automatic (most platforms).
+#                            -- this will be set to "m" on platforms where one must explicitly
+#                               pass the "-lm" linker flag.
+#
+# Copyright (c) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
+# Redistribution and use is allowed according to the terms of the 2-clause BSD license.
+
+
+include(CheckCXXSourceCompiles)
+
+# a little test program for c++ math functions.
+# notice the std:: is required on some platforms such as QNX
+
+set(find_standard_math_library_test_program
+"#include<cmath>
+int main() { std::sin(0.0); std::log(0.0f); }")
+
+# first try compiling/linking the test program without any linker flags
+
+set(CMAKE_REQUIRED_FLAGS "")
+set(CMAKE_REQUIRED_LIBRARIES "")
+CHECK_CXX_SOURCE_COMPILES(
+  "${find_standard_math_library_test_program}"
+  standard_math_library_linked_to_automatically
+)
+
+if(standard_math_library_linked_to_automatically)
+
+  # the test program linked successfully without any linker flag.
+  set(STANDARD_MATH_LIBRARY "")
+  set(STANDARD_MATH_LIBRARY_FOUND TRUE)
+
+else()
+
+  # the test program did not link successfully without any linker flag.
+  # This is a very uncommon case that so far we only saw on QNX. The next try is the
+  # standard name 'm' for the standard math library.
+
+  set(CMAKE_REQUIRED_LIBRARIES "m")
+  CHECK_CXX_SOURCE_COMPILES(
+    "${find_standard_math_library_test_program}"
+    standard_math_library_linked_to_as_m)
+
+  if(standard_math_library_linked_to_as_m)
+
+    # the test program linked successfully when linking to the 'm' library
+    set(STANDARD_MATH_LIBRARY "m")
+    set(STANDARD_MATH_LIBRARY_FOUND TRUE)
+
+  else()
+
+    # the test program still doesn't link successfully
+    set(STANDARD_MATH_LIBRARY_FOUND FALSE)
+
+  endif()
+
+endif()

doc/CMakeLists.txt

@@ -38,7 +38,6 @@ add_custom_target(
                                    ${CMAKE_CURRENT_BINARY_DIR}/html/
   COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/Eigen_Silly_Professor_64x64.png
                                    ${CMAKE_CURRENT_BINARY_DIR}/html/
-  COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/buildexamplelist.sh ${Eigen_SOURCE_DIR} > ${CMAKE_CURRENT_BINARY_DIR}/ExampleList.dox
   COMMAND doxygen
   COMMAND ${CMAKE_CURRENT_SOURCE_DIR}/cleanhierarchy.sh ${CMAKE_CURRENT_BINARY_DIR}/html/hierarchy.html
   WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}

doc/CustomizingEigen.dox

@@ -57,10 +57,10 @@ void makeFloor(const MatrixBase<OtherDerived>& other) { derived() = derived().cw
 template<typename OtherDerived>
 void makeCeil(const MatrixBase<OtherDerived>& other) { derived() = derived().cwise().max(other.derived()); }
 
-const typename Cwise<Derived>::ScalarAddReturnType
-operator+(const Scalar& scalar) const { return cwise() + scalar }
+const const CwiseUnaryOp<ei_scalar_add_op<Scalar>, Derived>
+operator+(const Scalar& scalar) const { return cwise() + scalar; }
 
-friend const typename Cwise<Derived>::ScalarAddReturnType
+friend const CwiseUnaryOp<ei_scalar_add_op<Scalar>, Derived>
 operator+(const Scalar& scalar, const MatrixBase<Derived>& mat) { return mat + scalar; }
 \endcode
 
@@ -81,7 +81,7 @@ In order to add support for a custom type \c T you need:
  3 - define a couple of math functions for your type such as: ei_sqrt, ei_abs, etc...
      (see the file Eigen/src/Core/MathFunctions.h)
 
-Here is a concrete example adding support for the Adolc's \c adouble type. <a href="http://www.math.tu-dresden.de/~adol-c/">Adolc</a> is an automatic differentiation library. The type \c adouble is basically a real value tracking the values of any number of partial derivatives.
+Here is a concrete example adding support for the Adolc's \c adouble type. <a href="https://projects.coin-or.org/ADOL-C">Adolc</a> is an automatic differentiation library. The type \c adouble is basically a real value tracking the values of any number of partial derivatives.
 
 \code
 #ifndef ADLOCSUPPORT_H

doc/D03_WrongStackAlignment.dox

@@ -4,6 +4,10 @@ namespace Eigen {
 
 This is an issue that, so far, we met only with GCC on Windows: for instance, MinGW and TDM-GCC.
 
+\htmlonly
+<b><big><big>It seems that this GCC bug has been <a href="http://gcc.gnu.org/gcc-4.5/changes.html">fixed in GCC 4.5</a>. We encourage all GCC/Windows users to upgrade to GCC 4.5.</big></big></b>
+\endhtmlonly
+
 By default, in a function like this,
 
 \code
@@ -16,9 +20,12 @@ void foo()
 
 GCC assumes that the stack is already 16-byte-aligned so that the object \a q will be created at a 16-byte-aligned location. For this reason, it doesn't take any special care to explicitly align the object \a q, as Eigen requires.
 
-The problem is that this assumption is wrong on Windows, where the stack is not guaranteed to have 16-byte alignment. This results in the object 'q' being created at an unaligned location, making your program crash with the \ref UnalignedArrayAssert "assertion on unaligned arrays".
+The problem is that, in some particular cases, this assumption can be wrong on Windows, where the stack is only guaranteed to have 4-byte alignment. Indeed, even though GCC takes care of aligning the stack in the main function and does its best to keep it aligned, when a function is called from another thread or from a binary compiled with another compiler, the stack alignment can be corrupted. This results in the object 'q' being created at an unaligned location, making your program crash with the \ref UnalignedArrayAssert "assertion on unaligned arrays". So far we found the three following solutions.
 
-A local solution is to mark your function with this attribute:
+
+\section sec_sol1 Local solution
+
+A local solution is to mark such a function with this attribute:
 \code
 __attribute__((force_align_arg_pointer)) void foo()
 {
@@ -26,17 +33,24 @@ __attribute__((force_align_arg_pointer)) void foo()
   //...
 }
 \endcode
-Read <a href="http://gcc.gnu.org/onlinedocs/gcc-4.4.0/gcc/Function-Attributes.html#Function-Attributes">this GCC documentation</a> to understand what this does. Of course this should only be done on GCC on Windows, so for portability you'll have to encapsulate this in a macro which you leave empty on other platforms. Also this needs to be done for every such function, which is inconvenient! So you may prefer the following global solution:
+Read <a href="http://gcc.gnu.org/onlinedocs/gcc-4.4.0/gcc/Function-Attributes.html#Function-Attributes">this GCC documentation</a> to understand what this does. Of course this should only be done on GCC on Windows, so for portability you'll have to encapsulate this in a macro which you leave empty on other platforms. The advantage of this solution is that you can finely select which function might have a corrupted stack alignment. Of course on the downside this has to be done for every such function, so you may prefer one of the following two global solutions.
+
+
+\section sec_sol2 Global solutions
 
 A global solution is to edit your project so that when compiling with GCC on Windows, you pass this option to GCC:
 \code
--mpreferred-stack-boundary=3
+-mincoming-stack-boundary=2
 \endcode
-or
+Explanation: this tells GCC that the stack is only required to be aligned to 2^2=4 bytes, so that GCC now knows that it really must take extra care to honor the 16 byte alignment of \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types" when needed.
+
+Another global solution is to pass this option to gcc:
 \code
--mpreferred-stack-boundary=2
+-mstackrealign
 \endcode
-Explanation: this tells GCC that the stack should only be required to be aligned to 2^3=8 byes (in the first version) or to 2^2=4 bytes (in the second version). In both cases, this is smaller than 16 bytes, so GCC now knows that it really must take extra care to honor the 16 byte alignment of \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types". However, you must make sure that you understand fully what this does, before applying this change to your project. Read the GCC manual page. A too low value of \c -mpreferred-stack-boundary can result in poor performance: for example, on various CPUs, double-precision numbers work much faster when aligned to 8 bytes boundary, which is guaranteed by \c -mpreferred-stack-boundary=3 but not by \c -mpreferred-stack-boundary=2, so passing \c -mpreferred-stack-boundary=2 can result in poor performance! On the other hand, notice that the higher the value of \c -mpreferred-stack-boundary, the bigger the code size.
+which has the same effect than adding the \c force_align_arg_pointer attribute to all functions.
+
+These global solutions are easy to use, but note that they may slowdown your program because they lead to extra prologue/epilogue instructions for every function.
 
 */
 

doc/Overview.dox

@@ -15,7 +15,7 @@ For a first contact with Eigen, the best place is to have a look at the \ref Tut
 
 Most of the API is available as methods in MatrixBase, so this is a good starting point for browsing. Also have a look at Matrix, as a few methods and the matrix constructors are there. Other notable classes for the Eigen API are Cwise, which contains the methods for doing certain coefficient-wise operations, and Part.
 
-In fact, except for advanced use, the only class that you'll have to explicitly name in your program, i.e. of which you'll explicitly contruct objects, is Matrix. For instance, vectors are handled as a special case of Matrix with one column. Typedefs are provided, e.g. Vector2f is a typedef for Matrix<float, 2, 1>. Finally, you might also have look at the \ref ExampleList "the list of selected examples".
+In fact, except for advanced use, the only class that you'll have to explicitly name in your program, i.e. of which you'll explicitly contruct objects, is Matrix. For instance, vectors are handled as a special case of Matrix with one column. Typedefs are provided, e.g. Vector2f is a typedef for Matrix<float, 2, 1>.
 
 Most of the other classes are just return types for MatrixBase methods.
 

doc/QuickStartGuide.dox

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

doc/StlContainers.dox

@@ -6,10 +6,12 @@ namespace Eigen {
   - \ref summary
   - \ref allocator
   - \ref vector
+  - \ref newvector
+
 
 \section summary Executive summary
 
-Using STL containers on \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types" requires taking the following two steps:
+Using STL containers on \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types", or classes having members of such types, requires taking the following two steps:
 
 \li A 16-byte-aligned allocator must be used. Eigen does provide one ready for use: aligned_allocator.
 \li If you want to use the std::vector container, you need to \#include <Eigen/StdVector>.
@@ -44,6 +46,18 @@ std::vector<Eigen::Vector4f>
 \endcode
 without having to worry about anything.
 
+\section newvector The new and improved workaround for std::vector
+
+Well, except that in Eigen 2.0 the <Eigen/StdVector> header causes some compatibility issues as it reimplements the std::vector<T> container in a not-fully-compatible way. If you want to avoid these issues, starting in Eigen 2.0.6 a new implementation is available, which will become default in the next major version of Eigen. You can enable it by defining EIGEN_USE_NEW_STDVECTOR:
+
+\code
+#define EIGEN_USE_NEW_STDVECTOR
+#include<Eigen/StdVector>
+\endcode
+
+This new implementation <a href="http://eigen.tuxfamily.org/dox-devel/StlContainers.html#vector">is documented here</a>. In particular, note that if you use it, you must specify Eigen::aligned_allocator<T> as the allocator type, otherwise it doesn't make any difference from the standard std::vector. This new std::vector implementation \b only overrides the standard one if used with this allocator, which guarantees that it doesn't break existing non-Eigen code.
+
+
 <span class="note">\b Explanation: The resize() method of std::vector takes a value_type argument (defaulting to value_type()). So with std::vector<Eigen::Vector4f>, some Eigen::Vector4f objects will be passed by value, which discards any alignment modifiers, so a Eigen::Vector4f can be created at an unaligned location. In order to avoid that, the only solution we saw was to specialize std::vector to make it work on a slight modification of, here, Eigen::Vector4f, that is able to deal properly with this situation.
 </span>
 

doc/StructHavingEigenMembers.dox

@@ -88,7 +88,7 @@ The solution is to let class Foo have an aligned "operator new", as we showed in
 
 \section movetotop Should I then put all the members of Eigen types at the beginning of my class?
 
-No, that's not needed. Since Eigen takes care of declaring 128-bit alignment, all members that need it are automatically 128-bit aligned relatively to the class. So when you have code like
+That's not required. Since Eigen takes care of declaring 128-bit alignment, all members that need it are automatically 128-bit aligned relatively to the class. So code like this works fine:
 
 \code
 class Foo
@@ -100,25 +100,13 @@ public:
 };
 \endcode
 
-it will work just fine. You do \b not need to rewrite it as
-
-\code
-class Foo
-{
-  Eigen::Vector2d v;
-  double x;
-public:
-  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
-};
-\endcode
-
 \section dynamicsize What about dynamic-size matrices and vectors?
 
-Dynamic-size matrices and vectors, such as Eigen::VectorXd, allocate dynamically their own array of coefficients, so they take care of requiring absolute alignment automatically. So they don't cause this issue. The issue discussed here is only with fixed-size matrices and vectors.
+Dynamic-size matrices and vectors, such as Eigen::VectorXd, allocate dynamically their own array of coefficients, so they take care of requiring absolute alignment automatically. So they don't cause this issue. The issue discussed here is only with \ref FixedSizeVectorizable  "fixed-size vectorizable matrices and vectors".
 
 \section bugineigen So is this a bug in Eigen?
 
-No, it's not our bug. It's more like an inherent problem of the C++ language -- though it must be said that any other existing language probably has the same problem. The problem is that there is no way that you can specify an aligned "operator new" that would propagate to classes having you as member data.
+No, it's not our bug. It's more like an inherent problem of the C++98 language specification, and seems to be taken care of in the upcoming language revision: <a href="http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2007/n2341.pdf">see this document</a>.
 
 \section conditional What if I want to do this conditionnally (depending on template parameters) ?
 

doc/TutorialGeometry.dox

@@ -152,7 +152,7 @@ OpenGL compatibility \b 3D </td><td>\code
 glLoadMatrixf(t.data());\endcode</td></tr>
 <tr><td>
 OpenGL compatibility \b 2D </td><td>\code
-Transform3f aux(Transform3f::Identity);
+Transform3f aux(Transform3f::Identity());
 aux.linear().corner<2,2>(TopLeft) = t.linear();
 aux.translation().start<2>() = t.translation();
 glLoadMatrixf(aux.data());\endcode</td></tr>

doc/UnalignedArrayAssert.dox

@@ -12,14 +12,27 @@ is explained here: http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html
 **** READ THIS WEB PAGE !!! ****"' failed.
 </pre>
 
-There are 3 known causes for this issue. Please read on to understand them and learn how to fix them.
+There are 4 known causes for this issue. Please read on to understand them and learn how to fix them.
 
 \b Table \b of \b contents
+ - \ref where
  - \ref c1
  - \ref c2
  - \ref c3
  - \ref c4
  - \ref explanation
+ - \ref getrid
+
+\section where Where in my own code is the cause of the problem?
+
+First of all, you need to find out where in your own code this assertion was triggered from. At first glance, the error message doesn't look helpful, as it refers to a file inside Eigen! However, since your program crashed, if you can reproduce the crash, you can get a backtrace using any debugger. For example, if you're using GCC, you can use the GDB debugger as follows:
+\code
+$ gdb ./my_program          # Start GDB on your program
+> run                       # Start running your program
+...                         # Now reproduce the crash!
+> bt                        # Obtain the backtrace
+\endcode
+Now that you know precisely where in your own code the problem is happening, read on to understand what you need to change.
 
 \section c1 Cause 1: Structures having Eigen objects as members
 
@@ -42,16 +55,17 @@ Note that here, Eigen::Vector2d is only used as an example, more generally the i
 
 \section c2 Cause 2: STL Containers
 
-If you use STL Containers such as std::vector, std::map, ..., with Eigen objects, like this,
+If you use STL Containers such as std::vector, std::map, ..., with Eigen objects, or with classes containing Eigen objects, like this,
 
 \code
 std::vector<Eigen::Matrix2f> my_vector;
-std::map<int, Eigen::Matrix2f> my_map;
+struct my_class { ... Eigen::Matrix2f m; ... };
+std::map<int, my_class> my_map;
 \endcode
 
 then you need to read this separate page: \ref StlContainers "Using STL Containers with Eigen".
 
-Note that here, Eigen::Matrix2f is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types".
+Note that here, Eigen::Matrix2f is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types" and \ref StructHavingEigenMembers "structures having such Eigen objects as member".
 
 \section c3 Cause 3: Passing Eigen objects by value
 
@@ -89,6 +103,16 @@ Eigen normally takes care of these alignment issues for you, by setting an align
 
 However there are a few corner cases where these alignment settings get overridden: they are the possible causes for this assertion.
 
+\section getrid I don't care about vectorization, how do I get rid of that stuff?
+
+Two possibilities:
+<ul>
+  <li>Define EIGEN_DONT_ALIGN. That disables all 128-bit alignment code, and in particular everything vectorization-related. But do note that this in particular breaks ABI compatibility with vectorized code. This requires Eigen 2.0.6 or later, and with Eigen 2.0.12 or later it automatically implies EIGEN_DONT_VECTORIZE (before 2.0.12 you had to define both).</li>
+  <li>Or define both EIGEN_DONT_VECTORIZE and EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT. This keeps the 128-bit alignment code and thus preserves ABI compatibility.</li>
+</ul>
+
+For more information, see <a href="http://eigen.tuxfamily.org/index.php?title=FAQ#I_disabled_vectorization.2C_but_I.27m_still_getting_annoyed_about_alignment_issues.21">this FAQ</a>.
+
 */
 
 }

doc/examples/CMakeLists.txt

@@ -5,6 +5,9 @@ ADD_CUSTOM_TARGET(all_examples)
 FOREACH(example_src ${examples_SRCS})
 GET_FILENAME_COMPONENT(example ${example_src} NAME_WE)
 ADD_EXECUTABLE(${example} ${example_src})
+if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
+  target_link_libraries(${example} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO})
+endif()
 GET_TARGET_PROPERTY(example_executable
                     ${example} LOCATION)
 ADD_CUSTOM_COMMAND(

doc/snippets/CMakeLists.txt

@@ -11,6 +11,9 @@ CONFIGURE_FILE(${CMAKE_CURRENT_SOURCE_DIR}/compile_snippet.cpp.in
                ${CMAKE_CURRENT_BINARY_DIR}/${compile_snippet_src})
 ADD_EXECUTABLE(${compile_snippet_target}
                ${CMAKE_CURRENT_BINARY_DIR}/${compile_snippet_src})
+if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
+  target_link_libraries(${compile_snippet_target} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO})
+endif()
 GET_TARGET_PROPERTY(compile_snippet_executable
                     ${compile_snippet_target} LOCATION)
 ADD_CUSTOM_COMMAND(

eigen2.pc.in

@@ -0,0 +1,7 @@
+
+Name: Eigen2
+Description: A C++ template library for linear algebra: vectors, matrices, and related algorithms
+Requires:
+Version: ${EIGEN_VERSION_NUMBER}
+Libs:
+Cflags: -I${INCLUDE_INSTALL_DIR}

scripts/eigen_gen_docs

@@ -0,0 +1,19 @@
+#!/bin/sh
+
+# configuration
+# You should call this script with USER set as you want, else some default
+# will be used
+USER=${USER:-'orzel'}
+
+# step 1 : build
+# todo if 'build is not there, create one:
+#mkdir build
+(cd build && cmake .. && make -j3 doc) || { echo "make failed"; exit 1; }
+#todo: n+1 where n = number of cpus
+
+#step 2 : upload
+# (the '/' at the end of path are very important, see rsync documentation)
+rsync -az build/doc/html/ $USER@ssh.tuxfamily.org:eigen/eigen.tuxfamily.org-web/htdocs/dox-2.0/ || { echo "upload failed"; exit 1; }
+
+echo "Uploaded successfully"
+

test/CMakeLists.txt

@@ -1,3 +1,7 @@
+option(EIGEN_DEFAULT_TO_ROW_MAJOR "Use row-major as default matrix storage order" OFF)
+if(EIGEN_DEFAULT_TO_ROW_MAJOR)
+  add_definitions("-DEIGEN_DEFAULT_TO_ROW_MAJOR")
+endif()
 
 find_package(GSL)
 if(GSL_FOUND AND GSL_VERSION_MINOR LESS 9)
@@ -34,7 +38,7 @@ else(CHOLMOD_FOUND)
   set(EIGEN_MISSING_BACKENDS ${EIGEN_MISSING_BACKENDS} Cholmod)
 endif(CHOLMOD_FOUND)
 
-option(EIGEN_TEST_NO_FORTRAN "Disable Fortran" OFF)
+option(EIGEN_TEST_NO_FORTRAN "Disable Fortran" ON)
 if(NOT MSVC AND NOT EIGEN_TEST_NO_FORTRAN)
   enable_language(Fortran OPTIONAL)
 endif(NOT MSVC AND NOT EIGEN_TEST_NO_FORTRAN)
@@ -93,12 +97,20 @@ else(CMAKE_COMPILER_IS_GNUCXX)
 endif(CMAKE_COMPILER_IS_GNUCXX)
 
 option(EIGEN_NO_ASSERTION_CHECKING "Disable checking of assertions" OFF)
+if(EIGEN_NO_ASSERTION_CHECKING)
+  add_definitions("-DEIGEN_NO_ASSERTION_CHECKING=1")
+endif()
+
 
 # similar to set_target_properties but append the property instead of overwriting it
 macro(ei_add_target_property target prop value)
 
   get_target_property(previous ${target} ${prop})
-  set_target_properties(${target} PROPERTIES ${prop} "${previous} ${value}")
+  if(previous MATCHES "NOTFOUND")
+    set_target_properties(${target} PROPERTIES ${prop} "${value}")
+  else()
+    set_target_properties(${target} PROPERTIES ${prop} "${previous} ${value}")
+  endif()
 
 endmacro(ei_add_target_property)
 
@@ -134,13 +146,9 @@ macro(ei_add_test testname)
 
     option(EIGEN_DEBUG_ASSERTS "Enable debuging of assertions" OFF)
     if(EIGEN_DEBUG_ASSERTS)
-      set_target_properties(${targetname} PROPERTIES COMPILE_DEFINITIONS "-DEIGEN_DEBUG_ASSERTS=1")
+      set_target_properties(${targetname} PROPERTIES COMPILE_DEFINITIONS "EIGEN_DEBUG_ASSERTS=1")
     endif(EIGEN_DEBUG_ASSERTS)
 
-  else(NOT EIGEN_NO_ASSERTION_CHECKING)
-
-    set_target_properties(${targetname} PROPERTIES COMPILE_DEFINITIONS "-DEIGEN_NO_ASSERTION_CHECKING=1")
-
   endif(NOT EIGEN_NO_ASSERTION_CHECKING)
 
   if(${ARGC} GREATER 1)
@@ -153,6 +161,10 @@ macro(ei_add_test testname)
     target_link_libraries(${targetname} Eigen2)
   endif(TEST_LIB)
 
+  if(EIGEN_STANDARD_LIBRARIES_TO_LINK_TO)
+    target_link_libraries(${targetname} ${EIGEN_STANDARD_LIBRARIES_TO_LINK_TO})
+  endif()
+
   target_link_libraries(${targetname} ${EXTERNAL_LIBS})
   if(${ARGC} GREATER 2)
     string(STRIP "${ARGV2}" ARGV2_stripped)
@@ -180,6 +192,7 @@ ei_add_test(meta)
 ei_add_test(sizeof)
 ei_add_test(dynalloc)
 ei_add_test(nomalloc)
+ei_add_test(first_aligned)
 ei_add_test(mixingtypes)
 ei_add_test(packetmath)
 ei_add_test(unalignedassert)
@@ -200,7 +213,7 @@ ei_add_test(array)
 ei_add_test(triangular)
 ei_add_test(cholesky " " "${GSL_LIBRARIES}")
 ei_add_test(lu ${EI_OFLAG})
-ei_add_test(determinant)
+ei_add_test(determinant ${EI_OFLAG})
 ei_add_test(inverse)
 ei_add_test(qr)
 ei_add_test(eigensolver " " "${GSL_LIBRARIES}")
@@ -211,13 +224,21 @@ ei_add_test(parametrizedline)
 ei_add_test(alignedbox)
 ei_add_test(regression)
 ei_add_test(stdvector)
+ei_add_test(newstdvector)
 if(QT4_FOUND)
-  ei_add_test(qtvector " " ${QT_QTCORE_LIBRARY})
+  ei_add_test(qtvector " " "${QT_QTCORE_LIBRARY}")
 endif(QT4_FOUND)
-ei_add_test(sparse_vector)
-ei_add_test(sparse_basic)
-ei_add_test(sparse_solvers " " "${SPARSE_LIBS}")
-ei_add_test(sparse_product)
+if(NOT EIGEN_DEFAULT_TO_ROW_MAJOR)
+  ei_add_test(sparse_vector)
+  ei_add_test(sparse_basic)
+  ei_add_test(sparse_solvers " " "${SPARSE_LIBS}")
+  ei_add_test(sparse_product)
+endif()
+ei_add_test(swap)
+ei_add_test(visitor)
+ei_add_test(bug_132)
+
+ei_add_test(prec_inverse_4x4 ${EI_OFLAG})
 
 # print a summary of the different options
 message("************************************************************")
@@ -262,11 +283,22 @@ else(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
   message("Explicit vec:      AUTO")
 endif(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
 
+if(EIGEN_DEFAULT_TO_ROW_MAJOR)
+  message("Default order:     Row-major")
+else()
+  message("Default order:     Column-major")
+endif()
+
+string(TOLOWER "${CMAKE_CXX_COMPILER}" cmake_cxx_compiler_tolower)
+if(cmake_cxx_compiler_tolower MATCHES "qcc")
+  set(CXX_IS_QCC "ON")
+endif()
+
   message("CXX:               ${CMAKE_CXX_COMPILER}")
-if(CMAKE_COMPILER_IS_GNUCXX)
+if(CMAKE_COMPILER_IS_GNUCXX AND NOT CXX_IS_QCC)
   execute_process(COMMAND ${CMAKE_CXX_COMPILER} --version COMMAND head -n 1 OUTPUT_VARIABLE EIGEN_CXX_VERSION_STRING OUTPUT_STRIP_TRAILING_WHITESPACE)
   message("CXX_VERSION:       ${EIGEN_CXX_VERSION_STRING}")
-endif(CMAKE_COMPILER_IS_GNUCXX)
+endif()
   message("CXX_FLAGS:         ${CMAKE_CXX_FLAGS}")
   message("Sparse lib flags:  ${SPARSE_LIBS}")
 

test/bug_132.cpp

@@ -0,0 +1,41 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+
+void test_bug_132() {
+  enum { size = 100 };
+  MatrixXd A(size, size);
+  VectorXd b(size), c(size);
+  {
+    VectorXd y = A.transpose() * (b-c);  // bug 132: infinite recursion in coeffRef
+    VectorXd z = (b-c).transpose() * A;  // bug 132: infinite recursion in coeffRef
+  }
+
+  // the following ones weren't failing, but let's include them for completeness:
+  {
+    VectorXd y = A * (b-c);
+    VectorXd z = (b-c).transpose() * A.transpose();
+  }
+}

test/cholesky.cpp

@@ -100,6 +100,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     VERIFY_IS_APPROX(symm * matX, matB);
   }
 
+#if 0 // cholesky is not rank-revealing anyway
   // test isPositiveDefinite on non definite matrix
   if (rows>4)
   {
@@ -109,6 +110,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m)
     LDLT<SquareMatrixType> cholnosqrt(symm);
     VERIFY(!cholnosqrt.isPositiveDefinite());
   }
+#endif
 }
 
 void test_cholesky()
@@ -122,4 +124,8 @@ void test_cholesky()
     CALL_SUBTEST( cholesky(MatrixXf(17,17)) );
     CALL_SUBTEST( cholesky(MatrixXd(33,33)) );
   }
+
+  MatrixXf m = MatrixXf::Zero(10,10);
+  VectorXf b = VectorXf::Zero(10);
+  VERIFY(!m.llt().isPositiveDefinite());
 }

test/dynalloc.cpp

@@ -35,7 +35,7 @@ void check_handmade_aligned_malloc()
   for(int i = 1; i < 1000; i++)
   {
     char *p = (char*)ei_handmade_aligned_malloc(i);
-    VERIFY(size_t(p)%ALIGNMENT==0);
+    VERIFY(std::size_t(p)%ALIGNMENT==0);
     // if the buffer is wrongly allocated this will give a bad write --> check with valgrind
     for(int j = 0; j < i; j++) p[j]=0;
     ei_handmade_aligned_free(p);
@@ -47,7 +47,7 @@ void check_aligned_malloc()
   for(int i = 1; i < 1000; i++)
   {
     char *p = (char*)ei_aligned_malloc(i);
-    VERIFY(size_t(p)%ALIGNMENT==0);
+    VERIFY(std::size_t(p)%ALIGNMENT==0);
     // if the buffer is wrongly allocated this will give a bad write --> check with valgrind
     for(int j = 0; j < i; j++) p[j]=0;
     ei_aligned_free(p);
@@ -59,7 +59,7 @@ void check_aligned_new()
   for(int i = 1; i < 1000; i++)
   {
     float *p = ei_aligned_new<float>(i);
-    VERIFY(size_t(p)%ALIGNMENT==0);
+    VERIFY(std::size_t(p)%ALIGNMENT==0);
     // if the buffer is wrongly allocated this will give a bad write --> check with valgrind
     for(int j = 0; j < i; j++) p[j]=0;
     ei_aligned_delete(p,i);
@@ -71,7 +71,7 @@ void check_aligned_stack_alloc()
   for(int i = 1; i < 1000; i++)
   {
     float *p = ei_aligned_stack_new(float,i);
-    VERIFY(size_t(p)%ALIGNMENT==0);
+    VERIFY(std::size_t(p)%ALIGNMENT==0);
     // if the buffer is wrongly allocated this will give a bad write --> check with valgrind
     for(int j = 0; j < i; j++) p[j]=0;
     ei_aligned_stack_delete(float,p,i);
@@ -98,7 +98,7 @@ class MyClassA
 template<typename T> void check_dynaligned()
 {
   T* obj = new T;
-  VERIFY(size_t(obj)%ALIGNMENT==0);
+  VERIFY(std::size_t(obj)%ALIGNMENT==0);
   delete obj;
 }
 
@@ -121,15 +121,15 @@ void test_dynalloc()
   
   // check static allocation, who knows ?
   {
-    MyStruct foo0;  VERIFY(size_t(foo0.avec.data())%ALIGNMENT==0);
-    MyClassA fooA;  VERIFY(size_t(fooA.avec.data())%ALIGNMENT==0);
+    MyStruct foo0;  VERIFY(std::size_t(foo0.avec.data())%ALIGNMENT==0);
+    MyClassA fooA;  VERIFY(std::size_t(fooA.avec.data())%ALIGNMENT==0);
   }
 
   // dynamic allocation, single object
   for (int i=0; i<g_repeat*100; ++i)
   {
-    MyStruct *foo0 = new MyStruct();  VERIFY(size_t(foo0->avec.data())%ALIGNMENT==0);
-    MyClassA *fooA = new MyClassA();  VERIFY(size_t(fooA->avec.data())%ALIGNMENT==0);
+    MyStruct *foo0 = new MyStruct();  VERIFY(std::size_t(foo0->avec.data())%ALIGNMENT==0);
+    MyClassA *fooA = new MyClassA();  VERIFY(std::size_t(fooA->avec.data())%ALIGNMENT==0);
     delete foo0;
     delete fooA;
   }
@@ -138,8 +138,8 @@ void test_dynalloc()
   const int N = 10;
   for (int i=0; i<g_repeat*100; ++i)
   {
-    MyStruct *foo0 = new MyStruct[N];  VERIFY(size_t(foo0->avec.data())%ALIGNMENT==0);
-    MyClassA *fooA = new MyClassA[N];  VERIFY(size_t(fooA->avec.data())%ALIGNMENT==0);
+    MyStruct *foo0 = new MyStruct[N];  VERIFY(std::size_t(foo0->avec.data())%ALIGNMENT==0);
+    MyClassA *fooA = new MyClassA[N];  VERIFY(std::size_t(fooA->avec.data())%ALIGNMENT==0);
     delete[] foo0;
     delete[] fooA;
   }

test/first_aligned.cpp

@@ -0,0 +1,64 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+
+template<typename Scalar>
+void test_first_aligned_helper(Scalar *array, int size)
+{
+  const int packet_size = sizeof(Scalar) * ei_packet_traits<Scalar>::size;
+  VERIFY(((std::size_t(array) + sizeof(Scalar) * ei_alignmentOffset(array, size)) % packet_size) == 0);
+}
+
+template<typename Scalar>
+void test_none_aligned_helper(Scalar *array, int size)
+{
+  VERIFY(ei_packet_traits<Scalar>::size == 1 || ei_alignmentOffset(array, size) == size);
+}
+
+struct some_non_vectorizable_type { float x; };
+
+void test_first_aligned()
+{
+  EIGEN_ALIGN_128 float array_float[100];
+  test_first_aligned_helper(array_float, 50);
+  test_first_aligned_helper(array_float+1, 50);
+  test_first_aligned_helper(array_float+2, 50);
+  test_first_aligned_helper(array_float+3, 50);
+  test_first_aligned_helper(array_float+4, 50);
+  test_first_aligned_helper(array_float+5, 50);
+  
+  EIGEN_ALIGN_128 double array_double[100];
+  test_first_aligned_helper(array_double, 50);
+  test_first_aligned_helper(array_double+1, 50);
+  test_first_aligned_helper(array_double+2, 50);
+  
+  double *array_double_plus_4_bytes = (double*)(std::size_t(array_double)+4);
+  test_none_aligned_helper(array_double_plus_4_bytes, 50);
+  test_none_aligned_helper(array_double_plus_4_bytes+1, 50);
+  
+  some_non_vectorizable_type array_nonvec[100];
+  test_first_aligned_helper(array_nonvec, 100);
+  test_none_aligned_helper(array_nonvec, 100);
+}

test/geometry.cpp

@@ -151,7 +151,13 @@ template<typename Scalar> void geometry(void)
     a = ei_random<Scalar>(-Scalar(0.4)*Scalar(M_PI), Scalar(0.4)*Scalar(M_PI));
   q1 = AngleAxisx(a, v0.normalized());
   Transform3 t0, t1, t2;
+  // first test setIdentity() and Identity()
   t0.setIdentity();
+  VERIFY_IS_APPROX(t0.matrix(), Transform3::MatrixType::Identity());
+  t0.matrix().setZero();
+  t0 = Transform3::Identity();
+  VERIFY_IS_APPROX(t0.matrix(), Transform3::MatrixType::Identity());
+
   t0.linear() = q1.toRotationMatrix();
   t1.setIdentity();
   t1.linear() = q1.toRotationMatrix();

test/hyperplane.cpp

@@ -121,7 +121,8 @@ template<typename Scalar> void lines()
     VERIFY_IS_APPROX(result, center);
 
     // check conversions between two types of lines
-    CoeffsType converted_coeffs(HLine(PLine(line_u)).coeffs());
+    PLine pl(line_u); // gcc 3.3 will commit suicide if we don't name this variable
+    CoeffsType converted_coeffs(HLine(pl).coeffs());
     converted_coeffs *= line_u.coeffs()(0)/converted_coeffs(0);
     VERIFY(line_u.coeffs().isApprox(converted_coeffs));
   }

test/inverse.cpp

@@ -43,11 +43,9 @@ template<typename MatrixType> void inverse(const MatrixType& m)
              mzero = MatrixType::Zero(rows, cols),
              identity = MatrixType::Identity(rows, rows);
 
-  if (ei_is_same_type<RealScalar,float>::ret)
+  while(ei_abs(m1.determinant()) < RealScalar(0.1) && rows <= 8)
   {
-    // let's build a more stable to inverse matrix
-    MatrixType a = MatrixType::Random(rows,cols);
-    m1 += m1 * m1.adjoint() + a * a.adjoint();
+    m1 = MatrixType::Random(rows, cols);
   }
 
   m2 = m1.inverse();

test/lu.cpp

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

test/main.h

@@ -142,7 +142,7 @@ namespace Eigen
 #define VERIFY(a) do { if (!(a)) { \
     std::cerr << "Test " << g_test_stack.back() << " failed in "EI_PP_MAKE_STRING(__FILE__) << " (" << EI_PP_MAKE_STRING(__LINE__) << ")" \
       << std::endl << "    " << EI_PP_MAKE_STRING(a) << std::endl << std::endl; \
-    exit(2); \
+    std::exit(2); \
   } } while (0)
 
 #define VERIFY_IS_APPROX(a, b) VERIFY(test_ei_isApprox(a, b))
@@ -257,7 +257,7 @@ int main(int argc, char *argv[])
           std::cout << "Argument " << argv[i] << " conflicting with a former argument" << std::endl;
           return 1;
         }
-        repeat = atoi(argv[i]+1);
+        repeat = std::atoi(argv[i]+1);
         has_set_repeat = true;
         if(repeat <= 0)
         {
@@ -272,7 +272,7 @@ int main(int argc, char *argv[])
           std::cout << "Argument " << argv[i] << " conflicting with a former argument" << std::endl;
           return 1;
         }
-        seed = int(strtoul(argv[i]+1, 0, 10));
+        seed = int(std::strtoul(argv[i]+1, 0, 10));
         has_set_seed = true;
         bool ok = seed!=0;
         if(!ok)
@@ -295,11 +295,11 @@ int main(int argc, char *argv[])
       return 1;
     }
 
-    if(!has_set_seed) seed = (unsigned int) time(NULL);
+    if(!has_set_seed) seed = (unsigned int) std::time(NULL);
     if(!has_set_repeat) repeat = DEFAULT_REPEAT;
 
     std::cout << "Initializing random number generator with seed " << seed << std::endl;
-    srand(seed);
+	std::srand(seed);
     std::cout << "Repeating each test " << repeat << " times" << std::endl;
 
     Eigen::g_repeat = repeat;

test/map.cpp

@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra. Eigen itself is part of the KDE project.
 //
-// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2007-2010 Benoit Jacob <jacob.benoit.1@gmail.com>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -24,7 +24,7 @@
 
 #include "main.h"
 
-template<typename VectorType> void map_class(const VectorType& m)
+template<typename VectorType> void map_class_vector(const VectorType& m)
 {
   typedef typename VectorType::Scalar Scalar;
 
@@ -34,7 +34,7 @@ template<typename VectorType> void map_class(const VectorType& m)
   Scalar* array1 = ei_aligned_new<Scalar>(size);
   Scalar* array2 = ei_aligned_new<Scalar>(size);
   Scalar* array3 = new Scalar[size+1];
-  Scalar* array3unaligned = size_t(array3)%16 == 0 ? array3+1 : array3;
+  Scalar* array3unaligned = std::size_t(array3)%16 == 0 ? array3+1 : array3;
   
   Map<VectorType, Aligned>(array1, size) = VectorType::Random(size);
   Map<VectorType>(array2, size) = Map<VectorType>(array1, size);
@@ -50,6 +50,34 @@ template<typename VectorType> void map_class(const VectorType& m)
   delete[] array3;
 }
 
+template<typename MatrixType> void map_class_matrix(const MatrixType& m)
+{
+  typedef typename MatrixType::Scalar Scalar;
+
+  int rows = m.rows(), cols = m.cols(), size = rows*cols;
+
+  // test Map.h
+  Scalar* array1 = ei_aligned_new<Scalar>(size);
+  for(int i = 0; i < size; i++) array1[i] = Scalar(1);
+  Scalar* array2 = ei_aligned_new<Scalar>(size);
+  for(int i = 0; i < size; i++) array2[i] = Scalar(1);
+  Scalar* array3 = new Scalar[size+1];
+  for(int i = 0; i < size+1; i++) array3[i] = Scalar(1);
+  Scalar* array3unaligned = std::size_t(array3)%16 == 0 ? array3+1 : array3;
+  Map<MatrixType, Aligned>(array1, rows, cols) = MatrixType::Ones(rows,cols);
+  Map<MatrixType>(array2, rows, cols) = Map<MatrixType>(array1, rows, cols);
+  Map<MatrixType>(array3unaligned, rows, cols) = Map<MatrixType>(array1, rows, cols);
+  MatrixType ma1 = Map<MatrixType>(array1, rows, cols);
+  MatrixType ma2 = Map<MatrixType, Aligned>(array2, rows, cols);
+  VERIFY_IS_APPROX(ma1, ma2);
+  MatrixType ma3 = Map<MatrixType>(array3unaligned, rows, cols);
+  VERIFY_IS_APPROX(ma1, ma3);
+  
+  ei_aligned_delete(array1, size);
+  ei_aligned_delete(array2, size);
+  delete[] array3;
+}
+
 template<typename VectorType> void map_static_methods(const VectorType& m)
 {
   typedef typename VectorType::Scalar Scalar;
@@ -60,7 +88,7 @@ template<typename VectorType> void map_static_methods(const VectorType& m)
   Scalar* array1 = ei_aligned_new<Scalar>(size);
   Scalar* array2 = ei_aligned_new<Scalar>(size);
   Scalar* array3 = new Scalar[size+1];
-  Scalar* array3unaligned = size_t(array3)%16 == 0 ? array3+1 : array3;
+  Scalar* array3unaligned = std::size_t(array3)%16 == 0 ? array3+1 : array3;
   
   VectorType::MapAligned(array1, size) = VectorType::Random(size);
   VectorType::Map(array2, size) = VectorType::Map(array1, size);
@@ -80,11 +108,17 @@ template<typename VectorType> void map_static_methods(const VectorType& m)
 void test_map()
 {
   for(int i = 0; i < g_repeat; i++) {
-    CALL_SUBTEST( map_class(Matrix<float, 1, 1>()) );
-    CALL_SUBTEST( map_class(Vector4d()) );
-    CALL_SUBTEST( map_class(RowVector4f()) );
-    CALL_SUBTEST( map_class(VectorXcf(8)) );
-    CALL_SUBTEST( map_class(VectorXi(12)) );
+    CALL_SUBTEST( map_class_vector(Matrix<float, 1, 1>()) );
+    CALL_SUBTEST( map_class_vector(Vector4d()) );
+    CALL_SUBTEST( map_class_vector(RowVector4f()) );
+    CALL_SUBTEST( map_class_vector(VectorXcf(8)) );
+    CALL_SUBTEST( map_class_vector(VectorXi(12)) );
+
+    CALL_SUBTEST( map_class_matrix(Matrix<float, 1, 1>()) );
+    CALL_SUBTEST( map_class_matrix(Matrix4d()) );
+    CALL_SUBTEST( map_class_matrix(Matrix<float,3,5>()) );
+    CALL_SUBTEST( map_class_matrix(MatrixXcf(ei_random<int>(1,10),ei_random<int>(1,10))) );
+    CALL_SUBTEST( map_class_matrix(MatrixXi(ei_random<int>(1,10),ei_random<int>(1,10))) );
 
     CALL_SUBTEST( map_static_methods(Matrix<double, 1, 1>()) );
     CALL_SUBTEST( map_static_methods(Vector3f()) );

test/newstdvector.cpp

@@ -0,0 +1,164 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#define EIGEN_USE_NEW_STDVECTOR
+#include "main.h"
+#include <Eigen/StdVector>
+#include <Eigen/Geometry>
+
+template<typename MatrixType>
+void check_stdvector_matrix(const MatrixType& m)
+{
+  int rows = m.rows();
+  int cols = m.cols();
+  MatrixType x = MatrixType::Random(rows,cols), y = MatrixType::Random(rows,cols);
+  std::vector<MatrixType,Eigen::aligned_allocator<MatrixType> > v(10, MatrixType(rows,cols)), w(20, y);
+  v[5] = x;
+  w[6] = v[5];
+  VERIFY_IS_APPROX(w[6], v[5]);
+  v = w;
+  for(int i = 0; i < 20; i++)
+  {
+    VERIFY_IS_APPROX(w[i], v[i]);
+  }
+
+  v.resize(21);
+  v[20] = x;
+  VERIFY_IS_APPROX(v[20], x);
+  v.resize(22,y);
+  VERIFY_IS_APPROX(v[21], y);
+  v.push_back(x);
+  VERIFY_IS_APPROX(v[22], x);
+  VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(MatrixType));
+
+  // do a lot of push_back such that the vector gets internally resized
+  // (with memory reallocation)
+  MatrixType* ref = &w[0];
+  for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i)
+    v.push_back(w[i%w.size()]);
+  for(unsigned int i=23; i<v.size(); ++i)
+  {
+    VERIFY(v[i]==w[(i-23)%w.size()]);
+  }
+}
+
+template<typename TransformType>
+void check_stdvector_transform(const TransformType&)
+{
+  typedef typename TransformType::MatrixType MatrixType;
+  TransformType x(MatrixType::Random()), y(MatrixType::Random());
+  std::vector<TransformType,Eigen::aligned_allocator<TransformType> > v(10), w(20, y);
+  v[5] = x;
+  w[6] = v[5];
+  VERIFY_IS_APPROX(w[6], v[5]);
+  v = w;
+  for(int i = 0; i < 20; i++)
+  {
+    VERIFY_IS_APPROX(w[i], v[i]);
+  }
+
+  v.resize(21);
+  v[20] = x;
+  VERIFY_IS_APPROX(v[20], x);
+  v.resize(22,y);
+  VERIFY_IS_APPROX(v[21], y);
+  v.push_back(x);
+  VERIFY_IS_APPROX(v[22], x);
+  VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(TransformType));
+
+  // do a lot of push_back such that the vector gets internally resized
+  // (with memory reallocation)
+  TransformType* ref = &w[0];
+  for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i)
+    v.push_back(w[i%w.size()]);
+  for(unsigned int i=23; i<v.size(); ++i)
+  {
+    VERIFY(v[i].matrix()==w[(i-23)%w.size()].matrix());
+  }
+}
+
+template<typename QuaternionType>
+void check_stdvector_quaternion(const QuaternionType&)
+{
+  typedef typename QuaternionType::Coefficients Coefficients;
+  QuaternionType x(Coefficients::Random()), y(Coefficients::Random());
+  std::vector<QuaternionType,Eigen::aligned_allocator<QuaternionType> > v(10), w(20, y);
+  v[5] = x;
+  w[6] = v[5];
+  VERIFY_IS_APPROX(w[6], v[5]);
+  v = w;
+  for(int i = 0; i < 20; i++)
+  {
+    VERIFY_IS_APPROX(w[i], v[i]);
+  }
+
+  v.resize(21);
+  v[20] = x;
+  VERIFY_IS_APPROX(v[20], x);
+  v.resize(22,y);
+  VERIFY_IS_APPROX(v[21], y);
+  v.push_back(x);
+  VERIFY_IS_APPROX(v[22], x);
+  VERIFY((std::size_t)&(v[22]) == (std::size_t)&(v[21]) + sizeof(QuaternionType));
+
+  // do a lot of push_back such that the vector gets internally resized
+  // (with memory reallocation)
+  QuaternionType* ref = &w[0];
+  for(int i=0; i<30 || ((ref==&w[0]) && i<300); ++i)
+    v.push_back(w[i%w.size()]);
+  for(unsigned int i=23; i<v.size(); ++i)
+  {
+    VERIFY(v[i].coeffs()==w[(i-23)%w.size()].coeffs());
+  }
+}
+
+void test_newstdvector()
+{
+  // some non vectorizable fixed sizes
+  CALL_SUBTEST(check_stdvector_matrix(Vector2f()));
+  CALL_SUBTEST(check_stdvector_matrix(Matrix3f()));
+  CALL_SUBTEST(check_stdvector_matrix(Matrix3d()));
+
+  // some vectorizable fixed sizes
+  CALL_SUBTEST(check_stdvector_matrix(Matrix2f()));
+  CALL_SUBTEST(check_stdvector_matrix(Vector4f()));
+  CALL_SUBTEST(check_stdvector_matrix(Matrix4f()));
+  CALL_SUBTEST(check_stdvector_matrix(Matrix4d()));
+
+  // some dynamic sizes
+  CALL_SUBTEST(check_stdvector_matrix(MatrixXd(1,1)));
+  CALL_SUBTEST(check_stdvector_matrix(VectorXd(20)));
+  CALL_SUBTEST(check_stdvector_matrix(RowVectorXf(20)));
+  CALL_SUBTEST(check_stdvector_matrix(MatrixXcf(10,10)));
+
+  // some Transform
+  CALL_SUBTEST(check_stdvector_transform(Transform2f()));
+  CALL_SUBTEST(check_stdvector_transform(Transform3f()));
+  CALL_SUBTEST(check_stdvector_transform(Transform3d()));
+  //CALL_SUBTEST(check_stdvector_transform(Transform4d()));
+
+  // some Quaternion
+  CALL_SUBTEST(check_stdvector_quaternion(Quaternionf()));
+  CALL_SUBTEST(check_stdvector_quaternion(Quaterniond()));
+}

test/prec_inverse_4x4.cpp

@@ -0,0 +1,99 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+#include <Eigen/LU>
+#include <algorithm>
+
+template<typename T> std::string type_name() { return "other"; }
+template<> std::string type_name<float>() { return "float"; }
+template<> std::string type_name<double>() { return "double"; }
+template<> std::string type_name<int>() { return "int"; }
+template<> std::string type_name<std::complex<float> >() { return "complex<float>"; }
+template<> std::string type_name<std::complex<double> >() { return "complex<double>"; }
+template<> std::string type_name<std::complex<int> >() { return "complex<int>"; }
+
+#define EIGEN_DEBUG_VAR(x) std::cerr << #x << " = " << x << std::endl;
+
+template<typename T> inline typename NumTraits<T>::Real epsilon()
+{
+ return std::numeric_limits<typename NumTraits<T>::Real>::epsilon();
+}
+
+template<typename MatrixType> void inverse_permutation_4x4()
+{
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  Vector4i indices(0,1,2,3);
+  for(int i = 0; i < 24; ++i)
+  {
+    MatrixType m = MatrixType::Zero();
+    m(indices(0),0) = 1;
+    m(indices(1),1) = 1;
+    m(indices(2),2) = 1;
+    m(indices(3),3) = 1;
+    MatrixType inv = m.inverse();
+    double error = double( (m*inv-MatrixType::Identity()).norm() / epsilon<Scalar>() );
+    VERIFY(error == 0.0);
+    std::next_permutation(indices.data(),indices.data()+4);
+  }
+}
+
+template<typename MatrixType> void inverse_general_4x4(int repeat)
+{
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  double error_sum = 0., error_max = 0.;
+  for(int i = 0; i < repeat; ++i)
+  {
+    MatrixType m;
+    RealScalar absdet;
+    do {
+      m = MatrixType::Random();
+      absdet = ei_abs(m.determinant());
+    } while(absdet < 10 * epsilon<Scalar>());
+    MatrixType inv = m.inverse();
+    double error = double( (m*inv-MatrixType::Identity()).norm() * absdet / epsilon<Scalar>() );
+    error_sum += error;
+    error_max = std::max(error_max, error);
+  }
+  std::cerr << "inverse_general_4x4, Scalar = " << type_name<Scalar>() << std::endl;
+  double error_avg = error_sum / repeat;
+  EIGEN_DEBUG_VAR(error_avg);
+  EIGEN_DEBUG_VAR(error_max);
+  VERIFY(error_avg < (NumTraits<Scalar>::IsComplex ? 8.0 : 1.0));
+  VERIFY(error_max < (NumTraits<Scalar>::IsComplex ? 64.0 : 20.0));
+}
+
+void test_prec_inverse_4x4()
+{
+  CALL_SUBTEST((inverse_permutation_4x4<Matrix4f>()));
+  CALL_SUBTEST(( inverse_general_4x4<Matrix4f>(200000 * g_repeat) ));
+
+  CALL_SUBTEST((inverse_permutation_4x4<Matrix<double,4,4,RowMajor> >()));
+  CALL_SUBTEST(( inverse_general_4x4<Matrix<double,4,4,RowMajor> >(200000 * g_repeat) ));
+
+  CALL_SUBTEST((inverse_permutation_4x4<Matrix4cf>()));
+  CALL_SUBTEST((inverse_general_4x4<Matrix4cf>(50000 * g_repeat)));
+}

test/product.h

@@ -46,7 +46,7 @@ template<typename MatrixType> void product(const MatrixType& m)
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> RowSquareMatrixType;
   typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> ColSquareMatrixType;
   typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
-                         MatrixType::Flags&RowMajorBit> OtherMajorMatrixType;
+                         MatrixType::Options^RowMajor> OtherMajorMatrixType;
 
   int rows = m.rows();
   int cols = m.cols();
@@ -77,6 +77,7 @@ template<typename MatrixType> void product(const MatrixType& m)
 
   // begin testing Product.h: only associativity for now
   // (we use Transpose.h but this doesn't count as a test for it)
+
   VERIFY_IS_APPROX((m1*m1.transpose())*m2,  m1*(m1.transpose()*m2));
   m3 = m1;
   m3 *= m1.transpose() * m2;
@@ -137,6 +138,7 @@ template<typename MatrixType> void product(const MatrixType& m)
   res2 = square2;
   res2 += (m1.transpose() * m2).lazy();
   VERIFY_IS_APPROX(res2, square2 + m1.transpose() * m2);
+
   if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1)
   {
     VERIFY(areNotApprox(res2,square2 + m2.transpose() * m1));

test/product_large.cpp

@@ -42,4 +42,17 @@ void test_product_large()
     m = (v+v).asDiagonal() * m;
     VERIFY_IS_APPROX(m, MatrixXf::Constant(N,3,2));
   }
+
+  {
+    // test deferred resizing in Matrix::operator=
+    MatrixXf a = MatrixXf::Random(10,4), b = MatrixXf::Random(4,10), c = a;
+    VERIFY_IS_APPROX((a = a * b), (c * b).eval());
+  }
+
+  {
+    MatrixXf mat1(10,10); mat1.setRandom();
+    MatrixXf mat2(32,10); mat2.setRandom();
+    MatrixXf result = mat1.row(2)*mat2.transpose();
+    VERIFY_IS_APPROX(result, (mat1.row(2)*mat2.transpose()).eval());
+  }
 }

test/qr.cpp

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

test/qtvector.cpp

@@ -26,10 +26,12 @@
 #define EIGEN_WORK_AROUND_QT_BUG_CALLING_WRONG_OPERATOR_NEW_FIXED_IN_QT_4_5
 
 #include "main.h"
-#include <QtCore/QVector>
+
 #include <Eigen/Geometry>
 #include <Eigen/QtAlignedMalloc>
 
+#include <QtCore/QVector>
+
 template<typename MatrixType>
 void check_qtvector_matrix(const MatrixType& m)
 {