Big rewrite in the Sparse module: SparseMatrixBase no longer inherits MatrixBase.

That means a lot of features which were available for sparse matrices
via the dense (and super slow) implemention are no longer available.
All features which make sense for sparse matrices (aka can be implemented efficiently) will be
implemented soon, but don't expect to see an API as rich as for the dense path.
Other changes:
* no block(), row(), col() anymore.
* instead use .innerVector() to get a col or row vector of a matrix.
* .segment(), start(), end() will be back soon, not sure for block()
* faster cwise product

31 files changed, 1919 insertions, 199 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index dd74b7674..f304401cb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -76,6 +76,12 @@ if(MSVC)
   endif(EIGEN_TEST_SSE2)
 endif(MSVC)
 
+option(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION "Disable explicit vectorization in tests/examples" OFF)
+if(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
+  add_definitions(EIGEN_DONT_VECTORIZE=1)
+  message("Disabling vectorization in tests/examples")
+endif(EIGEN_TEST_NO_EXPLICIT_VECTORIZATION)
+
 include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR})
 
 add_subdirectory(Eigen)
diff --git a/Eigen/Sparse b/Eigen/Sparse
index 168511c66..b73946c74 100644
--- a/Eigen/Sparse
+++ b/Eigen/Sparse
@@ -79,7 +79,15 @@ namespace Eigen {
 #include "src/Sparse/SparseMatrix.h"
 #include "src/Sparse/SparseVector.h"
 #include "src/Sparse/CoreIterators.h"
+#include "src/Sparse/SparseTranspose.h"
+#include "src/Sparse/SparseCwise.h"
+#include "src/Sparse/SparseCwiseUnaryOp.h"
+#include "src/Sparse/SparseCwiseBinaryOp.h"
+#include "src/Sparse/SparseDot.h"
+#include "src/Sparse/SparseAssign.h"
 #include "src/Sparse/SparseRedux.h"
+#include "src/Sparse/SparseFuzzy.h"
+#include "src/Sparse/SparseFlagged.h"
 #include "src/Sparse/SparseProduct.h"
 #include "src/Sparse/TriangularSolver.h"
 #include "src/Sparse/SparseLLT.h"
diff --git a/Eigen/src/Core/CwiseBinaryOp.h b/Eigen/src/Core/CwiseBinaryOp.h
index c5bead928..c4223e220 100644
--- a/Eigen/src/Core/CwiseBinaryOp.h
+++ b/Eigen/src/Core/CwiseBinaryOp.h
@@ -86,14 +86,12 @@ class CwiseBinaryOp : ei_no_assignment_operator,
     typedef typename ei_traits<CwiseBinaryOp>::LhsNested LhsNested;
     typedef typename ei_traits<CwiseBinaryOp>::RhsNested RhsNested;
 
-    class InnerIterator;
-
     EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp())
       : m_lhs(lhs), m_rhs(rhs), m_functor(func)
     {
       // we require Lhs and Rhs to have the same scalar type. Currently there is no example of a binary functor
       // that would take two operands of different types. If there were such an example, then this check should be
-      // moved to the BinaryOp functors, on a per-case basis. This would however require a change in the BinaryOp functors, as 
+      // moved to the BinaryOp functors, on a per-case basis. This would however require a change in the BinaryOp functors, as
       // currently they take only one typename Scalar template parameter.
       // It is tempting to always allow mixing different types but remember that this is often impossible in the vectorized paths.
       // So allowing mixing different types gives very unexpected errors when enabling vectorization, when the user tries to
diff --git a/Eigen/src/Core/CwiseUnaryOp.h b/Eigen/src/Core/CwiseUnaryOp.h
index 68be2ac54..076d568e0 100644
--- a/Eigen/src/Core/CwiseUnaryOp.h
+++ b/Eigen/src/Core/CwiseUnaryOp.h
@@ -64,8 +64,6 @@ class CwiseUnaryOp : ei_no_assignment_operator,
 
     EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseUnaryOp)
 
-    class InnerIterator;
-
     inline CwiseUnaryOp(const MatrixType& mat, const UnaryOp& func = UnaryOp())
       : m_matrix(mat), m_functor(func) {}
 
diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h
index 71e7030dc..b66fcbaae 100644
--- a/Eigen/src/Core/Dot.h
+++ b/Eigen/src/Core/Dot.h
@@ -143,13 +143,12 @@ struct ei_dot_vec_unroller<Derived1, Derived2, Index, Stop, true>
 
 template<typename Derived1, typename Derived2,
          int Vectorization = ei_dot_traits<Derived1, Derived2>::Vectorization,
-         int Unrolling = ei_dot_traits<Derived1, Derived2>::Unrolling,
-         int Storage = (ei_traits<Derived1>::Flags | ei_traits<Derived2>::Flags) & SparseBit
+         int Unrolling = ei_dot_traits<Derived1, Derived2>::Unrolling
 >
 struct ei_dot_impl;
 
 template<typename Derived1, typename Derived2>
-struct ei_dot_impl<Derived1, Derived2, NoVectorization, NoUnrolling, IsDense>
+struct ei_dot_impl<Derived1, Derived2, NoVectorization, NoUnrolling>
 {
   typedef typename Derived1::Scalar Scalar;
   static Scalar run(const Derived1& v1, const Derived2& v2)
@@ -164,12 +163,12 @@ struct ei_dot_impl<Derived1, Derived2, NoVectorization, NoUnrolling, IsDense>
 };
 
 template<typename Derived1, typename Derived2>
-struct ei_dot_impl<Derived1, Derived2, NoVectorization, CompleteUnrolling, IsDense>
+struct ei_dot_impl<Derived1, Derived2, NoVectorization, CompleteUnrolling>
   : public ei_dot_novec_unroller<Derived1, Derived2, 0, Derived1::SizeAtCompileTime>
 {};
 
 template<typename Derived1, typename Derived2>
-struct ei_dot_impl<Derived1, Derived2, LinearVectorization, NoUnrolling, IsDense>
+struct ei_dot_impl<Derived1, Derived2, LinearVectorization, NoUnrolling>
 {
   typedef typename Derived1::Scalar Scalar;
   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
@@ -222,7 +221,7 @@ struct ei_dot_impl<Derived1, Derived2, LinearVectorization, NoUnrolling, IsDense
 };
 
 template<typename Derived1, typename Derived2>
