various work on the Sparse module:

* added some glue to Eigen/Core (SparseBit, ei_eval, Matrix)
* add two new sparse matrix types:
   HashMatrix: based on std::map (for random writes)
   LinkedVectorMatrix: array of linked vectors
   (for outer coherent writes, e.g. to transpose a matrix)
* add a SparseSetter class to easily set/update any kind of matrices, e.g.:
   { SparseSetter<MatrixType,RandomAccessPattern> wrapper(mymatrix);
     for (...) wrapper->coeffRef(rand(),rand()) = rand(); }
* automatic shallow copy for RValue
* and a lot of mess !
plus:
* remove the remaining ArrayBit related stuff
* don't use alloca in product for very large memory allocation

16 files changed, 1106 insertions, 160 deletions

diff --git a/Eigen/Sparse b/Eigen/Sparse
index 6a5eaf6a0..962660cbb 100644
--- a/Eigen/Sparse
+++ b/Eigen/Sparse
@@ -4,12 +4,20 @@
 #include "Core"
 #include <vector>
 #include <map>
+#include <stdlib.h>
+#include <string.h>
 
 namespace Eigen {
 
+#include "src/Sparse/SparseUtil.h"
+#include "src/Sparse/SparseMatrixBase.h"
+#include "src/Sparse/SparseProduct.h"
 #include "src/Sparse/SparseArray.h"
 #include "src/Sparse/SparseMatrix.h"
+#include "src/Sparse/HashMatrix.h"
+#include "src/Sparse/LinkedVectorMatrix.h"
 #include "src/Sparse/CoreIterators.h"
+#include "src/Sparse/SparseSetter.h"
 
 } // namespace Eigen
 
diff --git a/Eigen/src/Core/CacheFriendlyProduct.h b/Eigen/src/Core/CacheFriendlyProduct.h
index 4a0e4e24a..0da84eeac 100644
--- a/Eigen/src/Core/CacheFriendlyProduct.h
+++ b/Eigen/src/Core/CacheFriendlyProduct.h
@@ -86,7 +86,12 @@ static void ei_cache_friendly_product(
   const int l2BlockRows = MaxL2BlockSize > rows ? rows : MaxL2BlockSize;
   const int l2BlockCols = MaxL2BlockSize > cols ? cols : MaxL2BlockSize;
   const int l2BlockSize = MaxL2BlockSize > size ? size : MaxL2BlockSize;
-  Scalar* __restrict__ block   = (Scalar*)alloca(sizeof(Scalar)*l2BlockRows*size);
+  Scalar* __restrict__ block = 0;
+  const int allocBlockSize = sizeof(Scalar)*l2BlockRows*size;
+  if (allocBlockSize>16000000)
+    block = (Scalar*)malloc(allocBlockSize);
+  else
+    block = (Scalar*)alloca(allocBlockSize);
   Scalar* __restrict__ rhsCopy = (Scalar*)alloca(sizeof(Scalar)*l2BlockSize);
 
   // loops on each L2 cache friendly blocks of the result
@@ -347,6 +352,9 @@ static void ei_cache_friendly_product(
         }
     }
   }
+
+  if (allocBlockSize>16000000)
+    free(block);
 }
 
 #endif // EIGEN_CACHE_FRIENDLY_PRODUCT_H
diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h
index 40947e997..f00a27b33 100644
--- a/Eigen/src/Core/Matrix.h
+++ b/Eigen/src/Core/Matrix.h
@@ -92,7 +92,8 @@ struct ei_traits<Matrix<_Scalar, _Rows, _Cols, _MaxRows, _MaxCols, _Flags> >
                 _Rows, _Cols, _MaxRows, _MaxCols,
                 _Flags
             >::ret,
-    CoeffReadCost = NumTraits<Scalar>::ReadCost
+    CoeffReadCost = NumTraits<Scalar>::ReadCost,
+    SupportedAccessPatterns = RandomAccessPattern
   };
 };
 
diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index a0b29449c..e3b4b18e2 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -258,7 +258,6 @@ template<typename OtherDerived>
 inline const typename ProductReturnType<Derived,OtherDerived>::Type
 MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
 {
-  assert( (Derived::Flags&ArrayBit) == (OtherDerived::Flags&ArrayBit) );
   return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
 }
 
diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h
index f19d4affa..c364e5f86 100644
--- a/Eigen/src/Core/Transpose.h
+++ b/Eigen/src/Core/Transpose.h
@@ -63,6 +63,8 @@ template<typename MatrixType> class Transpose
 
     EIGEN_GENERIC_PUBLIC_INTERFACE(Transpose)
 
