diff --git a/CMakeLists.txt b/CMakeLists.txt
index c96a1f0cd..70f92c215 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,5 +1,5 @@
 project(Eigen)
-set(EIGEN_VERSION_NUMBER "2.0.2")
+set(EIGEN_VERSION_NUMBER "2.0.3")
 
 #if the svnversion program is absent, this will leave the SVN_REVISION string empty,
 #but won't stop CMake.
diff --git a/Eigen/Sparse b/Eigen/Sparse
index 57e9af6a2..536c28454 100644
--- a/Eigen/Sparse
+++ b/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"
diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h
index 07eaf0747..01f01fdf2 100644
--- a/Eigen/src/Core/DiagonalMatrix.h
+++ b/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>
diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h
index 064badf1e..1c072c12b 100644
--- a/Eigen/src/Core/util/Macros.h
+++ b/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 3
 
 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
                                       (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \
diff --git a/Eigen/src/Core/util/StaticAssert.h b/Eigen/src/Core/util/StaticAssert.h
index 60055683e..2c13098a2 100644
--- a/Eigen/src/Core/util/StaticAssert.h
+++ b/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
       };
     };
 
diff --git a/Eigen/src/Sparse/CMakeLists.txt b/Eigen/src/Sparse/CMakeLists.txt
index b0ba6e140..aa1468812 100644
--- a/Eigen/src/Sparse/CMakeLists.txt
+++ b/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
   )
diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h
index d666fa1e9..dfd9c787a 100644
--- a/Eigen/src/Sparse/CholmodSupport.h
+++ b/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);
diff --git a/Eigen/src/Sparse/DynamicSparseMatrix.h b/Eigen/src/Sparse/DynamicSparseMatrix.h
index d59b81458..7119a84bd 100644
--- a/Eigen/src/Sparse/DynamicSparseMatrix.h
+++ b/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);
     }
diff --git a/Eigen/src/Sparse/MappedSparseMatrix.h b/Eigen/src/Sparse/MappedSparseMatrix.h
index b4a12cf0d..f4935d834 100644
--- a/Eigen/src/Sparse/MappedSparseMatrix.h
+++ b/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; }
 
diff --git a/Eigen/src/Sparse/SparseBlock.h b/Eigen/src/Sparse/SparseBlock.h
index cbec36e61..c39066676 100644
--- a/Eigen/src/Sparse/SparseBlock.h
+++ b/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>
diff --git a/Eigen/src/Sparse/SparseCwiseBinaryOp.h b/Eigen/src/Sparse/SparseCwiseBinaryOp.h
index 87fd429be..d19970efc 100644
--- a/Eigen/src/Sparse/SparseCwiseBinaryOp.h
+++ b/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>,
diff --git a/Eigen/src/Sparse/SparseCwiseUnaryOp.h b/Eigen/src/Sparse/SparseCwiseUnaryOp.h
index ff7c9c73f..b11c0f8a3 100644
--- a/Eigen/src/Sparse/SparseCwiseUnaryOp.h
+++ b/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>
diff --git a/Eigen/src/Sparse/SparseDiagonalProduct.h b/Eigen/src/Sparse/SparseDiagonalProduct.h
new file mode 100644
index 000000000..932daf220
--- /dev/null
+++ b/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
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index 23766d516..3f09596bc 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/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);
diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h
index 46c0c546a..468bc9e22 100644
--- a/Eigen/src/Sparse/SparseMatrixBase.h
+++ b/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
diff --git a/Eigen/src/Sparse/SparseProduct.h b/Eigen/src/Sparse/SparseProduct.h
index ec4961d9b..c98a71e99 100644
--- a/Eigen/src/Sparse/SparseProduct.h
+++ b/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);
     }
diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h
index 18d0ee238..393cdda6e 100644
--- a/Eigen/src/Sparse/SparseUtil.h
+++ b/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;
diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h
index 457984cad..8e5a6efed 100644
--- a/Eigen/src/Sparse/SparseVector.h
+++ b/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)
 //       {
diff --git a/Eigen/src/Sparse/SuperLUSupport.h b/Eigen/src/Sparse/SuperLUSupport.h
index ded145eeb..3c9a4fcce 100644
--- a/Eigen/src/Sparse/SuperLUSupport.h
+++ b/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)
 {
diff --git a/Eigen/src/Sparse/TaucsSupport.h b/Eigen/src/Sparse/TaucsSupport.h
index 375d29f32..4dddca7b6 100644
--- a/Eigen/src/Sparse/TaucsSupport.h
+++ b/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
diff --git a/Eigen/src/Sparse/UmfPackSupport.h b/Eigen/src/Sparse/UmfPackSupport.h
index 1bb40a7c3..b76ffb252 100644
--- a/Eigen/src/Sparse/UmfPackSupport.h
+++ b/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
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index f3c7454e3..70656de9a 100644
--- a/test/CMakeLists.txt
+++ b/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("************************************************************")
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp
index 439458128..410ef96a6 100644
--- a/test/sparse_basic.cpp
+++ b/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
   {
diff --git a/test/sparse_product.cpp b/test/sparse_product.cpp
new file mode 100644
index 000000000..dcfc58a14
--- /dev/null
+++ b/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)) );
+  }
+}
diff --git a/test/sparse_vector.cpp b/test/sparse_vector.cpp
index 8207e522a..934719f2c 100644
--- a/test/sparse_vector.cpp
+++ b/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));
 
 }