bump that too

that eigen2 owes all its performance to nonstandard restrict keyword.
well, this can also improve portability in case some compiler doesn't have __restrict.

.hgignore

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

CMakeLists.txt

@@ -1,22 +1,9 @@
 project(Eigen)
-set(EIGEN_VERSION_NUMBER "2.0-rc1")
-
-#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(EIGEN_VERSION_NUMBER "2.0.6")
+set(EIGEN_VERSION "${EIGEN_VERSION_NUMBER}")
+
 set(CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)
 
 option(EIGEN_BUILD_TESTS "Build tests" OFF)
@@ -85,6 +72,7 @@ endif(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
 include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR})
 
 add_subdirectory(Eigen)
+add_subdirectory(unsupported)
 
 if(EIGEN_BUILD_TESTS)
   include(CTest)

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/Cholesky

@@ -17,6 +17,9 @@
 namespace Eigen {
 
 /** \defgroup Cholesky_Module Cholesky module
+  *
+  * \nonstableyet
+  *
   * This module provides two variants of the Cholesky decomposition for selfadjoint (hermitian) matrices.
   * Those decompositions are accessible via the following MatrixBase methods:
   *  - MatrixBase::llt(),
@@ -35,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

@@ -7,11 +7,10 @@
 #ifdef _MSC_VER
   #include <malloc.h> // for _aligned_malloc -- need it regardless of whether vectorization is enabled
   #if (_MSC_VER >= 1500) // 2008 or later
-    // Remember that usage of defined() in a #define is undefined by the standard
-    #ifdef _M_IX86_FP
-      #if _M_IX86_FP >= 2
-        #define EIGEN_SSE2_ON_MSVC_2008_OR_LATER
-      #endif
+    // Remember that usage of defined() in a #define is undefined by the standard.
+    // a user reported that in 64-bit mode, MSVC doesn't care to define _M_IX86_FP.
+    #if (defined(_M_IX86_FP) && (_M_IX86_FP >= 2)) || defined(_M_X64)
+      #define EIGEN_SSE2_ON_MSVC_2008_OR_LATER
     #endif   
   #endif
 #endif

Eigen/Dense

@@ -0,0 +1,8 @@
+#include "Core"
+#include "Array"
+#include "LU"
+#include "Cholesky"
+#include "QR"
+#include "SVD"
+#include "Geometry"
+#include "LeastSquares"

Eigen/Eigen

@@ -0,0 +1,2 @@
+#include "Dense"
+#include "Sparse"

Eigen/LeastSquares

@@ -5,7 +5,6 @@
 
 #include "src/Core/util/DisableMSVCWarnings.h"
 
-#include "LU"
 #include "QR"
 #include "Geometry"
 

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 size_t    size_type;
+  typedef 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/Sparse

@@ -22,7 +22,6 @@
 #endif
 
 #ifdef EIGEN_TAUCS_SUPPORT
-
   // taucs.h declares a lot of mess
   #define isnan
   #define finite
@@ -40,7 +39,9 @@
   #ifdef max
     #undef max
   #endif
-
+  #ifdef complex
+    #undef complex
+  #endif
 #endif
 
 #ifdef EIGEN_SUPERLU_SUPPORT
@@ -102,6 +103,7 @@ namespace Eigen {
 #include "src/Sparse/SparseFuzzy.h"
 #include "src/Sparse/SparseFlagged.h"
 #include "src/Sparse/SparseProduct.h"
+#include "src/Sparse/SparseDiagonalProduct.h"
 #include "src/Sparse/TriangularSolver.h"
 #include "src/Sparse/SparseLLT.h"
 #include "src/Sparse/SparseLDLT.h"

Eigen/StdVector

@@ -1,15 +1,147 @@
+#ifdef EIGEN_USE_NEW_STDVECTOR
+#include "NewStdVector"
+#else
+
 #ifndef EIGEN_STDVECTOR_MODULE_H
 #define EIGEN_STDVECTOR_MODULE_H
 
-#include "Core"
-#include <vector>
+#if defined(_GLIBCXX_VECTOR) || defined(_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
+#ifdef __GNUC__
+  #define EIGEN_GNUC_AT_LEAST(x,y) ((__GNUC__>=x && __GNUC_MINOR__>=y) || __GNUC__>x)
+#else
+  #define EIGEN_GNUC_AT_LEAST(x,y) 0
+#endif
+#endif
+
+#define vector std_vector
+#include <vector>        
+#undef vector
 
 namespace Eigen {
-#include "src/StdVector/UnalignedType.h"
-} // namespace Eigen
+
+template<typename T> class aligned_allocator;
+
+// meta programming to determine if a class has a given member
+struct ei_does_not_have_aligned_operator_new_marker_sizeof {int a[1];};
+struct ei_has_aligned_operator_new_marker_sizeof {int a[2];};
+
+template<typename ClassType>
+struct ei_has_aligned_operator_new {
+    template<typename T>
+    static ei_has_aligned_operator_new_marker_sizeof
+    test(T const *, typename T::ei_operator_new_marker_type const * = 0);
+    static ei_does_not_have_aligned_operator_new_marker_sizeof
+    test(...);
+
+    // note that the following indirection is needed for gcc-3.3
+    enum {ret =  sizeof(test(static_cast<ClassType*>(0))) 
+              == sizeof(ei_has_aligned_operator_new_marker_sizeof) };
+};
+
+#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 {
-#include "src/StdVector/StdVector.h"
-} // 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;  \
+    explicit vector(const allocator_type& __a = allocator_type()) : vector_base(__a) {}  \
+    vector(const vector& c) : vector_base(c) {}  \
+    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,
+         typename AllocT = std::allocator<T>,
+         bool HasAlignedNew = Eigen::ei_has_aligned_operator_new<T>::ret>
+class vector : public std::std_vector<T,AllocT>
+{
+  typedef std_vector<T, AllocT> vector_base;
+  EIGEN_STD_VECTOR_SPECIALIZATION_BODY
+};
+
+template<typename T,typename DummyAlloc>
+class vector<T,DummyAlloc,true>
+  : public std::std_vector<_EIGEN_WORKAROUND_MSVC_STD_VECTOR(T),
+                           Eigen::aligned_allocator<_EIGEN_WORKAROUND_MSVC_STD_VECTOR(T)> >          
+{
+  typedef std_vector<_EIGEN_WORKAROUND_MSVC_STD_VECTOR(T),
+                     Eigen::aligned_allocator<_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());
+  }
+  #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)
+  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  
+};
+
+}
 
 #endif // EIGEN_STDVECTOR_MODULE_H
+
+#endif // EIGEN_USE_NEW_STDVECTOR

Eigen/src/Array/PartialRedux.h

@@ -61,7 +61,11 @@ struct ei_traits<PartialReduxExpr<MatrixType, MemberOp, Direction> >
     Flags = (unsigned int)_MatrixTypeNested::Flags & HereditaryBits,
     TraversalSize = Direction==Vertical ? RowsAtCompileTime : ColsAtCompileTime
   };
+  #if EIGEN_GNUC_AT_LEAST(3,4)
   typedef typename MemberOp::template Cost<InputScalar,int(TraversalSize)> CostOpType;
+  #else
+  typedef typename MemberOp::template Cost<InputScalar,TraversalSize> CostOpType;
+  #endif
   enum {
     CoeffReadCost = TraversalSize * ei_traits<_MatrixTypeNested>::CoeffReadCost + int(CostOpType::value)
   };
@@ -104,7 +108,7 @@ class PartialReduxExpr : ei_no_assignment_operator,
     { enum { value = COST }; };                                     \
     template<typename Derived>                                      \
     inline ResultType operator()(const MatrixBase<Derived>& mat) const     \
-    { return mat.MEMBER(); }                                        \
+    { return mat.MEMBER(); } \
   }
 
 EIGEN_MEMBER_FUNCTOR(squaredNorm, Size * NumTraits<Scalar>::MulCost + (Size-1)*NumTraits<Scalar>::AddCost);

Eigen/src/Array/Random.h

@@ -110,7 +110,7 @@ MatrixBase<Derived>::Random()
   * Example: \include MatrixBase_setRandom.cpp
   * Output: \verbinclude MatrixBase_setRandom.out
   *
-  * \sa class CwiseNullaryOp, MatrixBase::setRandom(int,int)
+  * \sa class CwiseNullaryOp, setRandom(int), setRandom(int,int)
   */
 template<typename Derived>
 inline Derived& MatrixBase<Derived>::setRandom()
@@ -118,4 +118,39 @@ inline Derived& MatrixBase<Derived>::setRandom()
   return *this = Random(rows(), cols());
 }
 
+/** Resizes to the given \a size, and sets all coefficients in this expression to random values.
+  *
+  * \only_for_vectors
+  *
+  * Example: \include Matrix_setRandom_int.cpp
+  * Output: \verbinclude Matrix_setRandom_int.out
+  *
+  * \sa MatrixBase::setRandom(), setRandom(int,int), class CwiseNullaryOp, MatrixBase::Random()
+  */
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+EIGEN_STRONG_INLINE Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>&
+Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setRandom(int size)
+{
+  resize(size);
+  return setRandom();
+}
+
+/** Resizes to the given size, and sets all coefficients in this expression to random values.
+  *
+  * \param rows the new number of rows
+  * \param cols the new number of columns
+  *
+  * Example: \include Matrix_setRandom_int_int.cpp
+  * Output: \verbinclude Matrix_setRandom_int_int.out
+  *
+  * \sa MatrixBase::setRandom(), setRandom(int), class CwiseNullaryOp, MatrixBase::Random()
+  */
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+EIGEN_STRONG_INLINE Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>&
+Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setRandom(int rows, int cols)
+{
+  resize(rows, cols);
+  return setRandom();
+}
+
 #endif // EIGEN_RANDOM_H

Eigen/src/CMakeLists.txt

@@ -7,4 +7,3 @@ ADD_SUBDIRECTORY(Array)
 ADD_SUBDIRECTORY(Geometry)
 ADD_SUBDIRECTORY(LeastSquares)
 ADD_SUBDIRECTORY(Sparse)
-ADD_SUBDIRECTORY(StdVector)

Eigen/src/Cholesky/LDLT.h

@@ -68,8 +68,8 @@ template<typename MatrixType> class LDLT
     /** \returns true if the matrix is positive definite */
     inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
 
-    template<typename RhsDerived, typename ResDerived>
-    bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
+    template<typename RhsDerived, typename ResultType>
+    bool solve(const MatrixBase<RhsDerived> &b, ResultType *result) const;
 
     template<typename Derived>
     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
@@ -152,9 +152,9 @@ void LDLT<MatrixType>::compute(const MatrixType& a)
   * \sa LDLT::solveInPlace(), MatrixBase::ldlt()
   */
 template<typename MatrixType>
-template<typename RhsDerived, typename ResDerived>
+template<typename RhsDerived, typename ResultType>
 bool LDLT<MatrixType>
-::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const
+::solve(const MatrixBase<RhsDerived> &b, ResultType *result) const
 {
   const int size = m_matrix.rows();
   ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b");

Eigen/src/Cholesky/LLT.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>
 //
@@ -41,11 +41,16 @@
   * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other
   * situations like generalised eigen problems with hermitian matrices.
   *
-  * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
-  * the strict lower part does not have to store correct values.
+  * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices,
+  * use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine whether a system of equations
+  * has a solution.
   *
   * \sa MatrixBase::llt(), class LDLT
   */
+ /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
+  * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
+  * the strict lower part does not have to store correct values.
+  */
 template<typename MatrixType> class LLT
 {
   private:
@@ -60,20 +65,33 @@ template<typename MatrixType> class LLT
 
   public:
 
+    /** 
+    * \brief Default Constructor.
+    *
+    * The default constructor is useful in cases in which the user intends to
+    * perform decompositions via LLT::compute(const MatrixType&).
+    */
+    LLT() : m_matrix(), m_isInitialized(false) {}
+
     LLT(const MatrixType& matrix)
-      : m_matrix(matrix.rows(), matrix.cols())
+      : m_matrix(matrix.rows(), matrix.cols()),
+        m_isInitialized(false)
     {
       compute(matrix);
     }
 
     /** \returns the lower triangular matrix L */
-    inline Part<MatrixType, LowerTriangular> matrixL(void) const { return m_matrix; }
+    inline Part<MatrixType, LowerTriangular> matrixL(void) const 
+    { 
+      ei_assert(m_isInitialized && "LLT is not initialized.");
+      return m_matrix; 
+    }
+    
+    /** \deprecated */
+    inline bool isPositiveDefinite(void) const { return m_isInitialized && m_isPositiveDefinite; }
 
-    /** \returns true if the matrix is positive definite */
-    inline bool isPositiveDefinite(void) const { return m_isPositiveDefinite; }
-
-    template<typename RhsDerived, typename ResDerived>
-    bool solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const;
+    template<typename RhsDerived, typename ResultType>
+    bool solve(const MatrixBase<RhsDerived> &b, ResultType *result) const;
 
     template<typename Derived>
     bool solveInPlace(MatrixBase<Derived> &bAndX) const;
@@ -86,6 +104,7 @@ template<typename MatrixType> class LLT
       * The strict upper part is not used and even not initialized.
       */
     MatrixType m_matrix;
+    bool m_isInitialized;
     bool m_isPositiveDefinite;
 };
 
@@ -95,24 +114,34 @@ template<typename MatrixType>
 void LLT<MatrixType>::compute(const MatrixType& a)
 {
   assert(a.rows()==a.cols());
+  m_isPositiveDefinite = true;
   const int size = a.rows();
   m_matrix.resize(size, size);
-  const RealScalar eps = ei_sqrt(precision<Scalar>());
-
+  // The biggest overall is the point of reference to which further diagonals
+  // are compared; if any diagonal is negligible compared
+  // to the largest overall, the algorithm bails.  This cutoff is suggested
+  // in "Analysis of the Cholesky Decomposition of a Semi-definite Matrix" by
+  // Nicholas J. Higham. Also see "Accuracy and Stability of Numerical
+  // Algorithms" page 217, also by Higham.
+  const RealScalar cutoff = machine_epsilon<Scalar>() * size * a.diagonal().cwise().abs().maxCoeff();
   RealScalar x;
   x = ei_real(a.coeff(0,0));
-  m_isPositiveDefinite = x > eps && ei_isMuchSmallerThan(ei_imag(a.coeff(0,0)), RealScalar(1));
   m_matrix.coeffRef(0,0) = ei_sqrt(x);
+  if(size==1)
+  {
+    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));
   for (int j = 1; j < size; ++j)
   {
-    Scalar tmp = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm();
-    x = ei_real(tmp);
-    if (x < eps || (!ei_isMuchSmallerThan(ei_imag(tmp), RealScalar(1))))
+    x = ei_real(a.coeff(j,j)) - m_matrix.row(j).start(j).squaredNorm();
+    if (x < cutoff)
     {
       m_isPositiveDefinite = false;
-      return;
+      continue;
     }
+
     m_matrix.coeffRef(j,j) = x = ei_sqrt(x);
 
     int endSize = size-j-1;
@@ -127,12 +156,14 @@ void LLT<MatrixType>::compute(const MatrixType& a)
         - m_matrix.col(j).end(endSize) ) / x;
     }
   }
+
+  m_isInitialized = true;
 }
 
 /** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
   * The result is stored in \a result
   *
-  * \returns true in case of success, false otherwise.
+  * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
   *
   * In other words, it computes \f$ b = A^{-1} b \f$ with
   * \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left.
@@ -143,9 +174,10 @@ void LLT<MatrixType>::compute(const MatrixType& a)
   * \sa LLT::solveInPlace(), MatrixBase::llt()
   */
 template<typename MatrixType>
