bump to 2.0.4

.

CMakeLists.txt

@@ -1,5 +1,5 @@
 project(Eigen)
-set(EIGEN_VERSION_NUMBER "2.0.2")
+set(EIGEN_VERSION_NUMBER "2.0.4")
 
 #if the svnversion program is absent, this will leave the SVN_REVISION string empty,
 #but won't stop CMake.
@@ -85,6 +85,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/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/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/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/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;

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

@@ -507,10 +507,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<(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.
       *

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();
 }
 
 /***************************************************************************

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/util/Macros.h

@@ -30,7 +30,7 @@
 
 #define EIGEN_WORLD_VERSION 2
 #define EIGEN_MAJOR_VERSION 0
-#define EIGEN_MINOR_VERSION 2
+#define EIGEN_MINOR_VERSION 4
 
 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
                                       (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \

Eigen/src/Core/util/Memory.h

@@ -238,16 +238,39 @@ inline static int ei_alignmentOffset(const Scalar* ptr, int maxOffset)
 
 
 #if EIGEN_ARCH_WANTS_ALIGNMENT
+  #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) throw() { \
+      void *operator new(size_t size) { \
         return Eigen::ei_conditional_aligned_malloc<NeedsToAlign>(size); \
       } \
-      void *operator new[](size_t size) throw() { \
+      void *operator new[](size_t size) { \
         return Eigen::ei_conditional_aligned_malloc<NeedsToAlign>(size); \
       } \
-      void operator delete(void * ptr) { Eigen::ei_conditional_aligned_free<NeedsToAlign>(ptr); } \
-      void operator delete[](void * ptr) { Eigen::ei_conditional_aligned_free<NeedsToAlign>(ptr); } \
-      void *operator new(size_t, void *ptr) throw() { return ptr; } \
+      void operator delete(void * ptr) throw() { Eigen::ei_conditional_aligned_free<NeedsToAlign>(ptr); } \
+      void operator delete[](void * ptr) throw() { Eigen::ei_conditional_aligned_free<NeedsToAlign>(ptr); } \
+      /* in-place new and delete. since (at least afaik) there is no actual   */ \
+      /* memory allocated we can safely let the default implementation handle */ \
+      /* this particular case. */ \
+      static void *operator new(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)

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/EulerAngles.h

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

Eigen/src/Geometry/Quaternion.h

@@ -346,7 +346,6 @@ inline Quaternion<Scalar>& Quaternion<Scalar>::setFromTwoVectors(const MatrixBas
 {
   Vector3 v0 = a.normalized();
   Vector3 v1 = b.normalized();
-  Vector3 axis = v0.cross(v1);
   Scalar c = v0.dot(v1);
 
   // if dot == 1, vectors are the same
@@ -354,7 +353,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/Transform.h

@@ -335,9 +335,9 @@ template<typename Scalar, int Dim>
 QMatrix Transform<Scalar,Dim>::toQMatrix(void) const
 {
   EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE)
-  return QMatrix(other.coeffRef(0,0), other.coeffRef(1,0),
-                 other.coeffRef(0,1), other.coeffRef(1,1),
-                 other.coeffRef(0,2), other.coeffRef(1,2));
+  return QMatrix(m_matrix.coeff(0,0), m_matrix.coeff(1,0),
+                 m_matrix.coeff(0,1), m_matrix.coeff(1,1),
+                 m_matrix.coeff(0,2), m_matrix.coeff(1,2));
 }
 
 /** Initialises \c *this from a QTransform assuming the dimension is 2.
@@ -369,12 +369,12 @@ Transform<Scalar,Dim>& Transform<Scalar,Dim>::operator=(const QTransform& other)
   * This function is available only if the token EIGEN_QT_SUPPORT is defined.
   */
 template<typename Scalar, int Dim>
-QMatrix Transform<Scalar,Dim>::toQTransform(void) const
+QTransform Transform<Scalar,Dim>::toQTransform(void) const
 {
   EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE)
-  return QTransform(other.coeffRef(0,0), other.coeffRef(1,0), other.coeffRef(2,0)
-                    other.coeffRef(0,1), other.coeffRef(1,1), other.coeffRef(2,1)
-                    other.coeffRef(0,2), other.coeffRef(1,2), other.coeffRef(2,2);
+  return QTransform(m_matrix.coeff(0,0), m_matrix.coeff(1,0), m_matrix.coeff(2,0),
+                    m_matrix.coeff(0,1), m_matrix.coeff(1,1), m_matrix.coeff(2,1),
+                    m_matrix.coeff(0,2), m_matrix.coeff(1,2), m_matrix.coeff(2,2));
 }
 #endif
 

Eigen/src/LU/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++)
@@ -506,7 +511,7 @@ bool LU<MatrixType>::solve(
     RealScalar biggest_in_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff();
     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/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

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/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.
+
 
 
 

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

@@ -15,11 +15,24 @@ is explained here: http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html
 There are 3 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
 
+\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
 
 If you have code like this,
@@ -64,6 +77,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.

test/CMakeLists.txt

@@ -217,6 +217,7 @@ 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)
 
 # print a summary of the different options
 message("************************************************************")

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/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/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 ")
 