-struct ei_dot_impl<Derived1, Derived2, LinearVectorization, CompleteUnrolling, IsDense>
+struct ei_dot_impl<Derived1, Derived2, LinearVectorization, CompleteUnrolling>
 {
   typedef typename Derived1::Scalar Scalar;
   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 42d3bece3..eecd24c85 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -251,8 +251,8 @@ template<typename Derived> class MatrixBase
     { return lazyAssign(other._expression()); }
 
     /** Overloaded for sparse product evaluation */
-    template<typename Derived1, typename Derived2>
-    Derived& lazyAssign(const Product<Derived1,Derived2,SparseProduct>& product);
+    /*template<typename Derived1, typename Derived2>
+    Derived& lazyAssign(const Product<Derived1,Derived2,SparseProduct>& product);*/
 
     CommaInitializer<Derived> operator<< (const Scalar& s);
 
diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index 2fd9d5875..305efc3dc 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -79,7 +79,6 @@ struct ProductReturnType<Lhs,Rhs,CacheFriendlyProduct>
  *    - NormalProduct
  *    - CacheFriendlyProduct
  *    - DiagonalProduct
- *    - SparseProduct
  */
 template<typename Lhs, typename Rhs> struct ei_product_mode
 {
@@ -87,8 +86,6 @@ template<typename Lhs, typename Rhs> struct ei_product_mode
 
     value = ((Rhs::Flags&Diagonal)==Diagonal) || ((Lhs::Flags&Diagonal)==Diagonal)
           ? DiagonalProduct
-          : (Rhs::Flags & Lhs::Flags & SparseBit)
-          ? SparseProduct
           : Lhs::MaxColsAtCompileTime == Dynamic
             && ( Lhs::MaxRowsAtCompileTime == Dynamic
               || Rhs::MaxColsAtCompileTime == Dynamic )
diff --git a/Eigen/src/Core/Sum.h b/Eigen/src/Core/Sum.h
index ab9cdad1f..e30392534 100644
--- a/Eigen/src/Core/Sum.h
+++ b/Eigen/src/Core/Sum.h
@@ -154,13 +154,12 @@ struct ei_sum_vec_unroller<Derived, Index, Stop, true>
 
 template<typename Derived,
          int Vectorization = ei_sum_traits<Derived>::Vectorization,
-         int Unrolling = ei_sum_traits<Derived>::Unrolling,
-         int Storage = ei_traits<Derived>::Flags & SparseBit
+         int Unrolling = ei_sum_traits<Derived>::Unrolling
 >
 struct ei_sum_impl;
 
 template<typename Derived>
-struct ei_sum_impl<Derived, NoVectorization, NoUnrolling, IsDense>
+struct ei_sum_impl<Derived, NoVectorization, NoUnrolling>
 {
   typedef typename Derived::Scalar Scalar;
   static Scalar run(const Derived& mat)
@@ -178,12 +177,12 @@ struct ei_sum_impl<Derived, NoVectorization, NoUnrolling, IsDense>
 };
 
 template<typename Derived>
-struct ei_sum_impl<Derived, NoVectorization, CompleteUnrolling, IsDense>
+struct ei_sum_impl<Derived, NoVectorization, CompleteUnrolling>
   : public ei_sum_novec_unroller<Derived, 0, Derived::SizeAtCompileTime>
 {};
 
 template<typename Derived>
-struct ei_sum_impl<Derived, LinearVectorization, NoUnrolling,IsDense>
+struct ei_sum_impl<Derived, LinearVectorization, NoUnrolling>
 {
   typedef typename Derived::Scalar Scalar;
   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
@@ -228,7 +227,7 @@ struct ei_sum_impl<Derived, LinearVectorization, NoUnrolling,IsDense>
 };
 
 template<typename Derived>
-struct ei_sum_impl<Derived, LinearVectorization, CompleteUnrolling, IsDense>
+struct ei_sum_impl<Derived, LinearVectorization, CompleteUnrolling>
 {
   typedef typename Derived::Scalar Scalar;
   typedef typename ei_packet_traits<Scalar>::type PacketScalar;
diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h
index 1137e42dd..93028ed49 100644
--- a/Eigen/src/Core/Transpose.h
+++ b/Eigen/src/Core/Transpose.h
@@ -63,8 +63,6 @@ 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)
@@ -185,7 +183,7 @@ struct ei_inplace_transpose_selector<MatrixType,false> { // non square matrix
   *
   * In most cases it is probably better to simply use the transposed expression
   * of a matrix. However, when transposing the matrix data itself is really needed,
-  * then this "in-place" version is probably the right choice because it provides 
+  * then this "in-place" version is probably the right choice because it provides
   * the following additional features:
   *  - less error prone: doing the same operation with .transpose() requires special care:
   *    \code m.set(m.transpose().eval()); \endcode
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h
index 10e50f43b..f2c76cc01 100644
--- a/Eigen/src/Core/util/Constants.h
+++ b/Eigen/src/Core/util/Constants.h
@@ -201,7 +201,7 @@ enum { ForceAligned, AsRequested };
 enum { ConditionalJumpCost = 5 };
 enum CornerType { TopLeft, TopRight, BottomLeft, BottomRight };
 enum DirectionType { Vertical, Horizontal };
-enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct, SparseProduct };
+enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct, DiagonalProduct };
 
 enum {
   /** \internal Equivalent to a slice vectorization for fixed-size matrices having good alignment
diff --git a/Eigen/src/Sparse/CoreIterators.h b/Eigen/src/Sparse/CoreIterators.h
index b1d113a95..d57eac287 100644
--- a/Eigen/src/Sparse/CoreIterators.h
+++ b/Eigen/src/Sparse/CoreIterators.h
@@ -28,19 +28,26 @@
 /* This file contains the respective InnerIterator definition of the expressions defined in Eigen/Core
  */
 
-template<typename Derived>
-class MatrixBase<Derived>::InnerIterator
+/** \class InnerIterator
+  * \brief An InnerIterator allows to loop over the element of a sparse (or dense) matrix or expression
+  *
+  * todo
+  */
+// template<typename XprType> InnerIterator;
+
+// generic version for dense matrix and expressions
+template<typename Derived> class MatrixBase<Derived>::InnerIterator
 {
     typedef typename Derived::Scalar Scalar;
   public:
-    EIGEN_STRONG_INLINE InnerIterator(const Derived& mat, int outer)
-      : m_matrix(mat), m_inner(0), m_outer(outer), m_end(mat.rows())
+    EIGEN_STRONG_INLINE InnerIterator(const Derived& expr, int outer)
+      : m_expression(expr), m_inner(0), m_outer(outer), m_end(expr.rows())
     {}
 
     EIGEN_STRONG_INLINE Scalar value() const
     {
-      return (Derived::Flags&RowMajorBit) ? m_matrix.coeff(m_outer, m_inner)
-                                          : m_matrix.coeff(m_inner, m_outer);
+      return (Derived::Flags&RowMajorBit) ? m_expression.coeff(m_outer, m_inner)
+                                          : m_expression.coeff(m_inner, m_outer);
     }
 
     EIGEN_STRONG_INLINE InnerIterator& operator++() { m_inner++; return *this; }
@@ -50,22 +57,14 @@ class MatrixBase<Derived>::InnerIterator
     EIGEN_STRONG_INLINE operator bool() const { return m_inner < m_end && m_inner>=0; }
 
   protected:
-    const Derived& m_matrix;
+    const Derived& m_expression;
     int m_inner;
     const int m_outer;
     const int m_end;
 };
 
-template<typename MatrixType>
-class Transpose<MatrixType>::InnerIterator : public MatrixType::InnerIterator
-{
-  public:
-
-    EIGEN_STRONG_INLINE InnerIterator(const Transpose& trans, int outer)
-      : MatrixType::InnerIterator(trans.m_matrix, outer)
-    {}
-};
 
+/*
 template<typename MatrixType, int BlockRows, int BlockCols, int PacketAccess, int _DirectAccessStatus>
 class Block<MatrixType, BlockRows, BlockCols, PacketAccess, _DirectAccessStatus>::InnerIterator
 {
@@ -138,50 +137,9 @@ class Block<MatrixType, BlockRows, BlockCols, PacketAccess, IsSparse>::InnerIter
     int m_start;
     int m_end;
     int m_offset;
-};
-
-template<typename UnaryOp, typename MatrixType>
-class CwiseUnaryOp<UnaryOp,MatrixType>::InnerIterator
-{
-    typedef typename CwiseUnaryOp::Scalar Scalar;
-    typedef typename ei_traits<CwiseUnaryOp>::_MatrixTypeNested _MatrixTypeNested;
-    typedef typename _MatrixTypeNested::InnerIterator MatrixTypeIterator;
-  public:
-
-    EIGEN_STRONG_INLINE InnerIterator(const CwiseUnaryOp& unaryOp, int outer)
-      : m_iter(unaryOp.m_matrix,outer), m_functor(unaryOp.m_functor), m_id(-1)
-    {
-      this->operator++();
-    }
-
-    EIGEN_STRONG_INLINE InnerIterator& operator++()
-    {
-      if (m_iter)
-      {
-        m_id = m_iter.index();
-        m_value = m_functor(m_iter.value());
-        ++m_iter;
-      }
-      else
-      {
-        m_id = -1;
-      }
-      return *this;
-    }
-
-    EIGEN_STRONG_INLINE Scalar value() const { return m_value; }
-
-    EIGEN_STRONG_INLINE int index() const { return m_id; }
-
-    EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; }
-
-  protected:
-    MatrixTypeIterator m_iter;
-    const UnaryOp& m_functor;
-    Scalar m_value;
-    int m_id;
-};
+};*/
 
+/*
 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 }; };
 
@@ -204,8 +162,8 @@ class CwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator
     EIGEN_STRONG_INLINE InnerIterator(const CwiseBinaryOp& binOp, int outer)
       : Base(binOp.m_lhs,binOp.m_rhs,binOp.m_functor,outer)
     {}
-};
-
+};*/
+/*
 template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived>
 class CwiseBinaryOpInnerIterator
 {
@@ -270,7 +228,7 @@ class CwiseBinaryOpInnerIterator
     Scalar m_value;
     int m_id;
 };
-/*
+
 template<typename T, typename Lhs, typename Rhs, typename Derived>
 class CwiseBinaryOpInnerIterator<ei_scalar_product_op<T>,Lhs,Rhs,Derived>
 {
diff --git a/Eigen/src/Sparse/SparseAssign.h b/Eigen/src/Sparse/SparseAssign.h
new file mode 100644
index 000000000..e69de29bb
--- /dev/null
+++ b/Eigen/src/Sparse/SparseAssign.h
diff --git a/Eigen/src/Sparse/SparseBlock.h b/Eigen/src/Sparse/SparseBlock.h
index 55590d9e2..48be3754f 100644
--- a/Eigen/src/Sparse/SparseBlock.h
+++ b/Eigen/src/Sparse/SparseBlock.h
@@ -23,9 +23,76 @@
 // License and a copy of the GNU General Public License along with
 // Eigen. If not, see <http://www.gnu.org/licenses/>.
 
-#ifndef EIGEN_SPARSEBLOCK_H
-#define EIGEN_SPARSEBLOCK_H
+#ifndef EIGEN_SPARSE_BLOCK_H
+#define EIGEN_SPARSE_BLOCK_H
 
+template<typename MatrixType>
+struct ei_traits<SparseInnerVector<MatrixType> >
+{
+  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,
+    CoeffReadCost = MatrixType::CoeffReadCost
+  };
+};
+
+template<typename MatrixType>
+class SparseInnerVector : ei_no_assignment_operator,
+  public SparseMatrixBase<SparseInnerVector<MatrixType> >
+{
+    enum {
+      IsRowMajor = ei_traits<SparseInnerVector>::IsRowMajor
+    };
+public:
+
+    EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseInnerVector)
+    class InnerIterator;
+
+    inline SparseInnerVector(const MatrixType& matrix, int outer)
+      : m_matrix(matrix), m_outer(outer)
+    {
+      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; }
+
+  protected:
+
+    const typename MatrixType::Nested m_matrix;
+    int m_outer;
+
+};
+
+template<typename MatrixType>
+class SparseInnerVector<MatrixType>::InnerIterator : public MatrixType::InnerIterator
+{
+public:
+  inline InnerIterator(const SparseInnerVector& xpr, int outer=0)
+    : MatrixType::InnerIterator(xpr.m_matrix, xpr.m_outer)
+  {
+    ei_assert(outer==0);
+  }
+};
+
+/** \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>
+SparseInnerVector<Derived> SparseMatrixBase<Derived>::innerVector(int outer)
+{ return SparseInnerVector<Derived>(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); }
+
+# if 0
 template<typename MatrixType, int BlockRows, int BlockCols, int PacketAccess>
 class Block<MatrixType,BlockRows,BlockCols,PacketAccess,IsSparse>
   : public SparseMatrixBase<Block<MatrixType,BlockRows,BlockCols,PacketAccess,IsSparse> >
@@ -117,6 +184,6 @@ public:
     const ei_int_if_dynamic<ColsAtCompileTime> m_blockCols;
 
 };
+#endif
 
-
-#endif // EIGEN_SPARSEBLOCK_H
+#endif // EIGEN_SPARSE_BLOCK_H
diff --git a/Eigen/src/Sparse/SparseCwise.h b/Eigen/src/Sparse/SparseCwise.h
new file mode 100644
index 000000000..0697439d3
--- /dev/null
+++ b/Eigen/src/Sparse/SparseCwise.h
@@ -0,0 +1,181 @@
+// 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 Benoit Jacob <jacob.benoit.1@gmail.com>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_SPARSE_CWISE_H
+#define EIGEN_SPARSE_CWISE_H
+
+/** \internal
+  * convenient macro to defined the return type of a cwise binary operation */
+#define EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(OP) \
+    CwiseBinaryOp<OP<typename ei_traits<ExpressionType>::Scalar>, ExpressionType, OtherDerived>
+
+#define EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE \
+    SparseCwiseBinaryOp< \
+      ei_scalar_product_op< \
+        typename ei_scalar_product_traits< \
+          typename ei_traits<ExpressionType>::Scalar, \
+          typename ei_traits<OtherDerived>::Scalar \
+        >::ReturnType \
+      >, \
+      ExpressionType, \
+      OtherDerived \
+    >
+
+/** \internal
+  * convenient macro to defined the return type of a cwise unary operation */
+#define EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(OP) \
+    SparseCwiseUnaryOp<OP<typename ei_traits<ExpressionType>::Scalar>, ExpressionType>
+
+/** \internal
+  * convenient macro to defined the return type of a cwise comparison to a scalar */
+/*#define EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(OP) \
+    CwiseBinaryOp<OP<typename ei_traits<ExpressionType>::Scalar>, ExpressionType, \
+        NestByValue<typename ExpressionType::ConstantReturnType> >*/
+
+template<typename ExpressionType> class SparseCwise
+{
+  public:
+
+    typedef typename ei_traits<ExpressionType>::Scalar Scalar;
+    typedef typename ei_meta_if<ei_must_nest_by_value<ExpressionType>::ret,
+        ExpressionType, const ExpressionType&>::ret ExpressionTypeNested;
+    typedef CwiseUnaryOp<ei_scalar_add_op<Scalar>, ExpressionType> ScalarAddReturnType;
+
+    inline SparseCwise(const ExpressionType& matrix) : m_matrix(matrix) {}
+
+    /** \internal */
+    inline const ExpressionType& _expression() const { return m_matrix; }
+
+    template<typename OtherDerived>
+    const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE
+    operator*(const SparseMatrixBase<OtherDerived> &other) const;
+
+    template<typename OtherDerived>
+    const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)
+    operator/(const SparseMatrixBase<OtherDerived> &other) const;
+
+    template<typename OtherDerived>
+    const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op)
+    min(const SparseMatrixBase<OtherDerived> &other) const;
+
+    template<typename OtherDerived>
+    const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op)
+    max(const SparseMatrixBase<OtherDerived> &other) const;
+
+    const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_abs_op)      abs() const;
+    const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_abs2_op)     abs2() const;
+//     const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_square_op)   square() const;
+//     const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_cube_op)     cube() const;
+//     const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_inverse_op)  inverse() const;
+//     const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_sqrt_op)     sqrt() const;
+//     const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_exp_op)      exp() const;
+//     const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_log_op)      log() const;
+//     const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_cos_op)      cos() const;
+//     const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_sin_op)      sin() const;
+//     const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_pow_op)      pow(const Scalar& exponent) const;
+
+    /*
+    const ScalarAddReturnType
+    operator+(const Scalar& scalar) const;
+
+    friend const ScalarAddReturnType
+    operator+(const Scalar& scalar, const Cwise& mat)
+    { return mat + scalar; }
+
+    ExpressionType& operator+=(const Scalar& scalar);
+
+    const ScalarAddReturnType
+    operator-(const Scalar& scalar) const;
+
+    ExpressionType& operator-=(const Scalar& scalar);
+
+    template<typename OtherDerived>
+    inline ExpressionType& operator*=(const MatrixBase<OtherDerived> &other);
+
+    template<typename OtherDerived>
+    inline ExpressionType& operator/=(const MatrixBase<OtherDerived> &other);
+
+    template<typename OtherDerived> const EIGEN_CWISE_BINOP_RETURN_TYPE(std::less)
+    operator<(const MatrixBase<OtherDerived>& other) const;
+
+    template<typename OtherDerived> const EIGEN_CWISE_BINOP_RETURN_TYPE(std::less_equal)
+    operator<=(const MatrixBase<OtherDerived>& other) const;
+
+    template<typename OtherDerived> const EIGEN_CWISE_BINOP_RETURN_TYPE(std::greater)
+    operator>(const MatrixBase<OtherDerived>& other) const;
+
+    template<typename OtherDerived> const EIGEN_CWISE_BINOP_RETURN_TYPE(std::greater_equal)
+    operator>=(const MatrixBase<OtherDerived>& other) const;
+
+    template<typename OtherDerived> const EIGEN_CWISE_BINOP_RETURN_TYPE(std::equal_to)
+    operator==(const MatrixBase<OtherDerived>& other) const;
+
+    template<typename OtherDerived> const EIGEN_CWISE_BINOP_RETURN_TYPE(std::not_equal_to)
+    operator!=(const MatrixBase<OtherDerived>& other) const;
+
+    // comparisons to a scalar value
+    const EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(std::less)
+    operator<(Scalar s) const;
+
+    const EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(std::less_equal)
+    operator<=(Scalar s) const;
+
+    const EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(std::greater)
+    operator>(Scalar s) const;
+
+    const EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(std::greater_equal)
+    operator>=(Scalar s) const;
+
+    const EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(std::equal_to)
+    operator==(Scalar s) const;
+
+    const EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(std::not_equal_to)
+    operator!=(Scalar s) const;
+    */
+
+    // allow to extend SparseCwise outside Eigen
+    #ifdef EIGEN_SPARSE_CWISE_PLUGIN
+    #include EIGEN_SPARSE_CWISE_PLUGIN
+    #endif
+
+  protected:
+    ExpressionTypeNested m_matrix;
+};
+
+template<typename Derived>
+inline const SparseCwise<Derived>
+SparseMatrixBase<Derived>::cwise() const
+{
+  return derived();
+}
+
+template<typename Derived>
+inline SparseCwise<Derived>
+SparseMatrixBase<Derived>::cwise()
+{
+  return derived();
+}
+
+#endif // EIGEN_SPARSE_CWISE_H
diff --git a/Eigen/src/Sparse/SparseCwiseBinaryOp.h b/Eigen/src/Sparse/SparseCwiseBinaryOp.h
new file mode 100644
index 000000000..99bbae0cd
--- /dev/null
+++ b/Eigen/src/Sparse/SparseCwiseBinaryOp.h
@@ -0,0 +1,560 @@
+// 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_SPARSE_CWISE_BINARY_OP_H
+#define EIGEN_SPARSE_CWISE_BINARY_OP_H
+
+// Here we have to handle 3 cases:
+//  1 - sparse op dense
+//  2 - dense op sparse
+//  3 - sparse op sparse
+// We also need to implement a 4th iterator for:
+//  4 - dense op dense
+// Finally, we also need to distinguish between the product and other operations :
+//                configuration      returned mode
+//  1 - sparse op dense    product      sparse
+//                         generic      dense
+//  2 - dense op sparse    product      sparse
+//                         generic      dense
+//  3 - sparse op sparse   product      sparse
+//                         generic      sparse
+//  4 - dense op dense     product      dense
+//                         generic      dense
+
+// template<typename BinaryOp, typename Lhs, typename Rhs> class ei_sparse_cwise_binary_op
+//   : ei_no_assignment_operator
+// {
+//     typedef CwiseBinaryOp<BinaryOp,Lhs,Rhs> CwiseBinaryXpr;
+//   public:
+//
+//     typedef typename ei_traits<CwiseBinaryXpr>::LhsNested LhsNested;
+//     typedef typename ei_traits<CwiseBinaryXpr>::RhsNested RhsNested;
+//
+//     EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp())
+//       : m_lhs(lhs), m_rhs(rhs), m_functor(func)
+//     {
+//       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)),
+//         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
+//       // require the sizes to match
+//       EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs)
+//       ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols());
+//     }
+//
+//     EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); }
+//     EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); }
+//
+//     EIGEN_STRONG_INLINE const LhsNested& _lhs() const { return m_lhs; }
+//     EIGEN_STRONG_INLINE const RhsNested& _rhs() const { return m_rhs; }
+//     EIGEN_STRONG_INLINE const UnaryOp& _functor() const { return m_functor; }
+//
+//   protected:
+//     const LhsNested m_lhs;
+//     const RhsNested m_rhs;
+//     const BinaryOp m_functor;
+// }
+
+template<typename BinaryOp, typename Lhs, typename Rhs>
+struct ei_traits<SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> >
+{
+  typedef typename ei_result_of<
+                     BinaryOp(
+                       typename Lhs::Scalar,
+                       typename Rhs::Scalar
+                     )
+                   >::type Scalar;
+  typedef typename Lhs::Nested LhsNested;
+  typedef typename Rhs::Nested RhsNested;
+  typedef typename ei_unref<LhsNested>::type _LhsNested;
+  typedef typename ei_unref<RhsNested>::type _RhsNested;
+  enum {
+    LhsCoeffReadCost = _LhsNested::CoeffReadCost,
+    RhsCoeffReadCost = _RhsNested::CoeffReadCost,
+    LhsFlags = _LhsNested::Flags,
+    RhsFlags = _RhsNested::Flags,
+    RowsAtCompileTime = Lhs::RowsAtCompileTime,
+    ColsAtCompileTime = Lhs::ColsAtCompileTime,
+    MaxRowsAtCompileTime = Lhs::MaxRowsAtCompileTime,
+    MaxColsAtCompileTime = Lhs::MaxColsAtCompileTime,
+    Flags = (int(LhsFlags) | int(RhsFlags)) & HereditaryBits,
+    CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + ei_functor_traits<BinaryOp>::Cost
+  };
+};
+
+template<typename BinaryOp, typename Lhs, typename Rhs>
+class SparseCwiseBinaryOp : ei_no_assignment_operator,
+  public SparseMatrixBase<SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> >
+{
+  public:
+
+    class InnerIterator;
+
+    EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseCwiseBinaryOp)
+    typedef typename ei_traits<SparseCwiseBinaryOp>::LhsNested LhsNested;
+    typedef typename ei_traits<SparseCwiseBinaryOp>::RhsNested RhsNested;
+    typedef typename ei_unref<LhsNested>::type _LhsNested;
+    typedef typename ei_unref<RhsNested>::type _RhsNested;
+
+    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((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)),
+        YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
+      // require the sizes to match
+      EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs)
+      ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols());
+    }
+
+    EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); }
+    EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); }
+
+    EIGEN_STRONG_INLINE const _LhsNested& _lhs() const { return m_lhs; }
+    EIGEN_STRONG_INLINE const _RhsNested& _rhs() const { return m_rhs; }
+    EIGEN_STRONG_INLINE const BinaryOp& _functor() const { return m_functor; }
+
+  protected:
+    const LhsNested m_lhs;
+    const RhsNested m_rhs;
+    const BinaryOp m_functor;
+};
+/*
+template<typename BinaryOp, typename Lhs, typename Rhs>
+class CwiseBinaryOp<BinaryOp, SparseMatrixBase<Lhs>, Rhs>
+  : public ei_sparse_cwise_binary_op<BinaryOp, Lhs, Rhs>,
+    public SparseMatrixBase<CwiseBinaryOp<BinaryOp, SparseMatrixBase<Lhs>, Rhs> >
+{
+  public:
+    EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp)
+    EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp())
+      : ei_sparse_cwise_binary_op<BinaryOp, Lhs, Rhs>(lhs, rhs, func)
+    {}
+};
+template<typename BinaryOp, typename Lhs, typename Rhs>
+class CwiseBinaryOp<BinaryOp, SparseMatrixBase<Lhs>, SparseMatrixBase<Rhs> >
+  : public ei_sparse_cwise_binary_op<BinaryOp, Lhs, Rhs>,
+    public SparseMatrixBase<CwiseBinaryOp<BinaryOp, Lhs, Rhs> >
+{
+  public:
+    EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp)
+    EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp())
+      : ei_sparse_cwise_binary_op<BinaryOp, Lhs, Rhs>(lhs, rhs, func)
+    {}
+};*/
+
+// template<typename BinaryOp, typename Lhs, typename RhsDerived>
+// class CwiseBinaryOp<BinaryOp, Lhs, SparseMatrixBase<Rhs> > : ei_no_assignment_operator,
+//   public SparseMatrixBase<CwiseBinaryOp<BinaryOp, SparseMatrixBase<LhsDerived>, SparseMatrixBase<Rhs> > >
+// {
+//   public:
+//
+//     EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp)
+//     typedef typename ei_traits<CwiseBinaryOp>::LhsNested LhsNested;
+//     typedef typename ei_traits<CwiseBinaryOp>::RhsNested RhsNested;
+//
+//     EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp())
+//       : m_lhs(lhs), m_rhs(rhs), m_functor(func)
+//     {
+//       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)),
+//         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
+//       // require the sizes to match
+//       EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs)
+//       ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols());
+//     }
+//
+//     EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); }
+//     EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); }
+//
+//   protected:
+//     const LhsNested m_lhs;
+//     const RhsNested m_rhs;
+//     const BinaryOp m_functor;
+// };
+//
+// template<typename BinaryOp, typename LhsDerived, typename RhsDerived>
+// class CwiseBinaryOp<BinaryOp, SparseMatrixBase<LhsDerived>, SparseMatrixBase<Rhs> > : ei_no_assignment_operator,
+//   public SparseMatrixBase<CwiseBinaryOp<BinaryOp, SparseMatrixBase<LhsDerived>, SparseMatrixBase<Rhs> > >
+// {
+//   public:
+//
+//     EIGEN_GENERIC_PUBLIC_INTERFACE(CwiseBinaryOp)
+//     typedef typename ei_traits<CwiseBinaryOp>::LhsNested LhsNested;
+//     typedef typename ei_traits<CwiseBinaryOp>::RhsNested RhsNested;
+//
+//     EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& lhs, const Rhs& rhs, const BinaryOp& func = BinaryOp())
+//       : m_lhs(lhs), m_rhs(rhs), m_functor(func)
+//     {
+//       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)),
+//         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
+//       // require the sizes to match
+//       EIGEN_STATIC_ASSERT_SAME_MATRIX_SIZE(Lhs, Rhs)
+//       ei_assert(lhs.rows() == rhs.rows() && lhs.cols() == rhs.cols());
+//     }
+//
+//     EIGEN_STRONG_INLINE int rows() const { return m_lhs.rows(); }
+//     EIGEN_STRONG_INLINE int cols() const { return m_lhs.cols(); }
+//
+//   protected:
+//     const LhsNested m_lhs;
+//     const RhsNested m_rhs;
+//     const BinaryOp m_functor;
+// };
+
+// 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 }; };
+
+template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived,
+  int _LhsStorageMode = int(Lhs::Flags) & SparseBit,
+  int _RhsStorageMode = int(Rhs::Flags) & SparseBit>
+class ei_sparse_cwise_binary_op_inner_iterator_selector;
+
+template<typename BinaryOp, typename Lhs, typename Rhs>
+class SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator
+  : public ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp,Lhs,Rhs, typename SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator>
+{
+    typedef ei_sparse_cwise_binary_op_inner_iterator_selector<
+      BinaryOp,Lhs,Rhs, typename SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator> Base;
+  public:
+    typedef typename SparseCwiseBinaryOp::Scalar Scalar;
+    typedef typename ei_traits<SparseCwiseBinaryOp>::_LhsNested _LhsNested;
+    typedef typename _LhsNested::InnerIterator LhsIterator;
+    typedef typename ei_traits<SparseCwiseBinaryOp>::_RhsNested _RhsNested;
+    typedef typename _RhsNested::InnerIterator RhsIterator;
+
+    EIGEN_STRONG_INLINE InnerIterator(const SparseCwiseBinaryOp& binOp, int outer)
+      : Base(binOp,outer)
+    {}
+};
+
+/***************************************************************************
+* 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
+
+
+// sparse - sparse  (generic)
+template<typename BinaryOp, typename Lhs, typename Rhs, typename Derived>
+class ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp, Lhs, Rhs, Derived, IsSparse, IsSparse>
+{
+    typedef SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> CwiseBinaryXpr;
+    typedef typename ei_traits<CwiseBinaryXpr>::Scalar Scalar;
+    typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
+    typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
+    typedef typename _LhsNested::InnerIterator LhsIterator;
+    typedef typename _RhsNested::InnerIterator RhsIterator;
+  public:
+
+    EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
+      : m_lhsIter(xpr._lhs(),outer), m_rhsIter(xpr._rhs()), m_functor(xpr._functor())
+    {
+      this->operator++();
+    }
+
+    EIGEN_STRONG_INLINE Derived& operator++()
+    {
+      if (m_lhsIter && m_rhsIter && (m_lhsIter.index() == m_rhsIter.index()))
+      {
+        m_id = m_lhsIter.index();
+        m_value = m_functor(m_lhsIter.value(), m_rhsIter.value());
+        ++m_lhsIter;
+        ++m_rhsIter;
+      }
+      else if (m_lhsIter && (!m_rhsIter || (m_lhsIter.index() < m_rhsIter.index())))
+      {
+        m_id = m_lhsIter.index();
+        m_value = m_functor(m_lhsIter.value(), Scalar(0));
+        ++m_lhsIter;
+      }
+      else if (m_rhsIter && (!m_lhsIter || (m_lhsIter.index() > m_rhsIter.index())))
+      {
+        m_id = m_rhsIter.index();
+        m_value = m_functor(Scalar(0), m_rhsIter.value());
+        ++m_rhsIter;
+      }
+      else
+      {
+        m_id = -1;
+      }
+      return *static_cast<Derived*>(this);
+    }
+
+    EIGEN_STRONG_INLINE Scalar value() const { return m_value; }
+
+    EIGEN_STRONG_INLINE int index() const { return m_id; }
+
+    EIGEN_STRONG_INLINE operator bool() const { return m_id>=0; }
+
+  protected:
+    LhsIterator m_lhsIter;
+    RhsIterator m_rhsIter;
+    const BinaryOp& m_functor;
+    Scalar m_value;
+    int m_id;
+};
+
+// template<typename BinaryOp, typename Lhs, typename Rhs>
+// class SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs>::InnerIterator
+//   : public ei_sparse_cwise_binary_op_inner_iterator_selector<BinaryOp,Lhs,Rhs, typename SparseCwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator>
+// {
+//     typedef CwiseBinaryOpInnerIterator<
+//     BinaryOp,Lhs,Rhs, typename CwiseBinaryOp<BinaryOp,Lhs,Rhs>::InnerIterator> Base;
+//   public:
+//     typedef typename CwiseBinaryOp::Scalar Scalar;
+//     typedef typename ei_traits<CwiseBinaryOp>::_LhsNested _LhsNested;
+//     typedef typename _LhsNested::InnerIterator LhsIterator;
+//     typedef typename ei_traits<CwiseBinaryOp>::_RhsNested _RhsNested;
+//     typedef typename _RhsNested::InnerIterator RhsIterator;
+// //   public:
+//     EIGEN_STRONG_INLINE InnerIterator(const CwiseBinaryOp& binOp, int outer)
+//       : Base(binOp.m_lhs,binOp.m_rhs,binOp.m_functor,outer)
+//     {}
+// };
+
+// sparse - sparse  (product)
+template<typename T, typename Lhs, typename Rhs, typename Derived>
+class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsSparse, IsSparse>
+{
+    typedef ei_scalar_product_op<T> BinaryOp;
+    typedef SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> CwiseBinaryXpr;
+    typedef typename CwiseBinaryXpr::Scalar Scalar;
+    typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
+    typedef typename _LhsNested::InnerIterator LhsIterator;
+    typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
+    typedef typename _RhsNested::InnerIterator RhsIterator;
+  public:
+
+    EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
+      : m_lhsIter(xpr._lhs(),outer), m_rhsIter(xpr._rhs()), m_functor(xpr._functor())
+    {
+      while (m_lhsIter && m_rhsIter && m_lhsIter.index() != m_rhsIter.index())
+      {
+        if (m_lhsIter.index() < m_rhsIter.index())
+          ++m_lhsIter;
+        else
+          ++m_rhsIter;
+      }
+    }
+
+    EIGEN_STRONG_INLINE Derived& operator++()
+    {
+      while (m_lhsIter && m_rhsIter)
+      {
+        if (m_lhsIter.index() == m_rhsIter.index())
+        {
+          ++m_lhsIter;
+          ++m_rhsIter;
+          break;
+        }
+        else if (m_lhsIter.index() < m_rhsIter.index())
+          ++m_lhsIter;
+        else
+          ++m_rhsIter;
+      }
+      return *static_cast<Derived*>(this);
+    }
+
+    EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_lhsIter.value(), m_rhsIter.value()); }
+
+    EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); }
+
+    EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter && m_rhsIter; }
+
+  protected:
+    LhsIterator m_lhsIter;
+    RhsIterator m_rhsIter;
+    const BinaryOp& m_functor;
+};
+
+// sparse - dense  (product)
+template<typename T, typename Lhs, typename Rhs, typename Derived>
+class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsSparse, IsDense>
+{
+    typedef ei_scalar_product_op<T> BinaryOp;
+    typedef SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> CwiseBinaryXpr;
+    typedef typename CwiseBinaryXpr::Scalar Scalar;
+    typedef typename ei_traits<CwiseBinaryXpr>::_LhsNested _LhsNested;
+    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)
+    {}
+
+    EIGEN_STRONG_INLINE Derived& operator++()
+    {
+      ++m_lhsIter;
+      return *static_cast<Derived*>(this);
+    }
+
+    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)); }
+
+    EIGEN_STRONG_INLINE int index() const { return m_lhsIter.index(); }
+
+    EIGEN_STRONG_INLINE operator bool() const { return m_lhsIter; }
+
+  protected:
+    const CwiseBinaryXpr& m_xpr;
+    LhsIterator m_lhsIter;
+    const BinaryOp& m_functor;
+    const int m_outer;
+};
+
+// sparse - dense  (product)
+template<typename T, typename Lhs, typename Rhs, typename Derived>
+class ei_sparse_cwise_binary_op_inner_iterator_selector<ei_scalar_product_op<T>, Lhs, Rhs, Derived, IsDense, IsSparse>
+{
+    typedef ei_scalar_product_op<T> BinaryOp;
+    typedef SparseCwiseBinaryOp<BinaryOp, Lhs, Rhs> CwiseBinaryXpr;
+    typedef typename CwiseBinaryXpr::Scalar Scalar;
+    typedef typename ei_traits<CwiseBinaryXpr>::_RhsNested _RhsNested;
+    typedef typename _RhsNested::InnerIterator RhsIterator;
+    enum { IsRowMajor = (int(Rhs::Flags)&RowMajorBit)==RowMajorBit };
+  public:
+
+    EIGEN_STRONG_INLINE ei_sparse_cwise_binary_op_inner_iterator_selector(const CwiseBinaryXpr& xpr, int outer)
+      : m_xpr(xpr), m_rhsIter(xpr.rhs(),outer), m_functor(xpr.functor()), m_outer(outer)
+    {}
+
+    EIGEN_STRONG_INLINE Derived& operator++()
+    {
+      ++m_rhsIter;
+      return *static_cast<Derived*>(this);
+    }
+
+    EIGEN_STRONG_INLINE Scalar value() const
+    { return m_functor(m_xpr.lhs().coeff(IsRowMajor?m_outer:m_rhsIter.index(),IsRowMajor?m_rhsIter.index():m_outer), m_rhsIter.value()); }
+
+    EIGEN_STRONG_INLINE int index() const { return m_rhsIter.index(); }
+
+    EIGEN_STRONG_INLINE operator bool() const { return m_rhsIter; }
+
+  protected:
+    const CwiseBinaryXpr& m_xpr;
+    RhsIterator m_rhsIter;
+    const BinaryOp& m_functor;
+    const int m_outer;
+};
+
+
+template<typename Derived>
+template<typename OtherDerived>
+EIGEN_STRONG_INLINE const SparseCwiseBinaryOp<ei_scalar_difference_op<typename ei_traits<Derived>::Scalar>,
+                                 Derived, OtherDerived>
+SparseMatrixBase<Derived>::operator-(const SparseMatrixBase<OtherDerived> &other) const
+{
+  return SparseCwiseBinaryOp<ei_scalar_difference_op<Scalar>,
+                       Derived, OtherDerived>(derived(), other.derived());
+}
+
+template<typename Derived>
+template<typename OtherDerived>
+EIGEN_STRONG_INLINE Derived &
+SparseMatrixBase<Derived>::operator-=(const SparseMatrixBase<OtherDerived> &other)
+{
+  return *this = *this - other;
+}
+
+template<typename Derived>
+template<typename OtherDerived>
+EIGEN_STRONG_INLINE const SparseCwiseBinaryOp<ei_scalar_sum_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
+SparseMatrixBase<Derived>::operator+(const SparseMatrixBase<OtherDerived> &other) const
+{
+  return SparseCwiseBinaryOp<ei_scalar_sum_op<Scalar>, Derived, OtherDerived>(derived(), other.derived());
+}
+
+template<typename Derived>
+template<typename OtherDerived>
+EIGEN_STRONG_INLINE Derived &
+SparseMatrixBase<Derived>::operator+=(const SparseMatrixBase<OtherDerived>& other)
+{
+  return *this = *this + other;
+}
+
+template<typename ExpressionType>
+template<typename OtherDerived>
+EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE
+SparseCwise<ExpressionType>::operator*(const SparseMatrixBase<OtherDerived> &other) const
+{
+  return EIGEN_SPARSE_CWISE_PRODUCT_RETURN_TYPE(_expression(), other.derived());
+}
+
+template<typename ExpressionType>
+template<typename OtherDerived>
+EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)
+SparseCwise<ExpressionType>::operator/(const SparseMatrixBase<OtherDerived> &other) const
+{
+  return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_quotient_op)(_expression(), other.derived());
+}
+
+// template<typename ExpressionType>
+// template<typename OtherDerived>
+// inline ExpressionType& Cwise<ExpressionType>::operator*=(const SparseMatrixBase<OtherDerived> &other)
+// {
+//   return m_matrix.const_cast_derived() = *this * other;
+// }
+
+// template<typename ExpressionType>
+// template<typename OtherDerived>
+// inline ExpressionType& Cwise<ExpressionType>::operator/=(const SparseMatrixBase<OtherDerived> &other)
+// {
+//   return m_matrix.const_cast_derived() = *this / other;
+// }
+
+template<typename ExpressionType>
+template<typename OtherDerived>
+EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op)
+SparseCwise<ExpressionType>::min(const SparseMatrixBase<OtherDerived> &other) const
+{
+  return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_min_op)(_expression(), other.derived());
+}
+
+template<typename ExpressionType>
+template<typename OtherDerived>
+EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op)
+SparseCwise<ExpressionType>::max(const SparseMatrixBase<OtherDerived> &other) const
+{
+  return EIGEN_SPARSE_CWISE_BINOP_RETURN_TYPE(ei_scalar_max_op)(_expression(), other.derived());
+}
+
+// template<typename Derived>
+// template<typename CustomBinaryOp, typename OtherDerived>
+// EIGEN_STRONG_INLINE const CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>
+// SparseMatrixBase<Derived>::binaryExpr(const SparseMatrixBase<OtherDerived> &other, const CustomBinaryOp& func) const
+// {
+//   return CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>(derived(), other.derived(), func);
+// }
+
+#endif // EIGEN_SPARSE_CWISE_BINARY_OP_H
diff --git a/Eigen/src/Sparse/SparseCwiseUnaryOp.h b/Eigen/src/Sparse/SparseCwiseUnaryOp.h
new file mode 100644
index 000000000..2b2050edd
--- /dev/null
+++ b/Eigen/src/Sparse/SparseCwiseUnaryOp.h
@@ -0,0 +1,175 @@
+// 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_SPARSE_CWISE_UNARY_OP_H
+#define EIGEN_SPARSE_CWISE_UNARY_OP_H
+
+template<typename UnaryOp, typename MatrixType>
+struct ei_traits<SparseCwiseUnaryOp<UnaryOp, MatrixType> > : ei_traits<MatrixType>
+{
+  typedef typename ei_result_of<
+                     UnaryOp(typename MatrixType::Scalar)
+                   >::type Scalar;
+  typedef typename MatrixType::Nested MatrixTypeNested;
+  typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested;
+  enum {
+    CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ei_functor_traits<UnaryOp>::Cost
+  };
+};
+
+template<typename UnaryOp, typename MatrixType>
+class SparseCwiseUnaryOp : ei_no_assignment_operator,
+  public SparseMatrixBase<SparseCwiseUnaryOp<UnaryOp, MatrixType> >
+{
+  public:
+
+    class InnerIterator;
+//     typedef typename ei_unref<LhsNested>::type _LhsNested;
+
+    EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseCwiseUnaryOp)
+
+    inline SparseCwiseUnaryOp(const MatrixType& mat, const UnaryOp& func = UnaryOp())
+      : m_matrix(mat), m_functor(func) {}
+
+    EIGEN_STRONG_INLINE int rows() const { return m_matrix.rows(); }
+    EIGEN_STRONG_INLINE int cols() const { return m_matrix.cols(); }
+
+//     EIGEN_STRONG_INLINE const typename MatrixType::Nested& _matrix() const { return m_matrix; }
+//     EIGEN_STRONG_INLINE const UnaryOp& _functor() const { return m_functor; }
+
+  protected:
+    const typename MatrixType::Nested m_matrix;
+    const UnaryOp m_functor;
+};
+
+
+template<typename UnaryOp, typename MatrixType>
+class SparseCwiseUnaryOp<UnaryOp,MatrixType>::InnerIterator
+{
+    typedef typename SparseCwiseUnaryOp::Scalar Scalar;
+    typedef typename ei_traits<SparseCwiseUnaryOp>::_MatrixTypeNested _MatrixTypeNested;
+    typedef typename _MatrixTypeNested::InnerIterator MatrixTypeIterator;
+  public:
+
+    EIGEN_STRONG_INLINE InnerIterator(const SparseCwiseUnaryOp& unaryOp, int outer)
+      : m_iter(unaryOp.m_matrix,outer), m_functor(unaryOp.m_functor)
+    {}
+
+    EIGEN_STRONG_INLINE InnerIterator& operator++()
+    { ++m_iter; return *this; }
+
+    EIGEN_STRONG_INLINE Scalar value() const { return m_functor(m_iter.value()); }
+
+    EIGEN_STRONG_INLINE int index() const { return m_iter.index(); }
+
+    EIGEN_STRONG_INLINE operator bool() const { return m_iter; }
+
+  protected:
+    MatrixTypeIterator m_iter;
+    const UnaryOp& m_functor;
+};
+
+template<typename Derived>
+template<typename CustomUnaryOp>
+EIGEN_STRONG_INLINE const SparseCwiseUnaryOp<CustomUnaryOp, Derived>
+SparseMatrixBase<Derived>::unaryExpr(const CustomUnaryOp& func) const
+{
+  return SparseCwiseUnaryOp<CustomUnaryOp, Derived>(derived(), func);
+}
+
+template<typename Derived>
+EIGEN_STRONG_INLINE const SparseCwiseUnaryOp<ei_scalar_opposite_op<typename ei_traits<Derived>::Scalar>,Derived>
+SparseMatrixBase<Derived>::operator-() const
+{
+  return derived();
+}
+
+template<typename ExpressionType>
+EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_abs_op)
+SparseCwise<ExpressionType>::abs() const
+{
+  return _expression();
+}
+
+template<typename ExpressionType>
+EIGEN_STRONG_INLINE const EIGEN_SPARSE_CWISE_UNOP_RETURN_TYPE(ei_scalar_abs2_op)
+SparseCwise<ExpressionType>::abs2() const
+{
+  return _expression();
+}
+
+template<typename Derived>
+EIGEN_STRONG_INLINE typename SparseMatrixBase<Derived>::ConjugateReturnType
+SparseMatrixBase<Derived>::conjugate() const
+{
+  return ConjugateReturnType(derived());
+}
+
+template<typename Derived>
+EIGEN_STRONG_INLINE const typename SparseMatrixBase<Derived>::RealReturnType
+SparseMatrixBase<Derived>::real() const { return derived(); }
+
+template<typename Derived>
+EIGEN_STRONG_INLINE const typename SparseMatrixBase<Derived>::ImagReturnType
+SparseMatrixBase<Derived>::imag() const { return derived(); }
+
+template<typename Derived>
+template<typename NewType>
+EIGEN_STRONG_INLINE const SparseCwiseUnaryOp<ei_scalar_cast_op<typename ei_traits<Derived>::Scalar, NewType>, Derived>
+SparseMatrixBase<Derived>::cast() const
+{
+  return derived();
+}
+
+template<typename Derived>
+EIGEN_STRONG_INLINE const SparseCwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<Derived>::Scalar>, Derived>
+SparseMatrixBase<Derived>::operator*(const Scalar& scalar) const
+{
+  return SparseCwiseUnaryOp<ei_scalar_multiple_op<Scalar>, Derived>
+    (derived(), ei_scalar_multiple_op<Scalar>(scalar));
+}
+
+template<typename Derived>
+EIGEN_STRONG_INLINE const SparseCwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<Derived>::Scalar>, Derived>
+SparseMatrixBase<Derived>::operator/(const Scalar& scalar) const
+{
+  return SparseCwiseUnaryOp<ei_scalar_quotient1_op<Scalar>, Derived>
+    (derived(), ei_scalar_quotient1_op<Scalar>(scalar));
+}
+
+template<typename Derived>
+EIGEN_STRONG_INLINE Derived&
+SparseMatrixBase<Derived>::operator*=(const Scalar& other)
+{
+  return *this = *this * other;
+}
+
+template<typename Derived>
+EIGEN_STRONG_INLINE Derived&
+SparseMatrixBase<Derived>::operator/=(const Scalar& other)
+{
+  return *this = *this / other;
+}
+
+#endif // EIGEN_SPARSE_CWISE_UNARY_OP_H
diff --git a/Eigen/src/Sparse/SparseDot.h b/Eigen/src/Sparse/SparseDot.h
new file mode 100644
index 000000000..7a26e0f4b
--- /dev/null
+++ b/Eigen/src/Sparse/SparseDot.h
@@ -0,0 +1,97 @@
+// 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_SPARSE_DOT_H
+#define EIGEN_SPARSE_DOT_H
+
+template<typename Derived>
+template<typename OtherDerived>
+typename ei_traits<Derived>::Scalar
+SparseMatrixBase<Derived>::dot(const MatrixBase<OtherDerived>& other) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
+  EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived)
+  EIGEN_STATIC_ASSERT((ei_is_same_type<Scalar, typename OtherDerived::Scalar>::ret),
+    YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
+
+  ei_assert(size() == other.size());
+  ei_assert(other.size()>0 && "you are using a non initialized vector");
+  
+  typename Derived::InnerIterator i(derived(),0);
+  Scalar res = 0;
+  while (i)
+  {
+    res += i.value() * ei_conj(other.coeff(i.index()));
+    ++i;
+  }
+  return res;
+}
+
+template<typename Derived>
+template<typename OtherDerived>
+typename ei_traits<Derived>::Scalar
+SparseMatrixBase<Derived>::dot(const SparseMatrixBase<OtherDerived>& other) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
+  EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(Derived,OtherDerived)
+  EIGEN_STATIC_ASSERT((ei_is_same_type<Scalar, typename OtherDerived::Scalar>::ret),
+    YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
+  
+  ei_assert(size() == other.size());
+  
+  typename Derived::InnerIterator i(derived(),0);
+  typename OtherDerived::InnerIterator j(other.derived(),0);
+  Scalar res = 0;
+  while (i && j)
+  {
+    if (i.index()==j.index())
+    {
+      res += i.value() * ei_conj(j.value());
+      ++i; ++j;
+    }
+    else if (i.index()<j.index())
+      ++i;
+    else
+      ++j;
+  }
+  return res;
+}
+
+template<typename Derived>
+inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
+SparseMatrixBase<Derived>::squaredNorm() const
+{
+  return ei_real((*this).cwise().abs2().sum());
+}
+
+template<typename Derived>
+inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real
+SparseMatrixBase<Derived>::norm() const
+{
+  return ei_sqrt(squaredNorm());
+}
+
+#endif // EIGEN_SPARSE_DOT_H
diff --git a/Eigen/src/Sparse/SparseFlagged.h b/Eigen/src/Sparse/SparseFlagged.h
new file mode 100644
index 000000000..2edac9c1e
--- /dev/null
+++ b/Eigen/src/Sparse/SparseFlagged.h
@@ -0,0 +1,98 @@
+// 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 Benoit Jacob <jacob.benoit.1@gmail.com>
+// 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_SPARSE_FLAGGED_H
+#define EIGEN_SPARSE_FLAGGED_H
+
+template<typename ExpressionType, unsigned int Added, unsigned int Removed>
+struct ei_traits<SparseFlagged<ExpressionType, Added, Removed> > : ei_traits<ExpressionType>
+{
+  enum { Flags = (ExpressionType::Flags | Added) & ~Removed };
+};
+
+template<typename ExpressionType, unsigned int Added, unsigned int Removed> class SparseFlagged
+  : public SparseMatrixBase<SparseFlagged<ExpressionType, Added, Removed> >
+{
+  public:
+
+    EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseFlagged)
+    class InnerIterator;
+    class ReverseInnerIterator;
+    
+    typedef typename ei_meta_if<ei_must_nest_by_value<ExpressionType>::ret,
+        ExpressionType, const ExpressionType&>::ret ExpressionTypeNested;
+    typedef typename ExpressionType::InnerIterator InnerIterator;
+
+    inline SparseFlagged(const ExpressionType& matrix) : m_matrix(matrix) {}
+
+    inline int rows() const { return m_matrix.rows(); }
+    inline int cols() const { return m_matrix.cols(); }
+    
+    // FIXME should be keep them ?
+    inline Scalar& coeffRef(int row, int col)
+    { return m_matrix.const_cast_derived().coeffRef(col, row); }
+
+    inline const Scalar coeff(int row, int col) const
+    { return m_matrix.coeff(col, row); }
+
+    inline const Scalar coeff(int index) const
+    { return m_matrix.coeff(index); }
+
+    inline Scalar& coeffRef(int index)
+    { return m_matrix.const_cast_derived().coeffRef(index); }
+
+  protected:
+    ExpressionTypeNested m_matrix;
+};
+
+template<typename ExpressionType, unsigned int Added, unsigned int Removed>
+  class SparseFlagged<ExpressionType,Added,Removed>::InnerIterator : public ExpressionType::InnerIterator
+{
+  public:
+
+    EIGEN_STRONG_INLINE InnerIterator(const SparseFlagged& xpr, int outer)
+      : ExpressionType::InnerIterator(xpr.m_matrix, outer)
+    {}
+};
+
+template<typename ExpressionType, unsigned int Added, unsigned int Removed>
+  class SparseFlagged<ExpressionType,Added,Removed>::ReverseInnerIterator : public ExpressionType::ReverseInnerIterator
+{
+  public:
+
+    EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseFlagged& xpr, int outer)
+      : ExpressionType::ReverseInnerIterator(xpr.m_matrix, outer)
+    {}
+};
+
+template<typename Derived>
+template<unsigned int Added>
+inline const SparseFlagged<Derived, Added, 0>
+SparseMatrixBase<Derived>::marked() const
+{
+  return derived();
+}
+
+#endif // EIGEN_SPARSE_FLAGGED_H
diff --git a/Eigen/src/Sparse/SparseFuzzy.h b/Eigen/src/Sparse/SparseFuzzy.h
new file mode 100644
index 000000000..355f4d52e
--- /dev/null
+++ b/Eigen/src/Sparse/SparseFuzzy.h
@@ -0,0 +1,41 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 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_SPARSE_FUZZY_H
+#define EIGEN_SPARSE_FUZZY_H
+
+// template<typename Derived>
+// template<typename OtherDerived>
+// bool SparseMatrixBase<Derived>::isApprox(
+//   const OtherDerived& other,
+//   typename NumTraits<Scalar>::Real prec
+// ) const
+// {
+//   const typename ei_nested<Derived,2>::type nested(derived());
+//   const typename ei_nested<OtherDerived,2>::type otherNested(other.derived());
+//   return    (nested - otherNested).cwise().abs2().sum()
+//          <= prec * prec * std::min(nested.cwise().abs2().sum(), otherNested.cwise().abs2().sum());
+// }
+
+#endif // EIGEN_SPARSE_FUZZY_H
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index 6820ae403..a732bdc31 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -56,14 +56,14 @@ class SparseMatrix
   : public SparseMatrixBase<SparseMatrix<_Scalar, _Flags> >
 {
   public:
-    EIGEN_GENERIC_PUBLIC_INTERFACE(SparseMatrix)
+    EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseMatrix)
 
   protected:
   public:
 
     typedef SparseMatrixBase<SparseMatrix> SparseBase;
     enum {
-      RowMajor = SparseBase::RowMajor
+      RowMajor = SparseBase::IsRowMajor
     };
     typedef SparseMatrix<Scalar,(Flags&~RowMajorBit)|(RowMajor?RowMajorBit:0)> TransposedSparseMatrix;
 
@@ -267,7 +267,7 @@ class SparseMatrix
     }
 
     template<typename OtherDerived>
-    inline SparseMatrix(const MatrixBase<OtherDerived>& other)
+    inline SparseMatrix(const SparseMatrixBase<OtherDerived>& other)
       : m_outerSize(0), m_innerSize(0), m_outerIndex(0)
     {
       *this = other.derived();
@@ -305,7 +305,7 @@ class SparseMatrix
     }
 
     template<typename OtherDerived>
-    inline SparseMatrix& operator=(const MatrixBase<OtherDerived>& other)
+    inline SparseMatrix& operator=(const SparseMatrixBase<OtherDerived>& other)
     {
       const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
       if (needToTranspose)
diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h
index d70c99259..d01fa1ec5 100644
--- a/Eigen/src/Sparse/SparseMatrixBase.h
+++ b/Eigen/src/Sparse/SparseMatrixBase.h
@@ -25,43 +25,131 @@
 #ifndef EIGEN_SPARSEMATRIXBASE_H
 #define EIGEN_SPARSEMATRIXBASE_H
 
-template<typename Derived>
-class SparseMatrixBase : public MatrixBase<Derived>
+template<typename Derived> class SparseMatrixBase
 {
   public:
 
-    typedef MatrixBase<Derived> Base;
-    typedef typename Base::Scalar Scalar;
-    typedef typename Base::RealScalar RealScalar;
+    typedef typename ei_traits<Derived>::Scalar Scalar;
+//     typedef typename Derived::InnerIterator InnerIterator;
+
     enum {
-      Flags = Base::Flags,
-      RowMajor = ei_traits<Derived>::Flags&RowMajorBit ? 1 : 0
+
+      RowsAtCompileTime = ei_traits<Derived>::RowsAtCompileTime,
+        /**< The number of rows at compile-time. This is just a copy of the value provided
+          * by the \a Derived type. If a value is not known at compile-time,
+          * it is set to the \a Dynamic constant.
+          * \sa MatrixBase::rows(), MatrixBase::cols(), ColsAtCompileTime, SizeAtCompileTime */
+
+      ColsAtCompileTime = ei_traits<Derived>::ColsAtCompileTime,
+        /**< The number of columns at compile-time. This is just a copy of the value provided
+          * by the \a Derived type. If a value is not known at compile-time,
+          * it is set to the \a Dynamic constant.
+          * \sa MatrixBase::rows(), MatrixBase::cols(), RowsAtCompileTime, SizeAtCompileTime */
+
+
+      SizeAtCompileTime = (ei_size_at_compile_time<ei_traits<Derived>::RowsAtCompileTime,
+                                                   ei_traits<Derived>::ColsAtCompileTime>::ret),
+        /**< This is equal to the number of coefficients, i.e. the number of
+          * rows times the number of columns, or to \a Dynamic if this is not
+          * known at compile-time. \sa RowsAtCompileTime, ColsAtCompileTime */
+
+      MaxRowsAtCompileTime = RowsAtCompileTime,
+      MaxColsAtCompileTime = ColsAtCompileTime,
+
+      MaxSizeAtCompileTime = (ei_size_at_compile_time<MaxRowsAtCompileTime,
+                                                      MaxColsAtCompileTime>::ret),
+
+      IsVectorAtCompileTime = RowsAtCompileTime == 1 || ColsAtCompileTime == 1,
+        /**< This is set to true if either the number of rows or the number of
+          * columns is known at compile-time to be equal to 1. Indeed, in that case,
+          * we are dealing with a column-vector (if there is only one column) or with
+          * a row-vector (if there is only one row). */
+
+      Flags = ei_traits<Derived>::Flags,
+        /**< This stores expression \ref flags flags which may or may not be inherited by new expressions
+          * constructed from this one. See the \ref flags "list of flags".
+          */
+
+      CoeffReadCost = ei_traits<Derived>::CoeffReadCost,
+        /**< This is a rough measure of how expensive it is to read one coefficient from
+          * this expression.
+          */
+
+      IsRowMajor = Flags&RowMajorBit ? 1 : 0
     };
 