-template<typename RhsDerived, typename ResDerived>
-bool LLT<MatrixType>::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDerived> *result) const
+template<typename RhsDerived, typename ResultType>
+bool LLT<MatrixType>::solve(const MatrixBase<RhsDerived> &b, ResultType *result) const
 {
+  ei_assert(m_isInitialized && "LLT is not initialized.");
   const int size = m_matrix.rows();
   ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b");
   return solveInPlace((*result) = b);
@@ -155,6 +187,8 @@ bool LLT<MatrixType>::solve(const MatrixBase<RhsDerived> &b, MatrixBase<ResDeriv
   *
   * \param bAndX represents both the right-hand side matrix b and result x.
   *
+  * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
+  *
   * This version avoids a copy when the right hand side matrix b is not
   * needed anymore.
   *
@@ -164,10 +198,9 @@ template<typename MatrixType>
 template<typename Derived>
 bool LLT<MatrixType>::solveInPlace(MatrixBase<Derived> &bAndX) const
 {
+  ei_assert(m_isInitialized && "LLT is not initialized.");
   const int size = m_matrix.rows();
   ei_assert(size==bAndX.rows());
-  if (!m_isPositiveDefinite)
-    return false;
   matrixL().solveTriangularInPlace(bAndX);
   m_matrix.adjoint().template part<UpperTriangular>().solveTriangularInPlace(bAndX);
   return true;

Eigen/src/Core/Assign.h

@@ -353,7 +353,7 @@ struct ei_assign_impl<Derived1, Derived2, SliceVectorization, NoUnrolling>
     const int outerSize = dst.outerSize();
     const int alignedStep = (packetSize - dst.stride() % packetSize) & packetAlignedMask;
     int alignedStart = ei_assign_traits<Derived1,Derived2>::DstIsAligned ? 0
-                     : ei_alignmentOffset(&dst.coeffRef(0), innerSize);
+                     : ei_alignmentOffset(&dst.coeffRef(0,0), innerSize);
 
     for(int i = 0; i < outerSize; ++i)
     {

Eigen/src/Core/Block.h

@@ -61,27 +61,28 @@
   *
   * \sa MatrixBase::block(int,int,int,int), MatrixBase::block(int,int), class VectorBlock
   */
+
 template<typename MatrixType, int BlockRows, int BlockCols, int _PacketAccess, int _DirectAccessStatus>
 struct ei_traits<Block<MatrixType, BlockRows, BlockCols, _PacketAccess, _DirectAccessStatus> >
 {
-  typedef typename MatrixType::Scalar Scalar;
-  typedef typename MatrixType::Nested MatrixTypeNested;
+  typedef typename ei_traits<MatrixType>::Scalar Scalar;
+  typedef typename ei_nested<MatrixType>::type MatrixTypeNested;
   typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
   enum{
-    RowsAtCompileTime = MatrixType::RowsAtCompileTime == 1 ? 1 : BlockRows,
-    ColsAtCompileTime = MatrixType::ColsAtCompileTime == 1 ? 1 : BlockCols,
+    RowsAtCompileTime = ei_traits<MatrixType>::RowsAtCompileTime == 1 ? 1 : BlockRows,
+    ColsAtCompileTime = ei_traits<MatrixType>::ColsAtCompileTime == 1 ? 1 : BlockCols,
     MaxRowsAtCompileTime = RowsAtCompileTime == 1 ? 1
-      : (BlockRows==Dynamic ? MatrixType::MaxRowsAtCompileTime : BlockRows),
+      : (BlockRows==Dynamic ? int(ei_traits<MatrixType>::MaxRowsAtCompileTime) : BlockRows),
     MaxColsAtCompileTime = ColsAtCompileTime == 1 ? 1
-      : (BlockCols==Dynamic ? MatrixType::MaxColsAtCompileTime : BlockCols),
-    RowMajor = int(MatrixType::Flags)&RowMajorBit,
-    InnerSize = RowMajor ? ColsAtCompileTime : RowsAtCompileTime,
-    InnerMaxSize = RowMajor ? MaxColsAtCompileTime : MaxRowsAtCompileTime,
+      : (BlockCols==Dynamic ? int(ei_traits<MatrixType>::MaxColsAtCompileTime) : BlockCols),
+    RowMajor = int(ei_traits<MatrixType>::Flags)&RowMajorBit,
+    InnerSize = RowMajor ? int(ColsAtCompileTime) : int(RowsAtCompileTime),
+    InnerMaxSize = RowMajor ? int(MaxColsAtCompileTime) : int(MaxRowsAtCompileTime),
     MaskPacketAccessBit = (InnerMaxSize == Dynamic || (InnerSize >= ei_packet_traits<Scalar>::size))
                         ? PacketAccessBit : 0,
     FlagsLinearAccessBit = (RowsAtCompileTime == 1 || ColsAtCompileTime == 1) ? LinearAccessBit : 0,
-    Flags = (MatrixType::Flags & (HereditaryBits | MaskPacketAccessBit | DirectAccessBit)) | FlagsLinearAccessBit,
-    CoeffReadCost = MatrixType::CoeffReadCost,
+    Flags = (ei_traits<MatrixType>::Flags & (HereditaryBits | MaskPacketAccessBit | DirectAccessBit)) | FlagsLinearAccessBit,
+    CoeffReadCost = ei_traits<MatrixType>::CoeffReadCost,
     PacketAccess = _PacketAccess
   };
   typedef typename ei_meta_if<int(PacketAccess)==ForceAligned,
@@ -122,7 +123,7 @@ template<typename MatrixType, int BlockRows, int BlockCols, int PacketAccess, in
       : m_matrix(matrix), m_startRow(startRow), m_startCol(startCol),
         m_blockRows(matrix.rows()), m_blockCols(matrix.cols())
     {
-      EIGEN_STATIC_ASSERT(RowsAtCompileTime!=Dynamic && RowsAtCompileTime!=Dynamic,THIS_METHOD_IS_ONLY_FOR_FIXED_SIZE)
+      EIGEN_STATIC_ASSERT(RowsAtCompileTime!=Dynamic && ColsAtCompileTime!=Dynamic,THIS_METHOD_IS_ONLY_FOR_FIXED_SIZE)
       ei_assert(startRow >= 0 && BlockRows >= 1 && startRow + BlockRows <= matrix.rows()
           && startCol >= 0 && BlockCols >= 1 && startCol + BlockCols <= matrix.cols());
     }
@@ -221,15 +222,13 @@ 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>;
 
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Block)
 
-    AlignedDerivedType forceAligned()
+    AlignedDerivedType _convertToForceAligned()
     {
-      if (PacketAccess==ForceAligned)
-        return *this;
-      else
-        return Block<MatrixType,BlockRows,BlockCols,ForceAligned,HasDirectAccess>
+      return Block<MatrixType,BlockRows,BlockCols,ForceAligned,HasDirectAccess>
                     (m_matrix, Base::m_data, Base::m_rows.value(), Base::m_cols.value());
     }
 
@@ -454,7 +453,7 @@ MatrixBase<Derived>::end(int size) const
   * \only_for_vectors
   *
   * The template parameter \a Size is the number of coefficients in the block
-  * 
+  *
   * \param start the index of the first element of the sub-vector
   *
   * Example: \include MatrixBase_template_int_segment.cpp

Eigen/src/Core/CacheFriendlyProduct.h

@@ -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;
@@ -361,13 +361,14 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
   #ifdef _EIGEN_ACCUMULATE_PACKETS
   #error _EIGEN_ACCUMULATE_PACKETS has already been defined
   #endif
-
-  #define _EIGEN_ACCUMULATE_PACKETS(A0,A13,A2,OFFSET) \
-    ei_pstore(&res[j OFFSET], \
-      ei_padd(ei_pload(&res[j OFFSET]), \
+  #define _EIGEN_ACCUMULATE_PACKETS(A0,A13,A2) \
+    ei_pstore(&res[j], \
+      ei_padd(ei_pload(&res[j]), \
         ei_padd( \
-          ei_padd(ei_pmul(ptmp0,ei_pload ## A0(&lhs0[j OFFSET])),ei_pmul(ptmp1,ei_pload ## A13(&lhs1[j OFFSET]))), \
-          ei_padd(ei_pmul(ptmp2,ei_pload ## A2(&lhs2[j OFFSET])),ei_pmul(ptmp3,ei_pload ## A13(&lhs3[j OFFSET]))) )))
+          ei_padd(ei_pmul(ptmp0,EIGEN_CAT(ei_ploa , A0)(&lhs0[j])), \
+                  ei_pmul(ptmp1,EIGEN_CAT(ei_ploa , A13)(&lhs1[j]))), \
+          ei_padd(ei_pmul(ptmp2,EIGEN_CAT(ei_ploa , A2)(&lhs2[j])), \
+                  ei_pmul(ptmp3,EIGEN_CAT(ei_ploa , A13)(&lhs3[j]))) )))
 
   typedef typename ei_packet_traits<Scalar>::type Packet;
   const int PacketSize = sizeof(Packet)/sizeof(Scalar);
@@ -397,7 +398,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
   if (PacketSize>1)
   {
     ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0 || size<PacketSize);
-    
+
     while (skipColumns<PacketSize &&
            alignedStart != ((lhsAlignmentOffset + alignmentStep*skipColumns)%PacketSize))
       ++skipColumns;
@@ -418,7 +419,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
 
   int offset1 = (FirstAligned && alignmentStep==1?3:1);
   int offset3 = (FirstAligned && alignmentStep==1?1:3);
-  
+
   int columnBound = ((rhs.size()-skipColumns)/columnsAtOnce)*columnsAtOnce + skipColumns;
   for (int i=skipColumns; i<columnBound; i+=columnsAtOnce)
   {
@@ -442,11 +443,11 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
         {
           case AllAligned:
             for (int j = alignedStart; j<alignedSize; j+=PacketSize)
-              _EIGEN_ACCUMULATE_PACKETS(,,,);
+              _EIGEN_ACCUMULATE_PACKETS(d,d,d);
             break;
           case EvenAligned:
             for (int j = alignedStart; j<alignedSize; j+=PacketSize)
-              _EIGEN_ACCUMULATE_PACKETS(,u,,);
+              _EIGEN_ACCUMULATE_PACKETS(d,du,d);
             break;
           case FirstAligned:
             if(peels>1)
@@ -482,11 +483,11 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
               }
             }
             for (int j = peeledSize; j<alignedSize; j+=PacketSize)
-              _EIGEN_ACCUMULATE_PACKETS(,u,u,);
+              _EIGEN_ACCUMULATE_PACKETS(d,du,du);
             break;
           default:
             for (int j = alignedStart; j<alignedSize; j+=PacketSize)
-              _EIGEN_ACCUMULATE_PACKETS(u,u,u,);
+              _EIGEN_ACCUMULATE_PACKETS(du,du,du);
             break;
         }
       }
@@ -494,7 +495,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_colmajor_times_vector(
 
     /* process remaining coeffs (or all if there is no explicit vectorization) */
     for (int j=alignedSize; j<size; ++j)
-	  res[j] += ei_pfirst(ptmp0)*lhs0[j] + ei_pfirst(ptmp1)*lhs1[j] + ei_pfirst(ptmp2)*lhs2[j] + ei_pfirst(ptmp3)*lhs3[j];
+      res[j] += ei_pfirst(ptmp0)*lhs0[j] + ei_pfirst(ptmp1)*lhs1[j] + ei_pfirst(ptmp2)*lhs2[j] + ei_pfirst(ptmp3)*lhs3[j];
   }
 
   // process remaining first and last columns (at most columnsAtOnce-1)
@@ -550,12 +551,12 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
   #error _EIGEN_ACCUMULATE_PACKETS has already been defined
   #endif
 
-  #define _EIGEN_ACCUMULATE_PACKETS(A0,A13,A2,OFFSET) {\
+  #define _EIGEN_ACCUMULATE_PACKETS(A0,A13,A2) {\
     Packet b = ei_pload(&rhs[j]); \
-    ptmp0 = ei_pmadd(b, ei_pload##A0 (&lhs0[j]), ptmp0); \
-    ptmp1 = ei_pmadd(b, ei_pload##A13(&lhs1[j]), ptmp1); \
-    ptmp2 = ei_pmadd(b, ei_pload##A2 (&lhs2[j]), ptmp2); \
-    ptmp3 = ei_pmadd(b, ei_pload##A13(&lhs3[j]), ptmp3); }
+    ptmp0 = ei_pmadd(b, EIGEN_CAT(ei_ploa,A0) (&lhs0[j]), ptmp0); \
+    ptmp1 = ei_pmadd(b, EIGEN_CAT(ei_ploa,A13)(&lhs1[j]), ptmp1); \
+    ptmp2 = ei_pmadd(b, EIGEN_CAT(ei_ploa,A2) (&lhs2[j]), ptmp2); \
+    ptmp3 = ei_pmadd(b, EIGEN_CAT(ei_ploa,A13)(&lhs3[j]), ptmp3); }
 
   typedef typename ei_packet_traits<Scalar>::type Packet;
   const int PacketSize = sizeof(Packet)/sizeof(Scalar);
@@ -580,13 +581,13 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
 
   // we cannot assume the first element is aligned because of sub-matrices
   const int lhsAlignmentOffset = ei_alignmentOffset(lhs,size);
-  
+
   // find how many rows do we have to skip to be aligned with rhs (if possible)
   int skipRows = 0;
   if (PacketSize>1)
   {
     ei_internal_assert(size_t(lhs+lhsAlignmentOffset)%sizeof(Packet)==0  || size<PacketSize);
-    
+
     while (skipRows<PacketSize &&
            alignedStart != ((lhsAlignmentOffset + alignmentStep*skipRows)%PacketSize))
       ++skipRows;
@@ -607,7 +608,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
 
   int offset1 = (FirstAligned && alignmentStep==1?3:1);
   int offset3 = (FirstAligned && alignmentStep==1?1:3);
-  
+
   int rowBound = ((res.size()-skipRows)/rowsAtOnce)*rowsAtOnce + skipRows;
   for (int i=skipRows; i<rowBound; i+=rowsAtOnce)
   {
@@ -621,7 +622,7 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
     {
       /* explicit vectorization */
       Packet ptmp0 = ei_pset1(Scalar(0)), ptmp1 = ei_pset1(Scalar(0)), ptmp2 = ei_pset1(Scalar(0)), ptmp3 = ei_pset1(Scalar(0));
-      
+
       // process initial unaligned coeffs
       // FIXME this loop get vectorized by the compiler !
       for (int j=0; j<alignedStart; ++j)
@@ -636,11 +637,11 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
         {
           case AllAligned:
             for (int j = alignedStart; j<alignedSize; j+=PacketSize)
-              _EIGEN_ACCUMULATE_PACKETS(,,,);
+              _EIGEN_ACCUMULATE_PACKETS(d,d,d);
             break;
           case EvenAligned:
             for (int j = alignedStart; j<alignedSize; j+=PacketSize)
-              _EIGEN_ACCUMULATE_PACKETS(,u,,);
+              _EIGEN_ACCUMULATE_PACKETS(d,du,d);
             break;
           case FirstAligned:
             if (peels>1)
@@ -679,11 +680,11 @@ static EIGEN_DONT_INLINE void ei_cache_friendly_product_rowmajor_times_vector(
               }
             }
             for (int j = peeledSize; j<alignedSize; j+=PacketSize)
-              _EIGEN_ACCUMULATE_PACKETS(,u,u,);
+              _EIGEN_ACCUMULATE_PACKETS(d,du,du);
             break;
           default:
             for (int j = alignedStart; j<alignedSize; j+=PacketSize)
-              _EIGEN_ACCUMULATE_PACKETS(u,u,u,);
+              _EIGEN_ACCUMULATE_PACKETS(du,du,du);
             break;
         }
         tmp0 += ei_predux(ptmp0);

Eigen/src/Core/CwiseNullaryOp.h

@@ -146,6 +146,7 @@ template<typename CustomNullaryOp>
 EIGEN_STRONG_INLINE const CwiseNullaryOp<CustomNullaryOp, Derived>
 MatrixBase<Derived>::NullaryExpr(int size, const CustomNullaryOp& func)
 {
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
   ei_assert(IsVectorAtCompileTime);
   if(RowsAtCompileTime == 1) return CwiseNullaryOp<CustomNullaryOp, Derived>(1, size, func);
   else return CwiseNullaryOp<CustomNullaryOp, Derived>(size, 1, func);
@@ -227,6 +228,7 @@ MatrixBase<Derived>::Constant(const Scalar& value)
   return NullaryExpr(RowsAtCompileTime, ColsAtCompileTime, ei_scalar_constant_op<Scalar>(value));
 }
 
+/** \returns true if all coefficients in this matrix are approximately equal to \a value, to within precision \a prec */
 template<typename Derived>
 bool MatrixBase<Derived>::isApproxToConstant
 (const Scalar& value, RealScalar prec) const
@@ -238,6 +240,16 @@ bool MatrixBase<Derived>::isApproxToConstant
   return true;
 }
 
+/** This is just an alias for isApproxToConstant().
+  *
+  * \returns true if all coefficients in this matrix are approximately equal to \a value, to within precision \a prec */
+template<typename Derived>
+bool MatrixBase<Derived>::isConstant
+(const Scalar& value, RealScalar prec) const
+{
+  return isApproxToConstant(value, prec);
+}
+
 /** Alias for setConstant(): sets all coefficients in this expression to \a value.
   *
   * \sa setConstant(), Constant(), class CwiseNullaryOp
@@ -250,7 +262,7 @@ EIGEN_STRONG_INLINE void MatrixBase<Derived>::fill(const Scalar& value)
 
 /** Sets all coefficients in this expression to \a value.
   *
-  * \sa fill(), Constant(), class CwiseNullaryOp, setZero(), setOnes()
+  * \sa fill(), setConstant(int,const Scalar&), setConstant(int,int,const Scalar&), setZero(), setOnes(), Constant(), class CwiseNullaryOp, setZero(), setOnes()
   */
 template<typename Derived>
 EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::setConstant(const Scalar& value)
@@ -258,6 +270,42 @@ EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::setConstant(const Scalar& valu
   return derived() = Constant(rows(), cols(), value);
 }
 
+/** Resizes to the given \a size, and sets all coefficients in this expression to the given \a value.
+  *
+  * \only_for_vectors
+  *
+  * Example: \include Matrix_set_int.cpp
+  * Output: \verbinclude Matrix_setConstant_int.out
+  *
+  * \sa MatrixBase::setConstant(const Scalar&), setConstant(int,int,const Scalar&), class CwiseNullaryOp, MatrixBase::Constant(const Scalar&)
+  */
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+EIGEN_STRONG_INLINE Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>&
+Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setConstant(int size, const Scalar& value)
+{
+  resize(size);
+  return setConstant(value);
+}
+
+/** Resizes to the given size, and sets all coefficients in this expression to the given \a value.
+  *
+  * \param rows the new number of rows
+  * \param cols the new number of columns
+  *
+  * Example: \include Matrix_setConstant_int_int.cpp
+  * Output: \verbinclude Matrix_setConstant_int_int.out
+  *
+  * \sa MatrixBase::setConstant(const Scalar&), setConstant(int,const Scalar&), class CwiseNullaryOp, MatrixBase::Constant(const Scalar&)
+  */
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+EIGEN_STRONG_INLINE Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>&
+Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setConstant(int rows, int cols, const Scalar& value)
+{
+  resize(rows, cols);
+  return setConstant(value);
+}
+
+
 // zero:
 
 /** \returns an expression of a zero matrix.
@@ -354,6 +402,41 @@ EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::setZero()
   return setConstant(Scalar(0));
 }
 
+/** Resizes to the given \a size, and sets all coefficients in this expression to zero.
+  *
+  * \only_for_vectors
+  *
+  * Example: \include Matrix_setZero_int.cpp
+  * Output: \verbinclude Matrix_setZero_int.out
+  *
+  * \sa MatrixBase::setZero(), setZero(int,int), class CwiseNullaryOp, MatrixBase::Zero()
+  */
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+EIGEN_STRONG_INLINE Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>&
+Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setZero(int size)
+{
+  resize(size);
+  return setConstant(Scalar(0));
+}
+
+/** Resizes to the given size, and sets all coefficients in this expression to zero.
+  *
+  * \param rows the new number of rows
+  * \param cols the new number of columns
+  *
+  * Example: \include Matrix_setZero_int_int.cpp
+  * Output: \verbinclude Matrix_setZero_int_int.out
+  *
+  * \sa MatrixBase::setZero(), setZero(int), class CwiseNullaryOp, MatrixBase::Zero()
+  */
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+EIGEN_STRONG_INLINE Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>&
+Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setZero(int rows, int cols)
+{
+  resize(rows, cols);
+  return setConstant(Scalar(0));
+}
+
 // ones:
 
 /** \returns an expression of a matrix where all coefficients equal one.
@@ -447,6 +530,41 @@ EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::setOnes()
   return setConstant(Scalar(1));
 }
 
+/** Resizes to the given \a size, and sets all coefficients in this expression to one.
+  *
+  * \only_for_vectors
+  *
+  * Example: \include Matrix_setOnes_int.cpp
+  * Output: \verbinclude Matrix_setOnes_int.out
+  *
+  * \sa MatrixBase::setOnes(), setOnes(int,int), class CwiseNullaryOp, MatrixBase::Ones()
+  */
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+EIGEN_STRONG_INLINE Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>&
+Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setOnes(int size)
+{
+  resize(size);
+  return setConstant(Scalar(1));
+}
+
+/** Resizes to the given size, and sets all coefficients in this expression to one.
+  *
+  * \param rows the new number of rows
+  * \param cols the new number of columns
+  *
+  * Example: \include Matrix_setOnes_int_int.cpp
+  * Output: \verbinclude Matrix_setOnes_int_int.out
+  *
+  * \sa MatrixBase::setOnes(), setOnes(int), class CwiseNullaryOp, MatrixBase::Ones()
+  */
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+EIGEN_STRONG_INLINE Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>&
+Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setOnes(int rows, int cols)
+{
+  resize(rows, cols);
+  return setConstant(Scalar(1));
+}
+
 // Identity:
 
 /** \returns an expression of the identity matrix (not necessarily square).
@@ -556,6 +674,24 @@ EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::setIdentity()
   return ei_setIdentity_impl<Derived>::run(derived());
 }
 
+/** Resizes to the given size, and writes the identity expression (not necessarily square) into *this.
+  *
+  * \param rows the new number of rows
+  * \param cols the new number of columns
+  *
+  * Example: \include Matrix_setIdentity_int_int.cpp
+  * Output: \verbinclude Matrix_setIdentity_int_int.out
+  *
+  * \sa MatrixBase::setIdentity(), class CwiseNullaryOp, MatrixBase::Identity()
+  */
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+EIGEN_STRONG_INLINE Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>&
+Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setIdentity(int rows, int cols)
+{
+  resize(rows, cols);
+  return setIdentity();
+}
+
 /** \returns an expression of the i-th unit (basis) vector.
   *
   * \only_for_vectors

Eigen/src/Core/DiagonalMatrix.h

@@ -62,6 +62,7 @@ class DiagonalMatrix : ei_no_assignment_operator,
   public:
 
     EIGEN_GENERIC_PUBLIC_INTERFACE(DiagonalMatrix)
+    typedef CoeffsVectorType _CoeffsVectorType;
 
     // needed to evaluate a DiagonalMatrix<Xpr> to a DiagonalMatrix<NestByValue<Vector> >
     template<typename OtherCoeffsVectorType>

Eigen/src/Core/Map.h

@@ -66,12 +66,9 @@ template<typename MatrixType, int PacketAccess> class Map
 
     inline int stride() const { return this->innerSize(); }
 
-    AlignedDerivedType forceAligned()
+    AlignedDerivedType _convertToForceAligned()
     {
-      if (PacketAccess==ForceAligned)
-        return *this;
-      else
-        return Map<MatrixType,ForceAligned>(Base::m_data, Base::m_rows.value(), Base::m_cols.value());
+      return Map<MatrixType,ForceAligned>(Base::m_data, Base::m_rows.value(), Base::m_cols.value());
     }
 
     inline Map(const Scalar* data) : Base(data) {}

Eigen/src/Core/MapBase.h

@@ -65,9 +65,20 @@ template<typename Derived> class MapBase
     inline int stride() const { return derived().stride(); }
     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(); }
+    };
+
+    template<typename Dummy> struct force_aligned_impl<false,Dummy> {
+      AlignedDerivedType static run(MapBase& a) { return a.derived()._convertToForceAligned(); }
+    };
+
     /** \returns an expression equivalent to \c *this but having the \c PacketAccess constant
       * set to \c ForceAligned. Must be reimplemented by the derived class. */
-    AlignedDerivedType forceAligned() { return derived().forceAligned(); }
+    AlignedDerivedType forceAligned()
+    {
+      return force_aligned_impl<int(PacketAccess)==int(ForceAligned),Derived>::run(*this);
+    }
 
     inline const Scalar& coeff(int row, int col) const
     {
@@ -96,7 +107,11 @@ template<typename Derived> class MapBase
 
     inline Scalar& coeffRef(int index)
     {
-      return *const_cast<Scalar*>(m_data + index);
+      ei_assert(Derived::IsVectorAtCompileTime || (ei_traits<Derived>::Flags & LinearAccessBit));
+      if ( ((RowsAtCompileTime == 1) == IsRowMajor) )
+        return const_cast<Scalar*>(m_data)[index];
+      else
+        return const_cast<Scalar*>(m_data)[index*stride()];
     }
 
     template<int LoadMode>
@@ -150,7 +165,20 @@ template<typename Derived> class MapBase
               || (   rows > 0 && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows)
                   && cols > 0 && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols)));
     }
+
+    Derived& operator=(const MapBase& other)
+    {
+      return Base::operator=(other);
+    }
+
+    template<typename OtherDerived>
+    Derived& operator=(const MatrixBase<OtherDerived>& other)
+    {
+      return Base::operator=(other);
+    }
     
+    using Base::operator*=;
+
     template<typename OtherDerived>
     Derived& operator+=(const MatrixBase<OtherDerived>& other)
     { return derived() = forceAligned() + other; }

Eigen/src/Core/MathFunctions.h

@@ -26,6 +26,7 @@
 #define EIGEN_MATHFUNCTIONS_H
 
 template<typename T> inline typename NumTraits<T>::Real precision();
+template<typename T> inline typename NumTraits<T>::Real machine_epsilon();
 template<typename T> inline T ei_random(T a, T b);
 template<typename T> inline T ei_random();
 template<typename T> inline T ei_random_amplitude()
@@ -49,6 +50,7 @@ template<typename T> inline T ei_hypot(T x, T y)
 **************/
 
 template<> inline int precision<int>() { return 0; }
+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; }
@@ -59,12 +61,8 @@ 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)
 {
@@ -93,6 +91,7 @@ inline bool ei_isApproxOrLessThan(int a, int b, int = precision<int>())
 **************/
 
 template<> inline float precision<float>() { return 1e-5f; }
+template<> inline float machine_epsilon<float>() { return 1.192e-07f; }
 inline float ei_real(float x)  { return x; }
 inline float ei_imag(float)    { return 0.f; }
 inline float ei_conj(float x)  { return x; }
@@ -103,6 +102,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)
@@ -138,6 +138,8 @@ inline bool ei_isApproxOrLessThan(float a, float b, float prec = precision<float
 **************/
 
 template<> inline double precision<double>() { return 1e-11; }
+template<> inline double machine_epsilon<double>() { return 2.220e-16; }
+
 inline double ei_real(double x)  { return x; }
 inline double ei_imag(double)    { return 0.; }
 inline double ei_conj(double x)  { return x; }
@@ -148,6 +150,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)
@@ -183,6 +186,7 @@ inline bool ei_isApproxOrLessThan(double a, double b, double prec = precision<do
 *********************/
 
 template<> inline float precision<std::complex<float> >() { return precision<float>(); }
+template<> inline float machine_epsilon<std::complex<float> >() { return machine_epsilon<float>(); }
 inline float ei_real(const std::complex<float>& x) { return std::real(x); }
 inline float ei_imag(const std::complex<float>& x) { return std::imag(x); }
 inline std::complex<float> ei_conj(const std::complex<float>& x) { return std::conj(x); }
@@ -191,6 +195,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()
 {
@@ -216,6 +221,7 @@ inline bool ei_isApprox(const std::complex<float>& a, const std::complex<float>&
 **********************/
 
 template<> inline double precision<std::complex<double> >() { return precision<double>(); }
+template<> inline double machine_epsilon<std::complex<double> >() { return machine_epsilon<double>(); }
 inline double ei_real(const std::complex<double>& x) { return std::real(x); }
 inline double ei_imag(const std::complex<double>& x) { return std::imag(x); }
 inline std::complex<double> ei_conj(const std::complex<double>& x) { return std::conj(x); }
@@ -224,6 +230,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()
 {
@@ -250,6 +257,7 @@ inline bool ei_isApprox(const std::complex<double>& a, const std::complex<double
 ******************/
 
 template<> inline long double precision<long double>() { return precision<double>(); }
+template<> inline long double machine_epsilon<long double>() { return 1.084e-19l; }
 inline long double ei_real(long double x)  { return x; }
 inline long double ei_imag(long double)    { return 0.; }
 inline long double ei_conj(long double x)  { return x; }
@@ -260,6 +268,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)

Eigen/src/Core/Matrix.h

@@ -137,6 +137,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(); }
@@ -226,7 +229,6 @@ class Matrix
       */
     inline void resize(int rows, int cols)
     {
-      ei_assert(rows > 0 && cols > 0 && "a matrix cannot be resized to 0 size");
       ei_assert((MaxRowsAtCompileTime == Dynamic || MaxRowsAtCompileTime >= rows)
              && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows)
              && (MaxColsAtCompileTime == Dynamic || MaxColsAtCompileTime >= cols)
@@ -240,7 +242,6 @@ class Matrix
       */
     inline void resize(int size)
     {
-      ei_assert(size>0 && "a vector cannot be resized to 0 length");
       EIGEN_STATIC_ASSERT_VECTOR_ONLY(Matrix)
       if(RowsAtCompileTime == 1)
         m_storage.resize(size, 1, size);
@@ -397,13 +398,8 @@ 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,
@@ -440,6 +436,25 @@ class Matrix
     { return AlignedMapType(data, rows, cols); }
     //@}
 
+    using Base::setConstant;
+    Matrix& setConstant(int size, const Scalar& value);
+    Matrix& setConstant(int rows, int cols, const Scalar& value);
+
+    using Base::setZero;
+    Matrix& setZero(int size);
+    Matrix& setZero(int rows, int cols);
+
+    using Base::setOnes;
+    Matrix& setOnes(int size);
+    Matrix& setOnes(int rows, int cols);
+
+    using Base::setRandom;
+    Matrix& setRandom(int size);
+    Matrix& setRandom(int rows, int cols);
+
+    using Base::setIdentity;
+    Matrix& setIdentity(int rows, int cols);
+
 /////////// Geometry module ///////////
 
     template<typename OtherDerived>
@@ -490,10 +505,16 @@ class Matrix
     template<typename OtherDerived>
     EIGEN_STRONG_INLINE Matrix& _set(const MatrixBase<OtherDerived>& other)
     {
-      _resize_to_match(other);
-      return Base::operator=(other);
+      _set_selector(other.derived(), typename ei_meta_if<bool(int(OtherDerived::Flags) & EvalBeforeAssigningBit), 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.
       *
@@ -508,17 +529,48 @@ class Matrix
       return ei_assign_selector<Matrix,OtherDerived,false>::run(*this, other.derived());
     }
 
-  static EIGEN_STRONG_INLINE void _check_template_params()
+    static EIGEN_STRONG_INLINE void _check_template_params()
+    {
+        EIGEN_STATIC_ASSERT((_Rows > 0
+                        && _Cols > 0
+                        && _MaxRows <= _Rows
+                        && _MaxCols <= _Cols
+                        && (_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)
   {
-      EIGEN_STATIC_ASSERT((_Rows > 0
-                       && _Cols > 0
-                       && _MaxRows <= _Rows
-                       && _MaxCols <= _Cols
-                       && (_Options & (AutoAlign|RowMajor)) == _Options),
-        INVALID_MATRIX_TEMPLATE_PARAMETERS)
+    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)
+{
+  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

@@ -463,6 +463,7 @@ template<typename Derived> class MatrixBase
                            RealScalar prec = precision<Scalar>()) const;
 
     bool isApproxToConstant(const Scalar& value, RealScalar prec = precision<Scalar>()) const;
+    bool isConstant(const Scalar& value, RealScalar prec = precision<Scalar>()) const;
     bool isZero(RealScalar prec = precision<Scalar>()) const;
     bool isOnes(RealScalar prec = precision<Scalar>()) const;
     bool isIdentity(RealScalar prec = precision<Scalar>()) const;
@@ -531,8 +532,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

Eigen/src/Core/MatrixStorage.h

@@ -2,7 +2,7 @@
 // for linear algebra. Eigen itself is part of the KDE project.
 //
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
-// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2006-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
@@ -41,7 +41,7 @@ template <typename T, int Size, int MatrixOptions,
   {
     #ifndef EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT
     ei_assert((reinterpret_cast<size_t>(array) & 0xf) == 0
-              && "this assertion is explained here: http://eigen.tuxfamily.org/api/UnalignedArrayAssert.html  **** READ THIS WEB PAGE !!! ****");
+              && "this assertion is explained here: http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html  **** READ THIS WEB PAGE !!! ****");
     #endif
   }
 
@@ -176,7 +176,10 @@ template<typename T, int _Options> class ei_matrix_storage<T, Dynamic, Dynamic,
       if(size != m_rows*m_cols)
       {
         ei_aligned_delete(m_data, m_rows*m_cols);
-        m_data = ei_aligned_new<T>(size);
+        if (size)
+          m_data = ei_aligned_new<T>(size);
+        else
+          m_data = 0;
       }
       m_rows = rows;
       m_cols = cols;
@@ -203,7 +206,10 @@ template<typename T, int _Rows, int _Options> class ei_matrix_storage<T, Dynamic
       if(size != _Rows*m_cols)
       {
         ei_aligned_delete(m_data, _Rows*m_cols);
-        m_data = ei_aligned_new<T>(size);
+        if (size)
+          m_data = ei_aligned_new<T>(size);
+        else
+          m_data = 0;
       }
       m_cols = cols;
     }
@@ -229,7 +235,10 @@ template<typename T, int _Cols, int _Options> class ei_matrix_storage<T, Dynamic
       if(size != m_rows*_Cols)
       {
         ei_aligned_delete(m_data, _Cols*m_rows);
-        m_data = ei_aligned_new<T>(size);
+        if (size)
+          m_data = ei_aligned_new<T>(size);
+        else
+          m_data = 0;
       }
       m_rows = rows;
     }

Eigen/src/Core/Part.h

@@ -26,7 +26,7 @@
 #ifndef EIGEN_PART_H
 #define EIGEN_PART_H
 
-/** \nonstableyet 
+/** \nonstableyet
   * \class Part
   *
   * \brief Expression of a triangular matrix extracted from a given matrix
@@ -117,10 +117,10 @@ template<typename MatrixType, unsigned int Mode> class Part
     const Block<Part, RowsAtCompileTime, 1> col(int i) { return Base::col(i); }
     const Block<Part, RowsAtCompileTime, 1> col(int i) const { return Base::col(i); }
 
-    template<typename OtherDerived/*, int OtherMode*/>
+    template<typename OtherDerived>
     void swap(const MatrixBase<OtherDerived>& other)
     {
-      Part<SwapWrapper<MatrixType>,Mode>(SwapWrapper<MatrixType>(const_cast<MatrixType&>(m_matrix))).lazyAssign(other.derived());
+      Part<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived());
     }
 
   protected:
@@ -128,7 +128,7 @@ template<typename MatrixType, unsigned int Mode> class Part
     const typename MatrixType::Nested m_matrix;
 };
 
-/** \nonstableyet 
+/** \nonstableyet
   * \returns an expression of a triangular matrix extracted from the current matrix
   *
   * The parameter \a Mode can have the following values: \c UpperTriangular, \c StrictlyUpperTriangular, \c UnitUpperTriangular,
@@ -280,7 +280,7 @@ void Part<MatrixType, Mode>::lazyAssign(const Other& other)
     >::run(m_matrix.const_cast_derived(), other.derived());
 }
 
-/** \nonstableyet 
+/** \nonstableyet
   * \returns a lvalue pseudo-expression allowing to perform special operations on \c *this.
   *
   * The \a Mode parameter can have the following values: \c UpperTriangular, \c StrictlyUpperTriangular, \c LowerTriangular,

Eigen/src/Core/Product.h

@@ -299,7 +299,7 @@ template<typename OtherDerived>
 inline Derived &
 MatrixBase<Derived>::operator*=(const MatrixBase<OtherDerived> &other)
 {
-  return *this = *this * other;
+  return derived() = derived() * other.derived();
 }
 
 /***************************************************************************
@@ -762,7 +762,7 @@ 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()
   );
 }
 

Eigen/src/Core/SolveTriangular.h

@@ -221,6 +221,8 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpLo,ColMajor|IsDense>
 };
 
 /** "in-place" version of MatrixBase::solveTriangular() where the result is written in \a other
+  *
+  * \nonstableyet
   *
   * The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here.
   * This function will const_cast it, so constness isn't honored here.
@@ -250,6 +252,8 @@ void MatrixBase<Derived>::solveTriangularInPlace(const MatrixBase<OtherDerived>&
 }
 
 /** \returns the product of the inverse of \c *this with \a other, \a *this being triangular.
+  *
+  * \nonstableyet
   *
   * This function computes the inverse-matrix matrix product inverse(\c *this) * \a other.
   * The matrix \c *this must be triangular and invertible (i.e., all the coefficients of the

Eigen/src/Core/Transpose.h

@@ -125,7 +125,20 @@ template<typename MatrixType> class Transpose
   * Example: \include MatrixBase_transpose.cpp
   * Output: \verbinclude MatrixBase_transpose.out
   *
-  * \sa adjoint(), class DiagonalCoeffs */
+  * \warning If you want to replace a matrix by its own transpose, do \b NOT do this:
+  * \code
+  * m = m.transpose(); // bug!!! caused by aliasing effect
+  * \endcode
+  * Instead, use the transposeInPlace() method:
+  * \code
+  * m.transposeInPlace();
+  * \endcode
+  * which gives Eigen good opportunities for optimization, or alternatively you can also do:
+  * \code
+  * m = m.transpose().eval();
+  * \endcode
+  *
+  * \sa transposeInPlace(), adjoint() */
 template<typename Derived>
 inline Transpose<Derived>
 MatrixBase<Derived>::transpose()
@@ -133,7 +146,11 @@ MatrixBase<Derived>::transpose()
   return derived();
 }
 
-/** This is the const version of transpose(). \sa adjoint() */
+/** This is the const version of transpose().
+  *
+  * Make sure you read the warning for transpose() !
+  *
+  * \sa transposeInPlace(), adjoint() */
 template<typename Derived>
 inline const Transpose<Derived>
 MatrixBase<Derived>::transpose() const
@@ -146,6 +163,15 @@ MatrixBase<Derived>::transpose() const
   * Example: \include MatrixBase_adjoint.cpp
   * Output: \verbinclude MatrixBase_adjoint.out
   *
+  * \warning If you want to replace a matrix by its own adjoint, do \b NOT do this:
+  * \code
+  * m = m.adjoint(); // bug!!! caused by aliasing effect
+  * \endcode
+  * Instead, do:
+  * \code
+  * m = m.adjoint().eval();
+  * \endcode
+  *
   * \sa transpose(), conjugate(), class Transpose, class ei_scalar_conjugate_op */
 template<typename Derived>
 inline const typename MatrixBase<Derived>::AdjointReturnType

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,16 @@ 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; }
+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
 template<> EIGEN_STRONG_INLINE float  ei_pfirst<__m128>(const __m128&  a) { return _mm_cvtss_f32(a); }
 template<> EIGEN_STRONG_INLINE double ei_pfirst<__m128d>(const __m128d& a) { return _mm_cvtsd_f64(a); }
 template<> EIGEN_STRONG_INLINE int    ei_pfirst<__m128i>(const __m128i& a) { return _mm_cvtsi128_si32(a); }
+#endif
 
 #ifdef __SSE3__
 // TODO implement SSE2 versions as well as integer versions
@@ -308,4 +315,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,12 +30,48 @@
 
 #define EIGEN_WORLD_VERSION 2
 #define EIGEN_MAJOR_VERSION 0
-#define EIGEN_MINOR_VERSION 0
+#define EIGEN_MINOR_VERSION 6
 
 #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))))
 
+// 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(__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)
+  #define EIGEN_ARCH_WANTS_ALIGNMENT 1
+#else
+  #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
@@ -147,18 +183,25 @@ 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 (defined __GNUC__)
+#if !EIGEN_ALIGN
+#define EIGEN_ALIGN_128
+#elif (defined __GNUC__)
 #define EIGEN_ALIGN_128 __attribute__((aligned(16)))
 #elif (defined _MSC_VER)
 #define EIGEN_ALIGN_128 __declspec(align(16))
 #else
-#define EIGEN_ALIGN_128
+#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
@@ -173,18 +216,18 @@ using Eigen::ei_cos;
 template<typename OtherDerived> \
 EIGEN_STRONG_INLINE Derived& operator Op(const Eigen::MatrixBase<OtherDerived>& other) \
 { \
-  return Eigen::MatrixBase<Derived>::operator Op(other.derived()); \
+  return Base::operator Op(other.derived()); \
 } \
 EIGEN_STRONG_INLINE Derived& operator Op(const Derived& other) \
 { \
-  return Eigen::MatrixBase<Derived>::operator Op(other); \
+  return Base::operator Op(other); \
 }
 
 #define EIGEN_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, Op) \
 template<typename Other> \
 EIGEN_STRONG_INLINE Derived& operator Op(const Other& scalar) \
 { \
-  return Eigen::MatrixBase<Derived>::operator Op(scalar); \
+  return Base::operator Op(scalar); \
 }
 
 #define EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Derived) \
@@ -216,4 +259,15 @@ _EIGEN_GENERIC_PUBLIC_INTERFACE(Derived, Eigen::MatrixBase<Derived>)
 #define EIGEN_ENUM_MIN(a,b) (((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 !
+#define EIGEN_EMPTY
+
+// concatenate two tokens
+#define EIGEN_CAT2(a,b) a ## b
+#define EIGEN_CAT(a,b) EIGEN_CAT2(a,b)
+
+// convert a token to a string
+#define EIGEN_MAKESTRING2(a) #a
+#define EIGEN_MAKESTRING(a) EIGEN_MAKESTRING2(a)
+
 #endif // EIGEN_MACROS_H

Eigen/src/Core/util/Memory.h

@@ -2,7 +2,7 @@
 // for linear algebra. Eigen itself is part of the KDE project.
 //
 // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
-// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2008-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
 // Copyright (C) 2009 Kenneth Riddile <kfriddile@yahoo.com>
 //
 // Eigen is free software; you can redistribute it and/or
@@ -27,14 +27,13 @@
 #ifndef EIGEN_MEMORY_H
 #define EIGEN_MEMORY_H
 
-// for NetBSD I didn't see any clear statement in the docs, but Mark Davies is confident about this.
-#if defined(__APPLE__) || defined(__FreeBSD__) || defined(__NetBSD__) || defined(_WIN64)
+#if defined(__APPLE__) || defined(_WIN64)
   #define EIGEN_MALLOC_ALREADY_ALIGNED 1
 #else
   #define EIGEN_MALLOC_ALREADY_ALIGNED 0
 #endif
 
-#if (defined _GNU_SOURCE) || ((defined _XOPEN_SOURCE) && (_XOPEN_SOURCE >= 600))
+#if ((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
@@ -74,28 +73,23 @@ inline void* ei_aligned_malloc(size_t size)
     ei_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)");
   #endif
 
-  void *result;
-  #if EIGEN_HAS_POSIX_MEMALIGN && !EIGEN_MALLOC_ALREADY_ALIGNED
-    #ifdef EIGEN_EXCEPTIONS
-      const int failed =
-    #endif
-    posix_memalign(&result, 16, size);
+  void *result;  
+  #if !EIGEN_ALIGN
+    result = malloc(size);
+  #elif EIGEN_MALLOC_ALREADY_ALIGNED
+    result = 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_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;
@@ -142,7 +136,9 @@ template<typename T, bool Align> inline T* ei_conditional_aligned_new(size_t siz
   */
 inline void ei_aligned_free(void *ptr)
 {
-  #if EIGEN_MALLOC_ALREADY_ALIGNED
+  #if !EIGEN_ALIGN
+    free(ptr);
+  #elif EIGEN_MALLOC_ALREADY_ALIGNED
     free(ptr);
   #elif EIGEN_HAS_POSIX_MEMALIGN
     free(ptr);
@@ -233,60 +229,51 @@ inline static int ei_alignmentOffset(const Scalar* ptr, int maxOffset)
 #define ei_aligned_stack_delete(TYPE,PTR,SIZE) do {ei_delete_elements_of_array<TYPE>(PTR, SIZE); \
                                                    ei_aligned_stack_free(PTR,sizeof(TYPE)*SIZE);} while(0)
 
-    
-/** \brief Overloads the operator new and delete of the class Type with operators that are aligned if NeedsToAlign is true
-  *
-  * When Eigen's explicit vectorization is enabled, Eigen assumes that some fixed sizes types are aligned
-  * on a 16 bytes boundary. Those include all Matrix types having a sizeof multiple of 16 bytes, e.g.:
-  *  - Vector2d, Vector4f, Vector4i, Vector4d,
-  *  - Matrix2d, Matrix4f, Matrix4i, Matrix4d,
-  *  - etc.
-  * When an object is statically allocated, the compiler will automatically and always enforces 16 bytes
-  * alignment of the data when needed. However some troubles might appear when data are dynamically allocated.
-  * Let's pick an example:
-  * \code
-  * struct Foo {
-  *   char dummy;
-  *   Vector4f some_vector;
-  * };
-  * Foo obj1;                           // static allocation
-  * obj1.some_vector = Vector4f(..);    // =>   OK
-  *
-  * Foo *pObj2 = new Foo;               // dynamic allocation
-  * pObj2->some_vector = Vector4f(..);  // =>  !! might segfault !!
-  * \endcode
-  * Here, the problem is that operator new is not aware of the compile time alignment requirement of the
-  * type Vector4f (and hence of the type Foo). Therefore "new Foo" does not necessarily returns a 16 bytes
-  * aligned pointer. The purpose of the class WithAlignedOperatorNew is exactly to overcome this issue by
-  * overloading the operator new to return aligned data when the vectorization is enabled.
-  * Here is a similar safe example:
-  * \code
-  * struct Foo {
-  *   EIGEN_MAKE_ALIGNED_OPERATOR_NEW
-  *   char dummy;
-  *   Vector4f some_vector;
-  * };
-  * Foo *pObj2 = new Foo;               // dynamic allocation
-  * pObj2->some_vector = Vector4f(..);  // =>  SAFE !
-  * \endcode
-  *
-  * \sa class ei_new_allocator
-  */
-#define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign) \
-    void *operator new(size_t size) throw() { \
-      return Eigen::ei_conditional_aligned_malloc<NeedsToAlign>(size); \
-    } \
-    void *operator new[](size_t size) throw() { \
-      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; }
+
+#if EIGEN_ALIGN
+  #ifdef EIGEN_EXCEPTIONS
+    #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \
+      void* operator new(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(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) { \
+        return Eigen::ei_conditional_aligned_malloc<NeedsToAlign>(size); \
+      } \
+      void *operator new[](size_t size) { \
+        return Eigen::ei_conditional_aligned_malloc<NeedsToAlign>(size); \
+      } \
+      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(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)
+#endif
 
 #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(true)
 #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(Scalar,Size) \
   EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(((Size)!=Eigen::Dynamic) && ((sizeof(Scalar)*(Size))%16==0))
 
+
 /** \class aligned_allocator
 *
 * \brief stl compatible allocator to use with with 16 byte aligned types
@@ -318,34 +305,34 @@ public:
         typedef aligned_allocator<U> other;
     };
 
-    pointer address( reference value ) const 
+    pointer address( reference value ) const
     {
         return &value;
     }
 
-    const_pointer address( const_reference value ) const 
+    const_pointer address( const_reference value ) const
     {
         return &value;
     }
 
-    aligned_allocator() throw() 
+    aligned_allocator() throw()
     {
     }
 
-    aligned_allocator( const aligned_allocator& ) throw() 
+    aligned_allocator( const aligned_allocator& ) throw()
     {
     }
 
     template<class U>
-    aligned_allocator( const aligned_allocator<U>& ) throw() 
+    aligned_allocator( const aligned_allocator<U>& ) throw()
     {
     }
 
-    ~aligned_allocator() throw() 
+    ~aligned_allocator() throw()
     {
     }
 
-    size_type max_size() const throw() 
+    size_type max_size() const throw()
     {
         return std::numeric_limits<size_type>::max();
     }
@@ -356,20 +343,26 @@ public:
         return static_cast<pointer>( ei_aligned_malloc( num * sizeof(T) ) );
     }
 
-    void construct( pointer p, const T& value ) 
+    void construct( pointer p, const T& value )
     {
         ::new( p ) T( value );
     }
 
-    void destroy( pointer p ) 
+    void destroy( pointer p )
     {
         p->~T();
     }
 
-    void deallocate( pointer p, size_type /*num*/ ) 
+    void deallocate( pointer p, size_type /*num*/ )
     {
         ei_aligned_free( p );
     }
+    
+    bool operator!=(const aligned_allocator<T>& other) const
+    { return false; }
+    
+    bool operator==(const aligned_allocator<T>& other) const
+    { return true; }
 };
 
 #endif // EIGEN_MEMORY_H

Eigen/src/Core/util/StaticAssert.h

@@ -73,7 +73,9 @@
         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY,
         THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES,
         THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES,
-        INVALID_MATRIX_TEMPLATE_PARAMETERS
+        INVALID_MATRIX_TEMPLATE_PARAMETERS,
+        BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER,
+        THIS_METHOD_IS_ONLY_FOR_DIAGONAL_MATRIX
       };
     };
 

Eigen/src/Geometry/AlignedBox.h

@@ -25,7 +25,7 @@
 #ifndef EIGEN_ALIGNEDBOX_H
 #define EIGEN_ALIGNEDBOX_H
 
-/** \geometry_module \ingroup GeometryModule
+/** \geometry_module \ingroup Geometry_Module
   * \nonstableyet
   *
   * \class AlignedBox

Eigen/src/Geometry/AngleAxis.h

@@ -25,7 +25,7 @@
 #ifndef EIGEN_ANGLEAXIS_H
 #define EIGEN_ANGLEAXIS_H
 
-/** \geometry_module \ingroup GeometryModule
+/** \geometry_module \ingroup Geometry_Module
   *
   * \class AngleAxis
   *
@@ -158,10 +158,10 @@ public:
   { return m_axis.isApprox(other.m_axis, prec) && ei_isApprox(m_angle,other.m_angle, prec); }
 };
 
-/** \ingroup GeometryModule
+/** \ingroup Geometry_Module
   * single precision angle-axis type */
 typedef AngleAxis<float> AngleAxisf;
-/** \ingroup GeometryModule
+/** \ingroup Geometry_Module
   * double precision angle-axis type */
 typedef AngleAxis<double> AngleAxisd;
 

Eigen/src/Geometry/EulerAngles.h

@@ -25,7 +25,7 @@
 #ifndef EIGEN_EULERANGLES_H
 #define EIGEN_EULERANGLES_H
 
-/** \geometry_module \ingroup GeometryModule
+/** \geometry_module \ingroup Geometry_Module
   * \nonstableyet
   *
   * \returns the Euler-angles of the rotation matrix \c *this using the convention defined by the triplet (\a a0,\a a1,\a a2)
@@ -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

@@ -26,7 +26,7 @@
 #ifndef EIGEN_HYPERPLANE_H
 #define EIGEN_HYPERPLANE_H
 
-/** \geometry_module \ingroup GeometryModule
+/** \geometry_module \ingroup Geometry_Module
   *
   * \class Hyperplane
   *

Eigen/src/Geometry/ParametrizedLine.h

@@ -26,7 +26,7 @@
 #ifndef EIGEN_PARAMETRIZEDLINE_H
 #define EIGEN_PARAMETRIZEDLINE_H
 
-/** \geometry_module \ingroup GeometryModule
+/** \geometry_module \ingroup Geometry_Module
   *
   * \class ParametrizedLine
   *

Eigen/src/Geometry/Quaternion.h

@@ -30,7 +30,7 @@ template<typename Other,
          int OtherCols=Other::ColsAtCompileTime>
 struct ei_quaternion_assign_impl;
 
-/** \geometry_module \ingroup GeometryModule
+/** \geometry_module \ingroup Geometry_Module
   *
   * \class Quaternion
   *
@@ -61,12 +61,12 @@ template<typename _Scalar>
 class Quaternion : public RotationBase<Quaternion<_Scalar>,3>
 {
   typedef RotationBase<Quaternion<_Scalar>,3> Base;
-  
+
 public:
   EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,4)
 
   using Base::operator*;
-  
+
   /** the scalar type of the coefficients */
   typedef _Scalar Scalar;
 
@@ -112,10 +112,6 @@ public:
   /** Default constructor leaving the quaternion uninitialized. */
   inline Quaternion() {}
 
-  inline Quaternion(ei_constructor_without_unaligned_array_assert)
-    : m_coeffs(ei_constructor_without_unaligned_array_assert()) {}
-
-
   /** Constructs and initializes the quaternion \f$ w+xi+yj+zk \f$ from
     * its four coefficients \a w, \a x, \a y and \a z.
     *
@@ -217,28 +213,56 @@ public:
   bool isApprox(const Quaternion& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const
   { return m_coeffs.isApprox(other.m_coeffs, prec); }
 
-protected: 
+protected:
   Coefficients m_coeffs;
 };
 
-/** \ingroup GeometryModule
+/** \ingroup Geometry_Module
   * single precision quaternion type */
 typedef Quaternion<float> Quaternionf;
-/** \ingroup GeometryModule
+/** \ingroup Geometry_Module
   * 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) */
@@ -350,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
@@ -358,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;

Eigen/src/Geometry/Rotation2D.h

@@ -25,7 +25,7 @@
 #ifndef EIGEN_ROTATION2D_H
 #define EIGEN_ROTATION2D_H
 
-/** \geometry_module \ingroup GeometryModule
+/** \geometry_module \ingroup Geometry_Module
   *
   * \class Rotation2D
   *
@@ -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
@@ -125,10 +125,10 @@ public:
   { return ei_isApprox(m_angle,other.m_angle, prec); }
 };
 
-/** \ingroup GeometryModule
+/** \ingroup Geometry_Module
   * single precision 2D rotation type */
 typedef Rotation2D<float> Rotation2Df;
-/** \ingroup GeometryModule
+/** \ingroup Geometry_Module
   * double precision 2D rotation type */
 typedef Rotation2D<double> Rotation2Dd;
 

Eigen/src/Geometry/Scaling.h

@@ -25,7 +25,7 @@
 #ifndef EIGEN_SCALING_H
 #define EIGEN_SCALING_H
 
-/** \geometry_module \ingroup GeometryModule
+/** \geometry_module \ingroup Geometry_Module
   *
   * \class Scaling
   *
@@ -149,7 +149,7 @@ public:
 
 };
 
-/** \addtogroup GeometryModule */
+/** \addtogroup Geometry_Module */
 //@{
 typedef Scaling<float, 2> Scaling2f;
 typedef Scaling<double,2> Scaling2d;

Eigen/src/Geometry/Transform.h

@@ -43,7 +43,7 @@ template< typename Other,
           int OtherCols=Other::ColsAtCompileTime>
 struct ei_transform_product_impl;
 
-/** \geometry_module \ingroup GeometryModule
+/** \geometry_module \ingroup Geometry_Module
   *
   * \class Transform
   *
@@ -95,11 +95,8 @@ public:
   /** Default constructor without initialization of the coefficients. */
   inline Transform() { }
 
-  inline Transform(ei_constructor_without_unaligned_array_assert)
-    : m_matrix(ei_constructor_without_unaligned_array_assert()) {}
-
   inline Transform(const Transform& other)
-  { 
+  {
     m_matrix = other.m_matrix;
   }
 
@@ -201,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);
@@ -286,17 +287,21 @@ 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:
 
 };
 
-/** \ingroup GeometryModule */
+/** \ingroup Geometry_Module */
 typedef Transform<float,2> Transform2f;
-/** \ingroup GeometryModule */
+/** \ingroup Geometry_Module */
 typedef Transform<float,3> Transform3f;
-/** \ingroup GeometryModule */
+/** \ingroup Geometry_Module */
 typedef Transform<double,2> Transform2d;
-/** \ingroup GeometryModule */
+/** \ingroup Geometry_Module */
 typedef Transform<double,3> Transform3d;
 
 /**************************
@@ -338,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.
@@ -372,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
 
@@ -648,7 +653,7 @@ template<typename Scalar, int Dim>
 template<typename ScalingMatrixType, typename RotationMatrixType>
 void Transform<Scalar,Dim>::computeScalingRotation(ScalingMatrixType *scaling, RotationMatrixType *rotation) const
 {
-  linear().svd().computeScalingRotation(scaling, rotation);  
+  linear().svd().computeScalingRotation(scaling, rotation);
 }
 
 /** Convenient method to set \c *this from a position, orientation and scale

Eigen/src/Geometry/Translation.h

@@ -25,7 +25,7 @@
 #ifndef EIGEN_TRANSLATION_H
 #define EIGEN_TRANSLATION_H
 
-/** \geometry_module \ingroup GeometryModule
+/** \geometry_module \ingroup Geometry_Module
   *
   * \class Translation
   *
@@ -152,7 +152,7 @@ public:
 
 };
 
-/** \addtogroup GeometryModule */
+/** \addtogroup Geometry_Module */
 //@{
 typedef Translation<float, 2> Translation2f;
 typedef Translation<double,2> Translation2d;

Eigen/src/LU/Inverse.h

@@ -132,21 +132,31 @@ void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* resu
     // 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(1).swap(m.row(2));
+    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 two
+      // 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(1).swap(result->col(2));
+      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.
-      m.row(1).swap(m.row(2));
+      // 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(3));
+      result->col(1).swap(result->col(swap1with));
+      result->col(0).swap(result->col(swap0with));
     }
   }
 }

Eigen/src/LU/LU.h

@@ -323,6 +323,7 @@ template<typename MatrixType> class LU
     IntRowVectorType m_q;
     int m_det_pq;
     int m_rank;
+    RealScalar m_precision;
 };
 
 template<typename MatrixType>
@@ -335,6 +336,10 @@ LU<MatrixType>::LU(const MatrixType& matrix)
   const int size = matrix.diagonal().size();
   const int rows = matrix.rows();
   const int cols = matrix.cols();
+  
+  // this formula comes from experimenting (see "LU precision tuning" thread on the list)
+  // and turns out to be identical to Higham's formula used already in LDLt.
+  m_precision = machine_epsilon<Scalar>() * size;
 
   IntColVectorType rows_transpositions(matrix.rows());
   IntRowVectorType cols_transpositions(matrix.cols());
@@ -355,7 +360,7 @@ LU<MatrixType>::LU(const MatrixType& matrix)
     if(k==0) biggest = biggest_in_corner;
 
     // if the corner is negligible, then we have less than full rank, and we can finish early
-    if(ei_isMuchSmallerThan(biggest_in_corner, biggest))
+    if(ei_isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
     {
       m_rank = k;
       for(int i = k; i < size; i++)
@@ -503,10 +508,10 @@ 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() : 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))
+        if(!ei_isMuchSmallerThan(c.coeff(row,col), biggest_in_c, m_precision))
           return false;
   }
   m_lu.corner(TopLeft, m_rank, m_rank)