unsupported/CMakeLists.txt

@@ -0,0 +1 @@
+add_subdirectory(Eigen)

unsupported/Eigen/CMakeLists.txt

@@ -0,0 +1,8 @@
+set(Eigen_HEADERS IterativeSolvers)
+
+install(FILES
+  ${Eigen_HEADERS}
+  DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen COMPONENT Devel
+  )
+
+add_subdirectory(src)

unsupported/Eigen/IterativeSolvers

@@ -0,0 +1,51 @@
+// 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-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_ITERATIVE_SOLVERS_MODULE_H
+#define EIGEN_ITERATIVE_SOLVERS_MODULE_H
+
+#include <Eigen/Core>
+
+namespace Eigen {
+
+/** \ingroup Unsupported_modules
+  * \defgroup IterativeSolvers_Module Iterative solvers module
+  * This module aims to provide various iterative linear and non linear solver algorithms.
+  * It currently provides:
+  *  - a constrained conjugate gradient
+  *
+  * \code
+  * #include <unsupported/Eigen/IterativeSolvers>
+  * \endcode
+  */
+//@{
+
+#include "src/IterativeSolvers/IterationController.h"
+#include "src/IterativeSolvers/ConstrainedConjGrad.h"
+
+//@}
+
+}
+
+#endif // EIGEN_ITERATIVE_SOLVERS_MODULE_H

unsupported/Eigen/src/CMakeLists.txt

@@ -0,0 +1 @@
+ADD_SUBDIRECTORY(IterativeSolvers)

unsupported/Eigen/src/IterativeSolvers/CMakeLists.txt

@@ -0,0 +1,6 @@
+FILE(GLOB Eigen_IterativeSolvers_SRCS "*.h")
+
+INSTALL(FILES
+  ${Eigen_IterativeSolvers_SRCS}
+  DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/IterativeSolvers COMPONENT Devel
+  )

unsupported/Eigen/src/IterativeSolvers/ConstrainedConjGrad.h