+    /** \internal the return type of MatrixBase::conjugate() */
+    typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex,
+                        const SparseCwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, Derived>,
+                        const Derived&
+                     >::ret ConjugateReturnType;
+    /** \internal the return type of MatrixBase::real() */
+    typedef CwiseUnaryOp<ei_scalar_real_op<Scalar>, Derived> RealReturnType;
+    /** \internal the return type of MatrixBase::imag() */
+    typedef CwiseUnaryOp<ei_scalar_imag_op<Scalar>, Derived> ImagReturnType;
+    /** \internal the return type of MatrixBase::adjoint() */
+    typedef Eigen::Transpose<NestByValue<typename ei_cleantype<ConjugateReturnType>::type> >
+            AdjointReturnType;
+
+#ifndef EIGEN_PARSED_BY_DOXYGEN
+    /** This is the "real scalar" type; if the \a Scalar type is already real numbers
+      * (e.g. int, float or double) then \a RealScalar is just the same as \a Scalar. If
+      * \a Scalar is \a std::complex<T> then RealScalar is \a T.
+      *
+      * \sa class NumTraits
+      */
+    typedef typename NumTraits<Scalar>::Real RealScalar;
+
+    /** type of the equivalent square matrix */
+    typedef Matrix<Scalar,EIGEN_ENUM_MAX(RowsAtCompileTime,ColsAtCompileTime),
+                          EIGEN_ENUM_MAX(RowsAtCompileTime,ColsAtCompileTime)> SquareMatrixType;
+
     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)); }