+    class InnerIterator;
+
     inline Transpose(const MatrixType& matrix) : m_matrix(matrix) {}
 
     EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Transpose)
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h
index 716d86243..eafcbca7f 100644
--- a/Eigen/src/Core/util/Constants.h
+++ b/Eigen/src/Core/util/Constants.h
@@ -138,10 +138,8 @@ const unsigned int LowerTriangularBit = 0x400;
 
 /** \ingroup flags
   *
-  * means the object is just an array of scalars, and operations on it are regarded as operations
-  * on every of these scalars taken separately.
-  */
-const unsigned int ArrayBit = 0x800;
+  * means the expression includes sparse matrices and the sparse path has to be taken. */
+const unsigned int SparseBit = 0x800;
 
 /** \ingroup flags
   *
@@ -155,7 +153,7 @@ const unsigned int HereditaryBits = RowMajorBit
                                   | EvalBeforeNestingBit
                                   | EvalBeforeAssigningBit
                                   | LargeBit
-                                  | ArrayBit;
+                                  | SparseBit;
 
 // Possible values for the Mode parameter of part() and of extract()
 const unsigned int Upper = UpperTriangularBit;
@@ -173,7 +171,7 @@ enum { Aligned=0, UnAligned=1 };
 enum { ConditionalJumpCost = 5 };
 enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight };
 enum DirectionType { Vertical, Horizontal };
-enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct };
+enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct, SparseProduct };
 
 enum {
   InnerVectorization,
@@ -188,5 +186,14 @@ enum {
   NoUnrolling
 };
 
+enum {
+  Dense   = 0,
+  Sparse  = SparseBit
+};
+
+const int FullyCoherentAccessPattern  = 0x1;
+const int InnerCoherentAccessPattern  = 0x2 | FullyCoherentAccessPattern;
+const int OuterCoherentAccessPattern  = 0x4 | InnerCoherentAccessPattern;
+const int RandomAccessPattern         = 0x8 | OuterCoherentAccessPattern;
 
 #endif // EIGEN_CONSTANTS_H
diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h
index 5d809f619..93d50441e 100644
--- a/Eigen/src/Core/util/Meta.h
+++ b/Eigen/src/Core/util/Meta.h
@@ -175,7 +175,9 @@ template<int _Rows, int _Cols> struct ei_size_at_compile_time
   enum { ret = (_Rows==Dynamic || _Cols==Dynamic) ? Dynamic : _Rows * _Cols };
 };
 
-template<typename T> class ei_eval
+template<typename T, int Sparseness = ei_traits<T>::Flags&SparseBit> class ei_eval;
+
+template<typename T> class ei_eval<T,Dense>
 {
     typedef typename ei_traits<T>::Scalar _Scalar;
     enum {_Rows = ei_traits<T>::RowsAtCompileTime,
diff --git a/Eigen/src/Sparse/CoreIterators.h b/Eigen/src/Sparse/CoreIterators.h
index 00f3ae781..f44db0c1f 100644
--- a/Eigen/src/Sparse/CoreIterators.h
+++ b/Eigen/src/Sparse/CoreIterators.h
@@ -25,30 +25,47 @@
 #ifndef EIGEN_COREITERATORS_H
 #define EIGEN_COREITERATORS_H
 
+/* This file contains the respective InnerIterator definition of the expressions defined in Eigen/Core
+ */
+
 template<typename Derived>
 class MatrixBase<Derived>::InnerIterator
 {
     typedef typename Derived::Scalar Scalar;
   public:
-    InnerIterator(const Derived& mat, int col)
-      : m_matrix(mat), m_row(0), m_col(col), m_end(mat.rows())
+    InnerIterator(const Derived& mat, int outer)
+      : m_matrix(mat), m_inner(0), m_outer(outer), m_end(mat.rows())
     {}
 
-    Scalar value() { return m_matrix.coeff(m_row, m_col); }
+    Scalar value()
+    {
+      return (Derived::Flags&RowMajorBit) ? m_matrix.coeff(m_outer, m_inner)
+                                          : m_matrix.coeff(m_inner, m_outer);
+    }
 
-    InnerIterator& operator++() { m_row++; return *this; }
+    InnerIterator& operator++() { m_inner++; return *this; }
 
-    int index() const { return m_row; }
+    int index() const { return m_inner; }
 
-    operator bool() const { return m_row < m_end && m_row>=0; }
+    operator bool() const { return m_inner < m_end && m_inner>=0; }
 
   protected:
     const Derived& m_matrix;
-    int m_row;
-    const int m_col;
+    int m_inner;
+    const int m_outer;
     const int m_end;
 };
 