@@ -0,0 +1,192 @@
+// 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>
+//
+// 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/>.
+
+/* NOTE The functions of this file have been adapted from the GMM++ library */
+
+//========================================================================
+//
+// Copyright (C) 2002-2007 Yves Renard
+//
+// This file is a part of GETFEM++
+//
+// Getfem++ 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; version 2.1 of the License.
+//
+// This program 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 for more details.
+// You should have received a copy of the GNU Lesser General Public
+// License along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301,
+// USA.
+//
+//========================================================================
+
+#ifndef EIGEN_CONSTRAINEDCG_H
+#define EIGEN_CONSTRAINEDCG_H
+
+#include <Eigen/Core>
+
+/** \ingroup IterativeSolvers_Module
+  * Compute the pseudo inverse of the non-square matrix C such that
+  * \f$ CINV = (C * C^T)^{-1} * C \f$ based on a conjugate gradient method.
+  *
+  * This function is internally used by ei_constrained_cg.
+  */
+template <typename CMatrix, typename CINVMatrix>
+void ei_pseudo_inverse(const CMatrix &C, CINVMatrix &CINV)
+{
+  // optimisable : copie de la ligne, precalcul de C * trans(C).
+  typedef typename CMatrix::Scalar Scalar;
+  // FIXME use sparse vectors ?
+  typedef Matrix<Scalar,Dynamic,1> TmpVec;
+
+  int rows = C.rows(), cols = C.cols();
+
+  TmpVec d(rows), e(rows), l(cols), p(rows), q(rows), r(rows);
+  Scalar rho, rho_1, alpha;
+  d.setZero();
+
+  CINV.startFill(); // FIXME estimate the number of non-zeros
+  for (int i = 0; i < rows; ++i)
+  {
+    d[i] = 1.0;
+    rho = 1.0;
+    e.setZero();
+    r = d;
+    p = d;
+
+    while (rho >= 1e-38)
+    { /* conjugate gradient to compute e             */
+      /* which is the i-th row of inv(C * trans(C))  */
+      l = C.transpose() * p;
+      q = C * l;
+      alpha = rho / p.dot(q);
+      e +=  alpha * p;
+      r += -alpha * q;
+      rho_1 = rho;
+      rho = r.dot(r);
+      p = (rho/rho_1) * p + r;
+    }
+
+    l = C.transpose() * e; // l is the i-th row of CINV
+    // FIXME add a generic "prune/filter" expression for both dense and sparse object to sparse
+    for (int j=0; j<l.size(); ++j)
+      if (l[j]<1e-15)
+        CINV.fill(i,j) = l[j];
+
+    d[i] = 0.0;
+  }
+  CINV.endFill();
+}
+
+
+
+/** \ingroup IterativeSolvers_Module
+  * Constrained conjugate gradient
+  *
+  * Computes the minimum of \f$ 1/2((Ax).x) - bx \f$ under the contraint \f$ Cx <= f @\$
+  */
+template<typename TMatrix, typename CMatrix,
+         typename VectorX, typename VectorB, typename VectorF>
+void ei_constrained_cg(const TMatrix& A, const CMatrix& C, VectorX& x,
+                       const VectorB& b, const VectorF& f, IterationController &iter)
+{
+  typedef typename TMatrix::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,1>  TmpVec;
+
+  Scalar rho = 1.0, rho_1, lambda, gamma;
+  int xSize = x.size();
+  TmpVec  p(xSize), q(xSize), q2(xSize),
+          r(xSize), old_z(xSize), z(xSize),
+          memox(xSize);
+  std::vector<bool> satured(C.rows());
+  p.setZero();
+  iter.setRhsNorm(ei_sqrt(b.dot(b))); // gael vect_sp(PS, b, b)
+  if (iter.rhsNorm() == 0.0) iter.setRhsNorm(1.0);
+
+  SparseMatrix<Scalar,RowMajor> CINV(C.rows(), C.cols());
+  ei_pseudo_inverse(C, CINV);
+
+  while(true)
+  {
+    // computation of residual
+    old_z = z;
+    memox = x;
+    r = b;
+    r += A * -x;
+    z = r;
+    bool transition = false;
+    for (int i = 0; i < C.rows(); ++i)
+    {
+      Scalar al = C.row(i).dot(x) - f.coeff(i);
+      if (al >= -1.0E-15)
+      {
+        if (!satured[i])
+        {
+          satured[i] = true;
+          transition = true;
+        }
+        Scalar bb = CINV.row(i).dot(z);
+        if (bb > 0.0)
+          // FIXME: we should allow that: z += -bb * C.row(i);
+          for (typename CMatrix::InnerIterator it(C,i); it; ++it)
+            z.coeffRef(it.index()) -= bb*it.value();
+      }
+      else
+        satured[i] = false;
+    }
+
+    // descent direction
+    rho_1 = rho;
+    rho = r.dot(z);
+
+    if (iter.finished(rho)) break;
+
+    if (iter.noiseLevel() > 0 && transition) std::cerr << "CCG: transition\n";
+    if (transition || iter.first()) gamma = 0.0;
+    else gamma = std::max(0.0, (rho - old_z.dot(z)) / rho_1);
+    p = z + gamma*p;
+
+    ++iter;
+    // one dimensionnal optimization
+    q = A * p;
+    lambda = rho / q.dot(p);
+    for (int i = 0; i < C.rows(); ++i)
+    {
+      if (!satured[i])
+      {
+        Scalar bb = C.row(i).dot(p) - f[i];
+        if (bb > 0.0)
+          lambda = std::min(lambda, (f.coeff(i)-C.row(i).dot(x)) / bb);
+      }
+    }
+    x += lambda * p;
+    memox -= x;
+  }
+}
+
+#endif // EIGEN_CONSTRAINEDCG_H

unsupported/Eigen/src/IterativeSolvers/IterationController.h