+#endif // not EIGEN_PARSED_BY_DOXYGEN
 
-    SparseMatrixBase()
-      : m_isRValue(false)
-    {}
+    /** \returns the number of rows. \sa cols(), RowsAtCompileTime */
+    inline int rows() const { return derived().rows(); }
+    /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/
+    inline int cols() const { return derived().cols(); }
+    /** \returns the number of coefficients, which is \a rows()*cols().
+      * \sa rows(), cols(), SizeAtCompileTime. */
+    inline int size() const { return rows() * cols(); }
+    /** \returns the number of nonzero coefficients which is in practice the number
+      * of stored coefficients. */
+    inline int nonZeros() const { return derived.nonZeros(); }
+    /** \returns true if either the number of rows or the number of columns is equal to 1.
+      * In other words, this function returns
+      * \code rows()==1 || cols()==1 \endcode
+      * \sa rows(), cols(), IsVectorAtCompileTime. */
+    inline bool isVector() const { return rows()==1 || cols()==1; }
+    /** \returns the size of the storage major dimension,
+      * i.e., the number of columns for a columns major matrix, and the number of rows otherwise */
+    int outerSize() const { return (int(Flags)&RowMajorBit) ? this->rows() : this->cols(); }
+    /** \returns the size of the inner dimension according to the storage order,
+      * i.e., the number of rows for a columns major matrix, and the number of cols otherwise */
+    int innerSize() const { return (int(Flags)&RowMajorBit) ? this->cols() : this->rows(); }
 
     bool isRValue() const { return m_isRValue; }
     Derived& markAsRValue() { m_isRValue = true; return derived(); }
 