+template<typename MatrixType>
+class Transpose<MatrixType>::InnerIterator : public MatrixType::InnerIterator
+{
+  public:
+
+    InnerIterator(const Transpose& trans, int outer)
+      : MatrixType::InnerIterator(trans.m_matrix, outer)
+    {}
+};
+
 template<typename UnaryOp, typename MatrixType>
 class CwiseUnaryOp<UnaryOp,MatrixType>::InnerIterator
 {
@@ -57,8 +74,8 @@ class CwiseUnaryOp<UnaryOp,MatrixType>::InnerIterator
     typedef typename _MatrixTypeNested::InnerIterator MatrixTypeIterator;
   public:
 
-    InnerIterator(const CwiseUnaryOp& unaryOp, int col)
-      : m_iter(unaryOp.m_matrix,col), m_functor(unaryOp.m_functor), m_id(-1)
+    InnerIterator(const CwiseUnaryOp& unaryOp, int outer)
+      : m_iter(unaryOp.m_matrix,outer), m_functor(unaryOp.m_functor), m_id(-1)
     {
       this->operator++();
     }
@@ -101,8 +118,8 @@ class CwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator
     typedef typename _RhsNested::InnerIterator RhsIterator;
   public:
 
-    InnerIterator(const CwiseBinaryOp& binOp, int col)
-      : m_lhsIter(binOp.m_lhs,col), m_rhsIter(binOp.m_rhs,col), m_functor(binOp.m_functor), m_id(-1)
+    InnerIterator(const CwiseBinaryOp& binOp, int outer)
+      : m_lhsIter(binOp.m_lhs,outer), m_rhsIter(binOp.m_rhs,outer), m_functor(binOp.m_functor), m_id(-1)
     {
       this->operator++();
     }
diff --git a/Eigen/src/Sparse/HashMatrix.h b/Eigen/src/Sparse/HashMatrix.h
new file mode 100644
index 000000000..753c7c5e9
--- /dev/null
+++ b/Eigen/src/Sparse/HashMatrix.h
@@ -0,0 +1,165 @@
+// 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/>.
+
+#ifndef EIGEN_HASHMATRIX_H
+#define EIGEN_HASHMATRIX_H
+
+template<typename _Scalar, int _Flags>
+struct ei_traits<HashMatrix<_Scalar, _Flags> >
+{
+  typedef _Scalar Scalar;
+  enum {
+    RowsAtCompileTime = Dynamic,
+    ColsAtCompileTime = Dynamic,
+    MaxRowsAtCompileTime = Dynamic,
+    MaxColsAtCompileTime = Dynamic,
+    Flags = _Flags,
+    CoeffReadCost = NumTraits<Scalar>::ReadCost,
+    SupportedAccessPatterns = RandomAccessPattern
+  };
+};
+
+// TODO reimplement this class using custom linked lists
+template<typename _Scalar, int _Flags>
+class HashMatrix : public SparseMatrixBase<HashMatrix<_Scalar, _Flags> >
+{
+  public:
+    EIGEN_GENERIC_PUBLIC_INTERFACE(HashMatrix)
+    class InnerIterator;
+  protected:
+
+    typedef typename std::map<int, Scalar>::iterator MapIterator;
+    typedef typename std::map<int, Scalar>::const_iterator ConstMapIterator;
+
+  public:
+    inline int rows() const { return m_innerSize; }
+    inline int cols() const { return m_data.size(); }
+
+    inline const Scalar& coeff(int row, int col) const
+    {
+      const MapIterator it = m_data[col].find(row);
+      if (it!=m_data[col].end())
+        return Scalar(0);
+      return it->second;
+    }
+
+    inline Scalar& coeffRef(int row, int col)
+    {
+      return m_data[col][row];
+    }
+
+  public:
+
+    inline void startFill(int reserveSize = 1000) {}
+
+    inline Scalar& fill(int row, int col) { return coeffRef(row, col); }
+
+    inline void endFill() {}
+
+    ~HashMatrix()
+    {}
+
+    inline void shallowCopy(const HashMatrix& other)
+    {
+      EIGEN_DBG_SPARSE(std::cout << "HashMatrix:: shallowCopy\n");
+      // FIXME implement a true shallow copy !!
+      resize(other.rows(), other.cols());
+      for (int j=0; j<this->outerSize(); ++j)
+        m_data[j] = other.m_data[j];
+    }
+
+    void resize(int _rows, int _cols)
+    {
+      if (cols() != _cols)
+      {
+        m_data.resize(_cols);
+      }
+      m_innerSize = _rows;
+    }
+
+    inline HashMatrix(int rows, int cols)
+      : m_innerSize(0)
+    {
+      resize(rows, cols);
+    }
+
+    template<typename OtherDerived>
+    inline HashMatrix(const MatrixBase<OtherDerived>& other)
+      : m_innerSize(0)
+    {
+      *this = other.derived();
+    }
+
+    inline HashMatrix& operator=(const HashMatrix& other)
+    {
+      if (other.isRValue())
+      {
+        shallowCopy(other);
+      }
+      else
+      {
+        resize(other.rows(), other.cols());
+        for (int col=0; col<cols(); ++col)
+          m_data[col] = other.m_data[col];
+      }
+      return *this;
+    }
+
+    template<typename OtherDerived>
+    inline HashMatrix& operator=(const MatrixBase<OtherDerived>& other)
+    {
+      return SparseMatrixBase<HashMatrix>::operator=(other);
+    }
+
+  protected:
+
+    std::vector<std::map<int, Scalar> > m_data;
+    int m_innerSize;
+
+};
+
+template<typename Scalar, int _Flags>
+class HashMatrix<Scalar,_Flags>::InnerIterator
+{
+  public:
+
+    InnerIterator(const HashMatrix& mat, int col)
+      : m_matrix(mat), m_it(mat.m_data[col].begin()), m_end(mat.m_data[col].end())
+    {}
+
+    InnerIterator& operator++() { m_it++; return *this; }
+
+    Scalar value() { return m_it->second; }
+
+    int index() const { return m_it->first; }
+
+    operator bool() const { return m_it!=m_end; }
+
+  protected:
+    const HashMatrix& m_matrix;
+    typename HashMatrix::ConstMapIterator m_it;
+    typename HashMatrix::ConstMapIterator m_end;
+};
+
+#endif // EIGEN_HASHMATRIX_H
diff --git a/Eigen/src/Sparse/LinkedVectorMatrix.h b/Eigen/src/Sparse/LinkedVectorMatrix.h
new file mode 100644
index 000000000..2e2618e26
--- /dev/null
+++ b/Eigen/src/Sparse/LinkedVectorMatrix.h
@@ -0,0 +1,288 @@
+// 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/>.
+
+#ifndef EIGEN_LINKEDVECTORMATRIX_H
+#define EIGEN_LINKEDVECTORMATRIX_H
+
+template<typename _Scalar, int _Flags>
+struct ei_traits<LinkedVectorMatrix<_Scalar,_Flags> >
+{
+  typedef _Scalar Scalar;
+  enum {
+    RowsAtCompileTime = Dynamic,
+    ColsAtCompileTime = Dynamic,
+    MaxRowsAtCompileTime = Dynamic,
+    MaxColsAtCompileTime = Dynamic,
+    Flags = _Flags,
+    CoeffReadCost = NumTraits<Scalar>::ReadCost,
+    SupportedAccessPatterns = InnerCoherentAccessPattern
+  };
+};
+
+template<typename Element, int BlockSize = 8>
+struct LinkedVector
+{
+  LinkedVector() : size(0), next(0), prev(0) {}
+  Element data[BlockSize];
+  LinkedVector* next;
+  LinkedVector* prev;
+  int size;
+  bool isFull() const { return size==BlockSize; }
+};
+
+template<typename _Scalar, int _Flags>
+class LinkedVectorMatrix : public SparseMatrixBase<LinkedVectorMatrix<_Scalar,_Flags> >
+{
+  public:
+    EIGEN_GENERIC_PUBLIC_INTERFACE(LinkedVectorMatrix)
+    class InnerIterator;
+  protected:
+
+    enum {
+      RowMajor = Flags&RowMajorBit ? 1 : 0
+    };
+
+    struct ValueIndex
+    {
+      ValueIndex() : value(0), index(0) {}
+      ValueIndex(Scalar v, int i) : value(v), index(i) {}
+      Scalar value;
+      int index;
+    };
+    typedef LinkedVector<ValueIndex,8> LinkedVectorBlock;
+
+    inline int find(LinkedVectorBlock** _el, int id)
+    {
+      LinkedVectorBlock* el = *_el;
+      while (el && el->data[el->size-1].index<id)
+        el = el->next;
+      *_el = el;
+      if (el)
+      {
+        // binary search
+        int maxI = el->size-1;
+        int minI = 0;
+        int i = el->size/2;
+        const ValueIndex* data = el->data;
+        while (data[i].index!=id)
+        {
+          if (data[i].index<id)
+          {
+            minI = i+1;
+            i = (maxI + minI)+2;
+          }
+          else
+          {
+            maxI = i-1;
+            i = (maxI + minI)+2;
+          }
+          if (minI>=maxI)
+            return -1;
+        }
+        if (data[i].index==id)
+          return i;
+      }
+      return -1;
+    }
+
+  public:
+    inline int rows() const { return RowMajor ? m_data.size() : m_innerSize; }
+    inline int cols() const { return RowMajor ? m_innerSize : m_data.size(); }
+
+    inline const Scalar& coeff(int row, int col) const
+    {
+      const int outer = RowMajor ? row : col;
+      const int inner = RowMajor ? col : row;
+
+      LinkedVectorBlock* el = m_data[outer];
+      int id = find(&el, inner);
+      if (id<0)
+        return Scalar(0);
+      return el->data[id].value;
+    }
+
+    inline Scalar& coeffRef(int row, int col)
+    {
+      const int outer = RowMajor ? row : col;
+      const int inner = RowMajor ? col : row;
+
+      LinkedVectorBlock* el = m_data[outer];
+      int id = find(&el, inner);
+      ei_assert(id>=0);
+//       if (id<0)
+//         return Scalar(0);
+      return el->data[id].value;
+    }
+
+  public:
+
+    inline void startFill(int reserveSize = 1000)
+    {
+      clear();
+      for (int i=0; i<m_data.size(); ++i)
+        m_ends[i] = m_data[i] = 0;
+    }
+
+    inline Scalar& fill(int row, int col)
+    {
+      const int outer = RowMajor ? row : col;
+      const int inner = RowMajor ? col : row;
+      if (m_ends[outer]==0)
+      {
+        m_data[outer] = m_ends[outer] = new LinkedVectorBlock();
+      }
+      else
+      {
+        ei_assert(m_ends[outer]->data[m_ends[outer]->size-1].index < inner);
+        if (m_ends[outer]->isFull())
+        {
+
+          LinkedVectorBlock* el = new LinkedVectorBlock();
+          m_ends[outer]->next = el;
+          el->prev = m_ends[outer];
+          m_ends[outer] = el;
+        }
+      }
+      m_ends[outer]->data[m_ends[outer]->size].index = inner;
+      return m_ends[outer]->data[m_ends[outer]->size++].value;
+    }
+
+    inline void endFill()
+    {
+    }
+
+    ~LinkedVectorMatrix()
+    {
+      if (this->isNotShared())
+        clear();
+    }
+
+    void clear()
+    {
+      for (int i=0; i<m_data.size(); ++i)
+      {
+        LinkedVectorBlock* el = m_data[i];
+        while (el)
+        {
+          LinkedVectorBlock* tmp = el;
+          el = el->next;
+          delete tmp;
+        }
+      }
+    }
+
+    void resize(int rows, int cols)
+    {
+      const int outers = RowMajor ? rows : cols;
+      const int inners = RowMajor ? cols : rows;
+
+      if (this->outerSize() != outers)
+      {
+        clear();
+        m_data.resize(outers);
+        m_ends.resize(outers);
+        for (int i=0; i<m_data.size(); ++i)
+          m_ends[i] = m_data[i] = 0;
+      }
+      m_innerSize = inners;
+    }
+
+    inline LinkedVectorMatrix(int rows, int cols)
+      : m_innerSize(0)
+    {
+      resize(rows, cols);
+    }
+
+    template<typename OtherDerived>
+    inline LinkedVectorMatrix(const MatrixBase<OtherDerived>& other)
+      : m_innerSize(0)
+    {
+      *this = other.derived();
+    }
+
+    inline void shallowCopy(const LinkedVectorMatrix& other)
+    {
+      EIGEN_DBG_SPARSE(std::cout << "LinkedVectorMatrix:: shallowCopy\n");
+      resize(other.rows(), other.cols());
+      for (int j=0; j<this->outerSize(); ++j)
+      {
+        m_data[j] = other.m_data[j];
+        m_ends[j] = other.m_ends[j];
+      }
+      other.markAsCopied();
+    }
+
+    inline LinkedVectorMatrix& operator=(const LinkedVectorMatrix& other)
+    {
+      if (other.isRValue())
+      {
+        shallowCopy(other);
+        return *this;
+      }
+      else
+      {
+        // TODO implement a specialized deep copy here
+        return operator=<LinkedVectorMatrix>(other);
+      }
+    }
+
+    template<typename OtherDerived>
+    inline LinkedVectorMatrix& operator=(const MatrixBase<OtherDerived>& other)
+    {
+      return SparseMatrixBase<LinkedVectorMatrix>::operator=(other.derived());
+    }
+
+  protected:
+
+    std::vector<LinkedVectorBlock*> m_data;
+    std::vector<LinkedVectorBlock*> m_ends;
+    int m_innerSize;
+
+};
+
+
+template<typename Scalar, int _Flags>
+class LinkedVectorMatrix<Scalar,_Flags>::InnerIterator
+{
+  public:
+
+    InnerIterator(const LinkedVectorMatrix& mat, int col)
+      : m_matrix(mat), m_el(mat.m_data[col]), m_it(0)
+    {}
+
+    InnerIterator& operator++() { if (m_it<m_el->size) m_it++; else {m_el = m_el->next; m_it=0;}; return *this; }
+
+    Scalar value() { return m_el->data[m_it].value; }
+
+    int index() const { return m_el->data[m_it].index; }
+
+    operator bool() const { return m_el && (m_el->next || m_it<m_el->size); }
+
+  protected:
+    const LinkedVectorMatrix& m_matrix;
+    LinkedVectorBlock* m_el;
+    int m_it;
+};
+
+#endif // EIGEN_LINKEDVECTORMATRIX_H
diff --git a/Eigen/src/Sparse/SparseArray.h b/Eigen/src/Sparse/SparseArray.h
index 2a87801a8..3ac181e04 100644
--- a/Eigen/src/Sparse/SparseArray.h
+++ b/Eigen/src/Sparse/SparseArray.h
@@ -32,11 +32,11 @@ template<typename Scalar> class SparseArray
 {
   public:
     SparseArray()
-      : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0)
+      : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0), m_isShared(0)
     {}
 
     SparseArray(int size)