Eigen/src/LeastSquares/LeastSquares.h

@@ -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) 2006-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
@@ -57,7 +57,7 @@
     Vector3d coeffs; // will store the coefficients a, b, c
     linearRegression(
       5,
-      points,
+      &points,
       &coeffs,
       1 // the coord to express as a function of
         // the other ones. 0 means x, 1 means y, 2 means z.
@@ -80,11 +80,11 @@
                   This vector must be of the same type and size as the
                   data points. The meaning of its coords is as follows.
                   For brevity, let \f$n=Size\f$,
-                  \f$r_i=retCoefficients[i]\f$,
+                  \f$r_i=result[i]\f$,
                   and \f$f=funcOfOthers\f$. Denote by
                   \f$x_0,\ldots,x_{n-1}\f$
                   the n coordinates in the n-dimensional space.
-                  Then the result equation is:
+                  Then the resulting equation is:
                   \f[ x_f = r_0 x_0 + \cdots + r_{f-1}x_{f-1}
                    + r_{f+1}x_{f+1} + \cdots + r_{n-1}x_{n-1} + r_n. \f]
   * @param funcOfOthers Determines which coord to express as a function of the
@@ -101,31 +101,15 @@ void linearRegression(int numPoints,
                       int funcOfOthers )
 {
   typedef typename VectorType::Scalar Scalar;
-  EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType)
-  ei_assert(numPoints >= 1);
-  int size = points[0]->size();
-  ei_assert(funcOfOthers >= 0 && funcOfOthers < size);
+  typedef Hyperplane<Scalar, VectorType::SizeAtCompileTime> HyperplaneType;
+  const int size = points[0]->size();
   result->resize(size);
-
-  Matrix<Scalar, Dynamic, VectorType::SizeAtCompileTime,
-         Dynamic, VectorType::MaxSizeAtCompileTime, RowMajorBit>
-    m(numPoints, size);
-  if(funcOfOthers>0)
-    for(int i = 0; i < numPoints; ++i)
-      m.row(i).start(funcOfOthers) = points[i]->start(funcOfOthers);
-  if(funcOfOthers<size-1)
-    for(int i = 0; i < numPoints; ++i)
-      m.row(i).block(funcOfOthers, size-funcOfOthers-1)
-        = points[i]->end(size-funcOfOthers-1);
-  for(int i = 0; i < numPoints; ++i)
-    m.row(i).coeffRef(size-1) = Scalar(1);
-
-  VectorType v(size);
-  v.setZero();
-  for(int i = 0; i < numPoints; ++i)
-    v += m.row(i).adjoint() * points[i]->coeff(funcOfOthers);
-
-  ei_assert((m.adjoint()*m).lu().solve(v, result));
+  HyperplaneType h(size);
+  fitHyperplane(numPoints, points, &h);
+  for(int i = 0; i < funcOfOthers; i++)
+    result->coeffRef(i) = - h.coeffs()[i] / h.coeffs()[funcOfOthers];
+  for(int i = funcOfOthers; i < size; i++)
+    result->coeffRef(i) = - h.coeffs()[i+1] / h.coeffs()[funcOfOthers];
 }
 
 /** \ingroup LeastSquares_Module

Eigen/src/QR/EigenSolver.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>
 //
@@ -53,9 +53,18 @@ template<typename _MatrixType> class EigenSolver
     typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType;
     typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
 
+    /** 
+    * \brief Default Constructor.
+    *
+    * The default constructor is useful in cases in which the user intends to
+    * perform decompositions via EigenSolver::compute(const MatrixType&).
+    */
+    EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) {}
+
     EigenSolver(const MatrixType& matrix)
       : m_eivec(matrix.rows(), matrix.cols()),