+    SparseMatrixBase() : m_isRValue(false) { /* TODO check flags */ }
+
     inline Derived& operator=(const Derived& other)
     {
 //       std::cout << "Derived& operator=(const Derived& other)\n";
-      if (other.isRValue())
-        derived().swap(other.const_cast_derived());
-      else
+//       if (other.isRValue())
+//         derived().swap(other.const_cast_derived());
+//       else
         this->operator=<Derived>(other);
       return derived();
     }
 
+
     template<typename OtherDerived>
-    inline Derived& operator=(const MatrixBase<OtherDerived>& other)
+    inline void assignGeneric(const OtherDerived& other)
     {
 //       std::cout << "Derived& operator=(const MatrixBase<OtherDerived>& other)\n";
       //const bool transpose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit);
@@ -87,9 +175,9 @@ class SparseMatrixBase : public MatrixBase<Derived>
       temp.endFill();
 
       derived() = temp.markAsRValue();
-      return derived();
     }
 
+
     template<typename OtherDerived>
     inline Derived& operator=(const SparseMatrixBase<OtherDerived>& other)
     {
@@ -110,7 +198,7 @@ class SparseMatrixBase : public MatrixBase<Derived>
             Scalar v = it.value();
             if (v!=Scalar(0))
             {
-              if (RowMajor) derived().fill(j,it.index()) = v;
+              if (IsRowMajor) derived().fill(j,it.index()) = v;
               else derived().fill(it.index(),j) = v;
             }
           }