-      : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0)
+      : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0), m_isShared(0)
     {
       resize(size);
     }
@@ -51,6 +51,28 @@ template<typename Scalar> class SparseArray
       resize(other.size());
       memcpy(m_values, other.m_values, m_size * sizeof(Scalar));
       memcpy(m_indices, other.m_indices, m_size * sizeof(int));
+      m_isShared = 0;
+    }
+
+    void shallowCopy(const SparseArray& other)
+    {
+      delete[] m_values;
+      delete[] m_indices;
+      m_values = other.m_values;
+      m_indices = other.m_indices;
+      m_size = other.m_size;
+      m_allocatedSize = other.m_allocatedSize;
+      m_isShared = false;
+      other.m_isShared = true;
+    }
+
+    ~SparseArray()
+    {
+      if (!m_isShared)
+      {
+        delete[] m_values;
+        delete[] m_indices;
+      }
     }
 
     void reserve(int size)
@@ -117,7 +139,11 @@ template<typename Scalar> class SparseArray
     Scalar* m_values;
     int* m_indices;
     int m_size;
-    int m_allocatedSize;
+    struct {
+      int m_allocatedSize:31;
+      mutable int m_isShared:1;
+    };
+
 };
 
 #endif // EIGEN_SPARSE_ARRAY_H
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index d3ee6e3b1..d9d404b8a 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -25,8 +25,6 @@
 #ifndef EIGEN_SPARSEMATRIX_H
 #define EIGEN_SPARSEMATRIX_H
 