-        m_eivalues(matrix.cols())
+        m_eivalues(matrix.cols()),
+        m_isInitialized(false)
     {
       compute(matrix);
     }
@@ -94,12 +103,20 @@ template<typename _MatrixType> class EigenSolver
       *
       * \sa pseudoEigenvalueMatrix()
       */
-    const MatrixType& pseudoEigenvectors() const { return m_eivec; }
+    const MatrixType& pseudoEigenvectors() const 
+    { 
+      ei_assert(m_isInitialized && "EigenSolver is not initialized.");
+      return m_eivec; 
+    }
 
     MatrixType pseudoEigenvalueMatrix() const;
 
     /** \returns the eigenvalues as a column vector */
-    EigenvalueType eigenvalues() const { return m_eivalues; }
+    EigenvalueType eigenvalues() const 
+    { 
+      ei_assert(m_isInitialized && "EigenSolver is not initialized.");
+      return m_eivalues; 
+    }
 
     void compute(const MatrixType& matrix);
 
@@ -111,6 +128,7 @@ template<typename _MatrixType> class EigenSolver
   protected:
     MatrixType m_eivec;
     EigenvalueType m_eivalues;
+    bool m_isInitialized;
 };
 
 /** \returns the real block diagonal matrix D of the eigenvalues.
@@ -120,6 +138,7 @@ template<typename _MatrixType> class EigenSolver
 template<typename MatrixType>
 MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
 {
+  ei_assert(m_isInitialized && "EigenSolver is not initialized.");
   int n = m_eivec.cols();
   MatrixType matD = MatrixType::Zero(n,n);
   for (int i=0; i<n; ++i)
@@ -143,6 +162,7 @@ MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
 template<typename MatrixType>
 typename EigenSolver<MatrixType>::EigenvectorType EigenSolver<MatrixType>::eigenvectors(void) const
 {
+  ei_assert(m_isInitialized && "EigenSolver is not initialized.");
   int n = m_eivec.cols();
   EigenvectorType matV(n,n);
   for (int j=0; j<n; ++j)
@@ -183,6 +203,8 @@ void EigenSolver<MatrixType>::compute(const MatrixType& matrix)
 
   // Reduce Hessenberg to real Schur form.
   hqr2(matH);
+
+  m_isInitialized = true;
 }
 
 // Nonsymmetric reduction to Hessenberg form.

Eigen/src/QR/QR.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>
 //
@@ -49,51 +49,146 @@ template<typename MatrixType> class QR
     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixTypeR;
     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
 
+    /** 
+    * \brief Default Constructor.
+    *
+    * The default constructor is useful in cases in which the user intends to
+    * perform decompositions via QR::compute(const MatrixType&).
+    */
+    QR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {}
+
     QR(const MatrixType& matrix)
       : m_qr(matrix.rows(), matrix.cols()),
-        m_hCoeffs(matrix.cols())
+        m_hCoeffs(matrix.cols()),
+        m_isInitialized(false)
     {
-      _compute(matrix);
+      compute(matrix);
+    }
+    
+    /** \deprecated use isInjective()
+      * \returns whether or not the matrix is of full rank
+      *
+      * \note Since the rank is computed only once, i.e. the first time it is needed, this
+      *       method almost does not perform any further computation.
+      */
+    EIGEN_DEPRECATED bool isFullRank() const 
+    { 
+      ei_assert(m_isInitialized && "QR is not initialized.");
+      return rank() == m_qr.cols(); 
+    }
+    
+    /** \returns the rank of the matrix of which *this is the QR decomposition.
+      *
+      * \note Since the rank is computed only once, i.e. the first time it is needed, this
+      *       method almost does not perform any further computation.
+      */
+    int rank() const;
+    
+    /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
+      *
+      * \note Since the rank is computed only once, i.e. the first time it is needed, this
+      *       method almost does not perform any further computation.
+      */
+    inline int dimensionOfKernel() const
+    {
+      ei_assert(m_isInitialized && "QR is not initialized.");
+      return m_qr.cols() - rank();
+    }
+    
+    /** \returns true if the matrix of which *this is the QR decomposition represents an injective
+      *          linear map, i.e. has trivial kernel; false otherwise.
+      *
+      * \note Since the rank is computed only once, i.e. the first time it is needed, this
+      *       method almost does not perform any further computation.
+      */
+    inline bool isInjective() const
+    {
+      ei_assert(m_isInitialized && "QR is not initialized.");
+      return rank() == m_qr.cols();
+    }
+    
+    /** \returns true if the matrix of which *this is the QR decomposition represents a surjective
+      *          linear map; false otherwise.
+      *
+      * \note Since the rank is computed only once, i.e. the first time it is needed, this
+      *       method almost does not perform any further computation.
+      */
+    inline bool isSurjective() const
+    {
+      ei_assert(m_isInitialized && "QR is not initialized.");
+      return rank() == m_qr.rows();
     }
 
-    /** \returns whether or not the matrix is of full rank */
-    bool isFullRank() const { return rank() == std::min(m_qr.rows(),m_qr.cols()); }
+    /** \returns true if the matrix of which *this is the QR decomposition is invertible.
+      *
+      * \note Since the rank is computed only once, i.e. the first time it is needed, this
+      *       method almost does not perform any further computation.
+      */
+    inline bool isInvertible() const
+    {
+      ei_assert(m_isInitialized && "QR is not initialized.");
+      return isInjective() && isSurjective();
+    }
     
-    int rank() const;
-
     /** \returns a read-only expression of the matrix R of the actual the QR decomposition */
     const Part<NestByValue<MatrixRBlockType>, UpperTriangular>
     matrixR(void) const
     {
+      ei_assert(m_isInitialized && "QR is not initialized.");
       int cols = m_qr.cols();
       return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part<UpperTriangular>();
     }
 
+    /** This method finds a solution x to the equation Ax=b, where A is the matrix of which
+      * *this is the QR decomposition, if any exists.
+      *
+      * \param b the right-hand-side of the equation to solve.
+      *
+      * \param result a pointer to the vector/matrix in which to store the solution, if any exists.
+      *          Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols().
+      *          If no solution exists, *result is left with undefined coefficients.
+      *
+      * \returns true if any solution exists, false if no solution exists.
+      *
+      * \note If there exist more than one solution, this method will arbitrarily choose one.
+      *       If you need a complete analysis of the space of solutions, take the one solution obtained
+      *       by this method and add to it elements of the kernel, as determined by kernel().
+      *
+      * \note The case where b is a matrix is not yet implemented. Also, this
+      *       code is space inefficient.
+      *
+      * Example: \include QR_solve.cpp
+      * Output: \verbinclude QR_solve.out
+      *
+      * \sa MatrixBase::solveTriangular(), kernel(), computeKernel(), inverse(), computeInverse()
+      */
+    template<typename OtherDerived, typename ResultType>
+    bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
+
     MatrixType matrixQ(void) const;
 
-  private:
-
-    void _compute(const MatrixType& matrix);
+    void compute(const MatrixType& matrix);
 
   protected:
     MatrixType m_qr;
     VectorType m_hCoeffs;
     mutable int m_rank;
     mutable bool m_rankIsUptodate;
+    bool m_isInitialized;
 };
 
 /** \returns the rank of the matrix of which *this is the QR decomposition. */
 template<typename MatrixType>
 int QR<MatrixType>::rank() const
 {
+  ei_assert(m_isInitialized && "QR is not initialized.");
   if (!m_rankIsUptodate)
   {
-    RealScalar maxCoeff = m_qr.diagonal().maxCoeff();
-    int n = std::min(m_qr.rows(),m_qr.cols());
-    m_rank = n;
-    for (int i=0; i<n; ++i)
-      if (ei_isMuchSmallerThan(m_qr.diagonal().coeff(i), maxCoeff))
-        --m_rank;
+    RealScalar maxCoeff = m_qr.diagonal().cwise().abs().maxCoeff();
+    int n = m_qr.cols();
+    m_rank = 0;
+    while(m_rank<n && !ei_isMuchSmallerThan(m_qr.diagonal().coeff(m_rank), maxCoeff))
+      ++m_rank;
     m_rankIsUptodate = true;
   }
   return m_rank;
@@ -102,12 +197,15 @@ int QR<MatrixType>::rank() const
 #ifndef EIGEN_HIDE_HEAVY_CODE
 
 template<typename MatrixType>
-void QR<MatrixType>::_compute(const MatrixType& matrix)
-{
+void QR<MatrixType>::compute(const MatrixType& matrix)
+{ 
   m_rankIsUptodate = false;
   m_qr = matrix;
+  m_hCoeffs.resize(matrix.cols());
+
   int rows = matrix.rows();
   int cols = matrix.cols();
+  RealScalar eps2 = precision<RealScalar>()*precision<RealScalar>();
 
   for (int k = 0; k < cols; ++k)
   {
@@ -132,7 +230,8 @@ void QR<MatrixType>::_compute(const MatrixType& matrix)
         m_hCoeffs.coeffRef(k) = 0;
       }
     }
-    else if ( (!ei_isMuchSmallerThan(beta=m_qr.col(k).end(remainingSize-1).squaredNorm(),static_cast<Scalar>(1))) || ei_imag(v0)==0 )
+    else if ((beta=m_qr.col(k).end(remainingSize-1).squaredNorm())>eps2)
+    // FIXME what about ei_imag(v0) ??
     {
       // form k-th Householder vector
       beta = ei_sqrt(ei_abs2(v0)+beta);
@@ -158,12 +257,46 @@ void QR<MatrixType>::_compute(const MatrixType& matrix)
       m_hCoeffs.coeffRef(k) = 0;
     }
   }
+  m_isInitialized = true;
+}
+
+template<typename MatrixType>
+template<typename OtherDerived, typename ResultType>
+bool QR<MatrixType>::solve(
+  const MatrixBase<OtherDerived>& b,
+  ResultType *result
+) const
+{
+  ei_assert(m_isInitialized && "QR is not initialized.");
+  const int rows = m_qr.rows();
+  ei_assert(b.rows() == rows);
+  result->resize(rows, 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(!isSurjective())
+  {
+    // is result is in the image of R ?
+    RealScalar biggest_in_res = result->corner(TopLeft, m_rank, result->cols()).cwise().abs().maxCoeff();
+    for(int col = 0; col < result->cols(); ++col)
+      for(int row = m_rank; row < result->rows(); ++row)
+        if(!ei_isMuchSmallerThan(result->coeff(row,col), biggest_in_res))
+          return false;
+  }
+  m_qr.corner(TopLeft, m_rank, m_rank)
+      .template marked<UpperTriangular>()
+      .solveTriangularInPlace(result->corner(TopLeft, m_rank, result->cols()));
+
+  return true;
 }
 
 /** \returns the matrix Q */
 template<typename MatrixType>
-MatrixType QR<MatrixType>::matrixQ(void) const
+MatrixType QR<MatrixType>::matrixQ() const
 {
+  ei_assert(m_isInitialized && "QR is not initialized.");
   // compute the product Q_0 Q_1 ... Q_n-1,
   // where Q_k is the k-th Householder transformation I - h_k v_k v_k'
   // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]

Eigen/src/QR/SelfAdjointEigenSolver.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>
 //
@@ -52,8 +52,8 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     typedef Tridiagonalization<MatrixType> TridiagonalizationType;
 
     SelfAdjointEigenSolver()
-        : m_eivec(Size, Size),
-          m_eivalues(Size)
+        : m_eivec(int(Size), int(Size)),
+          m_eivalues(int(Size))
     {
       ei_assert(Size!=Dynamic);
     }
@@ -189,6 +189,14 @@ void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool
   assert(matrix.cols() == matrix.rows());
   int n = matrix.cols();
   m_eivalues.resize(n,1);
+
+  if(n==1)
+  {
+    m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0));
+    m_eivec.setOnes();
+    return;
+  }
+
   m_eivec = matrix;
 
   // FIXME, should tridiag be a local variable of this function or an attribute of SelfAdjointEigenSolver ?

Eigen/src/QR/Tridiagonalization.h

@@ -201,6 +201,7 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
     // squared norm of the vector v skipping the first element
     RealScalar v1norm2 = matA.col(i).end(n-(i+2)).squaredNorm();
 
+    // FIXME comparing against 1
     if (ei_isMuchSmallerThan(v1norm2,static_cast<Scalar>(1)))
     {
       hCoeffs.coeffRef(i) = 0.;
@@ -331,7 +332,8 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
     if (ei_real(v0)>=0.)
       beta = -beta;
     matA.col(i).coeffRef(i+1) = beta;
-    hCoeffs.coeffRef(i) = (beta - v0) / beta;
+    if(ei_isMuchSmallerThan(beta, Scalar(1))) hCoeffs.coeffRef(i) = Scalar(0);
+    else hCoeffs.coeffRef(i) = (beta - v0) / beta;
   }
   else
   {

Eigen/src/Sparse/CMakeLists.txt

@@ -2,5 +2,5 @@ FILE(GLOB Eigen_Sparse_SRCS "*.h")
 
 INSTALL(FILES 
   ${Eigen_Sparse_SRCS}
-  DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Sparse
+  DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Sparse COMPONENT Devel
   )

Eigen/src/Sparse/CholmodSupport.h

@@ -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) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -54,16 +54,17 @@ void ei_cholmod_configure_matrix(CholmodType& mat)
   }
 }
 
-template<typename Scalar, int Flags>
-cholmod_sparse SparseMatrixBase<Scalar,Flags>::asCholmodMatrix()
+template<typename Derived>
+cholmod_sparse SparseMatrixBase<Derived>::asCholmodMatrix()
 {
+  typedef typename Derived::Scalar Scalar;
   cholmod_sparse res;
   res.nzmax   = nonZeros();
   res.nrow    = rows();;
   res.ncol    = cols();
-  res.p       = _outerIndexPtr();
-  res.i       = _innerIndexPtr();
-  res.x       = _valuePtr();
+  res.p       = derived()._outerIndexPtr();
+  res.i       = derived()._innerIndexPtr();
+  res.x       = derived()._valuePtr();
   res.xtype = CHOLMOD_REAL;
   res.itype = CHOLMOD_INT;
   res.sorted = 1;
@@ -73,11 +74,11 @@ cholmod_sparse SparseMatrixBase<Scalar,Flags>::asCholmodMatrix()
 
   ei_cholmod_configure_matrix<Scalar>(res);
 
-  if (Flags & SelfAdjoint)
+  if (Derived::Flags & SelfAdjoint)
   {
-    if (Flags & UpperTriangular)
+    if (Derived::Flags & UpperTriangular)
       res.stype = 1;
-    else if (Flags & LowerTriangular)
+    else if (Derived::Flags & LowerTriangular)
       res.stype = -1;
     else
       res.stype = 0;
@@ -108,14 +109,14 @@ cholmod_dense ei_cholmod_map_eigen_to_dense(MatrixBase<Derived>& mat)
 }
 
 template<typename Scalar, int Flags>
-MappedSparseMatrix<Scalar,Flags>::MappedSparseMatrix(taucs_ccs_matrix& taucsMat)
+MappedSparseMatrix<Scalar,Flags>::MappedSparseMatrix(cholmod_sparse& cm)
 {
   m_innerSize = cm.nrow;
   m_outerSize = cm.ncol;
   m_outerIndex = reinterpret_cast<int*>(cm.p);
   m_innerIndices = reinterpret_cast<int*>(cm.i);
   m_values = reinterpret_cast<Scalar*>(cm.x);
-  m_nnz = res.m_outerIndex[cm.ncol]);
+  m_nnz = m_outerIndex[cm.ncol];
 }
 
 template<typename MatrixType>