@@ -118,12 +206,14 @@ class SparseMatrixBase : public MatrixBase<Derived>
         derived().endFill();
       }
       else
-        this->operator=<OtherDerived>(static_cast<const MatrixBase<OtherDerived>&>(other));
+      {
+        assignGeneric(other.derived());
+      }
       return derived();
     }
 
     template<typename Lhs, typename Rhs>
-    inline Derived& operator=(const Product<Lhs,Rhs,SparseProduct>& product);
+    inline Derived& operator=(const SparseProduct<Lhs,Rhs>& product);
 
     friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m)
     {
@@ -167,6 +257,292 @@ class SparseMatrixBase : public MatrixBase<Derived>
       return s;
     }
 
+    const SparseCwiseUnaryOp<ei_scalar_opposite_op<typename ei_traits<Derived>::Scalar>,Derived> operator-() const;
+
+    template<typename OtherDerived>
+    const SparseCwiseBinaryOp<ei_scalar_sum_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
+    operator+(const SparseMatrixBase<OtherDerived> &other) const;
+
+    template<typename OtherDerived>
+    const SparseCwiseBinaryOp<ei_scalar_difference_op<typename ei_traits<Derived>::Scalar>, Derived, OtherDerived>
+    operator-(const SparseMatrixBase<OtherDerived> &other) const;
+
+    template<typename OtherDerived>
+    Derived& operator+=(const SparseMatrixBase<OtherDerived>& other);
+    template<typename OtherDerived>
+    Derived& operator-=(const SparseMatrixBase<OtherDerived>& other);
+
+//     template<typename Lhs,typename Rhs>
+//     Derived& operator+=(const Flagged<Product<Lhs,Rhs,CacheFriendlyProduct>, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit>& other);
+
+    Derived& operator*=(const Scalar& other);
+    Derived& operator/=(const Scalar& other);
+
+    const SparseCwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<Derived>::Scalar>, Derived>
+    operator*(const Scalar& scalar) const;
+    const SparseCwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<Derived>::Scalar>, Derived>
+    operator/(const Scalar& scalar) const;
+
+    inline friend const SparseCwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<Derived>::Scalar>, Derived>
+    operator*(const Scalar& scalar, const SparseMatrixBase& matrix)
+    { return matrix*scalar; }
+
+
+    template<typename OtherDerived>
+    const typename SparseProductReturnType<Derived,OtherDerived>::Type
+    operator*(const SparseMatrixBase<OtherDerived> &other) const;
+
+    template<typename OtherDerived>
+    Derived& operator*=(const SparseMatrixBase<OtherDerived>& other);
+
+    template<typename OtherDerived>
+    typename ei_plain_matrix_type_column_major<OtherDerived>::type
+    solveTriangular(const MatrixBase<OtherDerived>& other) const;
+
+    template<typename OtherDerived>
+    void solveTriangularInPlace(MatrixBase<OtherDerived>& other) const;
+
+    template<typename OtherDerived> Scalar dot(const MatrixBase<OtherDerived>& other) const;
+    template<typename OtherDerived> Scalar dot(const SparseMatrixBase<OtherDerived>& other) const;
+    RealScalar squaredNorm() const;
+    RealScalar norm()  const;
+//     const PlainMatrixType normalized() const;
+//     void normalize();
+
+    SparseTranspose<Derived> transpose() { return derived(); }
+    const SparseTranspose<Derived> transpose() const { return derived(); }
+    // void transposeInPlace();
+    // const AdjointReturnType adjoint() 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;
+
+//     typename BlockReturnType<Derived>::Type block(int startRow, int startCol, int blockRows, int blockCols);
+//     const typename BlockReturnType<Derived>::Type
+//     block(int startRow, int startCol, int blockRows, int blockCols) const;
+//
+//     typename BlockReturnType<Derived>::SubVectorType segment(int start, int size);
+//     const typename BlockReturnType<Derived>::SubVectorType segment(int start, int size) const;
+//
+//     typename BlockReturnType<Derived,Dynamic>::SubVectorType start(int size);
+//     const typename BlockReturnType<Derived,Dynamic>::SubVectorType start(int size) const;
+//
+//     typename BlockReturnType<Derived,Dynamic>::SubVectorType end(int size);
+//     const typename BlockReturnType<Derived,Dynamic>::SubVectorType end(int size) const;
+//
+//     typename BlockReturnType<Derived>::Type corner(CornerType type, int cRows, int cCols);
+//     const typename BlockReturnType<Derived>::Type corner(CornerType type, int cRows, int cCols) const;
+//
+//     template<int BlockRows, int BlockCols>
+//     typename BlockReturnType<Derived, BlockRows, BlockCols>::Type block(int startRow, int startCol);
+//     template<int BlockRows, int BlockCols>
+//     const typename BlockReturnType<Derived, BlockRows, BlockCols>::Type block(int startRow, int startCol) const;
+
+//     template<int CRows, int CCols>
+//     typename BlockReturnType<Derived, CRows, CCols>::Type corner(CornerType type);
+//     template<int CRows, int CCols>
+//     const typename BlockReturnType<Derived, CRows, CCols>::Type corner(CornerType type) const;
+
+//     template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType start(void);
+//     template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType start() const;
+
+//     template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType end();
+//     template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType end() const;
+
+//     template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType segment(int start);
+//     template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType segment(int start) const;
+
+//     DiagonalCoeffs<Derived> diagonal();
+//     const DiagonalCoeffs<Derived> diagonal() const;
+
+//     template<unsigned int Mode> Part<Derived, Mode> part();
+//     template<unsigned int Mode> const Part<Derived, Mode> part() const;
+
+
+//     static const ConstantReturnType Constant(int rows, int cols, const Scalar& value);
+//     static const ConstantReturnType Constant(int size, const Scalar& value);
+//     static const ConstantReturnType Constant(const Scalar& value);
+
+//     template<typename CustomNullaryOp>
+//     static const CwiseNullaryOp<CustomNullaryOp, Derived> NullaryExpr(int rows, int cols, const CustomNullaryOp& func);
+//     template<typename CustomNullaryOp>
+//     static const CwiseNullaryOp<CustomNullaryOp, Derived> NullaryExpr(int size, const CustomNullaryOp& func);
+//     template<typename CustomNullaryOp>
+//     static const CwiseNullaryOp<CustomNullaryOp, Derived> NullaryExpr(const CustomNullaryOp& func);
+
+//     static const ConstantReturnType Zero(int rows, int cols);
+//     static const ConstantReturnType Zero(int size);
+//     static const ConstantReturnType Zero();
+//     static const ConstantReturnType Ones(int rows, int cols);
+//     static const ConstantReturnType Ones(int size);
+//     static const ConstantReturnType Ones();
+//     static const IdentityReturnType Identity();
+//     static const IdentityReturnType Identity(int rows, int cols);
+//     static const BasisReturnType Unit(int size, int i);
+//     static const BasisReturnType Unit(int i);
+//     static const BasisReturnType UnitX();
+//     static const BasisReturnType UnitY();
+//     static const BasisReturnType UnitZ();
+//     static const BasisReturnType UnitW();
+
+//     const DiagonalMatrix<Derived> asDiagonal() const;
+
+//     Derived& setConstant(const Scalar& value);
+//     Derived& setZero();
+//     Derived& setOnes();
+//     Derived& setRandom();
+//     Derived& setIdentity();
+
+      Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> toDense() const
+      {
+        Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> res(rows(),cols());
+        res.setZero();
+        for (int j=0; j<outerSize(); ++j)
+        {
+          for (typename Derived::InnerIterator i(derived(),j); i; ++i)
+            if(IsRowMajor)
+              res.coeffRef(j,i.index()) = i.value();
+            else
+              res.coeffRef(i.index(),j) = i.value();
+        }
+        return res;
+      }
+
+    template<typename OtherDerived>
+    bool isApprox(const SparseMatrixBase<OtherDerived>& other,
+                  RealScalar prec = precision<Scalar>()) const
+    { return toDense().isApprox(other.toDense(),prec); }
+
+    template<typename OtherDerived>
+    bool isApprox(const MatrixBase<OtherDerived>& other,
+                  RealScalar prec = precision<Scalar>()) const
+    { return toDense().isApprox(other,prec); }
+//     bool isMuchSmallerThan(const RealScalar& other,
+//                            RealScalar prec = precision<Scalar>()) const;
+//     template<typename OtherDerived>
+//     bool isMuchSmallerThan(const MatrixBase<OtherDerived>& other,
+//                            RealScalar prec = precision<Scalar>()) const;
+
+//     bool isApproxToConstant(const Scalar& value, RealScalar prec = precision<Scalar>()) const;
+//     bool isZero(RealScalar prec = precision<Scalar>()) const;
+//     bool isOnes(RealScalar prec = precision<Scalar>()) const;
+//     bool isIdentity(RealScalar prec = precision<Scalar>()) const;
+//     bool isDiagonal(RealScalar prec = precision<Scalar>()) const;
+
+//     bool isUpperTriangular(RealScalar prec = precision<Scalar>()) const;
+//     bool isLowerTriangular(RealScalar prec = precision<Scalar>()) const;
+
+//     template<typename OtherDerived>
+//     bool isOrthogonal(const MatrixBase<OtherDerived>& other,
+//                       RealScalar prec = precision<Scalar>()) const;
+//     bool isUnitary(RealScalar prec = precision<Scalar>()) const;
+
+//     template<typename OtherDerived>
+//     inline bool operator==(const MatrixBase<OtherDerived>& other) const
+//     { return (cwise() == other).all(); }
+
+//     template<typename OtherDerived>
+//     inline bool operator!=(const MatrixBase<OtherDerived>& other) const
+//     { return (cwise() != other).any(); }
+
+
+    template<typename NewType>
+    const SparseCwiseUnaryOp<ei_scalar_cast_op<typename ei_traits<Derived>::Scalar, NewType>, Derived> cast() const;
+
+    /** \returns the matrix or vector obtained by evaluating this expression.
+      *
+      * Notice that in the case of a plain matrix or vector (not an expression) this function just returns
+      * a const reference, in order to avoid a useless copy.
+      */
+    EIGEN_STRONG_INLINE const typename ei_eval<Derived>::type eval() const
+    { return typename ei_eval<Derived>::type(derived()); }
+
+//     template<typename OtherDerived>
+//     void swap(const MatrixBase<OtherDerived>& other);
+
+    template<unsigned int Added>
+    const SparseFlagged<Derived, Added, 0> marked() const;
+//     const Flagged<Derived, 0, EvalBeforeNestingBit | EvalBeforeAssigningBit> lazy() const;
+
+    /** \returns number of elements to skip to pass from one row (resp. column) to another
+      * for a row-major (resp. column-major) matrix.
+      * Combined with coeffRef() and the \ref flags flags, it allows a direct access to the data
+      * of the underlying matrix.
+      */
+//     inline int stride(void) const { return derived().stride(); }
+
+//     inline const NestByValue<Derived> nestByValue() const;
+
+
+    ConjugateReturnType conjugate() const;
+    const RealReturnType real() const;
+    const ImagReturnType imag() const;
+
+    template<typename CustomUnaryOp>
+    const SparseCwiseUnaryOp<CustomUnaryOp, Derived> unaryExpr(const CustomUnaryOp& func = CustomUnaryOp()) const;
+
+//     template<typename CustomBinaryOp, typename OtherDerived>
+//     const CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived>
+//     binaryExpr(const MatrixBase<OtherDerived> &other, const CustomBinaryOp& func = CustomBinaryOp()) const;
+
+
+    Scalar sum() const;
+//     Scalar trace() const;
+
+//     typename ei_traits<Derived>::Scalar minCoeff() const;
+//     typename ei_traits<Derived>::Scalar maxCoeff() const;
+
+//     typename ei_traits<Derived>::Scalar minCoeff(int* row, int* col = 0) const;
+//     typename ei_traits<Derived>::Scalar maxCoeff(int* row, int* col = 0) const;
+
+//     template<typename BinaryOp>
+//     typename ei_result_of<BinaryOp(typename ei_traits<Derived>::Scalar)>::type
+//     redux(const BinaryOp& func) const;
+
+//     template<typename Visitor>
+//     void visit(Visitor& func) const;
+
+
+    const SparseCwise<Derived> cwise() const;
+    SparseCwise<Derived> cwise();
+
+//     inline const WithFormat<Derived> format(const IOFormat& fmt) const;
+
+/////////// Array module ///////////
+    /*
+    bool all(void) const;
+    bool any(void) const;
+
+    const PartialRedux<Derived,Horizontal> rowwise() const;
+    const PartialRedux<Derived,Vertical> colwise() const;
+
+    static const CwiseNullaryOp<ei_scalar_random_op<Scalar>,Derived> Random(int rows, int cols);
+    static const CwiseNullaryOp<ei_scalar_random_op<Scalar>,Derived> Random(int size);
+    static const CwiseNullaryOp<ei_scalar_random_op<Scalar>,Derived> Random();
+
+    template<typename ThenDerived,typename ElseDerived>
+    const Select<Derived,ThenDerived,ElseDerived>
+    select(const MatrixBase<ThenDerived>& thenMatrix,
+           const MatrixBase<ElseDerived>& elseMatrix) const;
+
+    template<typename ThenDerived>
+    inline const Select<Derived,ThenDerived, NestByValue<typename ThenDerived::ConstantReturnType> >
+    select(const MatrixBase<ThenDerived>& thenMatrix, typename ThenDerived::Scalar elseScalar) const;
+
+    template<typename ElseDerived>
+    inline const Select<Derived, NestByValue<typename ElseDerived::ConstantReturnType>, ElseDerived >
+    select(typename ElseDerived::Scalar thenScalar, const MatrixBase<ElseDerived>& elseMatrix) const;
+
+    template<int p> RealScalar lpNorm() const;
+    */
+
+
 //     template<typename OtherDerived>
 //     Scalar dot(const MatrixBase<OtherDerived>& other) const
 //     {