-template<typename _Scalar> class SparseMatrix;
-
 /** \class SparseMatrix
   *
   * \brief Sparse matrix
@@ -36,8 +34,8 @@ template<typename _Scalar> class SparseMatrix;
   * See http://www.netlib.org/linalg/html_templates/node91.html for details on the storage scheme.
   *
   */
-template<typename _Scalar>
-struct ei_traits<SparseMatrix<_Scalar> >
+template<typename _Scalar, int _Flags>
+struct ei_traits<SparseMatrix<_Scalar, _Flags> >
 {
   typedef _Scalar Scalar;
   enum {
@@ -45,13 +43,16 @@ struct ei_traits<SparseMatrix<_Scalar> >
     ColsAtCompileTime = Dynamic,
     MaxRowsAtCompileTime = Dynamic,
     MaxColsAtCompileTime = Dynamic,
-    Flags = 0,
-    CoeffReadCost = NumTraits<Scalar>::ReadCost
+    Flags = _Flags,
+    CoeffReadCost = NumTraits<Scalar>::ReadCost,
+    SupportedAccessPatterns = FullyCoherentAccessPattern
   };
 };
 
-template<typename _Scalar>
-class SparseMatrix : public MatrixBase<SparseMatrix<_Scalar> >
+
+
+template<typename _Scalar, int _Flags>
+class SparseMatrix : public SparseMatrixBase<SparseMatrix<_Scalar, _Flags> >
 {
   public:
     EIGEN_GENERIC_PUBLIC_INTERFACE(SparseMatrix)
@@ -63,6 +64,8 @@ class SparseMatrix : public MatrixBase<SparseMatrix<_Scalar> >
     int m_rows;
     int m_cols;
 
+  public:
+
     inline int rows() const { return m_rows; }
     inline int cols() const { return m_cols; }
 
@@ -81,8 +84,8 @@ class SparseMatrix : public MatrixBase<SparseMatrix<_Scalar> >
 
     inline Scalar& coeffRef(int row, int col)
     {
-      int id = m_colPtrs[cols];
-      int end = m_colPtrs[cols+1];
+      int id = m_colPtrs[col];
+      int end = m_colPtrs[col+1];
       while (id<end && m_data.index(id)!=row)
       {
         ++id;
@@ -95,21 +98,9 @@ class SparseMatrix : public MatrixBase<SparseMatrix<_Scalar> >
 
     class InnerIterator;
 
-    inline int rows() const { return rows(); }
-    inline int cols() const { return cols(); }
     /** \returns the number of non zero coefficients */
     inline int nonZeros() const  { return m_data.size(); }
 
-    inline const Scalar& operator() (int row, int col) const
-    {
-      return coeff(row, col);
-    }
-
-    inline Scalar& operator() (int row, int col)
-    {
-      return coeffRef(row, col);
-    }
-
     inline void startFill(int reserveSize = 1000)
     {
       m_data.clear();
@@ -170,143 +161,80 @@ class SparseMatrix : public MatrixBase<SparseMatrix<_Scalar> >
       resize(rows, cols);
     }
 
-    inline SparseMatrix& operator=(const SparseMatrix& other)
+    inline void shallowCopy(const SparseMatrix& other)
     {
-      resize(other.rows(), other.cols());
+      EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: shallowCopy\n");
+      delete[] m_colPtrs;
+      m_colPtrs = 0;
+      m_rows = other.rows();
+      m_cols = other.cols();
       m_colPtrs = other.m_colPtrs;
-      for (int col=0; col<=cols(); ++col)
-        m_colPtrs[col] = other.m_colPtrs[col];
-      m_data = other.m_data;
-      return *this;
+      m_data.shallowCopy(other.m_data);
+      other.markAsCopied();
     }
 
-    template<typename OtherDerived>
-    inline SparseMatrix& operator=(const MatrixBase<OtherDerived>& other)
+    inline SparseMatrix& operator=(const SparseMatrix& other)
     {
-      resize(other.rows(), other.cols());
-      startFill(std::max(m_rows,m_cols)*2);
-      for (int col=0; col<cols(); ++col)
+      if (other.isRValue())
       {
-        for (typename OtherDerived::InnerIterator it(other.derived(), col); it; ++it)
-        {
-          Scalar v = it.value();
-          if (v!=Scalar(0))
-            fill(it.index(),col) = v;
-        }
+        shallowCopy(other);
+      }
+      else
+      {
+        resize(other.rows(), other.cols());
+        for (int col=0; col<=cols(); ++col)
+          m_colPtrs[col] = other.m_colPtrs[col];
+        m_data = other.m_data;
+        return *this;
       }
-      endFill();
-      return *this;
     }
 
+    template<typename OtherDerived>
+    inline SparseMatrix& operator=(const MatrixBase<OtherDerived>& other)
+    {
+      return SparseMatrixBase<SparseMatrix>::operator=(other);
+    }
+
+    template<typename OtherDerived>
+    SparseMatrix<Scalar> operator*(const MatrixBase<OtherDerived>& other)
+    {
+      SparseMatrix<Scalar> res(rows(), other.cols());
+      ei_sparse_product<SparseMatrix,OtherDerived>(*this,other.derived(),res);
+      return res;
+    }
 
-    // old explicit operator+
-//     template<typename Other>
-//     SparseMatrix operator+(const Other& other)
-//     {
-//       SparseMatrix res(rows(), cols());
-//       res.startFill(nonZeros()*3);
-//       for (int col=0; col<cols(); ++col)
-//       {
-//         InnerIterator row0(*this,col);
-//         typename Other::InnerIterator row1(other,col);
-//         while (row0 && row1)
-//         {
-//           if (row0.index()==row1.index())
-//           {
-//             std::cout << "both " << col << " " << row0.index() << "\n";
-//             Scalar v = row0.value() + row1.value();
-//             if (v!=Scalar(0))
-//               res.fill(row0.index(),col) = v;
-//             ++row0;
-//             ++row1;
-//           }
-//           else if (row0.index()<row1.index())
-//           {
-//             std::cout << "row0 " << col << " " << row0.index() << "\n";
-//             Scalar v = row0.value();
-//             if (v!=Scalar(0))
-//               res.fill(row0.index(),col) = v;
-//             ++row0;
-//           }
-//           else if (row1)
-//           {
-//             std::cout << "row1 " << col << " " << row0.index() << "\n";
-//             Scalar v = row1.value();
-//             if (v!=Scalar(0))
-//               res.fill(row1.index(),col) = v;
-//             ++row1;
-//           }
-//         }
-//         while (row0)
-//         {
-//           std::cout << "row0 " << col << " " << row0.index() << "\n";
-//           Scalar v = row0.value();
-//           if (v!=Scalar(0))
-//             res.fill(row0.index(),col) = v;
-//           ++row0;
-//         }
-//         while (row1)
-//         {
-//           std::cout << "row1 " << col << " " << row1.index() << "\n";
-//           Scalar v = row1.value();
-//           if (v!=Scalar(0))
-//             res.fill(row1.index(),col) = v;
-//           ++row1;
-//         }
-//       }
-//       res.endFill();
-//       return res;
-// //       return binaryOp(other, ei_scalar_sum_op<Scalar>());
-//     }
-
-
-    // WARNING for efficiency reason it currently outputs the transposed matrix
     friend std::ostream & operator << (std::ostream & s, const SparseMatrix& m)
     {
-      s << "Nonzero entries:\n";
-      for (uint i=0; i<m.nonZeros(); ++i)
-      {
-        s << "(" << m.m_data.value(i) << "," << m.m_data.index(i) << ") ";
-      }
-      s << std::endl;
-      s << std::endl;
-      s << "Column pointers:\n";
-      for (uint i=0; i<m.cols(); ++i)
-      {
-        s << m.m_colPtrs[i] << " ";
-      }
-      s << std::endl;
-      s << std::endl;
-      s << "Matrix (transposed):\n";
-      for (int j=0; j<m.cols(); j++ )
-      {
-        int end = m.m_colPtrs[j+1];
-        int i=0;
-        for (int id=m.m_colPtrs[j]; id<end; id++)
+      EIGEN_DBG_SPARSE(
+        s << "Nonzero entries:\n";
+        for (uint i=0; i<m.nonZeros(); ++i)
         {
-          int row = m.m_data.index(id);
-          // fill with zeros
-          for (int k=i; k<row; ++k)
-            s << "0 ";
-          i = row+1;
-          s << m.m_data.value(id) << " ";
+          s << "(" << m.m_data.value(i) << "," << m.m_data.index(i) << ") ";
         }
-        for (int k=i; k<m.rows(); ++k)
-          s << "0 ";
         s << std::endl;
-      }
+        s << std::endl;
+        s << "Column pointers:\n";
+        for (uint i=0; i<m.cols(); ++i)
+        {
+          s << m.m_colPtrs[i] << " ";
+        }
+        s << std::endl;
+        s << std::endl;
+      );
+      s << static_cast<const SparseMatrixBase<SparseMatrix>&>(m);
       return s;
     }
 
     /** Destructor */
     inline ~SparseMatrix()
     {
-      delete[] m_colPtrs;
+      if (this->isNotShared())
+        delete[] m_colPtrs;
     }
 };
 
-template<typename Scalar>
-class SparseMatrix<Scalar>::InnerIterator
+template<typename Scalar, int _Flags>
+class SparseMatrix<Scalar,_Flags>::InnerIterator
 {
   public:
     InnerIterator(const SparseMatrix& mat, int col)
diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h
new file mode 100644
index 000000000..6357cbeff
--- /dev/null
+++ b/Eigen/src/Sparse/SparseMatrixBase.h
@@ -0,0 +1,163 @@
+// 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/>.
+
+#ifndef EIGEN_SPARSEMATRIXBASE_H
+#define EIGEN_SPARSEMATRIXBASE_H
+
+template<typename Derived>
+class SparseMatrixBase : public MatrixBase<Derived>
+{
+  public:
+
+    typedef MatrixBase<Derived> Base;
+    typedef typename Base::Scalar Scalar;
+    enum {
+      Flags = Base::Flags,
+      RowMajor = ei_traits<Derived>::Flags&RowMajorBit ? 1 : 0
+    };
+
+    inline const Derived& derived() const { return *static_cast<const Derived*>(this); }
+    inline Derived& derived() { return *static_cast<Derived*>(this); }
+    inline Derived& const_cast_derived() const
+    { return *static_cast<Derived*>(const_cast<SparseMatrixBase*>(this)); }
+
+    SparseMatrixBase()
+      : m_isRValue(false), m_hasBeenCopied(false)
+    {}
+
+    bool isRValue() const { return m_isRValue; }
+    Derived& temporary() { m_isRValue = true; return derived(); }
+
+    int outerSize() const { return RowMajor ? this->rows() : this->cols(); }
+    int innerSize() const { return RowMajor ? this->cols() : this->rows(); }
+
+    inline Derived& operator=(const Derived& other)
+    {
+      if (other.isRValue())
+      {
+        m_hasBeenCopied = true;
+        derived().shallowCopy(other);
+      }
+      else
+        this->operator=<Derived>(other);
+      return derived();
+    }
+
+
+
+    template<typename OtherDerived>
+    inline Derived& operator=(const MatrixBase<OtherDerived>& other)
+    {
+      const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
+      const int outerSize = (int(OtherDerived::Flags) & RowMajorBit) ? other.rows() : other.cols();
+      // eval to a temporary and then do a shallow copy
+      typename ei_meta_if<transpose, LinkedVectorMatrix<Scalar,Flags&RowMajorBit>, Derived>::ret temp(other.rows(), other.cols());
+
+      temp.startFill(std::max(this->rows(),this->cols())*2);
+//       std::cout << other.rows() << " xm " << other.cols() << "\n";
+      for (int j=0; j<outerSize; ++j)
+      {
+        for (typename OtherDerived::InnerIterator it(other.derived(), j); it; ++it)
+        {
+//           std::cout << other.rows() << " x " << other.cols() << "\n";
+//           std::cout << it.m_matrix.rows() << "\n";
+          Scalar v = it.value();
+          if (v!=Scalar(0))
+            if (RowMajor) temp.fill(j,it.index()) = v;
+            else temp.fill(it.index(),j) = v;
+        }
+      }
+      temp.endFill();
+      derived() = temp.temporary();
+      return derived();
+    }
+
+    template<typename OtherDerived>
+    inline Derived& operator=(const SparseMatrixBase<OtherDerived>& other)
+    {
+      const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
+      const int outerSize = (int(OtherDerived::Flags) & RowMajorBit) ? other.rows() : other.cols();
+      if ((!transpose) && other.isRValue())
+      {
+        // eval without temporary
+        derived().resize(other.rows(), other.cols());
+        derived().startFill(std::max(this->rows(),this->cols())*2);
+        for (int j=0; j<outerSize; ++j)
+        {
+          for (typename OtherDerived::InnerIterator it(other.derived(), j); it; ++it)
+          {
+            Scalar v = it.value();
+            if (v!=Scalar(0))
+              if (RowMajor) derived().fill(j,it.index()) = v;
+              else derived().fill(it.index(),j) = v;
+          }
+        }
+        derived().endFill();
+      }
+      else
+        this->operator=<OtherDerived>(static_cast<const MatrixBase<OtherDerived>&>(other));
+      return derived();
+    }
+
+    friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m)
+    {
+
+      if (Flags&RowMajorBit)
+      {
+        for (int row=0; row<m.rows(); ++row)
+        {
+          int col = 0;
+          for (typename Derived::InnerIterator it(m.derived(), row); it; ++it)
+          {
+            for ( ; col<it.index(); ++col)
+              s << "0 ";
+            std::cout << it.value() << " ";
+          }
+          for ( ; col<m.cols(); ++col)
+            s << "0 ";
+          s << std::endl;
+        }
+      }
+      else
+      {
+        LinkedVectorMatrix<Scalar, RowMajorBit> trans = m.derived();
+        s << trans;
+      }
+      return s;
+    }
+
+  protected:
+
+    bool isNotShared() { return !m_hasBeenCopied; }
+    void markAsCopied() const { m_hasBeenCopied = true; }
+
+  protected:
+
+    bool m_isRValue;
+    mutable bool m_hasBeenCopied;
+};
+
+
+
+#endif // EIGEN_SPARSEMATRIXBASE_H
diff --git a/Eigen/src/Sparse/SparseProduct.h b/Eigen/src/Sparse/SparseProduct.h
new file mode 100644
index 000000000..9abddbd27
--- /dev/null
+++ b/Eigen/src/Sparse/SparseProduct.h
@@ -0,0 +1,169 @@
+// 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/>.
+
+#ifndef EIGEN_SPARSEPRODUCT_H
+#define EIGEN_SPARSEPRODUCT_H
+
+#define DENSE_TMP 1
+#define MAP_TMP 2
+#define LIST_TMP 3
+
+#define TMP_TMP 3
+
+template<typename Scalar>
+struct ListEl
+{
+  int next;
+  int index;
+  Scalar value;
+};
+
+template<typename Lhs, typename Rhs>
+static void ei_sparse_product(const Lhs& lhs, const Rhs& rhs, SparseMatrix<typename ei_traits<Lhs>::Scalar>& res)
+{
+  int rows = lhs.rows();
+  int cols = rhs.rows();
+  int size = lhs.cols();
+
+  float ratio = std::max(float(lhs.nonZeros())/float(lhs.rows()*lhs.cols()), float(rhs.nonZeros())/float(rhs.rows()*rhs.cols()));
+  std::cout << ratio << "\n";
+
+  ei_assert(size == rhs.rows());
+  typedef typename ei_traits<Lhs>::Scalar Scalar;
+  #if (TMP_TMP == MAP_TMP)
+  std::map<int,Scalar> tmp;
+  #elif (TMP_TMP == LIST_TMP)
+  std::vector<ListEl<Scalar> > tmp(2*rows);
+  #else
+  std::vector<Scalar> tmp(rows);
+  #endif
+  res.resize(rows, cols);
+  res.startFill(2*std::max(rows, cols));
+  for (int j=0; j<cols; ++j)
+  {
+    #if (TMP_TMP == MAP_TMP)
+    tmp.clear();
+    #elif (TMP_TMP == LIST_TMP)
+    int tmp_size = 0;
+    int tmp_start = -1;
+    #else
+    for (int k=0; k<rows; ++k)
+      tmp[k] = 0;
+    #endif
+    for (typename Rhs::InnerIterator rhsIt(rhs, j); rhsIt; ++rhsIt)
+    {
+      #if (TMP_TMP == MAP_TMP)
+      typename std::map<int,Scalar>::iterator hint = tmp.begin();
+      typename std::map<int,Scalar>::iterator r;
+      #elif (TMP_TMP == LIST_TMP)
+      int tmp_el = tmp_start;
+      #endif
+      for (typename Lhs::InnerIterator lhsIt(lhs, rhsIt.index()); lhsIt; ++lhsIt)
+      {
+        #if (TMP_TMP == MAP_TMP)
+        r = hint;
+        Scalar v = lhsIt.value() * rhsIt.value();
+        int id = lhsIt.index();
+        while (r!=tmp.end() && r->first < id)
+          ++r;
+        if (r!=tmp.end() && r->first==id)
+        {
+          r->second += v;
+          hint = r;
+        }
+        else
+          hint = tmp.insert(r, std::pair<int,Scalar>(id, v));
+        ++hint;
+        #elif (TMP_TMP == LIST_TMP)
+        Scalar v = lhsIt.value() * rhsIt.value();
+        int id = lhsIt.index();
+        if (tmp_size==0)
+        {
+          tmp_start = 0;
+          tmp_el = 0;
+          tmp_size++;
+          tmp[0].value = v;
+          tmp[0].index = id;
+          tmp[0].next = -1;
+        }
+        else if (id<tmp[tmp_start].index)
+        {
+          tmp[tmp_size].value = v;
+          tmp[tmp_size].index = id;
+          tmp[tmp_size].next = tmp_start;
+          tmp_start = tmp_size;
+          tmp_size++;
+        }
+        else
+        {
+          int nextel = tmp[tmp_el].next;
+          while (nextel >= 0 && tmp[nextel].index<=id)
+          {
+            tmp_el = nextel;
+            nextel = tmp[nextel].next;
+          }
+
+          if (tmp[tmp_el].index==id)
+          {
+            tmp[tmp_el].value += v;
+          }
+          else
+          {
+            tmp[tmp_size].value = v;
+            tmp[tmp_size].index = id;
+            tmp[tmp_size].next = tmp[tmp_el].next;
+            tmp[tmp_el].next = tmp_size;
+            tmp_size++;
+          }
+        }
+        #else
+        tmp[lhsIt.index()] += lhsIt.value() * rhsIt.value();
+        #endif
+        //res.coeffRef(lhsIt.index(), j) += lhsIt.value() * rhsIt.value();
+      }
+    }
+    #if (TMP_TMP == MAP_TMP)
+    for (typename std::map<int,Scalar>::const_iterator k=tmp.begin(); k!=tmp.end(); ++k)
+      if (k->second!=0)
+        res.fill(k->first, j) = k->second;
+    #elif (TMP_TMP == LIST_TMP)
+    int k = tmp_start;
+    while (k>=0)
+    {
+      if (tmp[k].value!=0)
+        res.fill(tmp[k].index, j) = tmp[k].value;
+      k = tmp[k].next;
+    }
+    #else
+    for (int k=0; k<rows; ++k)
+      if (tmp[k]!=0)
+        res.fill(k, j) = tmp[k];
+    #endif
+  }
+  res.endFill();
+
+  std::cout << "  => " << float(res.nonZeros())/float(res.rows()*res.cols()) << "\n";
+}
+
+#endif // EIGEN_SPARSEPRODUCT_H
diff --git a/Eigen/src/Sparse/SparseSetter.h b/Eigen/src/Sparse/SparseSetter.h
new file mode 100644
index 000000000..93b634322
--- /dev/null
+++ b/Eigen/src/Sparse/SparseSetter.h
@@ -0,0 +1,102 @@
+// 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/>.
+
+#ifndef EIGEN_SPARSESETTER_H
+#define EIGEN_SPARSESETTER_H
+
+template<typename MatrixType, int AccessPattern,
+  int IsSupported = ei_support_access_pattern<MatrixType,AccessPattern>::ret>
+struct ei_sparse_setter_selector;
+
+template<typename MatrixType,
+         int AccessPattern,
+         typename WrapperType = typename ei_sparse_setter_selector<MatrixType,AccessPattern>::type>
+class SparseSetter
+{
+    typedef typename ei_unref<WrapperType>::type _WrapperType;
+  public:
+
+    inline SparseSetter(MatrixType& matrix) : m_wrapper(matrix), mp_matrix(&matrix) {}
+
+    ~SparseSetter()
+    { *mp_matrix = m_wrapper; }
+
+    inline _WrapperType* operator->() { return &m_wrapper; }
+
+    inline _WrapperType& operator*() { return m_wrapper; }
+
+  protected:
+
+    WrapperType m_wrapper;
+    MatrixType* mp_matrix;
+};
+
+template<typename MatrixType, int AccessPattern>
+struct ei_sparse_setter_selector<MatrixType, AccessPattern, AccessPatternSupported>
+{
+  typedef MatrixType& type;
+};
+
+// forward each derived of SparseMatrixBase to the generic SparseMatrixBase specializations
+template<typename Scalar, int Flags, int AccessPattern>
+struct ei_sparse_setter_selector<SparseMatrix<Scalar,Flags>, AccessPattern, AccessPatternNotSupported>
+: public ei_sparse_setter_selector<SparseMatrixBase<SparseMatrix<Scalar,Flags> >,AccessPattern, AccessPatternNotSupported>
+{};
+
+template<typename Scalar, int Flags, int AccessPattern>
+struct ei_sparse_setter_selector<LinkedVectorMatrix<Scalar,Flags>, AccessPattern, AccessPatternNotSupported>
+: public ei_sparse_setter_selector<LinkedVectorMatrix<SparseMatrix<Scalar,Flags> >,AccessPattern, AccessPatternNotSupported>
+{};
+
+template<typename Scalar, int Flags, int AccessPattern>
+struct ei_sparse_setter_selector<HashMatrix<Scalar,Flags>, AccessPattern, AccessPatternNotSupported>
+: public ei_sparse_setter_selector<HashMatrix<SparseMatrix<Scalar,Flags> >,AccessPattern, AccessPatternNotSupported>
+{};
+
+// generic SparseMatrixBase specializations
+template<typename Derived>
+struct ei_sparse_setter_selector<SparseMatrixBase<Derived>, RandomAccessPattern, AccessPatternNotSupported>
+{
+  typedef HashMatrix<typename Derived::Scalar, Derived::Flags> type;
+};
+
+template<typename Derived>
+struct ei_sparse_setter_selector<SparseMatrixBase<Derived>, OuterCoherentAccessPattern, AccessPatternNotSupported>
+{
+  typedef HashMatrix<typename Derived::Scalar, Derived::Flags> type;
+};
+
+template<typename Derived>
+struct ei_sparse_setter_selector<SparseMatrixBase<Derived>, InnerCoherentAccessPattern, AccessPatternNotSupported>
+{
+  typedef LinkedVectorMatrix<typename Derived::Scalar, Derived::Flags> type;
+};
+
+template<typename Derived>
+struct ei_sparse_setter_selector<SparseMatrixBase<Derived>, FullyCoherentAccessPattern, AccessPatternNotSupported>
+{
+  typedef SparseMatrix<typename Derived::Scalar, Derived::Flags> type;
+};
+
+#endif // EIGEN_SPARSESETTER_H
diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h
new file mode 100644
index 000000000..e89a7c322
--- /dev/null
+++ b/Eigen/src/Sparse/SparseUtil.h
@@ -0,0 +1,61 @@
+// 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/>.
+
+#ifndef EIGEN_SPARSEUTIL_H
+#define EIGEN_SPARSEUTIL_H
+
+#ifdef NDEBUG
+#define EIGEN_DBG_SPARSE(X)
+#else
+#define EIGEN_DBG_SPARSE(X) X
+#endif
+
+template<typename _Scalar, int _Flags = 0> class SparseMatrix;
+template<typename _Scalar, int _Flags = 0> class HashMatrix;
+template<typename _Scalar, int _Flags = 0> class LinkedVectorMatrix;
+
+const int AccessPatternNotSupported = 0x0;
+const int AccessPatternSupported    = 0x1;
+
+
+template<typename MatrixType, int AccessPattern> struct ei_support_access_pattern
+{
+  enum { ret = (int(ei_traits<MatrixType>::SupportedAccessPatterns) & AccessPattern) == AccessPattern
+             ? AccessPatternSupported
+             : AccessPatternNotSupported
+  };
+};
+
+template<typename T> class ei_eval<T,Sparse>
+{
+    typedef typename ei_traits<T>::Scalar _Scalar;
+    enum {
+          _Flags = ei_traits<T>::Flags
+    };
+
+  public:
+    typedef SparseMatrix<_Scalar, _Flags> type;
+};
+
+#endif // EIGEN_SPARSEUTIL_H