@@ -123,8 +124,8 @@ class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType>
 {
   protected:
     typedef SparseLLT<MatrixType> Base;
-    using typename Base::Scalar;
-    using Base::RealScalar;
+    typedef typename Base::Scalar Scalar;
+    typedef typename Base::RealScalar RealScalar;
     using Base::MatrixLIsDirty;
     using Base::SupernodalFactorIsDirty;
     using Base::m_flags;
@@ -205,7 +206,7 @@ SparseLLT<MatrixType,Cholmod>::matrixL() const
     ei_assert(!(m_status & SupernodalFactorIsDirty));
 
     cholmod_sparse* cmRes = cholmod_factor_to_sparse(m_cholmodFactor, &m_cholmod);
-    const_cast<typename Base::CholMatrixType&>(m_matrix) = Base::CholMatrixType::Map(*cmRes);
+    const_cast<typename Base::CholMatrixType&>(m_matrix) = MappedSparseMatrix<Scalar>(*cmRes);
     free(cmRes);
 
     m_status = (m_status & ~MatrixLIsDirty);

Eigen/src/Sparse/DynamicSparseMatrix.h

@@ -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) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -83,6 +83,9 @@ class DynamicSparseMatrix
     inline int innerSize() const { return m_innerSize; }
     inline int outerSize() const { return m_data.size(); }
     inline int innerNonZeros(int j) const { return m_data[j].size(); }
+    
+    std::vector<CompressedStorage<Scalar> >& _data() { return m_data; }
+    const std::vector<CompressedStorage<Scalar> >& _data() const { return m_data; }
 
     /** \returns the coefficient value at given position \a row, \a col
       * This operation involes a log(rho*outer_size) binary search.
@@ -125,11 +128,14 @@ class DynamicSparseMatrix
     /** Set the matrix to zero and reserve the memory for \a reserveSize nonzero coefficients. */
     inline void startFill(int reserveSize = 1000)
     {
-      int reserveSizePerVector = std::max(reserveSize/outerSize(),4);
-      for (int j=0; j<outerSize(); ++j)
+      if (outerSize()>0)
       {
-        m_data[j].clear();
-        m_data[j].reserve(reserveSizePerVector);
+        int reserveSizePerVector = std::max(reserveSize/outerSize(),4);
+        for (int j=0; j<outerSize(); ++j)
+        {
+          m_data[j].clear();
+          m_data[j].reserve(reserveSizePerVector);
+        }
       }
     }
 
@@ -215,7 +221,7 @@ class DynamicSparseMatrix
     }
 
     inline DynamicSparseMatrix()
-      : m_innerSize(0)
+      : m_innerSize(0), m_data(0)
     {
       ei_assert(innerSize()==0 && outerSize()==0);
     }

Eigen/src/Sparse/MappedSparseMatrix.h

@@ -65,10 +65,10 @@ class MappedSparseMatrix
 
     //----------------------------------------
     // direct access interface
-    inline const Scalar* _valuePtr() const { return &m_values; }
-    inline Scalar* _valuePtr() { return &m_values; }
+    inline const Scalar* _valuePtr() const { return m_values; }
+    inline Scalar* _valuePtr() { return m_values; }
 
-    inline const int* _innerIndexPtr() const { return &m_innerIndices; }
+    inline const int* _innerIndexPtr() const { return m_innerIndices; }
     inline int* _innerIndexPtr() { return m_innerIndices; }
 
     inline const int* _outerIndexPtr() const { return m_outerIndex; }
@@ -108,7 +108,7 @@ class MappedSparseMatrix
       ei_assert((*r==inner) && (id<end) && "coeffRef cannot be called on a zero coefficient");
       return m_values[id];
     }
-    
+
     class InnerIterator;
 
     /** \returns the number of non zero coefficients */
@@ -140,21 +140,25 @@ class MappedSparseMatrix<Scalar,_Flags>::InnerIterator
 {
   public:
     InnerIterator(const MappedSparseMatrix& mat, int outer)
-      : m_matrix(mat), m_outer(outer), m_id(mat._outerIndexPtr[outer]), m_start(m_id), m_end(mat._outerIndexPtr[outer+1])
+      : m_matrix(mat),
+        m_outer(outer),
+        m_id(mat._outerIndexPtr()[outer]),
+        m_start(m_id),
+        m_end(mat._outerIndexPtr()[outer+1])
     {}
 
     template<unsigned int Added, unsigned int Removed>
     InnerIterator(const Flagged<MappedSparseMatrix,Added,Removed>& mat, int outer)
-      : m_matrix(mat._expression()), m_id(m_matrix._outerIndexPtr[outer]),
-        m_start(m_id), m_end(m_matrix._outerIndexPtr[outer+1])
+      : m_matrix(mat._expression()), m_id(m_matrix._outerIndexPtr()[outer]),
+        m_start(m_id), m_end(m_matrix._outerIndexPtr()[outer+1])
     {}
 
     inline InnerIterator& operator++() { m_id++; return *this; }
 
-    inline Scalar value() const { return m_matrix.m_valuePtr[m_id]; }
-    inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix._valuePtr[m_id]); }
+    inline Scalar value() const { return m_matrix._valuePtr()[m_id]; }
+    inline Scalar& valueRef() { return const_cast<Scalar&>(m_matrix._valuePtr()[m_id]); }
 
-    inline int index() const { return m_matrix._innerIndexPtr(m_id); }
+    inline int index() const { return m_matrix._innerIndexPtr()[m_id]; }
     inline int row() const { return IsRowMajor ? m_outer : index(); }
     inline int col() const { return IsRowMajor ? index() : m_outer; }
 

Eigen/src/Sparse/SparseBlock.h

@@ -26,70 +26,246 @@
 #ifndef EIGEN_SPARSE_BLOCK_H
 #define EIGEN_SPARSE_BLOCK_H
 
-template<typename MatrixType>
-struct ei_traits<SparseInnerVector<MatrixType> >
+template<typename MatrixType, int Size>
+struct ei_traits<SparseInnerVectorSet<MatrixType, Size> >
 {
   typedef typename ei_traits<MatrixType>::Scalar Scalar;
   enum {
     IsRowMajor = (int(MatrixType::Flags)&RowMajorBit)==RowMajorBit,
     Flags = MatrixType::Flags,
-    RowsAtCompileTime = IsRowMajor ? 1 : MatrixType::RowsAtCompileTime,
-    ColsAtCompileTime = IsRowMajor ? MatrixType::ColsAtCompileTime : 1,
+    RowsAtCompileTime = IsRowMajor ? Size : MatrixType::RowsAtCompileTime,
+    ColsAtCompileTime = IsRowMajor ? MatrixType::ColsAtCompileTime : Size,
     CoeffReadCost = MatrixType::CoeffReadCost
   };
 };
 
-template<typename MatrixType>
-class SparseInnerVector : ei_no_assignment_operator,
-  public SparseMatrixBase<SparseInnerVector<MatrixType> >
+template<typename MatrixType, int Size>
+class SparseInnerVectorSet : ei_no_assignment_operator,
+  public SparseMatrixBase<SparseInnerVectorSet<MatrixType, Size> >
 {
-    enum {
-      IsRowMajor = ei_traits<SparseInnerVector>::IsRowMajor
-    };
-public:
+    enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor };
+  public:
 
-    EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVector)
-    class InnerIterator;
-
-    inline SparseInnerVector(const MatrixType& matrix, int outer)
-      : m_matrix(matrix), m_outer(outer)
+    EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet)
+    class InnerIterator: public MatrixType::InnerIterator
     {
+      public:
+        inline InnerIterator(const SparseInnerVectorSet& xpr, int outer)
+          : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer)
+        {}
+    };
+
+    inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize)
+      : m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize)
+    {
+      ei_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) );
+    }
+
+    inline SparseInnerVectorSet(const MatrixType& matrix, int outer)
+      : m_matrix(matrix), m_outerStart(outer), m_outerSize(Size)
+    {
+      ei_assert(Size!=Dynamic);
       ei_assert( (outer>=0) && (outer<matrix.outerSize()) );
     }
 
-    EIGEN_STRONG_INLINE int rows() const { return IsRowMajor ? 1 : m_matrix.rows(); }
-    EIGEN_STRONG_INLINE int cols() const { return IsRowMajor ? m_matrix.cols() : 1; }
+//     template<typename OtherDerived>
+//     inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
+//     {
+//       return *this;
+//     }
+
+//     template<typename Sparse>
+//     inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
+//     {
+//       return *this;
+//     }
+
+    EIGEN_STRONG_INLINE int rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
+    EIGEN_STRONG_INLINE int cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
 
   protected:
 
     const typename MatrixType::Nested m_matrix;
-    int m_outer;
+    int m_outerStart;
+    const ei_int_if_dynamic<Size> m_outerSize;
 
 };
 
-template<typename MatrixType>
-class SparseInnerVector<MatrixType>::InnerIterator : public MatrixType::InnerIterator
+/***************************************************************************
+* specialisation for DynamicSparseMatrix
+***************************************************************************/
+
+template<typename _Scalar, int _Options, int Size>
+class SparseInnerVectorSet<DynamicSparseMatrix<_Scalar, _Options>, Size>
+  : public SparseMatrixBase<SparseInnerVectorSet<DynamicSparseMatrix<_Scalar, _Options>, Size> >
 {
-public:
-  inline InnerIterator(const SparseInnerVector& xpr, int outer=0)
-    : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outer)
-  {
-    ei_assert(outer==0);
-  }
+    typedef DynamicSparseMatrix<_Scalar, _Options> MatrixType;
+    enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor };
+  public:
+
+    EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet)
+    class InnerIterator: public MatrixType::InnerIterator
+    {
+      public:
+        inline InnerIterator(const SparseInnerVectorSet& xpr, int outer)
+          : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer)
+        {}
+    };
+
+    inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize)
+      : m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize)
+    {
+      ei_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) );
+    }
+
+    inline SparseInnerVectorSet(const MatrixType& matrix, int outer)
+      : m_matrix(matrix), m_outerStart(outer), m_outerSize(Size)
+    {
+      ei_assert(Size!=Dynamic);
+      ei_assert( (outer>=0) && (outer<matrix.outerSize()) );
+    }
+
+    template<typename OtherDerived>
+    inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
+    {
+      if (IsRowMajor != ((OtherDerived::Flags&RowMajorBit)==RowMajorBit))
+      {
+        // need to transpose => perform a block evaluation followed by a big swap
+        DynamicSparseMatrix<Scalar,IsRowMajor?RowMajorBit:0> aux(other);
+        *this = aux.markAsRValue();
+      }
+      else
+      {
+        // evaluate/copy vector per vector
+        for (int j=0; j<m_outerSize.value(); ++j)
+        {
+          SparseVector<Scalar,IsRowMajor ? RowMajorBit : 0> aux(other.innerVector(j));
+          m_matrix.const_cast_derived()._data()[m_outerStart+j].swap(aux._data());
+        }
+      }
+      return *this;
+    }
+
+    inline SparseInnerVectorSet& operator=(const SparseInnerVectorSet& other)
+    {
+      return operator=<SparseInnerVectorSet>(other);
+    }
+
+//     template<typename Sparse>
+//     inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
+//     {
+//       return *this;
+//     }
+
+    EIGEN_STRONG_INLINE int rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
+    EIGEN_STRONG_INLINE int cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
+
+  protected:
+
+    const typename MatrixType::Nested m_matrix;
+    int m_outerStart;
+    const ei_int_if_dynamic<Size> m_outerSize;
+
 };
 
+
+/***************************************************************************
+* specialisation for SparseMatrix
+***************************************************************************/
+/*
+template<typename _Scalar, int _Options, int Size>
+class SparseInnerVectorSet<SparseMatrix<_Scalar, _Options>, Size>
+  : public SparseMatrixBase<SparseInnerVectorSet<SparseMatrix<_Scalar, _Options>, Size> >
+{
+    typedef DynamicSparseMatrix<_Scalar, _Options> MatrixType;
+    enum { IsRowMajor = ei_traits<SparseInnerVectorSet>::IsRowMajor };
+  public:
+
+    EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVectorSet)
+    class InnerIterator: public MatrixType::InnerIterator
+    {
+      public:
+        inline InnerIterator(const SparseInnerVectorSet& xpr, int outer)
+          : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer)
+        {}
+    };
+
+    inline SparseInnerVectorSet(const MatrixType& matrix, int outerStart, int outerSize)
+      : m_matrix(matrix), m_outerStart(outerStart), m_outerSize(outerSize)
+    {
+      ei_assert( (outerStart>=0) && ((outerStart+outerSize)<=matrix.outerSize()) );
+    }
+
+    inline SparseInnerVectorSet(const MatrixType& matrix, int outer)
+      : m_matrix(matrix), m_outerStart(outer)
+    {
+      ei_assert(Size==1);
+      ei_assert( (outer>=0) && (outer<matrix.outerSize()) );
+    }
+
+    template<typename OtherDerived>
+    inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
+    {
+      if (IsRowMajor != ((OtherDerived::Flags&RowMajorBit)==RowMajorBit))
+      {
+        // need to transpose => perform a block evaluation followed by a big swap
+        DynamicSparseMatrix<Scalar,IsRowMajor?RowMajorBit:0> aux(other);
+        *this = aux.markAsRValue();
+      }
+      else
+      {
+        // evaluate/copy vector per vector
+        for (int j=0; j<m_outerSize.value(); ++j)
+        {
+          SparseVector<Scalar,IsRowMajor ? RowMajorBit : 0> aux(other.innerVector(j));
+          m_matrix.const_cast_derived()._data()[m_outerStart+j].swap(aux._data());
+        }
+      }
+      return *this;
+    }
+
+    inline SparseInnerVectorSet& operator=(const SparseInnerVectorSet& other)
+    {
+      return operator=<SparseInnerVectorSet>(other);
+    }
+
+    inline const Scalar* _valuePtr() const
+    { return m_matrix._valuePtr() + m_matrix._outerIndexPtr()[m_outerStart]; }
+    inline const int* _innerIndexPtr() const
+    { return m_matrix._innerIndexPtr() + m_matrix._outerIndexPtr()[m_outerStart]; }
+    inline const int* _outerIndexPtr() const { return m_matrix._outerIndexPtr() + m_outerStart; }
+
+//     template<typename Sparse>
+//     inline SparseInnerVectorSet& operator=(const SparseMatrixBase<OtherDerived>& other)
+//     {
+//       return *this;
+//     }
+
+    EIGEN_STRONG_INLINE int rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); }
+    EIGEN_STRONG_INLINE int cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); }
+
+  protected:
+
+    const typename MatrixType::Nested m_matrix;
+    int m_outerStart;
+    const ei_int_if_dynamic<Size> m_outerSize;
+
+};
+*/
+//----------
+
 /** \returns the i-th row of the matrix \c *this. For row-major matrix only. */
 template<typename Derived>
-SparseInnerVector<Derived> SparseMatrixBase<Derived>::row(int i)
+SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::row(int i)
 {
   EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
   return innerVector(i);
 }
 
-/** \returns the i-th row of the matrix \c *this. For row-major matrix only. 
+/** \returns the i-th row of the matrix \c *this. For row-major matrix only.
   * (read-only version) */
 template<typename Derived>
-const SparseInnerVector<Derived> SparseMatrixBase<Derived>::row(int i) const
+const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::row(int i) const
 {
   EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
   return innerVector(i);
@@ -97,18 +273,18 @@ const SparseInnerVector<Derived> SparseMatrixBase<Derived>::row(int i) const
 
 /** \returns the i-th column of the matrix \c *this. For column-major matrix only. */
 template<typename Derived>
-SparseInnerVector<Derived> SparseMatrixBase<Derived>::col(int i)
+SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::col(int i)
 {
-  EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
+  EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
   return innerVector(i);
 }
 
-/** \returns the i-th column of the matrix \c *this. For column-major matrix only. 
+/** \returns the i-th column of the matrix \c *this. For column-major matrix only.
   * (read-only version) */
 template<typename Derived>
-const SparseInnerVector<Derived> SparseMatrixBase<Derived>::col(int i) const
+const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::col(int i) const
 {
-  EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
+  EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
   return innerVector(i);
 }
 
@@ -116,15 +292,65 @@ const SparseInnerVector<Derived> SparseMatrixBase<Derived>::col(int i) const
   * is col-major (resp. row-major).
   */
 template<typename Derived>
-SparseInnerVector<Derived> SparseMatrixBase<Derived>::innerVector(int outer)
-{ return SparseInnerVector<Derived>(derived(), outer); }
+SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::innerVector(int outer)
+{ return SparseInnerVectorSet<Derived,1>(derived(), outer); }
 
 /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this
   * is col-major (resp. row-major). Read-only.
   */
 template<typename Derived>
-const SparseInnerVector<Derived> SparseMatrixBase<Derived>::innerVector(int outer) const
-{ return SparseInnerVector<Derived>(derived(), outer); }
+const SparseInnerVectorSet<Derived,1> SparseMatrixBase<Derived>::innerVector(int outer) const
+{ return SparseInnerVectorSet<Derived,1>(derived(), outer); }
+
+//----------
+
+/** \returns the i-th row of the matrix \c *this. For row-major matrix only. */
+template<typename Derived>
+SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::subrows(int start, int size)
+{
+  EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
+  return innerVectors(start, size);
+}
+
+/** \returns the i-th row of the matrix \c *this. For row-major matrix only.
+  * (read-only version) */
+template<typename Derived>
+const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::subrows(int start, int size) const
+{
+  EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
+  return innerVectors(start, size);
+}
+
+/** \returns the i-th column of the matrix \c *this. For column-major matrix only. */
+template<typename Derived>
+SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::subcols(int start, int size)
+{
+  EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
+  return innerVectors(start, size);
+}
+
+/** \returns the i-th column of the matrix \c *this. For column-major matrix only.
+  * (read-only version) */
+template<typename Derived>
+const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::subcols(int start, int size) const
+{
+  EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
+  return innerVectors(start, size);
+}
+
+/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this
+  * is col-major (resp. row-major).
+  */
+template<typename Derived>
+SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::innerVectors(int outerStart, int outerSize)
+{ return SparseInnerVectorSet<Derived,Dynamic>(derived(), outerStart, outerSize); }
+
+/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this
+  * is col-major (resp. row-major). Read-only.
+  */
+template<typename Derived>
+const SparseInnerVectorSet<Derived,Dynamic> SparseMatrixBase<Derived>::innerVectors(int outerStart, int outerSize) const
+{ return SparseInnerVectorSet<Derived,Dynamic>(derived(), outerStart, outerSize); }
 
 # if 0
 template<typename MatrixType, int BlockRows, int BlockCols, int PacketAccess>

Eigen/src/Sparse/SparseCwiseBinaryOp.h

@@ -86,6 +86,8 @@ class SparseCwiseBinaryOp : ei_no_assignment_operator,
     EIGEN_STRONG_INLINE SparseCwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp())
       : m_lhs(lhs), m_rhs(rhs), m_functor(func)
     {
+      EIGEN_STATIC_ASSERT((_LhsNested::Flags&RowMajorBit)==(_RhsNested::Flags&RowMajorBit),
+        BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER)
       EIGEN_STATIC_ASSERT((ei_functor_allows_mixing_real_and_complex<BinaryOp>::ret
                            ? int(ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret)
                            : int(ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret)),
@@ -130,11 +132,10 @@ class SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator
 * Implementation of inner-iterators
 ***************************************************************************/
 
-// template<typename T> struct ei_is_scalar_product { enum { ret = false }; };
-// template<typename T> struct ei_is_scalar_product<ei_scalar_product_op<T> > { enum { ret = true }; };
-
-// helper class
+// template<typename T> struct ei_func_is_conjunction { enum { ret = false }; };
+// template<typename T> struct ei_func_is_conjunction<ei_scalar_product_op<T> > { enum { ret = true }; };
 
+// TODO generalize the ei_scalar_product_op specialization to all conjunctions if any !
 
 // sparse - sparse  (generic)
 template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived>
@@ -259,12 +260,13 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>,
     typedef SparseCwiseBinaryOp<BinaryFunc, Lhs, Rhs> CwiseBinaryXpr;
     typedef typename CwiseBinaryXpr::Scalar Scalar;
     typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
+    typedef typename ei_traits<CwiseBinaryXpr>::RhsNested RhsNested;
     typedef typename _LhsNested::InnerIterator LhsIterator;
     enum { IsRowMajor = (int(Lhs::Flags)&RowMajorBit)==RowMajorBit };
   public:
 
     EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
-      : m_xpr(xpr), m_lhsIter(xpr.lhs(),outer), m_functor(xpr.functor()), m_outer(outer)
+      : m_rhs(xpr.rhs()), m_lhsIter(xpr.lhs(),outer), m_functor(xpr.functor()), m_outer(outer)
     {}
 
     EIGEN_STRONG_INLINE Derived& operator++()
@@ -275,7 +277,7 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>,
 
     EIGEN_STRONG_INLINE Scalar value() const
     { return m_functor(m_lhsIter.value(),
-                       m_xpr.rhs().coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); }
+                       m_rhs.coeff(IsRowMajor?m_outer:m_lhsIter.index(),IsRowMajor?m_lhsIter.index():m_outer)); }
 
     EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); }
     EIGEN_STRONG_INLINE int row() const { return m_lhsIter.row(); }
@@ -284,9 +286,9 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>,
     EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter; }
 
   protected:
-    const CwiseBinaryXpr& m_xpr;
+    const RhsNested m_rhs;
     LhsIterator m_lhsIter;
-    const BinaryFunc& m_functor;
+    const BinaryFunc m_functor;
     const int m_outer;
 };
 
@@ -329,6 +331,10 @@ class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>,
 };
 
 
+/***************************************************************************
+* Implementation of SparseMatrixBase and SparseCwise functions/operators
+***************************************************************************/
+
 template<typename Derived>
 template<typename OtherDerived>
 EIGEN_STRONG_INLINE const SparseCwiseBinaryOp<ei_scalar_difference_op<typename ei_traits<Derived>::Scalar>,

Eigen/src/Sparse/SparseCwiseUnaryOp.h

@@ -89,7 +89,7 @@ class SparseCwiseUnaryOp<UnaryOp,MatrixType>::InnerIterator
 
   protected:
     MatrixTypeIterator m_iter;
-    const UnaryOp& m_functor;
+    const UnaryOp m_functor;
 };
 
 template<typename Derived>

Eigen/src/Sparse/SparseDiagonalProduct.h