diff --git a/Eigen/src/Sparse/SparseProduct.h b/Eigen/src/Sparse/SparseProduct.h
index 36dfb31ac..b4ba2ee6f 100644
--- a/Eigen/src/Sparse/SparseProduct.h
+++ b/Eigen/src/Sparse/SparseProduct.h
@@ -27,7 +27,7 @@
 
 // sparse product return type specialization
 template<typename Lhs, typename Rhs>
-struct ProductReturnType<Lhs,Rhs,SparseProduct>
+struct SparseProductReturnType
 {
   typedef typename ei_traits<Lhs>::Scalar Scalar;
   enum {
@@ -47,12 +47,11 @@ struct ProductReturnType<Lhs,Rhs,SparseProduct>
     SparseMatrix<Scalar,0>,
     const typename ei_nested<Rhs,Lhs::RowsAtCompileTime>::type>::ret RhsNested;
 
-  typedef Product<LhsNested,
-                  RhsNested, SparseProduct> Type;
+  typedef SparseProduct<LhsNested, RhsNested> Type;
 };
 
 template<typename LhsNested, typename RhsNested>
-struct ei_traits<Product<LhsNested, RhsNested, SparseProduct> >
+struct ei_traits<SparseProduct<LhsNested, RhsNested> >
 {
   // clean the nested types:
   typedef typename ei_cleantype<LhsNested>::type _LhsNested;
@@ -87,30 +86,28 @@ struct ei_traits<Product<LhsNested, RhsNested, SparseProduct> >
   };
 };
 
-template<typename LhsNested, typename RhsNested> class Product<LhsNested,RhsNested,SparseProduct> : ei_no_assignment_operator,
-  public MatrixBase<Product<LhsNested, RhsNested, SparseProduct> >
+template<typename LhsNested, typename RhsNested>
+class SparseProduct : ei_no_assignment_operator,
+  public SparseMatrixBase<SparseProduct<LhsNested, RhsNested> >
 {
   public:
 
-    EIGEN_GENERIC_PUBLIC_INTERFACE(Product)
+    EIGEN_GENERIC_PUBLIC_INTERFACE(SparseProduct)
 
   private:
 
-    typedef typename ei_traits<Product>::_LhsNested _LhsNested;
-    typedef typename ei_traits<Product>::_RhsNested _RhsNested;
+    typedef typename ei_traits<SparseProduct>::_LhsNested _LhsNested;
+    typedef typename ei_traits<SparseProduct>::_RhsNested _RhsNested;
 
   public:
 
     template<typename Lhs, typename Rhs>
-    inline Product(const Lhs& lhs, const Rhs& rhs)
+    inline SparseProduct(const Lhs& lhs, const Rhs& rhs)
       : m_lhs(lhs), m_rhs(rhs)
     {
       ei_assert(lhs.cols() == rhs.rows());
     }
 
-    Scalar coeff(int, int) const { ei_assert(false && "eigen internal error"); }
-    Scalar& coeffRef(int, int) { ei_assert(false && "eigen internal error"); }
-
     inline int rows() const { return m_lhs.rows(); }
     inline int cols() const { return m_rhs.cols(); }
 
@@ -231,28 +228,51 @@ struct ei_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,RowMajor,ColMajor>
 // by ProductReturnType which transform it to (col col *) by evaluating rhs.
 
 
+// template<typename Derived>
+// template<typename Lhs, typename Rhs>
+// inline Derived& SparseMatrixBase<Derived>::lazyAssign(const SparseProduct<Lhs,Rhs>& product)
+// {
+// //   std::cout << "sparse product to dense\n";
+//   ei_sparse_product_selector<
+//     typename ei_cleantype<Lhs>::type,
+//     typename ei_cleantype<Rhs>::type,
+//     typename ei_cleantype<Derived>::type>::run(product.lhs(),product.rhs(),derived());
+//   return derived();
+// }
+
 template<typename Derived>
 template<typename Lhs, typename Rhs>
-inline Derived& MatrixBase<Derived>::lazyAssign(const Product<Lhs,Rhs,SparseProduct>& product)
+inline Derived& SparseMatrixBase<Derived>::operator=(const SparseProduct<Lhs,Rhs>& product)
 {
-//   std::cout << "sparse product to dense\n";
+//   std::cout << "sparse product to sparse\n";
   ei_sparse_product_selector<
     typename ei_cleantype<Lhs>::type,
     typename ei_cleantype<Rhs>::type,
-    typename ei_cleantype<Derived>::type>::run(product.lhs(),product.rhs(),derived());
+    Derived>::run(product.lhs(),product.rhs(),derived());
   return derived();
 }
 
 template<typename Derived>
-template<typename Lhs, typename Rhs>
-inline Derived& SparseMatrixBase<Derived>::operator=(const Product<Lhs,Rhs,SparseProduct>& product)
+template<typename OtherDerived>
+inline const typename SparseProductReturnType<Derived,OtherDerived>::Type
+SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other) const
 {
-//   std::cout << "sparse product to sparse\n";
-  ei_sparse_product_selector<
-    typename ei_cleantype<Lhs>::type,
-    typename ei_cleantype<Rhs>::type,
-    Derived>::run(product.lhs(),product.rhs(),derived());
-  return derived();
+  enum {
+    ProductIsValid =  Derived::ColsAtCompileTime==Dynamic
+                   || OtherDerived::RowsAtCompileTime==Dynamic
+                   || int(Derived::ColsAtCompileTime)==int(OtherDerived::RowsAtCompileTime),
+    AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
+    SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived,OtherDerived)
+  };
+  // note to the lost user:
+  //    * for a dot product use: v1.dot(v2)
+  //    * for a coeff-wise product use: v1.cwise()*v2
+  EIGEN_STATIC_ASSERT(ProductIsValid || !(AreVectors && SameSizes),
+    INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
+  EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
+    INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
+  EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
+  return typename SparseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
 }
 
 #endif // EIGEN_SPARSEPRODUCT_H
diff --git a/Eigen/src/Sparse/SparseRedux.h b/Eigen/src/Sparse/SparseRedux.h
index 746ae2f8a..f0d370548 100644
--- a/Eigen/src/Sparse/SparseRedux.h
+++ b/Eigen/src/Sparse/SparseRedux.h
@@ -25,46 +25,16 @@
 #ifndef EIGEN_SPARSEREDUX_H
 #define EIGEN_SPARSEREDUX_H
 
-
-template<typename Derived, int Vectorization, int Unrolling>
-struct ei_sum_impl<Derived, Vectorization, Unrolling, IsSparse>
-{
-  typedef typename Derived::Scalar Scalar;
-  static Scalar run(const Derived& mat)
-  {
-    ei_assert(mat.rows()>0 && mat.cols()>0 && "you are using a non initialized matrix");
-    Scalar res = 0;
-    for (int j=0; j<mat.outerSize(); ++j)
-      for (typename Derived::InnerIterator iter(mat,j); iter; ++iter)
-        res += iter.value();
-    return res;
-  }
-};
-
-template<typename Derived1, typename Derived2, int Vectorization, int Unrolling>
-struct ei_dot_impl<Derived1, Derived2, Vectorization, Unrolling, IsSparse>
+template<typename Derived>
+typename ei_traits<Derived>::Scalar
+SparseMatrixBase<Derived>::sum() const
 {
-  typedef typename Derived1::Scalar Scalar;
-  static Scalar run(const Derived1& v1, const Derived2& v2)
-  {
-    ei_assert(v1.size()>0 && "you are using a non initialized vector");
-    typename Derived1::InnerIterator i(v1,0);
-    typename Derived2::InnerIterator j(v2,0);
-    Scalar res = 0;
-    while (i && j)
-    {
-      if (i.index()==j.index())
-      {
-        res += i.value() * ei_conj(j.value());
-        ++i; ++j;
-      }
-      else if (i.index()<j.index())
-        ++i;
-      else
-        ++j;
-    }
-    return res;
-  }
-};
+  ei_assert(rows()>0 && cols()>0 && "you are using a non initialized matrix");
+  Scalar res = 0;
+  for (int j=0; j<outerSize(); ++j)
+    for (typename Derived::InnerIterator iter(derived(),j); iter; ++iter)
+      res += iter.value();
+  return res;
+}
 
 #endif // EIGEN_SPARSEREDUX_H