@@ -0,0 +1,166 @@
+// 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-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/>.
+
+/* NOTE The class IterationController has been adapted from the iteration
+ *      class of the GMM++ and ITL libraries.
+ */
+
+//=======================================================================
+// Copyright (C) 1997-2001
+// Authors: Andrew Lumsdaine <lums@osl.iu.edu> 
+//          Lie-Quan Lee     <llee@osl.iu.edu>
+//
+// This file is part of the Iterative Template Library
+//
+// You should have received a copy of the License Agreement for the
+// Iterative Template Library along with the software;  see the
+// file LICENSE.  
+//
+// Permission to modify the code and to distribute modified code is
+// granted, provided the text of this NOTICE is retained, a notice that
+// the code was modified is included with the above COPYRIGHT NOTICE and
+// with the COPYRIGHT NOTICE in the LICENSE file, and that the LICENSE
+// file is distributed with the modified code.
+//
+// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
+// By way of example, but not limitation, Licensor MAKES NO
+// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
+// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
+// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
+// OR OTHER RIGHTS.
+//=======================================================================
+
+//========================================================================
+//
+// Copyright (C) 2002-2007 Yves Renard
+//
+// This file is a part of GETFEM++
+//
+// Getfem++ 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; version 2.1 of the License.
+//
+// This program 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 for more details.
+// You should have received a copy of the GNU Lesser General Public
+// License along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301,
+// USA.
+//
+//========================================================================
+
+#ifndef EIGEN_ITERATION_CONTROLLER_H
+#define EIGEN_ITERATION_CONTROLLER_H
+
+/** \ingroup IterativeSolvers_Module
+  * \class IterationController
+  *
+  * \brief Controls the iterations of the iterative solvers
+  *
+  * This class has been adapted from the iteration class of GMM++ and ITL libraries.
+  *
+  */
+class IterationController
+{
+  protected :
+    double m_rhsn;        ///< Right hand side norm
+    size_t m_maxiter;     ///< Max. number of iterations
+    int m_noise;          ///< if noise > 0 iterations are printed
+    double m_resmax;      ///< maximum residual
+    double m_resminreach, m_resadd;
+    size_t m_nit;         ///< iteration number
+    double m_res;         ///< last computed residual
+    bool m_written;
+    void (*m_callback)(const IterationController&);
+  public :
+
+    void init()
+    {
+      m_nit = 0; m_res = 0.0; m_written = false;
+      m_resminreach = 1E50; m_resadd = 0.0;
+      m_callback = 0;
+    }
+
+    IterationController(double r = 1.0E-8, int noi = 0, size_t mit = size_t(-1))
+      : m_rhsn(1.0), m_maxiter(mit), m_noise(noi), m_resmax(r) { init(); }
+
+    void operator ++(int) { m_nit++; m_written = false; m_resadd += m_res; }
+    void operator ++() { (*this)++; }
+
+    bool first() { return m_nit == 0; }
+
+    /* get/set the "noisyness" (verbosity) of the solvers */
+    int noiseLevel() const { return m_noise; }
+    void setNoiseLevel(int n) { m_noise = n; }
+    void reduceNoiseLevel() { if (m_noise > 0) m_noise--; }
+
+    double maxResidual() const { return m_resmax; }
+    void setMaxResidual(double r) { m_resmax = r; }
+
+    double residual() const { return m_res; }
+
+    /* change the user-definable callback, called after each iteration */
+    void setCallback(void (*t)(const IterationController&))
+    {
+      m_callback = t;
+    }
+
+    size_t iteration() const { return m_nit; }
+    void setIteration(size_t i) { m_nit = i; }
+
+    size_t maxIterarions() const { return m_maxiter; }
+    void setMaxIterations(size_t i) { m_maxiter = i; }
+
+    double rhsNorm() const { return m_rhsn; }
+    void setRhsNorm(double r) { m_rhsn = r; }
+
+    bool converged() const { return m_res <= m_rhsn * m_resmax; }
+    bool converged(double nr)
+    {
+      m_res = ei_abs(nr); 
+      m_resminreach = std::min(m_resminreach, m_res);
+      return converged();
+    }
+    template<typename VectorType> bool converged(const VectorType &v)
+    { return converged(v.squaredNorm()); }
+
+    bool finished(double nr)
+    {
+      if (m_callback) m_callback(*this);
+      if (m_noise > 0 && !m_written)
+      {
+        converged(nr);
+        m_written = true;
+      }
+      return (m_nit >= m_maxiter || converged(nr));
+    }
+    template <typename VectorType>
+    bool finished(const MatrixBase<VectorType> &v)
+    { return finished(double(v.squaredNorm())); }
+
+};
+
+#endif // EIGEN_ITERATION_CONTROLLER_H