@@ -0,0 +1,157 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
+//
+// 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_SPARSE_DIAGONAL_PRODUCT_H
+#define EIGEN_SPARSE_DIAGONAL_PRODUCT_H
+
+// the product a diagonal matrix with a sparse matrix can be easily
+// implemented using expression template. We have two very different cases:
+// 1 - diag * row-major sparse
+//     => each inner vector <=> scalar * sparse vector product
+//     => so we can reuse CwiseUnaryOp::InnerIterator
+// 2 - diag * col-major sparse
+//     => each inner vector <=> densevector * sparse vector cwise product
+//     => again, we can reuse specialization of CwiseBinaryOp::InnerIterator
+//        for that particular case
+// The two other cases are symmetric.
+
+template<typename Lhs, typename Rhs>
+struct ei_traits<SparseDiagonalProduct<Lhs, Rhs> > : ei_traits<SparseProduct<Lhs, Rhs, DiagonalProduct> >
+{
+  typedef typename ei_cleantype<Lhs>::type _Lhs;
+  typedef typename ei_cleantype<Rhs>::type _Rhs;
+  enum {
+    SparseFlags = ((int(_Lhs::Flags)&Diagonal)==Diagonal) ? int(_Rhs::Flags) : int(_Lhs::Flags),
+    Flags = SparseBit | (SparseFlags&RowMajorBit)
+  };
+};
+
+enum {SDP_IsDiagonal, SDP_IsSparseRowMajor, SDP_IsSparseColMajor};
+template<typename Lhs, typename Rhs, typename SparseDiagonalProductType, int RhsMode, int LhsMode> 
+class ei_sparse_diagonal_product_inner_iterator_selector;
+
+template<typename LhsNested, typename RhsNested>
+class SparseDiagonalProduct : public SparseMatrixBase<SparseDiagonalProduct<LhsNested,RhsNested> >, ei_no_assignment_operator
+{
+    typedef typename ei_traits<SparseDiagonalProduct>::_LhsNested _LhsNested;
+    typedef typename ei_traits<SparseDiagonalProduct>::_RhsNested _RhsNested;
+    
+    enum {
+      LhsMode = (_LhsNested::Flags&Diagonal)==Diagonal ? SDP_IsDiagonal
+              : (_LhsNested::Flags&RowMajorBit) ? SDP_IsSparseRowMajor : SDP_IsSparseColMajor,
+      RhsMode = (_RhsNested::Flags&Diagonal)==Diagonal ? SDP_IsDiagonal
+              : (_RhsNested::Flags&RowMajorBit) ? SDP_IsSparseRowMajor : SDP_IsSparseColMajor
+    };
+
+  public:
+
+    EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseDiagonalProduct)
+    
+    typedef ei_sparse_diagonal_product_inner_iterator_selector
+                <_LhsNested,_RhsNested,SparseDiagonalProduct,LhsMode,RhsMode> InnerIterator;
+
+    template<typename Lhs, typename Rhs>
+    EIGEN_STRONG_INLINE SparseDiagonalProduct(const Lhs& lhs, const Rhs& rhs)
+      : m_lhs(lhs), m_rhs(rhs)
+    {
+      ei_assert(lhs.cols() == rhs.rows() && "invalid sparse matrix * diagonal matrix product");
+    }
+
+    EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); }
+    EIGEN_STRONG_INLINE int cols() const { return m_rhs.cols(); }
+
+    EIGEN_STRONG_INLINE const _LhsNested& lhs() const { return m_lhs; }
+    EIGEN_STRONG_INLINE const _RhsNested& rhs() const { return m_rhs; }
+
+  protected:
+    LhsNested m_lhs;
+    RhsNested m_rhs;
+};
+
+
+template<typename Lhs, typename Rhs, typename SparseDiagonalProductType>
+class ei_sparse_diagonal_product_inner_iterator_selector
+<Lhs,Rhs,SparseDiagonalProductType,SDP_IsDiagonal,SDP_IsSparseRowMajor>
+  : public SparseCwiseUnaryOp<ei_scalar_multiple_op<typename Lhs::Scalar>,Rhs>::InnerIterator
+{
+    typedef typename SparseCwiseUnaryOp<ei_scalar_multiple_op<typename Lhs::Scalar>,Rhs>::InnerIterator Base;
+  public:
+    inline ei_sparse_diagonal_product_inner_iterator_selector(
+              const SparseDiagonalProductType& expr, int outer)
+      : Base(expr.rhs()*(expr.lhs().diagonal().coeff(outer)), outer)
+    {}
+};
+
+template<typename Lhs, typename Rhs, typename SparseDiagonalProductType>
+class ei_sparse_diagonal_product_inner_iterator_selector
+<Lhs,Rhs,SparseDiagonalProductType,SDP_IsDiagonal,SDP_IsSparseColMajor>
+  : public SparseCwiseBinaryOp<
+      ei_scalar_product_op<typename Lhs::Scalar>,
+      SparseInnerVectorSet<Rhs,1>,
+      typename Lhs::_CoeffsVectorType>::InnerIterator
+{
+    typedef typename SparseCwiseBinaryOp<
+      ei_scalar_product_op<typename Lhs::Scalar>,
+      SparseInnerVectorSet<Rhs,1>,
+      typename Lhs::_CoeffsVectorType>::InnerIterator Base;
+  public:
+    inline ei_sparse_diagonal_product_inner_iterator_selector(
+              const SparseDiagonalProductType& expr, int outer)
+      : Base(expr.rhs().innerVector(outer) .cwise()* expr.lhs().diagonal(), 0)
+    {}
+};
+
+template<typename Lhs, typename Rhs, typename SparseDiagonalProductType>
+class ei_sparse_diagonal_product_inner_iterator_selector
+<Lhs,Rhs,SparseDiagonalProductType,SDP_IsSparseColMajor,SDP_IsDiagonal>
+  : public SparseCwiseUnaryOp<ei_scalar_multiple_op<typename Rhs::Scalar>,Lhs>::InnerIterator
+{
+    typedef typename SparseCwiseUnaryOp<ei_scalar_multiple_op<typename Rhs::Scalar>,Lhs>::InnerIterator Base;
+  public:
+    inline ei_sparse_diagonal_product_inner_iterator_selector(
+              const SparseDiagonalProductType& expr, int outer)
+      : Base(expr.lhs()*expr.rhs().diagonal().coeff(outer), outer)
+    {}
+};
+
+template<typename Lhs, typename Rhs, typename SparseDiagonalProductType>
+class ei_sparse_diagonal_product_inner_iterator_selector
+<Lhs,Rhs,SparseDiagonalProductType,SDP_IsSparseRowMajor,SDP_IsDiagonal>
+  : public SparseCwiseBinaryOp<
+      ei_scalar_product_op<typename Rhs::Scalar>,
+      SparseInnerVectorSet<Lhs,1>,
+      NestByValue<Transpose<typename Rhs::_CoeffsVectorType> > >::InnerIterator
+{
+    typedef typename SparseCwiseBinaryOp<
+      ei_scalar_product_op<typename Rhs::Scalar>,
+      SparseInnerVectorSet<Lhs,1>,
+      NestByValue<Transpose<typename Rhs::_CoeffsVectorType> > >::InnerIterator Base;
+  public:
+    inline ei_sparse_diagonal_product_inner_iterator_selector(
+              const SparseDiagonalProductType& expr, int outer)
+      : Base(expr.lhs().innerVector(outer) .cwise()* expr.rhs().diagonal().transpose().nestByValue(), 0)
+    {}
+};
+
+#endif // EIGEN_SPARSE_DIAGONAL_PRODUCT_H

Eigen/src/Sparse/SparseMatrix.h

@@ -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) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -62,7 +62,7 @@ class SparseMatrix
     // FIXME: why are these operator already alvailable ???
     // EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(SparseMatrix, *=)
     // EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(SparseMatrix, /=)
-    
+
     typedef MappedSparseMatrix<Scalar,Flags> Map;
 
   protected:
@@ -79,7 +79,7 @@ class SparseMatrix
 
     inline int rows() const { return IsRowMajor ? m_outerSize : m_innerSize; }
     inline int cols() const { return IsRowMajor ? m_innerSize : m_outerSize; }
-    
+
     inline int innerSize() const { return m_innerSize; }
     inline int outerSize() const { return m_outerSize; }
     inline int innerNonZeros(int j) const { return m_outerIndex[j+1]-m_outerIndex[j]; }
@@ -138,7 +138,6 @@ class SparseMatrix
       */
     inline void startFill(int reserveSize = 1000)
     {
-//       std::cerr << this << " startFill\n";
       setZero();
       m_data.reserve(reserveSize);
     }
@@ -161,6 +160,10 @@ class SparseMatrix
         }
         m_outerIndex[outer+1] = m_outerIndex[outer];
       }
+      else
+      {
+        ei_assert(m_data.index(m_data.size()-1)<inner && "wrong sorted insertion");
+      }
       assert(size_t(m_outerIndex[outer+1]) == m_data.size());
       int id = m_outerIndex[outer+1];
       ++m_outerIndex[outer+1];
@@ -192,7 +195,7 @@ class SparseMatrix
       // FIXME let's make sure sizeof(long int) == sizeof(size_t)
       size_t id = m_outerIndex[outer+1];
       ++m_outerIndex[outer+1];
-      
+
       float reallocRatio = 1;
       if (m_data.allocatedSize()<id+1)
       {
@@ -214,7 +217,7 @@ class SparseMatrix
         m_data.value(id) = m_data.value(id-1);
         --id;
       }
-      
+
       m_data.index(id) = inner;
       return (m_data.value(id) = 0);
     }
@@ -233,7 +236,7 @@ class SparseMatrix
         ++i;
       }
     }
-    
+
     void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>())
     {
       int k = 0;
@@ -390,11 +393,11 @@ class SparseMatrix
         s << std::endl;
         s << std::endl;
         s << "Column pointers:\n";
-        for (int i=0; i<m.cols(); ++i)
+        for (int i=0; i<m.outerSize(); ++i)
         {
           s << m.m_outerIndex[i] << " ";
         }
-        s << std::endl;
+        s << " $" << std::endl;
         s << std::endl;
       );
       s << static_cast<const SparseMatrixBase<SparseMatrix>&>(m);

Eigen/src/Sparse/SparseMatrixBase.h

@@ -327,18 +327,21 @@ template<typename Derived> class SparseMatrixBase
     // void transposeInPlace();
     const AdjointReturnType adjoint() const { return conjugate()/*.nestByValue()*/; }
 
-    SparseInnerVector<Derived> row(int i);
-    const SparseInnerVector<Derived> row(int i) const;
-    SparseInnerVector<Derived> col(int j);
-    const SparseInnerVector<Derived> col(int j) const;
-    SparseInnerVector<Derived> innerVector(int outer);
-    const SparseInnerVector<Derived> innerVector(int outer) const;
-
-//     RowXpr row(int i);
-//     const RowXpr row(int i) const;
-
-//     ColXpr col(int i);
-//     const ColXpr col(int i) const;
+    // sub-vector
+    SparseInnerVectorSet<Derived,1> row(int i);
+    const SparseInnerVectorSet<Derived,1> row(int i) const;
+    SparseInnerVectorSet<Derived,1> col(int j);
+    const SparseInnerVectorSet<Derived,1> col(int j) const;
+    SparseInnerVectorSet<Derived,1> innerVector(int outer);
+    const SparseInnerVectorSet<Derived,1> innerVector(int outer) const;
+    
+    // set of sub-vectors
+    SparseInnerVectorSet<Derived,Dynamic> subrows(int start, int size);
+    const SparseInnerVectorSet<Derived,Dynamic> subrows(int start, int size) const;
+    SparseInnerVectorSet<Derived,Dynamic> subcols(int start, int size);
+    const SparseInnerVectorSet<Derived,Dynamic> subcols(int start, int size) const;
+    SparseInnerVectorSet<Derived,Dynamic> innerVectors(int outerStart, int outerSize);
+    const SparseInnerVectorSet<Derived,Dynamic> innerVectors(int outerStart, int outerSize) const;
 
 //     typename BlockReturnType<Derived>::Type block(int startRow, int startCol, int blockRows, int blockCols);
 //     const typename BlockReturnType<Derived>::Type

Eigen/src/Sparse/SparseProduct.h

@@ -29,7 +29,9 @@ template<typename Lhs, typename Rhs> struct ei_sparse_product_mode
 {
   enum {
 
-    value = (Rhs::Flags&Lhs::Flags&SparseBit)==SparseBit
+    value = ((Lhs::Flags&Diagonal)==Diagonal || (Rhs::Flags&Diagonal)==Diagonal)
+          ? DiagonalProduct
+          :  (Rhs::Flags&Lhs::Flags&SparseBit)==SparseBit
           ? SparseTimeSparseProduct
           : (Lhs::Flags&SparseBit)==SparseBit
           ? SparseTimeDenseProduct
@@ -45,6 +47,15 @@ struct SparseProductReturnType
   typedef SparseProduct<LhsNested, RhsNested, ProductMode> Type;
 };
 
+template<typename Lhs, typename Rhs>
+struct SparseProductReturnType<Lhs,Rhs,DiagonalProduct>
+{
+  typedef const typename ei_nested<Lhs,Rhs::RowsAtCompileTime>::type LhsNested;
+  typedef const typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type RhsNested;
+
+  typedef SparseDiagonalProduct<LhsNested, RhsNested> Type;
+};
+
 // sparse product return type specialization
 template<typename Lhs, typename Rhs>
 struct SparseProductReturnType<Lhs,Rhs,SparseTimeSparseProduct>
@@ -95,7 +106,7 @@ struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >
 //     RhsIsRowMajor = (RhsFlags & RowMajorBit)==RowMajorBit,
 
     EvalToRowMajor = (RhsFlags & LhsFlags & RowMajorBit),
-    ResultIsSparse = ProductMode==SparseTimeSparseProduct,
+    ResultIsSparse = ProductMode==SparseTimeSparseProduct || ProductMode==DiagonalProduct,
 
     RemovedBits = ~( (EvalToRowMajor ? 0 : RowMajorBit) | (ResultIsSparse ? 0 : SparseBit) ),
 
@@ -105,14 +116,15 @@ struct ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >
 
     CoeffReadCost = Dynamic
   };
-  
+
   typedef typename ei_meta_if<ResultIsSparse,
     SparseMatrixBase<SparseProduct<LhsNested, RhsNested, ProductMode> >,
     MatrixBase<SparseProduct<LhsNested, RhsNested, ProductMode> > >::ret Base;
 };
 
 template<typename LhsNested, typename RhsNested, int ProductMode>
-class SparseProduct : ei_no_assignment_operator, public ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >::Base
+class SparseProduct : ei_no_assignment_operator,
+  public ei_traits<SparseProduct<LhsNested, RhsNested, ProductMode> >::Base
 {
   public:
 
@@ -130,7 +142,7 @@ class SparseProduct : ei_no_assignment_operator, public ei_traits<SparseProduct<
       : m_lhs(lhs), m_rhs(rhs)
     {
       ei_assert(lhs.cols() == rhs.rows());
-      
+
       enum {
         ProductIsValid = _LhsNested::ColsAtCompileTime==Dynamic
                       || _RhsNested::RowsAtCompileTime==Dynamic
@@ -159,6 +171,55 @@ class SparseProduct : ei_no_assignment_operator, public ei_traits<SparseProduct<
     RhsNested m_rhs;
 };
 
+// perform a pseudo in-place sparse * sparse product assuming all matrices are col major
+template<typename Lhs, typename Rhs, typename ResultType>
+static void ei_sparse_product_impl(const Lhs& lhs, const Rhs& rhs, ResultType& res)
+{
+  typedef typename ei_traits<typename ei_cleantype<Lhs>::type>::Scalar Scalar;
+
+  // make sure to call innerSize/outerSize since we fake the storage order.
+  int rows = lhs.innerSize();
+  int cols = rhs.outerSize();
+  //int size = lhs.outerSize();
+  ei_assert(lhs.outerSize() == rhs.innerSize());
+
+  // allocate a temporary buffer
+  AmbiVector<Scalar> tempVector(rows);
+
+  // estimate the number of non zero entries
+  float ratioLhs = float(lhs.nonZeros())/(float(lhs.rows())*float(lhs.cols()));
+  float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
+  float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f);
+
+  res.resize(rows, cols);
+  res.startFill(int(ratioRes*rows*cols));
+  for (int j=0; j<cols; ++j)
+  {
+    // let's do a more accurate determination of the nnz ratio for the current column j of res
+    //float ratioColRes = std::min(ratioLhs * rhs.innerNonZeros(j), 1.f);
+    // FIXME find a nice way to get the number of nonzeros of a sub matrix (here an inner vector)
+    float ratioColRes = ratioRes;
+    tempVector.init(ratioColRes);
+    tempVector.setZero();
+    for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
+    {
+      // FIXME should be written like this: tmp += rhsIt.value() * lhs.col(rhsIt.index())
+      tempVector.restart();
+      Scalar x = rhsIt.value();
+      for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
+      {
+        tempVector.coeffRef(lhsIt.index()) += lhsIt.value() * x;
+      }
+    }
+    for (typename AmbiVector<Scalar>::Iterator it(tempVector); it; ++it)
+      if (ResultType::Flags&RowMajorBit)
+        res.fill(j,it.index()) = it.value();
+      else
+        res.fill(it.index(), j) = it.value();
+  }
+  res.endFill();
+}
+
 template<typename Lhs, typename Rhs, typename ResultType,
   int LhsStorageOrder = ei_traits<Lhs>::Flags&RowMajorBit,
   int RhsStorageOrder = ei_traits<Rhs>::Flags&RowMajorBit,
@@ -172,58 +233,21 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,ColMajor>
 
   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
   {
-    // make sure to call innerSize/outerSize since we fake the storage order.
-    int rows = lhs.innerSize();
-    int cols = rhs.outerSize();
-    //int size = lhs.outerSize();
-    ei_assert(lhs.outerSize() == rhs.innerSize());
-
-    // allocate a temporary buffer
-    AmbiVector<Scalar> tempVector(rows);
-
-    // estimate the number of non zero entries
-    float ratioLhs = float(lhs.nonZeros())/float(lhs.rows()*lhs.cols());
-    float avgNnzPerRhsColumn = float(rhs.nonZeros())/float(cols);
-    float ratioRes = std::min(ratioLhs * avgNnzPerRhsColumn, 1.f);
-
-    res.resize(rows, cols);
-    res.startFill(int(ratioRes*rows*cols));
-    for (int j=0; j<cols; ++j)
-    {
-      // let's do a more accurate determination of the nnz ratio for the current column j of res
-      //float ratioColRes = std::min(ratioLhs * rhs.innerNonZeros(j), 1.f);
-      // FIXME find a nice way to get the number of nonzeros of a sub matrix (here an inner vector)
-      float ratioColRes = ratioRes;
-      tempVector.init(ratioColRes);
-      tempVector.setZero();
-      for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
-      {
-        // FIXME should be written like this: tmp += rhsIt.value() * lhs.col(rhsIt.index())
-        tempVector.restart();
-        Scalar x = rhsIt.value();
-        for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
-        {
-          tempVector.coeffRef(lhsIt.index()) += lhsIt.value() * x;
-        }
-      }
-      for (typename AmbiVector<Scalar>::Iterator it(tempVector); it; ++it)
-        if (ResultType::Flags&RowMajorBit)
-          res.fill(j,it.index()) = it.value();
-        else
-          res.fill(it.index(), j) = it.value();
-    }
-    res.endFill();
+    typename ei_cleantype<ResultType>::type _res(res.rows(), res.cols());
+    ei_sparse_product_impl<Lhs,Rhs,ResultType>(lhs, rhs, _res);
+    res.swap(_res);
   }
 };
 
 template<typename Lhs, typename Rhs, typename ResultType>
 struct ei_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,ColMajor,RowMajor>
 {
-  typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
   {
+    // we need a col-major matrix to hold the result
+    typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
     SparseTemporaryType _res(res.rows(), res.cols());
-    ei_sparse_product_selector<Lhs,Rhs,SparseTemporaryType,ColMajor,ColMajor,ColMajor>::run(lhs, rhs, _res);
+    ei_sparse_product_impl<Lhs,Rhs,SparseTemporaryType>(lhs, rhs, _res);
     res = _res;
   }
 };
@@ -234,20 +258,21 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,RowMajor>
   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
   {
     // let's transpose the product to get a column x column product
-    ei_sparse_product_selector<Rhs,Lhs,ResultType,ColMajor,ColMajor,ColMajor>::run(rhs, lhs, res);
+    typename ei_cleantype<ResultType>::type _res(res.rows(), res.cols());
+    ei_sparse_product_impl<Rhs,Lhs,ResultType>(rhs, lhs, _res);
+    res.swap(_res);
   }
 };
 
 template<typename Lhs, typename Rhs, typename ResultType>
 struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
 {
-  typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
   static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res)
   {
     // let's transpose the product to get a column x column product
+    typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType;
     SparseTemporaryType _res(res.cols(), res.rows());
-    ei_sparse_product_selector<Rhs,Lhs,SparseTemporaryType,ColMajor,ColMajor,ColMajor>
-      ::run(rhs, lhs, _res);
+    ei_sparse_product_impl<Rhs,Lhs,SparseTemporaryType>(rhs, lhs, _res);
     res = _res.transpose();
   }
 };
@@ -285,7 +310,6 @@ template<typename Derived>
 template<typename Lhs, typename Rhs>
 inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs,SparseTimeSparseProduct>& product)
 {
-//   std::cout << "sparse product to sparse\n";
   ei_sparse_product_selector<
     typename ei_cleantype<Lhs>::type,
     typename ei_cleantype<Rhs>::type,
@@ -333,7 +357,7 @@ Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeD
       derived().row(j) += i.value() * product.rhs().row(j);
       ++i;
     }
-    Block<Derived,1,Derived::ColsAtCompileTime> foo = derived().row(j);
+    Block<Derived,1,Derived::ColsAtCompileTime> res(derived().row(LhsIsRowMajor ? j : 0));
     for (; (ProcessFirstHalf ? i && i.index() < j : i) ; ++i)
     {
       if (LhsIsSelfAdjoint)
@@ -345,7 +369,7 @@ Derived& MatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs,SparseTimeD
         derived().row(b) += ei_conj(v) * product.rhs().row(a);
       }
       else if (LhsIsRowMajor)
-        foo += i.value() * product.rhs().row(i.index());
+        res += i.value() * product.rhs().row(i.index());
       else
         derived().row(i.index()) += i.value() * product.rhs().row(j);
     }

Eigen/src/Sparse/SparseUtil.h

@@ -107,12 +107,13 @@ template<typename _Scalar, int _Flags = 0> class SparseVector;
 template<typename _Scalar, int _Flags = 0> class MappedSparseMatrix;
 
 template<typename MatrixType>                            class SparseTranspose;
-template<typename MatrixType>                            class SparseInnerVector;
+template<typename MatrixType, int Size>                  class SparseInnerVectorSet;
 template<typename Derived>                               class SparseCwise;
 template<typename UnaryOp,   typename MatrixType>        class SparseCwiseUnaryOp;
 template<typename BinaryOp,  typename Lhs, typename Rhs> class SparseCwiseBinaryOp;
 template<typename ExpressionType,
          unsigned int Added, unsigned int Removed>       class SparseFlagged;
+template<typename Lhs, typename Rhs>                     class SparseDiagonalProduct;
 
 template<typename Lhs, typename Rhs> struct ei_sparse_product_mode;
 template<typename Lhs, typename Rhs, int ProductMode = ei_sparse_product_mode<Lhs,Rhs>::value> struct SparseProductReturnType;

Eigen/src/Sparse/SparseVector.h

@@ -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) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -59,6 +59,7 @@ class SparseVector
     EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseVector)
     EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseVector, +=)
     EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseVector, -=)
+//     EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(SparseVector, =)
 
   protected:
   public:
@@ -68,6 +69,9 @@ class SparseVector
 
     CompressedStorage<Scalar> m_data;
     int m_size;
+    
+    CompressedStorage<Scalar>& _data() { return m_data; }
+    CompressedStorage<Scalar>& _data() const { return m_data; }
 
   public:
 
@@ -198,6 +202,13 @@ class SparseVector
     {
       *this = other.derived();
     }
+    
+    template<typename OtherDerived>
+    inline SparseVector(const SparseMatrixBase<OtherDerived>& other)
+      : m_size(0)
+    {
+      *this = other.derived();
+    }
 
     inline SparseVector(const SparseVector& other)
       : m_size(0)
@@ -225,9 +236,12 @@ class SparseVector
       return *this;
     }
 
-//     template<typename OtherDerived>
-//     inline SparseVector& operator=(const MatrixBase<OtherDerived>& other)
-//     {
+    template<typename OtherDerived>
+    inline SparseVector& operator=(const SparseMatrixBase<OtherDerived>& other)
+    {
+      return Base::operator=(other);
+    }
+    
 //       const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
 //       if (needToTranspose)
 //       {

Eigen/src/Sparse/SuperLUSupport.h