diff --git a/Eigen/src/Sparse/SparseTranspose.h b/Eigen/src/Sparse/SparseTranspose.h
new file mode 100644
index 000000000..89a14d707
--- /dev/null
+++ b/Eigen/src/Sparse/SparseTranspose.h
@@ -0,0 +1,85 @@
+// 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_SPARSETRANSPOSE_H
+#define EIGEN_SPARSETRANSPOSE_H
+
+template<typename MatrixType>
+struct ei_traits<SparseTranspose<MatrixType> > : ei_traits<Transpose<MatrixType> >
+{};
+
+template<typename MatrixType> class SparseTranspose
+  : public SparseMatrixBase<SparseTranspose<MatrixType> >
+{
+  public:
+
+    EIGEN_GENERIC_PUBLIC_INTERFACE(SparseTranspose)
+
+    class InnerIterator;
+    class ReverseInnerIterator;
+
+    inline SparseTranspose(const MatrixType& matrix) : m_matrix(matrix) {}
+
+    //EIGEN_INHERIT_ASSIGNMENT_OPERATORS(SparseTranspose)
+
+    inline int rows() const { return m_matrix.cols(); }
+    inline int cols() const { return m_matrix.rows(); }
+    inline int nonZeros() const { return m_matrix.nonZeros(); }
+
+    // FIXME should be keep them ?
+    inline Scalar& coeffRef(int row, int col)
+    { return m_matrix.const_cast_derived().coeffRef(col, row); }
+
+    inline const Scalar coeff(int row, int col) const
+    { return m_matrix.coeff(col, row); }
+
+    inline const Scalar coeff(int index) const
+    { return m_matrix.coeff(index); }
+
+    inline Scalar& coeffRef(int index)
+    { return m_matrix.const_cast_derived().coeffRef(index); }
+
+  protected:
+    const typename MatrixType::Nested m_matrix;
+};
+
+template<typename MatrixType> class SparseTranspose<MatrixType>::InnerIterator : public MatrixType::InnerIterator
+{
+  public:
+
+    EIGEN_STRONG_INLINE InnerIterator(const SparseTranspose& trans, int outer)
+      : MatrixType::InnerIterator(trans.m_matrix, outer)
+    {}
+};
+
+template<typename MatrixType> class SparseTranspose<MatrixType>::ReverseInnerIterator : public MatrixType::ReverseInnerIterator
+{
+  public:
+
+    EIGEN_STRONG_INLINE ReverseInnerIterator(const SparseTranspose& xpr, int outer)
+      : MatrixType::ReverseInnerIterator(xpr.m_matrix, outer)
+    {}
+};
+
+#endif // EIGEN_SPARSETRANSPOSE_H
diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h
index dcafcc008..724fb9efb 100644
--- a/Eigen/src/Sparse/SparseUtil.h
+++ b/Eigen/src/Sparse/SparseUtil.h
@@ -31,6 +31,46 @@
 #define EIGEN_DBG_SPARSE(X) X
 #endif
 
+#define EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(Derived, Op) \
+template<typename OtherDerived> \
+EIGEN_STRONG_INLINE Derived& operator Op(const Eigen::SparseMatrixBase<OtherDerived>& other) \
+{ \
+  return Base::operator Op(other.derived()); \
+} \
+EIGEN_STRONG_INLINE Derived& operator Op(const Derived& other) \
+{ \
+  return Base::operator Op(other); \
+}
+
+#define EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, Op) \
+template<typename Other> \
+EIGEN_STRONG_INLINE Derived& operator Op(const Other& scalar) \
+{ \
+  return Base::operator Op(scalar); \
+}
+
+#define EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATORS(Derived) \
+EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(Derived, =) \
+EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(Derived, +=) \
+EIGEN_SPARSE_INHERIT_ASSIGNMENT_OPERATOR(Derived, -=) \
+EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, *=) \
+EIGEN_SPARSE_INHERIT_SCALAR_ASSIGNMENT_OPERATOR(Derived, /=)
+
+#define _EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(Derived, BaseClass) \
+typedef BaseClass Base; \
+typedef typename Eigen::ei_traits<Derived>::Scalar Scalar; \
+typedef typename Eigen::NumTraits<Scalar>::Real RealScalar; \
+typedef typename Eigen::ei_nested<Derived>::type Nested; \
+enum { RowsAtCompileTime = Eigen::ei_traits<Derived>::RowsAtCompileTime, \
+       ColsAtCompileTime = Eigen::ei_traits<Derived>::ColsAtCompileTime, \
+       Flags = Eigen::ei_traits<Derived>::Flags, \
+       CoeffReadCost = Eigen::ei_traits<Derived>::CoeffReadCost, \
+       SizeAtCompileTime = Base::SizeAtCompileTime, \
+       IsVectorAtCompileTime = Base::IsVectorAtCompileTime };
+
+#define EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(Derived) \
+_EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(Derived, Eigen::SparseMatrixBase<Derived>)
+
 enum SparseBackend {
   DefaultBackend,
   Taucs,
@@ -64,6 +104,16 @@ template<typename Derived> class SparseMatrixBase;
 template<typename _Scalar, int _Flags = 0> class SparseMatrix;
 template<typename _Scalar, int _Flags = 0> class SparseVector;
 
+template<typename MatrixType> class SparseTranspose;
+template<typename MatrixType> class SparseInnerVector;
+template<typename Derived> class SparseCwise;
+template<typename UnaryOp,   typename MatrixType>         class SparseCwiseUnaryOp;
+template<typename BinaryOp,  typename Lhs, typename Rhs>  class SparseCwiseBinaryOp;
+template<typename Lhs, typename Rhs>                      class SparseProduct;
+template<typename ExpressionType, unsigned int Added, unsigned int Removed> class SparseFlagged;
+
+template<typename Lhs, typename Rhs> struct SparseProductReturnType;
+
 const int AccessPatternNotSupported = 0x0;
 const int AccessPatternSupported    = 0x1;
 
diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h
index 7a2ce70e5..b189cb5ee 100644
--- a/Eigen/src/Sparse/SparseVector.h
+++ b/Eigen/src/Sparse/SparseVector.h
@@ -58,7 +58,7 @@ class SparseVector
   : public SparseMatrixBase<SparseVector<_Scalar, _Flags> >
 {
   public:
-    EIGEN_GENERIC_PUBLIC_INTERFACE(SparseVector)
+    EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseVector)
 
   protected:
   public:
diff --git a/Eigen/src/Sparse/TriangularSolver.h b/Eigen/src/Sparse/TriangularSolver.h
index d44938273..8948ae45e 100644
--- a/Eigen/src/Sparse/TriangularSolver.h
+++ b/Eigen/src/Sparse/TriangularSolver.h
@@ -144,4 +144,35 @@ struct ei_solve_triangular_selector<Lhs,Rhs,UpperTriangular,ColMajor|IsSparse>
   }
 };
 
+template<typename Derived>
+template<typename OtherDerived>
+void SparseMatrixBase<Derived>::solveTriangularInPlace(MatrixBase<OtherDerived>& other) const
+{
+  ei_assert(derived().cols() == derived().rows());
+  ei_assert(derived().cols() == other.rows());
+  ei_assert(!(Flags & ZeroDiagBit));
+  ei_assert(Flags & (UpperTriangularBit|LowerTriangularBit));
+
+  enum { copy = ei_traits<OtherDerived>::Flags & RowMajorBit };
+
+  typedef typename ei_meta_if<copy,
+    typename ei_plain_matrix_type_column_major<OtherDerived>::type, OtherDerived&>::ret OtherCopy;
+  OtherCopy otherCopy(other.derived());
+
+  ei_solve_triangular_selector<Derived, typename ei_unref<OtherCopy>::type>::run(derived(), otherCopy);
+
+  if (copy)
+    other = otherCopy;
+}
+
+template<typename Derived>
+template<typename OtherDerived>
+typename ei_plain_matrix_type_column_major<OtherDerived>::type
+SparseMatrixBase<Derived>::solveTriangular(const MatrixBase<OtherDerived>& other) const
+{
+  typename ei_plain_matrix_type_column_major<OtherDerived>::type res(other);
+  solveTriangularInPlace(res);
+  return res;
+}
+
 #endif // EIGEN_SPARSETRIANGULARSOLVER_H
diff --git a/cmake/FindSuperLU.cmake b/cmake/FindSuperLU.cmake
index e10561b60..fbefc943a 100644
--- a/cmake/FindSuperLU.cmake
+++ b/cmake/FindSuperLU.cmake
@@ -16,10 +16,6 @@ if(BLAS_FOUND)
   )
 
   find_library(SUPERLU_LIBRARIES superlu PATHS $ENV{SUPERLUDIR} ${LIB_INSTALL_DIR})
-
-  if(SUPERLU_LIBRARIES AND CMAKE_COMPILER_IS_GNUCXX)
-    set(SUPERLU_LIBRARIES ${SUPERLU_LIBRARIES} -lgfortran)
-  endif(SUPERLU_LIBRARIES AND CMAKE_COMPILER_IS_GNUCXX)
   
   if(SUPERLU_LIBRARIES)
     set(SUPERLU_LIBRARIES ${SUPERLU_LIBRARIES} ${BLAS_LIBRARIES})
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index bc9a5ca73..021799f5e 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -195,8 +195,8 @@ ei_add_test(parametrizedline)
 ei_add_test(alignedbox)
 ei_add_test(regression)
 ei_add_test(stdvector)
-ei_add_test(sparse_basic)
 ei_add_test(sparse_vector)
+ei_add_test(sparse_basic)
 ei_add_test(sparse_solvers " " "${SPARSE_LIBS}")
 
 # print a summary of the different options
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp
index c50682810..54272d871 100644
--- a/test/sparse_basic.cpp
+++ b/test/sparse_basic.cpp
@@ -72,7 +72,7 @@ template<typename Scalar> void sparse_basic(int rows, int cols)
   refMat.coeffRef(nonzeroCoords[0].x(), nonzeroCoords[0].y()) = Scalar(5);
 
   VERIFY_IS_APPROX(m, refMat);
-/*
+  /*
   // test InnerIterators and Block expressions
   for (int t=0; t<10; ++t)
   {
@@ -81,23 +81,23 @@ template<typename Scalar> void sparse_basic(int rows, int cols)
     int w = ei_random<int>(1,cols-j-1);
     int h = ei_random<int>(1,rows-i-1);
 
-    VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w));
+//     VERIFY_IS_APPROX(m.block(i,j,h,w), refMat.block(i,j,h,w));
     for(int c=0; c<w; c++)
     {
       VERIFY_IS_APPROX(m.block(i,j,h,w).col(c), refMat.block(i,j,h,w).col(c));
       for(int r=0; r<h; r++)
       {
-        VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r));
-      }
-    }
-    for(int r=0; r<h; r++)
-    {
-      VERIFY_IS_APPROX(m.block(i,j,h,w).row(r), refMat.block(i,j,h,w).row(r));
-      for(int c=0; c<w; c++)
-      {
-        VERIFY_IS_APPROX(m.block(i,j,h,w).row(r).coeff(c), refMat.block(i,j,h,w).row(r).coeff(c));
+//         VERIFY_IS_APPROX(m.block(i,j,h,w).col(c).coeff(r), refMat.block(i,j,h,w).col(c).coeff(r));
       }
     }
+//     for(int r=0; r<h; r++)
+//     {
+//       VERIFY_IS_APPROX(m.block(i,j,h,w).row(r), refMat.block(i,j,h,w).row(r));
+//       for(int c=0; c<w; c++)
+//       {
+//         VERIFY_IS_APPROX(m.block(i,j,h,w).row(r).coeff(c), refMat.block(i,j,h,w).row(r).coeff(c));
+//       }
+//     }
   }
 
   for(int c=0; c<cols; c++)
@@ -171,7 +171,7 @@ template<typename Scalar> void sparse_basic(int rows, int cols)
       }
       m2.endFill();
       std::cerr << m1 << "\n\n" << m2 << "\n";
-      VERIFY_IS_APPROX(m1,m2);
+      VERIFY_IS_APPROX(m2,m1);
     }
 //   {
 //     m.setZero();
@@ -191,6 +191,17 @@ template<typename Scalar> void sparse_basic(int rows, int cols)
 //   std::cerr << m.transpose() << "\n\n"  << refMat.transpose() << "\n\n";
 //   VERIFY_IS_APPROX(m, refMat);
 
+  // test innerVector()
+  {
+    DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
+    SparseMatrix<Scalar> m2(rows, rows);
+    initSparse<Scalar>(density, refMat2, m2);
+    int j0 = ei_random(0,rows-1);
+    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));
+  }
+
   // test transpose
   {
     DenseMatrix refMat2 = DenseMatrix::Zero(rows, rows);
diff --git a/test/sparse_vector.cpp b/test/sparse_vector.cpp
index 0a25884ca..18ef74833 100644
--- a/test/sparse_vector.cpp
+++ b/test/sparse_vector.cpp
@@ -62,7 +62,8 @@ template<typename Scalar> void sparse_vector(int rows, int cols)
     for (typename SparseVectorType::InnerIterator it(v1); it; ++it,++j)
     {
       VERIFY(nonzerocoords[j]==it.index());
-      VERIFY(it.value()==v1[it.index()]);
+      VERIFY(it.value()==v1.coeff(it.index()));
+      VERIFY(it.value()==refV1.coeff(it.index()));
     }
   }
   VERIFY_IS_APPROX(v1, refV1);
@@ -76,7 +77,7 @@ template<typename Scalar> void sparse_vector(int rows, int cols)
 
   VERIFY_IS_APPROX(v1*s1-v2, refV1*s1-refV2);
 
-  std::cerr << v1.dot(v2) << " == " << refV1.dot(refV2) << "\n";
+//   std::cerr << v1.dot(v2) << " == " << refV1.dot(refV2) << "\n";
   VERIFY_IS_APPROX(v1.dot(v2), refV1.dot(refV2));
 
 }
@@ -85,7 +86,7 @@ void test_sparse_vector()
 {
   for(int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST( sparse_vector<double>(8, 8) );
-//     CALL_SUBTEST( sparse_vector<std::complex<double> >(16, 16) );
+    CALL_SUBTEST( sparse_vector<std::complex<double> >(16, 16) );
     CALL_SUBTEST( sparse_vector<double>(299, 535) );
   }
 }