@@ -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) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -35,7 +35,8 @@
          FLOATTYPE *recip_pivot_growth,                      \
          FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr, \
          SuperLUStat_t *stats, int *info, KEYTYPE) {         \
-    NAMESPACE::mem_usage_t mem_usage;                                    \
+    using namespace NAMESPACE; \
+    mem_usage_t mem_usage;                                    \
     NAMESPACE::FNAME(options, A, perm_c, perm_r, etree, equed, R, C, L,  \
          U, work, lwork, B, X, recip_pivot_growth, rcond,    \
          ferr, berr, &mem_usage, stats, info);               \
@@ -59,7 +60,10 @@ struct SluMatrixMapHelper;
   */
 struct SluMatrix : SuperMatrix
 {
-  SluMatrix() {}
+  SluMatrix()
+  {
+    Store = &storage;
+  }
 
   SluMatrix(const SluMatrix& other)
     : SuperMatrix(other)
@@ -67,6 +71,14 @@ struct SluMatrix : SuperMatrix
     Store = &storage;
     storage = other.storage;
   }
+  
+  SluMatrix& operator=(const SluMatrix& other)
+  {
+    SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
+    Store = &storage;
+    storage = other.storage;
+    return *this;
+  }
 
   struct
   {
@@ -104,7 +116,7 @@ struct SluMatrix : SuperMatrix
       ei_assert(false && "Scalar type not supported by SuperLU");
     }
   }
-  
+
   template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
   static SluMatrix Map(Matrix<Scalar,Rows,Cols,Options,MRows,MCols>& mat)
   {
@@ -223,6 +235,7 @@ SluMatrix SparseMatrixBase<Derived>::asSluMatrix()
   return SluMatrix::Map(derived());
 }
 
+/** View a Super LU matrix as an Eigen expression */
 template<typename Scalar, int Flags>
 MappedSparseMatrix<Scalar,Flags>::MappedSparseMatrix(SluMatrix& sluMat)
 {

Eigen/src/Sparse/TaucsSupport.h

@@ -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) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -32,9 +32,9 @@ taucs_ccs_matrix SparseMatrixBase<Derived>::asTaucsMatrix()
   res.n         = cols();
   res.m         = rows();
   res.flags     = 0;
-  res.colptr    = _outerIndexPtr();
-  res.rowind    = _innerIndexPtr();
-  res.values.v  = _valuePtr();
+  res.colptr    = derived()._outerIndexPtr();
+  res.rowind    = derived()._innerIndexPtr();
+  res.values.v  = derived()._valuePtr();
   if (ei_is_same_type<Scalar,int>::ret)
     res.flags |= TAUCS_INT;
   else if (ei_is_same_type<Scalar,float>::ret)
@@ -78,8 +78,8 @@ class SparseLLT<MatrixType,Taucs> : public SparseLLT<MatrixType>
 {
   protected:
     typedef SparseLLT<MatrixType> Base;
-    using Base::Scalar;
-    using Base::RealScalar;
+    typedef typename Base::Scalar Scalar;
+    typedef typename Base::RealScalar RealScalar;
     using Base::MatrixLIsDirty;
     using Base::SupernodalFactorIsDirty;
     using Base::m_flags;
@@ -129,7 +129,10 @@ void SparseLLT<MatrixType,Taucs>::compute(const MatrixType& a)
   {
     taucs_ccs_matrix taucsMatA = const_cast<MatrixType&>(a).asTaucsMatrix();
     taucs_ccs_matrix* taucsRes = taucs_ccs_factor_llt(&taucsMatA, Base::m_precision, 0);
-    m_matrix = Base::CholMatrixType::Map(*taucsRes);
+    // the matrix returned by Taucs is not necessarily sorted,
+    // so let's copy it in two steps
+    DynamicSparseMatrix<Scalar,RowMajor> tmp = MappedSparseMatrix<Scalar>(*taucsRes);
+    m_matrix = tmp;
     free(taucsRes);
     m_status = (m_status & ~(CompleteFactorization|MatrixLIsDirty))
              | IncompleteFactorization
@@ -161,7 +164,11 @@ SparseLLT<MatrixType,Taucs>::matrixL() const
     ei_assert(!(m_status & SupernodalFactorIsDirty));
 
     taucs_ccs_matrix* taucsL = taucs_supernodal_factor_to_ccs(m_taucsSupernodalFactor);
-    const_cast<typename Base::CholMatrixType&>(m_matrix) = Base::CholMatrixType::Map(*taucsL);
+
+    // the matrix returned by Taucs is not necessarily sorted,
+    // so let's copy it in two steps
+    DynamicSparseMatrix<Scalar,RowMajor> tmp = MappedSparseMatrix<Scalar>(*taucsL);
+    const_cast<typename Base::CholMatrixType&>(m_matrix) = tmp;
     free(taucsL);
     m_status = (m_status & ~MatrixLIsDirty);
   }
@@ -172,22 +179,32 @@ template<typename MatrixType>
 template<typename Derived>
 void SparseLLT<MatrixType,Taucs>::solveInPlace(MatrixBase<Derived> &b) const
 {
-  if (m_status & MatrixLIsDirty)
+  bool inputIsCompatibleWithTaucs = (Derived::Flags&RowMajorBit)==0;
+
+  if (!inputIsCompatibleWithTaucs)
   {
-    // TODO use taucs's supernodal solver, in particular check types, storage order, etc.
-    // VectorXb x(b.rows());
-    // for (int j=0; j<b.cols(); ++j)
-    // {
-    //   taucs_supernodal_solve_llt(m_taucsSupernodalFactor,x.data(),&b.col(j).coeffRef(0));
-    //   b.col(j) = x;
-    // }
     matrixL();
-  }
-
-
-  {
     Base::solveInPlace(b);
   }
+  else if (m_flags & IncompleteFactorization)
+  {
+    taucs_ccs_matrix taucsLLT = const_cast<typename Base::CholMatrixType&>(m_matrix).asTaucsMatrix();
+    typename ei_plain_matrix_type<Derived>::type x(b.rows());
+    for (int j=0; j<b.cols(); ++j)
+    {
+      taucs_ccs_solve_llt(&taucsLLT,x.data(),&b.col(j).coeffRef(0));
+      b.col(j) = x;
+    }
+  }
+  else
+  {
+    typename ei_plain_matrix_type<Derived>::type x(b.rows());
+    for (int j=0; j<b.cols(); ++j)
+    {
+      taucs_supernodal_solve_llt(m_taucsSupernodalFactor,x.data(),&b.col(j).coeffRef(0));
+      b.col(j) = x;
+    }
+  }
 }
 
 #endif // EIGEN_TAUCSSUPPORT_H

Eigen/src/Sparse/UmfPackSupport.h

@@ -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) 2008 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public

Eigen/src/StdVector/CMakeLists.txt

@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_StdVector_SRCS "*.h")
-
-INSTALL(FILES
-  ${Eigen_StdVector_SRCS}
-  DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/StdVector
-  )

Eigen/src/StdVector/StdVector.h

@@ -1,73 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra. Eigen itself is part of the KDE project.
-//
-// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
-// Copyright (C) 2009 Alex Stapleton <alex.stapleton@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/>.
-
-#ifndef EIGEN_STDVECTOR_H
-#define EIGEN_STDVECTOR_H
-
-#define EIGEN_STD_VECTOR_SPECIALIZATION_BODY \
-  typedef Eigen::aligned_allocator<value_type> allocator_type; \
-  typedef vector<value_type, allocator_type > unaligned_base; \
-  typedef typename unaligned_base::size_type size_type; \
-  typedef typename unaligned_base::iterator iterator; \
-  explicit vector(const allocator_type& __a = allocator_type()) : unaligned_base(__a) {} \
-  vector(const vector& c) : unaligned_base(c) {} \
-  vector(size_type num, const value_type& val = value_type()) : unaligned_base(num, val) {}\
-  vector(iterator start, iterator end) : unaligned_base(start, end) {} \
-  vector& operator=(const vector& __x) { \
-    unaligned_base::operator=(__x); \
-    return *this; \
-  }
-
-template <typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols, typename _Alloc>
-class vector<Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>, _Alloc>
-  : public vector<Eigen::ei_unaligned_type<Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> >,
-                  Eigen::aligned_allocator<Eigen::ei_unaligned_type<Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> > > >
-{
-public:
-  typedef Eigen::ei_unaligned_type<Eigen::Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> > value_type;
-  EIGEN_STD_VECTOR_SPECIALIZATION_BODY
-};
-
-template <typename _Scalar, int _Dim, typename _Alloc>
-class vector<Eigen::Transform<_Scalar,_Dim>, _Alloc>
-  : public vector<Eigen::ei_unaligned_type<Eigen::Transform<_Scalar,_Dim> >,
-                  Eigen::aligned_allocator<Eigen::ei_unaligned_type<Eigen::Transform<_Scalar,_Dim> > > >
-{
-public:
-  typedef Eigen::ei_unaligned_type<Eigen::Transform<_Scalar,_Dim> > value_type;
-  EIGEN_STD_VECTOR_SPECIALIZATION_BODY
-};
-
-template <typename _Scalar, typename _Alloc>
-class vector<Eigen::Quaternion<_Scalar>, _Alloc>
-  : public vector<Eigen::ei_unaligned_type<Eigen::Quaternion<_Scalar> >,
-                  Eigen::aligned_allocator<Eigen::ei_unaligned_type<Eigen::Quaternion<_Scalar> > > >
-{
-public:
-  typedef Eigen::ei_unaligned_type<Eigen::Quaternion<_Scalar> > value_type;
-  EIGEN_STD_VECTOR_SPECIALIZATION_BODY
-};
-
-#endif // EIGEN_STDVECTOR_H

Eigen/src/StdVector/UnalignedType.h

@@ -1,105 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra. Eigen itself is part of the KDE project.
-//
-// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
-// Copyright (C) 2009 Alex Stapleton <alex.stapleton@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/>.
-
-#ifndef EIGEN_UNALIGNEDTYPE_H
-#define EIGEN_UNALIGNEDTYPE_H
-
-template<typename aligned_type> class ei_unaligned_type;
-
-template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
-class ei_unaligned_type<Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> >
-  : public Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols>
-{
-  private:
-    template<typename Other> void _unaligned_copy(const Other& other)
-    {
-      if(other.size() == 0) return;
-      resize(other.rows(), other.cols());
-      ei_assign_impl<ei_unaligned_type,aligned_base,NoVectorization>::run(*this, other);
-    }
-
-  public:
-    typedef Matrix<_Scalar,_Rows,_Cols,_Options,_MaxRows,_MaxCols> aligned_base;
-    ei_unaligned_type() : aligned_base(ei_constructor_without_unaligned_array_assert()) {}
-    ei_unaligned_type(const aligned_base& other) : aligned_base(ei_constructor_without_unaligned_array_assert())
-    {
-      _unaligned_copy(other);
-    }
-    ei_unaligned_type(const ei_unaligned_type& other) : aligned_base(ei_constructor_without_unaligned_array_assert())
-    {
-      _unaligned_copy(other);
-    }
-};
-
-template<typename _Scalar, int _Dim>
-class ei_unaligned_type<Transform<_Scalar,_Dim> >
-  : public Transform<_Scalar,_Dim>
-{
-  private:
-    template<typename Other> void _unaligned_copy(const Other& other)
-    {
-      // no resizing here, it's fixed-size anyway
-      ei_assign_impl<MatrixType,MatrixType,NoVectorization>::run(this->matrix(), other.matrix());
-    }
-  public:
-    typedef Transform<_Scalar,_Dim> aligned_base;
-    typedef typename aligned_base::MatrixType MatrixType;
-    ei_unaligned_type() : aligned_base(ei_constructor_without_unaligned_array_assert()) {}
-    ei_unaligned_type(const aligned_base& other) : aligned_base(ei_constructor_without_unaligned_array_assert())
-    {
-      _unaligned_copy(other);
-    }
-    ei_unaligned_type(const ei_unaligned_type& other) : aligned_base(ei_constructor_without_unaligned_array_assert())
-    {
-      _unaligned_copy(other);
-    }
-};
-
-template<typename _Scalar>
-class ei_unaligned_type<Quaternion<_Scalar> >
-  : public Quaternion<_Scalar>
-{
-  private:
-    template<typename Other> void _unaligned_copy(const Other& other)
-    {
-      // no resizing here, it's fixed-size anyway
-      ei_assign_impl<Coefficients,Coefficients,NoVectorization>::run(this->coeffs(), other.coeffs());
-    }
-  public:
-    typedef Quaternion<_Scalar> aligned_base;
-    typedef typename aligned_base::Coefficients Coefficients;
-    ei_unaligned_type() : aligned_base(ei_constructor_without_unaligned_array_assert()) {}
-    ei_unaligned_type(const aligned_base& other) : aligned_base(ei_constructor_without_unaligned_array_assert())
-    {
-      _unaligned_copy(other);
-    }
-    ei_unaligned_type(const ei_unaligned_type& other) : aligned_base(ei_constructor_without_unaligned_array_assert())
-    {
-      _unaligned_copy(other);
-    }
-};
-
-
-#endif // EIGEN_UNALIGNEDTYPE_H

Mainpage.dox

@@ -120,6 +120,6 @@ In order to generate the documentation of Eigen, please follow these steps:
 After doing that, you will find the HTML documentation in the doc/html/ subdirectory of the build directory.
 
 <h2>Note however that the documentation is available online here:
-<a href="http://eigen.tuxfamily.org/api/">http://eigen.tuxfamily.org/api</a></h2>
+<a href="http://eigen.tuxfamily.org/dox/">http://eigen.tuxfamily.org/dox</a></h2>
 
 */

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/D03_WrongStackAlignment.dox

@@ -0,0 +1,53 @@
+namespace Eigen {
+
+/** \page WrongStackAlignment Troubleshooting - Compiler making a wrong assumption on stack alignment
+
+This is an issue that, so far, we met only with GCC on Windows: for instance, MinGW and TDM-GCC.
+
+By default, in a function like this,
+
+\code
+void foo()
+{
+  Eigen::Quaternionf q;
+  //...
+}
+\endcode
+
+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, 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.
+
+
+\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()
+{
+  Eigen::Quaternionf q;
+  //...
+}
+\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. 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
+-mincoming-stack-boundary=2
+\endcode
+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
+-mstackrealign
+\endcode
+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/Experimental.dox

@@ -9,6 +9,12 @@ namespace Eigen {
 
 \section summary Summary
 
+With the 2.0 release, Eigen's API is, to a large extent, stable. However, we wish to retain the freedom to make API incompatible changes. To that effect, we call many parts of Eigen "experimental" which means that they are not subject to API stability guarantee.
+
+Our goal is that for the 2.1 release (expected in July 2009) most of these parts become API-stable too.
+
+We are aware that API stability is a major concern for our users. That's why it's a priority for us to reach it, but at the same time we're being serious about not calling Eigen API-stable too early.
+
 Experimental features may at any time:
 \li be removed;
 \li be subject to an API incompatible change;
@@ -16,11 +22,12 @@ Experimental features may at any time:
 
 \section modules Experimental modules
 
-The following modules are considered entirely experimental:
+The following modules are considered entirely experimental, and we make no firm API stability guarantee about them for the time being:
 \li SVD
 \li QR
+\li Cholesky
 \li Sparse
-\li Geometry
+\li Geometry (this one should be mostly stable, but it's a little too early to make a formal guarantee)
 
 \section core Experimental parts of the Core module
 
@@ -37,7 +44,7 @@ The only classes subject to (even partial) API stability guarantee (meaning that
 
 All other classes offer no direct API guarantee, e.g. their methods can be changed; however notice that most classes inherit MatrixBase and that this is where most of their API comes from -- so in practice most of the API is stable.
 
-Here are the MatrixBase methods that are considered experimental, hence not part of any API stability guarantee:
+A few MatrixBase methods are considered experimental, hence not part of any API stability guarantee:
 \li all methods documented as internal
 \li all methods hidden in the Doxygen documentation
 \li all methods marked as experimental

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

@@ -41,8 +41,6 @@ There is no library to link to. For good performance, add the \c -O2 compile-fla
 
 On the x86 architecture, the SSE2 instruction set is not enabled by default. Use \c -msse2 to enable it, and Eigen will then automatically enable its vectorized paths. On x86-64 and AltiVec-based architectures, vectorization is enabled by default.
 
-<a name="warningarraymodule"></a>
-\warning \redstar In most cases it is enough to include the \c Eigen/Core header only to get started with Eigen. However, some features presented in this tutorial require the Array module to be included (\c \#include \c <Eigen/Array>). Those features are highlighted with a red star \redstar.
 
 \section TutorialCoreSimpleExampleFixedSize Simple example with fixed-size matrices and vectors
 
@@ -69,6 +67,13 @@ output:
 </td></tr></table>
 
 
+<a name="warningarraymodule"></a>
+\warning \redstar In most cases it is enough to include the \c Eigen/Core header only to get started with Eigen. However, some features presented in this tutorial require the Array module to be included (\c \#include \c <Eigen/Array>). Those features are highlighted with a red star \redstar. Notice that if you want to include all Eigen functionality at once, you can do:
+\code
+#include <Eigen/Eigen>
+\endcode
+This slows compilation down but at least you don't have to worry anymore about including the correct files! There also is the Eigen/Dense header including all dense functionality i.e. leaving out the Sparse module.
+
 
 
 
@@ -357,8 +362,8 @@ absolute value (\link Cwise::abs() abs \endlink, \link Cwise::abs2() abs2 \endli
 </td><td>\code
 mat3 = mat1.cwise().min(mat2);
 mat3 = mat1.cwise().max(mat2);
-mat3 = mat1.cwise().abs(mat2);
-mat3 = mat1.cwise().abs2(mat2);
+mat3 = mat1.cwise().abs();
+mat3 = mat1.cwise().abs2();
 \endcode</td></tr>
 </table>
 </td></tr></table>

doc/StlContainers.dox

@@ -6,6 +6,8 @@ namespace Eigen {
   - \ref summary
   - \ref allocator
   - \ref vector
+  - \ref newvector
+
 
 \section summary Executive summary
 
@@ -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

@@ -8,17 +8,31 @@ my_program: path/to/eigen2/Eigen/src/Core/MatrixStorage.h:44:
 Eigen::ei_matrix_array<T, Size, MatrixOptions, Align>::ei_matrix_array()
 [with T = double, int Size = 2, int MatrixOptions = 2, bool Align = true]:
 Assertion `(reinterpret_cast<size_t>(array) & 0xf) == 0 && "this assertion
-is explained here: http://eigen.tuxfamily.org/api/UnalignedArrayAssert.html
+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
 
@@ -64,6 +78,22 @@ then you need to read this separate page: \ref PassingByValue "Passing Eigen obj
 
 Note that here, Eigen::Vector4d is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types".
 
+\section c4 Cause 4: Compiler making a wrong assumption on stack alignment (for instance GCC on Windows)
+
+This is a must-read for people using GCC on Windows (like MinGW or TDM-GCC). If you have this assertion failure in an innocent function declaring a local variable like this:
+
+\code
+void foo()
+{
+  Eigen::Quaternionf q;
+  //...
+}
+\endcode
+
+then you need to read this separate page: \ref WrongStackAlignment "Compiler making a wrong assumption on stack alignment".
+
+Note that here, Eigen::Quaternionf is only used as an example, more generally the issue arises for all \ref FixedSizeVectorizable "fixed-size vectorizable Eigen types".
+
 \section explanation General explanation of this assertion
 
 \ref FixedSizeVectorizable "fixed-size vectorizable Eigen objects" must absolutely be created at 16-byte-aligned locations, otherwise SIMD instructions adressing them will crash.
@@ -72,6 +102,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.</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/snippets/MatrixBase_part.cpp

@@ -1,4 +1,4 @@
-Matrix3d m = Matrix3i::Zero();
+Matrix3d m = Matrix3d::Zero();
 m.part<Eigen::StrictlyUpperTriangular>().setOnes();
 cout << "Here is the matrix m:" << endl << m << endl;
 cout << "And let us now compute m*m.adjoint() in a very optimized way" << endl

scripts/eigen_gen_docs

@@ -0,0 +1,27 @@
+#!/bin/sh
+
+# TODO : actually exit on exit, currently it only exit from the ()
+# TODO : display error msg on stderr instead of stdout
+
+# configuration
+USER='orzel'
+
+# step 1 : update
+hg pull -u || (echo "update failed"; exit 1)
+
+# step 2 : 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 3 : upload
+BRANCH=`hg branch`
+if [ $BRANCH == "default" ]
+then
+    BRANCH='devel'
+fi
+# (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-$BRANCH/ ||  (echo "upload failed"; exit 1)
+
+

test/CMakeLists.txt

@@ -34,7 +34,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)
@@ -211,12 +211,16 @@ 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)
+ei_add_test(swap)
+ei_add_test(visitor)
 
 # print a summary of the different options
 message("************************************************************")

test/basicstuff.cpp

@@ -103,7 +103,10 @@ template<typename MatrixType> void basicStuff(const MatrixType& m)
   m3 = m1;
   m1.swap(m2);
   VERIFY_IS_APPROX(m3, m2);
-  VERIFY_IS_NOT_APPROX(m3, m1);
+  if(rows*cols>=3)
+  {
+    VERIFY_IS_NOT_APPROX(m3, m1);
+  }
 }
 
 void test_basicstuff()

test/cwiseop.cpp

@@ -54,7 +54,9 @@ template<typename MatrixType> void cwiseops(const MatrixType& m)
              square = Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime>::Random(rows, rows);
   VectorType v1 = VectorType::Random(rows),
              v2 = VectorType::Random(rows),
-             vzero = VectorType::Zero(rows);
+             vzero = VectorType::Zero(rows),
+             vones = VectorType::Ones(rows),
+             v3(rows);
 
   int r = ei_random<int>(0, rows-1),
       c = ei_random<int>(0, cols-1);
@@ -70,12 +72,22 @@ template<typename MatrixType> void cwiseops(const MatrixType& m)
       VERIFY_IS_APPROX(mones(i,j), Scalar(1));
       VERIFY_IS_APPROX(m3(i,j), s1);
     }
+  VERIFY(mzero.isZero());
+  VERIFY(mones.isOnes());
+  VERIFY(m3.isConstant(s1));
+  VERIFY(identity.isIdentity());
   VERIFY_IS_APPROX(m4.setConstant(s1), m3);
+  VERIFY_IS_APPROX(m4.setConstant(rows,cols,s1), m3);
   VERIFY_IS_APPROX(m4.setZero(), mzero);
+  VERIFY_IS_APPROX(m4.setZero(rows,cols), mzero);
   VERIFY_IS_APPROX(m4.setOnes(), mones);
+  VERIFY_IS_APPROX(m4.setOnes(rows,cols), mones);
   m4.fill(s1);
   VERIFY_IS_APPROX(m4, m3);
   
+  VERIFY_IS_APPROX(v3.setConstant(rows, s1), VectorType::Constant(rows,s1));
+  VERIFY_IS_APPROX(v3.setZero(rows), vzero);
+  VERIFY_IS_APPROX(v3.setOnes(rows), vones);
 
   m2 = m2.template binaryExpr<AddIfNull<Scalar> >(mones);
 

test/dynalloc.cpp

@@ -24,12 +24,18 @@
 
 #include "main.h"
 
+#if EIGEN_ARCH_WANTS_ALIGNMENT
+#define ALIGNMENT 16
+#else
+#define ALIGNMENT 1
+#endif
+
 void check_handmade_aligned_malloc()
 {
   for(int i = 1; i < 1000; i++)
   {
     char *p = (char*)ei_handmade_aligned_malloc(i);
-    VERIFY(size_t(p)%16==0);
+    VERIFY(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);
@@ -41,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)%16==0);
+    VERIFY(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);
@@ -53,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)%16==0);
+    VERIFY(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);
@@ -65,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)%16==0);
+    VERIFY(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);
@@ -92,7 +98,7 @@ class MyClassA
 template<typename T> void check_dynaligned()
 {
   T* obj = new T;
-  VERIFY(size_t(obj)%16==0);
+  VERIFY(size_t(obj)%ALIGNMENT==0);
   delete obj;
 }
 
@@ -115,15 +121,15 @@ void test_dynalloc()
   
   // check static allocation, who knows ?
   {
-    MyStruct foo0;  VERIFY(size_t(foo0.avec.data())%16==0);
-    MyClassA fooA;  VERIFY(size_t(fooA.avec.data())%16==0);
+    MyStruct foo0;  VERIFY(size_t(foo0.avec.data())%ALIGNMENT==0);
+    MyClassA fooA;  VERIFY(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())%16==0);
-    MyClassA *fooA = new MyClassA();  VERIFY(size_t(fooA->avec.data())%16==0);
+    MyStruct *foo0 = new MyStruct();  VERIFY(size_t(foo0->avec.data())%ALIGNMENT==0);
+    MyClassA *fooA = new MyClassA();  VERIFY(size_t(fooA->avec.data())%ALIGNMENT==0);
     delete foo0;
     delete fooA;
   }
@@ -132,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())%16==0);
-    MyClassA *fooA = new MyClassA[N];  VERIFY(size_t(fooA->avec.data())%16==0);
+    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);
     delete[] foo0;
     delete[] fooA;
   }

test/geometry.cpp

@@ -50,7 +50,7 @@ template<typename Scalar> void geometry(void)
 
   Scalar largeEps = test_precision<Scalar>();
   if (ei_is_same_type<Scalar,float>::ret)
-    largeEps = 1e-3f;
+    largeEps = 1e-2f;
 
   Vector3 v0 = Vector3::Random(),
     v1 = Vector3::Random(),
@@ -96,14 +96,18 @@ template<typename Scalar> void geometry(void)
   Scalar refangle = ei_abs(AngleAxisx(q1.inverse()*q2).angle());
   if (refangle>Scalar(M_PI))
     refangle = Scalar(2)*Scalar(M_PI) - refangle;
-  VERIFY(ei_isApprox(q1.angularDistance(q2), refangle, largeEps));
+  
+  if((q1.coeffs()-q2.coeffs()).norm() > 10*largeEps)
+  {
+    VERIFY(ei_isApprox(q1.angularDistance(q2), refangle, largeEps));
+  }
 
   // rotation matrix conversion
   VERIFY_IS_APPROX(q1 * v2, q1.toRotationMatrix() * v2);
   VERIFY_IS_APPROX(q1 * q2 * v2,
     q1.toRotationMatrix() * q2.toRotationMatrix() * v2);
 
-  VERIFY( (q2*q1).isApprox(q1*q2) || !(q2 * q1 * v2).isApprox(
+  VERIFY( (q2*q1).isApprox(q1*q2, largeEps) || !(q2 * q1 * v2).isApprox(
     q1.toRotationMatrix() * q2.toRotationMatrix() * v2));
 
   q2 = q1.toRotationMatrix();
@@ -147,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();
@@ -157,7 +167,7 @@ template<typename Scalar> void geometry(void)
   t1.prescale(v0);
 
   VERIFY_IS_APPROX( (t0 * Vector3(1,0,0)).norm(), v0.x());
-  VERIFY(!ei_isApprox((t1 * Vector3(1,0,0)).norm(), v0.x()));
+  //VERIFY(!ei_isApprox((t1 * Vector3(1,0,0)).norm(), v0.x()));
 
   t0.setIdentity();
   t1.setIdentity();

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((size_t)&(v[22]) == (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((size_t)&(v[22]) == (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((size_t)&(v[22]) == (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/product_large.cpp

@@ -42,4 +42,10 @@ 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());
+  }
 }

test/regression.cpp

@@ -62,6 +62,20 @@ void makeNoisyCohyperplanarPoints(int numPoints,
     *(points[i]) += noiseAmplitude * VectorType::Random(size);
 }
 
+template<typename VectorType>
+void check_linearRegression(int numPoints,
+                            VectorType **points,
+                            const VectorType& original,
+                            typename VectorType::Scalar tolerance)
+{
+  int size = points[0]->size();
+  assert(size==2);
+  VectorType result(size);
+  linearRegression(numPoints, points, &result, 1);
+  typename VectorType::Scalar error = (result - original).norm() / original.norm();
+  VERIFY(ei_abs(error) < ei_abs(tolerance));
+}
+
 template<typename VectorType,
          typename HyperplaneType>
 void check_fitHyperplane(int numPoints,
@@ -81,6 +95,20 @@ void test_regression()
 {
   for(int i = 0; i < g_repeat; i++)
   {
+    {
+      Vector2f points2f [1000];
+      Vector2f *points2f_ptrs [1000];
+      for(int i = 0; i < 1000; i++) points2f_ptrs[i] = &(points2f[i]);
+      Vector2f coeffs2f;
+      Hyperplane<float,2> coeffs3f;
+      makeNoisyCohyperplanarPoints(1000, points2f_ptrs, &coeffs3f, 0.01f);
+      coeffs2f[0] = -coeffs3f.coeffs()[0]/coeffs3f.coeffs()[1];
+      coeffs2f[1] = -coeffs3f.coeffs()[2]/coeffs3f.coeffs()[1];
+      CALL_SUBTEST(check_linearRegression(10, points2f_ptrs, coeffs2f, 0.05f));
+      CALL_SUBTEST(check_linearRegression(100, points2f_ptrs, coeffs2f, 0.01f));
+      CALL_SUBTEST(check_linearRegression(1000, points2f_ptrs, coeffs2f, 0.002f));
+    }
+    
     {
       Vector2f points2f [1000];
       Vector2f *points2f_ptrs [1000];

test/sparse_basic.cpp

@@ -193,7 +193,6 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
         }
       }
       m2.endFill();
-      //std::cerr << m1 << "\n\n" << m2 << "\n";
       VERIFY_IS_APPROX(m2,m1);
     }
   
@@ -239,6 +238,8 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
     VERIFY_IS_APPROX(m1+=m2, refM1+=refM2);
     VERIFY_IS_APPROX(m1-=m2, refM1-=refM2);
     
+    VERIFY_IS_APPROX(m1.col(0).dot(refM2.row(0)), refM1.col(0).dot(refM2.row(0)));
+    
     refM4.setRandom();
     // sparse cwise* dense
     VERIFY_IS_APPROX(m3.cwise()*refM4, refM3.cwise()*refM4);
@@ -254,6 +255,24 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
     int j1 = ei_random(0,rows-1);
     VERIFY_IS_APPROX(m2.innerVector(j0), refMat2.col(j0));
     VERIFY_IS_APPROX(m2.innerVector(j0)+m2.innerVector(j1), refMat2.col(j0)+refMat2.col(j1));
+    //m2.innerVector(j0) = 2*m2.innerVector(j1);
+    //refMat2.col(j0) = 2*refMat2.col(j1);
+    //VERIFY_IS_APPROX(m2, refMat2);
+  }
+  
+  // test innerVectors()
+  {
+    DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
+    SparseMatrixType m2(rows, rows);
+    initSparse<Scalar>(density, refMat2, m2);
+    int j0 = ei_random(0,rows-2);
+    int j1 = ei_random(0,rows-2);
+    int n0 = ei_random<int>(1,rows-std::max(j0,j1));
+    VERIFY_IS_APPROX(m2.innerVectors(j0,n0), refMat2.block(0,j0,rows,n0));
+    VERIFY_IS_APPROX(m2.innerVectors(j0,n0)+m2.innerVectors(j1,n0),
+                     refMat2.block(0,j0,rows,n0)+refMat2.block(0,j1,rows,n0));
+    //m2.innerVectors(j0,n0) = m2.innerVectors(j0,n0) + m2.innerVectors(j1,n0);
+    //refMat2.block(0,j0,rows,n0) = refMat2.block(0,j0,rows,n0) + refMat2.block(0,j1,rows,n0);
   }
 
   // test transpose
@@ -264,69 +283,6 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
     VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval());
     VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose());
   }
-
-  // test matrix product
-  {
-    DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
-    DenseMatrix refMat3 = DenseMatrix::Zero(rows, rows);
-    DenseMatrix refMat4 = DenseMatrix::Zero(rows, rows);
-    DenseMatrix dm4 = DenseMatrix::Zero(rows, rows);
-    SparseMatrixType m2(rows, rows);
-    SparseMatrixType m3(rows, rows);
-    SparseMatrixType m4(rows, rows);
-    initSparse<Scalar>(density, refMat2, m2);
-    initSparse<Scalar>(density, refMat3, m3);
-    initSparse<Scalar>(density, refMat4, m4);
-    VERIFY_IS_APPROX(m4=m2*m3, refMat4=refMat2*refMat3);
-    VERIFY_IS_APPROX(m4=m2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
-    VERIFY_IS_APPROX(m4=m2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
-    VERIFY_IS_APPROX(m4=m2*m3.transpose(), refMat4=refMat2*refMat3.transpose());
-
-    // sparse * dense
-    VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3);
-    VERIFY_IS_APPROX(dm4=m2*refMat3.transpose(), refMat4=refMat2*refMat3.transpose());
-    VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3, refMat4=refMat2.transpose()*refMat3);
-    VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
-
-    // dense * sparse
-    VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3);
-    VERIFY_IS_APPROX(dm4=refMat2*m3.transpose(), refMat4=refMat2*refMat3.transpose());
-    VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
-    VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
-  }
-  
-  // test self adjoint products
-  {
-    DenseMatrix b = DenseMatrix::Random(rows, rows);
-    DenseMatrix x = DenseMatrix::Random(rows, rows);
-    DenseMatrix refX = DenseMatrix::Random(rows, rows);
-    DenseMatrix refUp = DenseMatrix::Zero(rows, rows);
-    DenseMatrix refLo = DenseMatrix::Zero(rows, rows);
-    DenseMatrix refS = DenseMatrix::Zero(rows, rows);
-    SparseMatrixType mUp(rows, rows);
-    SparseMatrixType mLo(rows, rows);
-    SparseMatrixType mS(rows, rows);
-    do {
-      initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular);
-    } while (refUp.isZero());
-    refLo = refUp.transpose().conjugate();
-    mLo = mUp.transpose().conjugate();
-    refS = refUp + refLo;
-    refS.diagonal() *= 0.5;
-    mS = mUp + mLo;
-    for (int k=0; k<mS.outerSize(); ++k)
-      for (typename SparseMatrixType::InnerIterator it(mS,k); it; ++it)
-        if (it.index() == k)
-          it.valueRef() *= 0.5;
-    
-    VERIFY_IS_APPROX(refS.adjoint(), refS);
-    VERIFY_IS_APPROX(mS.transpose().conjugate(), mS);
-    VERIFY_IS_APPROX(mS, refS);
-    VERIFY_IS_APPROX(x=mS*b, refX=refS*b);
-    VERIFY_IS_APPROX(x=mUp.template marked<UpperTriangular|SelfAdjoint>()*b, refX=refS*b);
-    VERIFY_IS_APPROX(x=mLo.template marked<LowerTriangular|SelfAdjoint>()*b, refX=refS*b);
-    VERIFY_IS_APPROX(x=mS.template marked<SelfAdjoint>()*b, refX=refS*b);
-  }
   
   // test prune
   {

test/sparse_product.cpp

@@ -0,0 +1,130 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@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 "sparse.h"
+
+template<typename SparseMatrixType> void sparse_product(const SparseMatrixType& ref)
+{
+  const int rows = ref.rows();
+  const int cols = ref.cols();
+  typedef typename SparseMatrixType::Scalar Scalar;
+  enum { Flags = SparseMatrixType::Flags };
+
+  double density = std::max(8./(rows*cols), 0.01);
+  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+  typedef Matrix<Scalar,Dynamic,1> DenseVector;
+
+  // test matrix-matrix product
+  {
+    DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
+    DenseMatrix refMat3 = DenseMatrix::Zero(rows, rows);
+    DenseMatrix refMat4 = DenseMatrix::Zero(rows, rows);
+    DenseMatrix dm4 = DenseMatrix::Zero(rows, rows);
+    SparseMatrixType m2(rows, rows);
+    SparseMatrixType m3(rows, rows);
+    SparseMatrixType m4(rows, rows);
+    initSparse<Scalar>(density, refMat2, m2);
+    initSparse<Scalar>(density, refMat3, m3);
+    initSparse<Scalar>(density, refMat4, m4);
+    VERIFY_IS_APPROX(m4=m2*m3, refMat4=refMat2*refMat3);
+    VERIFY_IS_APPROX(m4=m2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
+    VERIFY_IS_APPROX(m4=m2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
+    VERIFY_IS_APPROX(m4=m2*m3.transpose(), refMat4=refMat2*refMat3.transpose());
+
+    // sparse * dense
+    VERIFY_IS_APPROX(dm4=m2*refMat3, refMat4=refMat2*refMat3);
+    VERIFY_IS_APPROX(dm4=m2*refMat3.transpose(), refMat4=refMat2*refMat3.transpose());
+    VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3, refMat4=refMat2.transpose()*refMat3);
+    VERIFY_IS_APPROX(dm4=m2.transpose()*refMat3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
+
+    // dense * sparse
+    VERIFY_IS_APPROX(dm4=refMat2*m3, refMat4=refMat2*refMat3);
+    VERIFY_IS_APPROX(dm4=refMat2*m3.transpose(), refMat4=refMat2*refMat3.transpose());
+    VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3, refMat4=refMat2.transpose()*refMat3);
+    VERIFY_IS_APPROX(dm4=refMat2.transpose()*m3.transpose(), refMat4=refMat2.transpose()*refMat3.transpose());
+
+    VERIFY_IS_APPROX(m3=m3*m3, refMat3=refMat3*refMat3);
+  }
+
+  // test matrix - diagonal product
+  if(false) // it compiles, but the precision is terrible. probably doesn't matter in this branch....
+  {
+    DenseMatrix refM2 = DenseMatrix::Zero(rows, rows);
+    DenseMatrix refM3 = DenseMatrix::Zero(rows, rows);
+    DiagonalMatrix<DenseVector> d1(DenseVector::Random(rows));
+    SparseMatrixType m2(rows, rows);
+    SparseMatrixType m3(rows, rows);
+    initSparse<Scalar>(density, refM2, m2);
+    initSparse<Scalar>(density, refM3, m3);
+    VERIFY_IS_APPROX(m3=m2*d1, refM3=refM2*d1);
+    VERIFY_IS_APPROX(m3=m2.transpose()*d1, refM3=refM2.transpose()*d1);
+    VERIFY_IS_APPROX(m3=d1*m2, refM3=d1*refM2);
+    VERIFY_IS_APPROX(m3=d1*m2.transpose(), refM3=d1 * refM2.transpose());
+  }
+
+  // test self adjoint products
+  {
+    DenseMatrix b = DenseMatrix::Random(rows, rows);
+    DenseMatrix x = DenseMatrix::Random(rows, rows);
+    DenseMatrix refX = DenseMatrix::Random(rows, rows);
+    DenseMatrix refUp = DenseMatrix::Zero(rows, rows);
+    DenseMatrix refLo = DenseMatrix::Zero(rows, rows);
+    DenseMatrix refS = DenseMatrix::Zero(rows, rows);
+    SparseMatrixType mUp(rows, rows);
+    SparseMatrixType mLo(rows, rows);
+    SparseMatrixType mS(rows, rows);
+    do {
+      initSparse<Scalar>(density, refUp, mUp, ForceRealDiag|/*ForceNonZeroDiag|*/MakeUpperTriangular);
+    } while (refUp.isZero());
+    refLo = refUp.transpose().conjugate();
+    mLo = mUp.transpose().conjugate();
+    refS = refUp + refLo;
+    refS.diagonal() *= 0.5;
+    mS = mUp + mLo;
+    for (int k=0; k<mS.outerSize(); ++k)
+      for (typename SparseMatrixType::InnerIterator it(mS,k); it; ++it)
+        if (it.index() == k)
+          it.valueRef() *= 0.5;
+
+    VERIFY_IS_APPROX(refS.adjoint(), refS);
+    VERIFY_IS_APPROX(mS.transpose().conjugate(), mS);
+    VERIFY_IS_APPROX(mS, refS);
+    VERIFY_IS_APPROX(x=mS*b, refX=refS*b);
+    VERIFY_IS_APPROX(x=mUp.template marked<UpperTriangular|SelfAdjoint>()*b, refX=refS*b);
+    VERIFY_IS_APPROX(x=mLo.template marked<LowerTriangular|SelfAdjoint>()*b, refX=refS*b);
+    VERIFY_IS_APPROX(x=mS.template marked<SelfAdjoint>()*b, refX=refS*b);
+  }
+
+}
+
+void test_sparse_product()
+{
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST( sparse_product(SparseMatrix<double>(8, 8)) );
+    CALL_SUBTEST( sparse_product(SparseMatrix<std::complex<double> >(16, 16)) );
+    CALL_SUBTEST( sparse_product(SparseMatrix<double>(33, 33)) );
+
+    CALL_SUBTEST( sparse_product(DynamicSparseMatrix<double>(8, 8)) );
+  }
+}

test/sparse_vector.cpp

@@ -84,6 +84,7 @@ template<typename Scalar> void sparse_vector(int rows, int cols)
   VERIFY_IS_APPROX(v1-=v2, refV1-=refV2);
 
   VERIFY_IS_APPROX(v1.dot(v2), refV1.dot(refV2));
+  VERIFY_IS_APPROX(v1.dot(refV2), refV1.dot(refV2));
 
 }
 

test/stdvector.cpp

@@ -22,8 +22,8 @@
 // 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/StdVector>
+#include "main.h"
 #include <Eigen/Geometry>
 
 template<typename MatrixType>

test/swap.cpp

@@ -0,0 +1,98 @@
+// 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/>.
+
+#define EIGEN_NO_STATIC_ASSERT
+#include "main.h"
+
+template<typename T>
+struct other_matrix_type
+{
+  typedef int type;
+};
+
+template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols>
+struct other_matrix_type<Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> >
+{
+  typedef Matrix<_Scalar, _Rows, _Cols, _Options^RowMajor, _MaxRows, _MaxCols> type;
+};
+
+template<typename MatrixType> void swap(const MatrixType& m)
+{
+  typedef typename other_matrix_type<MatrixType>::type OtherMatrixType;
+  typedef typename MatrixType::Scalar Scalar;
+
+  ei_assert((!ei_is_same_type<MatrixType,OtherMatrixType>::ret));
+  int rows = m.rows();
+  int cols = m.cols();
+  
+  // construct 3 matrix guaranteed to be distinct
+  MatrixType m1 = MatrixType::Random(rows,cols);
+  MatrixType m2 = MatrixType::Random(rows,cols) + Scalar(100) * MatrixType::Identity(rows,cols);
+  OtherMatrixType m3 = OtherMatrixType::Random(rows,cols) + Scalar(200) * OtherMatrixType::Identity(rows,cols);
+  
+  MatrixType m1_copy = m1;
+  MatrixType m2_copy = m2;
+  OtherMatrixType m3_copy = m3;
+  
+  // test swapping 2 matrices of same type
+  m1.swap(m2);
+  VERIFY_IS_APPROX(m1,m2_copy);
+  VERIFY_IS_APPROX(m2,m1_copy);
+  m1 = m1_copy;
+  m2 = m2_copy;
+  
+  // test swapping 2 matrices of different types
+  m1.swap(m3);
+  VERIFY_IS_APPROX(m1,m3_copy);
+  VERIFY_IS_APPROX(m3,m1_copy);
+  m1 = m1_copy;
+  m3 = m3_copy;
+  
+  // test swapping matrix with expression
+  m1.swap(m2.block(0,0,rows,cols));
+  VERIFY_IS_APPROX(m1,m2_copy);
+  VERIFY_IS_APPROX(m2,m1_copy);
+  m1 = m1_copy;
+  m2 = m2_copy;
+
+  // test swapping two expressions of different types
+  m1.transpose().swap(m3.transpose());
+  VERIFY_IS_APPROX(m1,m3_copy);
+  VERIFY_IS_APPROX(m3,m1_copy);
+  m1 = m1_copy;
+  m3 = m3_copy;
+  
+  // test assertion on mismatching size -- matrix case
+  VERIFY_RAISES_ASSERT(m1.swap(m1.row(0)));
+  // test assertion on mismatching size -- xpr case
+  VERIFY_RAISES_ASSERT(m1.row(0).swap(m1));
+}
+
+void test_swap()
+{
+  CALL_SUBTEST( swap(Matrix3f()) ); // fixed size, no vectorization 
+  CALL_SUBTEST( swap(Matrix4d()) ); // fixed size, possible vectorization 
+  CALL_SUBTEST( swap(MatrixXd(3,3)) ); // dyn size, no vectorization 
+  CALL_SUBTEST( swap(MatrixXf(30,30)) ); // dyn size, possible vectorization 
+}

test/testsuite.cmake

@@ -39,7 +39,7 @@
 # VERSION=opensuse-11.1
 # WORK_DIR=/home/gael/Coding/eigen2/cdash
 # # get the last version of the script
-# svn cat svn://anonsvn.kde.org/home/kde/trunk/kdesupport/eigen2/test/testsuite.cmake > $WORK_DIR/testsuite.cmake
+# wget http://bitbucket.org/eigen/eigen2/raw/tip/test/testsuite.cmake -o $WORK_DIR/testsuite.cmake
 # COMMON="ctest -S $WORK_DIR/testsuite.cmake,EIGEN_WORK_DIR=$WORK_DIR,EIGEN_SITE=$SITE,EIGEN_MODE=$1,EIGEN_BUILD_STRING=$OS_VERSION-$ARCH"
 # $COMMON-gcc-3.4.6,EIGEN_CXX=g++-3.4
 # $COMMON-gcc-4.0.1,EIGEN_CXX=g++-4.0.1
@@ -132,11 +132,12 @@ endif(NOT EIGEN_MODE)
 
 ## mandatory variables (the default should be ok in most cases):
 
-SET (CTEST_CVS_COMMAND "svn")
-SET (CTEST_CVS_CHECKOUT "${CTEST_CVS_COMMAND} co svn://anonsvn.kde.org/home/kde/trunk/kdesupport/eigen2 \"${CTEST_SOURCE_DIRECTORY}\"")
+SET (CTEST_CVS_COMMAND "hg")
+SET (CTEST_CVS_CHECKOUT "${CTEST_CVS_COMMAND} clone -r 2.0 http://bitbucket.org/eigen/eigen2 \"${CTEST_SOURCE_DIRECTORY}\"")
 
 # which ctest command to use for running the dashboard
 SET (CTEST_COMMAND "${EIGEN_CMAKE_DIR}ctest -D ${EIGEN_MODE}")
+
 # what cmake command to use for configuring this dashboard
 SET (CTEST_CMAKE_COMMAND "${EIGEN_CMAKE_DIR}cmake -DEIGEN_BUILD_TESTS=on ")
 

test/triangular.cpp

@@ -114,6 +114,14 @@ template<typename MatrixType> void triangular(const MatrixType& m)
 
   VERIFY((m1.template part<Eigen::UpperTriangular>() * m2.template part<Eigen::UpperTriangular>()).isUpperTriangular());
 
+  // test swap
+  m1.setOnes();
+  m2.setZero();
+  m2.template part<Eigen::UpperTriangular>().swap(m1);
+  m3.setZero();
+  m3.template part<Eigen::UpperTriangular>().setOnes();
+  VERIFY_IS_APPROX(m2,m3);
+
 }
 
 void test_triangular()