* added an in-place version of inverseProduct which

  might be twice faster fot small fixed size matrix
* added a sparse triangular solver (sparse version
  of inverseProduct)
* various other improvements in the Sparse module

12 files changed, 572 insertions, 199 deletions

diff --git a/Eigen/Sparse b/Eigen/Sparse
index 089a77af5..12ffc0c66 100644
--- a/Eigen/Sparse
+++ b/Eigen/Sparse
@@ -18,6 +18,7 @@ namespace Eigen {
 #include "src/Sparse/CoreIterators.h"
 #include "src/Sparse/SparseSetter.h"
 #include "src/Sparse/SparseProduct.h"
+#include "src/Sparse/TriangularSolver.h"
 
 } // namespace Eigen
 
diff --git a/Eigen/src/Core/InverseProduct.h b/Eigen/src/Core/InverseProduct.h
index 40496e01d..57dbf7509 100755
--- a/Eigen/src/Core/InverseProduct.h
+++ b/Eigen/src/Core/InverseProduct.h
@@ -25,66 +25,73 @@
 #ifndef EIGEN_INVERSEPRODUCT_H
 #define EIGEN_INVERSEPRODUCT_H
 
-/** \returns the product of the inverse of \c *this with \a other.
-  *
-  * This function computes the inverse-matrix matrix product inverse(\c*this) * \a other
-  * It works as a forward (resp. backward) substitution if \c *this is an upper (resp. lower)
-  * triangular matrix.
-  *
-  * It is required that \c *this be marked as either an upper or a lower triangular matrix, as
-  * can be done by marked(), and as is automatically the case with expressions such as those returned
-  * by extract().
-  * Example: \include MatrixBase_marked.cpp
-  * Output: \verbinclude MatrixBase_marked.out
+
+/** "in-place" version of MatrixBase::inverseProduct() where the result is written in \a other
   *
-  * \sa marked(), extract()
+  * \sa inverseProduct()
   */
 template<typename Derived>
 template<typename OtherDerived>
-typename OtherDerived::Eval MatrixBase<Derived>::inverseProduct(const MatrixBase<OtherDerived>& other) const
+void MatrixBase<Derived>::inverseProductInPlace(MatrixBase<OtherDerived>& other) const
 {
-  assert(cols() == other.rows());
-  assert(!(Flags & ZeroDiagBit));
-  assert(Flags & (UpperTriangularBit|LowerTriangularBit));
-
-  typename OtherDerived::Eval res(other.rows(), other.cols());
+  ei_assert(cols() == other.rows());
+  ei_assert(!(Flags & ZeroDiagBit));
+  ei_assert(Flags & (UpperTriangularBit|LowerTriangularBit));
 
   for(int c=0 ; c<other.cols() ; ++c)
   {
     if(Flags & LowerTriangularBit)
     {
       // forward substitution
-      if(Flags & UnitDiagBit)
-        res.coeffRef(0,c) = other.coeff(0,c);
-      else
-        res.coeffRef(0,c) = other.coeff(0,c)/coeff(0, 0);
+      if(!(Flags & UnitDiagBit))
+        other.coeffRef(0,c) = other.coeff(0,c)/coeff(0, 0);
       for(int i=1; i<rows(); ++i)
       {
-        Scalar tmp = other.coeff(i,c) - ((this->row(i).start(i)) * res.col(c).start(i)).coeff(0,0);
+        Scalar tmp = other.coeff(i,c) - ((this->row(i).start(i)) * other.col(c).start(i)).coeff(0,0);
         if (Flags & UnitDiagBit)
-          res.coeffRef(i,c) = tmp;
+          other.coeffRef(i,c) = tmp;
         else
-          res.coeffRef(i,c) = tmp/coeff(i,i);
+          other.coeffRef(i,c) = tmp/coeff(i,i);
       }
     }
     else
     {
       // backward substitution
-      if(Flags & UnitDiagBit)
-        res.coeffRef(cols()-1,c) = other.coeff(cols()-1,c);
-      else
-        res.coeffRef(cols()-1,c) = other.coeff(cols()-1, c)/coeff(rows()-1, cols()-1);
+      if(!(Flags & UnitDiagBit))
+        other.coeffRef(cols()-1,c) = other.coeff(cols()-1, c)/coeff(rows()-1, cols()-1);
       for(int i=rows()-2 ; i>=0 ; --i)
       {
         Scalar tmp = other.coeff(i,c)
-                   - ((this->row(i).end(cols()-i-1)) * res.col(c).end(cols()-i-1)).coeff(0,0);
+                   - ((this->row(i).end(cols()-i-1)) * other.col(c).end(cols()-i-1)).coeff(0,0);
         if (Flags & UnitDiagBit)
-          res.coeffRef(i,c) = tmp;
+          other.coeffRef(i,c) = tmp;
         else
-          res.coeffRef(i,c) = tmp/coeff(i,i);
+          other.coeffRef(i,c) = tmp/coeff(i,i);
       }
     }
   }
+}
+
+/** \returns the product of the inverse of \c *this with \a other.
+  *
+  * This function computes the inverse-matrix matrix product inverse(\c*this) * \a other
+  * It works as a forward (resp. backward) substitution if \c *this is an upper (resp. lower)
+  * triangular matrix.
+  *
+  * It is required that \c *this be marked as either an upper or a lower triangular matrix, as
+  * can be done by marked(), and as is automatically the case with expressions such as those returned
+  * by extract().
+  * Example: \include MatrixBase_marked.cpp
+  * Output: \verbinclude MatrixBase_marked.out
+  *
+  * \sa marked(), extract()
+  */
+template<typename Derived>
+template<typename OtherDerived>
+typename OtherDerived::Eval MatrixBase<Derived>::inverseProduct(const MatrixBase<OtherDerived>& other) const
+{
+  typename OtherDerived::Eval res(other);
+  inverseProductInPlace(res);
   return res;
 }
 
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index eeb40c2f5..63817102f 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -296,6 +296,9 @@ template<typename Derived> class MatrixBase
     template<typename OtherDerived>
     typename OtherDerived::Eval inverseProduct(const MatrixBase<OtherDerived>& other) const;
 
+    template<typename OtherDerived>
+    void inverseProductInPlace(MatrixBase<OtherDerived>& other) const;
+
 
     template<typename OtherDerived>
     Scalar dot(const MatrixBase<OtherDerived>& other) const;
diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h
index 7e0c5650d..613cb053c 100644
--- a/Eigen/src/Core/util/Constants.h
+++ b/Eigen/src/Core/util/Constants.h
@@ -191,6 +191,11 @@ enum {
   Sparse  = SparseBit
 };
 
+enum {
+  ColMajor = 0,
+  RowMajor = RowMajorBit
+};
+
 const int FullyCoherentAccessPattern  = 0x1;
 const int InnerCoherentAccessPattern  = 0x2 | FullyCoherentAccessPattern;
 const int OuterCoherentAccessPattern  = 0x4 | InnerCoherentAccessPattern;
diff --git a/Eigen/src/Sparse/LinkedVectorMatrix.h b/Eigen/src/Sparse/LinkedVectorMatrix.h
index 76e3a4c12..969aee57b 100644
--- a/Eigen/src/Sparse/LinkedVectorMatrix.h
+++ b/Eigen/src/Sparse/LinkedVectorMatrix.h
@@ -40,15 +40,15 @@ struct ei_traits<LinkedVectorMatrix<_Scalar,_Flags> >
   };
 };
 
-template<typename Element, int BlockSize = 8>
-struct LinkedVector
+template<typename Element, int ChunkSize = 8>
+struct LinkedVectorChunk
 {
-  LinkedVector() : size(0), next(0), prev(0) {}
-  Element data[BlockSize];
-  LinkedVector* next;
-  LinkedVector* prev;
+  LinkedVectorChunk() : size(0), next(0), prev(0) {}
+  Element data[ChunkSize];
+  LinkedVectorChunk* next;
+  LinkedVectorChunk* prev;
   int size;
-  bool isFull() const { return size==BlockSize; }
+  bool isFull() const { return size==ChunkSize; }
 };
 
 template<typename _Scalar, int _Flags>
@@ -70,11 +70,11 @@ class LinkedVectorMatrix : public SparseMatrixBase<LinkedVectorMatrix<_Scalar,_F
       Scalar value;
       int index;
     };
-    typedef LinkedVector<ValueIndex,8> LinkedVectorBlock;
+    typedef LinkedVectorChunk<ValueIndex,8> VectorChunk;
 
-    inline int find(LinkedVectorBlock** _el, int id)
+    inline int find(VectorChunk** _el, int id)
     {
-      LinkedVectorBlock* el = *_el;
+      VectorChunk* el = *_el;
       while (el && el->data[el->size-1].index<id)
         el = el->next;
       *_el = el;
@@ -115,7 +115,7 @@ class LinkedVectorMatrix : public SparseMatrixBase<LinkedVectorMatrix<_Scalar,_F
       const int outer = RowMajor ? row : col;
       const int inner = RowMajor ? col : row;
 
-      LinkedVectorBlock* el = m_data[outer];
+      VectorChunk* el = m_data[outer];
       int id = find(&el, inner);
       if (id<0)
         return Scalar(0);
@@ -127,7 +127,7 @@ class LinkedVectorMatrix : public SparseMatrixBase<LinkedVectorMatrix<_Scalar,_F
       const int outer = RowMajor ? row : col;
       const int inner = RowMajor ? col : row;
 
-      LinkedVectorBlock* el = m_data[outer];
+      VectorChunk* el = m_data[outer];
       int id = find(&el, inner);
       ei_assert(id>=0);
 //       if (id<0)
@@ -150,7 +150,7 @@ class LinkedVectorMatrix : public SparseMatrixBase<LinkedVectorMatrix<_Scalar,_F
       const int inner = RowMajor ? col : row;
       if (m_ends[outer]==0)
       {
-        m_data[outer] = m_ends[outer] = new LinkedVectorBlock();
+        m_data[outer] = m_ends[outer] = new VectorChunk();
       }
       else
       {
@@ -158,7 +158,7 @@ class LinkedVectorMatrix : public SparseMatrixBase<LinkedVectorMatrix<_Scalar,_F
         if (m_ends[outer]->isFull())
         {
 
-          LinkedVectorBlock* el = new LinkedVectorBlock();
+          VectorChunk* el = new VectorChunk();
           m_ends[outer]->next = el;
           el->prev = m_ends[outer];
           m_ends[outer] = el;
@@ -168,24 +168,21 @@ class LinkedVectorMatrix : public SparseMatrixBase<LinkedVectorMatrix<_Scalar,_F
       return m_ends[outer]->data[m_ends[outer]->size++].value;
     }
 
-    inline void endFill()
-    {
-    }
+    inline void endFill() { }
 
     ~LinkedVectorMatrix()
     {
-      if (this->isNotShared())
-        clear();
+      clear();
     }
 
     void clear()
     {
       for (int i=0; i<m_data.size(); ++i)
       {
-        LinkedVectorBlock* el = m_data[i];
+        VectorChunk* el = m_data[i];
         while (el)
         {
-          LinkedVectorBlock* tmp = el;
+          VectorChunk* tmp = el;
           el = el->next;
           delete tmp;
         }
@@ -221,30 +218,26 @@ class LinkedVectorMatrix : public SparseMatrixBase<LinkedVectorMatrix<_Scalar,_F
       *this = other.derived();
     }
 
-    inline void shallowCopy(const LinkedVectorMatrix& other)
+    inline void swap(LinkedVectorMatrix& other)
     {
-      EIGEN_DBG_SPARSE(std::cout << "LinkedVectorMatrix:: shallowCopy\n");
+      EIGEN_DBG_SPARSE(std::cout << "LinkedVectorMatrix:: swap\n");
       resize(other.rows(), other.cols());
-      for (int j=0; j<this->outerSize(); ++j)
-      {
-        m_data[j] = other.m_data[j];
-        m_ends[j] = other.m_ends[j];
-      }
-      other.markAsCopied();
+      m_data.swap(other.m_data);
+      m_ends.swap(other.m_ends);
     }
 
     inline LinkedVectorMatrix& operator=(const LinkedVectorMatrix& other)
     {
       if (other.isRValue())
       {
-        shallowCopy(other);
-        return *this;
+        swap(other.const_cast_derived());
       }
       else
       {
         // TODO implement a specialized deep copy here
         return operator=<LinkedVectorMatrix>(other);
       }
+      return *this;
     }
 
     template<typename OtherDerived>
@@ -255,8 +248,10 @@ class LinkedVectorMatrix : public SparseMatrixBase<LinkedVectorMatrix<_Scalar,_F
 
   protected:
 
-    std::vector<LinkedVectorBlock*> m_data;
-    std::vector<LinkedVectorBlock*> m_ends;
+    // outer vector of inner linked vector chunks
+    std::vector<VectorChunk*> m_data;
+    // stores a reference to the last vector chunk for efficient filling
+    std::vector<VectorChunk*> m_ends;
     int m_innerSize;
 
 };
@@ -281,7 +276,7 @@ class LinkedVectorMatrix<Scalar,_Flags>::InnerIterator
 
   protected:
     const LinkedVectorMatrix& m_matrix;
-    LinkedVectorBlock* m_el;
+    VectorChunk* m_el;
     int m_it;
 };
 
diff --git a/Eigen/src/Sparse/SparseArray.h b/Eigen/src/Sparse/SparseArray.h
index 3ac181e04..bd82e6e9a 100644
--- a/Eigen/src/Sparse/SparseArray.h
+++ b/Eigen/src/Sparse/SparseArray.h
@@ -32,11 +32,11 @@ template<typename Scalar> class SparseArray
 {
   public:
     SparseArray()
-      : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0), m_isShared(0)
+      : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0)
     {}
 
     SparseArray(int size)
-      : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0), m_isShared(0)
+      : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0)
     {
       resize(size);
     }
@@ -51,66 +51,39 @@ template<typename Scalar> class SparseArray
       resize(other.size());
       memcpy(m_values, other.m_values, m_size * sizeof(Scalar));
       memcpy(m_indices, other.m_indices, m_size * sizeof(int));
-      m_isShared = 0;
     }
 
-    void shallowCopy(const SparseArray& other)
+    void swap(SparseArray& other)
     {
-      delete[] m_values;
-      delete[] m_indices;
-      m_values = other.m_values;
-      m_indices = other.m_indices;
-      m_size = other.m_size;
-      m_allocatedSize = other.m_allocatedSize;
-      m_isShared = false;
-      other.m_isShared = true;
+      std::swap(m_values, other.m_values);
+      std::swap(m_indices, other.m_indices);
+      std::swap(m_size, other.m_size);
+      std::swap(m_allocatedSize, other.m_allocatedSize);
     }
 
     ~SparseArray()
     {
-      if (!m_isShared)
-      {
-        delete[] m_values;
-        delete[] m_indices;
-      }
+      delete[] m_values;
+      delete[] m_indices;
     }
 
     void reserve(int size)
     {
       int newAllocatedSize = m_size + size;
       if (newAllocatedSize > m_allocatedSize)
-      {
-        Scalar* newValues  = new Scalar[newAllocatedSize];
-        int* newIndices = new int[newAllocatedSize];
-        // copy
-        memcpy(newValues,  m_values,  m_size * sizeof(Scalar));
-        memcpy(newIndices, m_indices, m_size * sizeof(int));
-        // delete old stuff
-        delete[] m_values;
-        delete[] m_indices;
-        m_values = newValues;
-        m_indices = newIndices;
-        m_allocatedSize = newAllocatedSize;
-      }
+        reallocate(newAllocatedSize);
+    }
+
+    void squeeze()
+    {
+      if (m_allocatedSize>m_size)
+        reallocate(m_size);
     }
 
     void resize(int size, int reserveSizeFactor = 0)
     {
       if (m_allocatedSize<size)
-      {
-        int newAllocatedSize = size + reserveSizeFactor*size;
-        Scalar* newValues  = new Scalar[newAllocatedSize];
-        int* newIndices = new int[newAllocatedSize];
-        // copy
-        memcpy(newValues,  m_values,  m_size * sizeof(Scalar));
-        memcpy(newIndices, m_indices, m_size * sizeof(int));
-        // delete old stuff
-        delete[] m_values;
-        delete[] m_indices;
-        m_values = newValues;
-        m_indices = newIndices;
-        m_allocatedSize = newAllocatedSize;
-      }
+        reallocate(size + reserveSizeFactor*size);
       m_size = size;
     }
 
@@ -123,26 +96,37 @@ template<typename Scalar> class SparseArray
     }
 
     int size() const { return m_size; }
-
-    void clear()
-    {
-      m_size = 0;
-    }
+    void clear() { m_size = 0; }
 
     Scalar& value(int i) { return m_values[i]; }
-    Scalar value(int i) const { return m_values[i]; }
+    const Scalar& value(int i) const { return m_values[i]; }
 
     int& index(int i) { return m_indices[i]; }
-    int index(int i) const { return m_indices[i]; }
+    const int& index(int i) const { return m_indices[i]; }
+
+  protected:
+
+    void reallocate(int size)
+    {
+      Scalar* newValues  = new Scalar[size];
+      int* newIndices = new int[size];
+      int copySize = std::min(size, m_size);
+      // copy
+      memcpy(newValues,  m_values,  copySize * sizeof(Scalar));
+      memcpy(newIndices, m_indices, copySize * sizeof(int));
+      // delete old stuff
+      delete[] m_values;
+      delete[] m_indices;
+      m_values = newValues;
+      m_indices = newIndices;
+      m_allocatedSize = size;
+    }
 
   protected:
     Scalar* m_values;
     int* m_indices;
     int m_size;
-    struct {
-      int m_allocatedSize:31;
-      mutable int m_isShared:1;
-    };
+    int m_allocatedSize;
 
 };
 
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index cbc504592..cc13e21c7 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -58,41 +58,50 @@ class SparseMatrix : public SparseMatrixBase<SparseMatrix<_Scalar, _Flags> >
     EIGEN_GENERIC_PUBLIC_INTERFACE(SparseMatrix)
 
   protected:
+  public:
+
+    typedef SparseMatrixBase<SparseMatrix> SparseBase;
+    enum {
+      RowMajor = SparseBase::RowMajor
+    };
 
-    int* m_colPtrs;
+    int m_outerSize;
+    int m_innerSize;
+    int* m_outerIndex;
     SparseArray<Scalar> m_data;
-    int m_rows;
-    int m_cols;
+
 
   public:
 
-    inline int rows() const { return m_rows; }
-    inline int cols() const { return m_cols; }
-    inline int innerNonZeros(int j) const { return m_colPtrs[j+1]-m_colPtrs[j]; }
+    inline int rows() const { return RowMajor ? m_outerSize : m_innerSize; }
+    inline int cols() const { return RowMajor ? m_innerSize : m_outerSize; }
+    inline int innerSize() const { return m_innerSize; }
+    inline int outerSize() const { return m_outerSize; }
+    inline int innerNonZeros(int j) const { return m_outerIndex[j+1]-m_outerIndex[j]; }
 
     inline Scalar coeff(int row, int col) const
     {
-      int id = m_colPtrs[col];
-      int end = m_colPtrs[col+1];
-      while (id<end && m_data.index(id)!=row)
-      {
-        ++id;
-      }
-      if (id==end)
-        return 0;
-      return m_data.value(id);
+      const int outer = RowMajor ? row : col;
+      const int inner = RowMajor ? col : row;
+
+      int id = m_outerIndex[outer];
+      int end = m_outerIndex[outer+1]-1;
+      if (m_data.index(end)==inner)
+        return m_data.value(end);
+      const int* r = std::lower_bound(&m_data.index(id),&m_data.index(end),inner);
+      return (*r==inner) ? m_data.value(*r) : Scalar(0);
     }
 
     inline Scalar& coeffRef(int row, int col)
     {
-      int id = m_colPtrs[col];
-      int end = m_colPtrs[col+1];
-      while (id<end && m_data.index(id)!=row)
-      {
-        ++id;
-      }
-      ei_assert(id!=end);
-      return m_data.value(id);
+      const int outer = RowMajor ? row : col;
+      const int inner = RowMajor ? col : row;
+
+      int id = m_outerIndex[outer];
+      int end = m_outerIndex[outer+1];
+      int* r = std::lower_bound(&m_data.index(id),&m_data.index(end),inner);
+      ei_assert(*r==inner);
+      return m_data.value(*r);
     }
 
   public:
@@ -106,92 +115,94 @@ class SparseMatrix : public SparseMatrixBase<SparseMatrix<_Scalar, _Flags> >
     {
       m_data.clear();
       m_data.reserve(reserveSize);
-      for (int i=0; i<=m_cols; ++i)
-        m_colPtrs[i] = 0;
+      for (int i=0; i<=m_outerSize; ++i)
+        m_outerIndex[i] = 0;
     }
 
     inline Scalar& fill(int row, int col)
     {
-      if (m_colPtrs[col+1]==0)
+      const int outer = RowMajor ? row : col;
+      const int inner = RowMajor ? col : row;
+
+      if (m_outerIndex[outer+1]==0)
       {
         int i=col;
-        while (i>=0 && m_colPtrs[i]==0)
+        while (i>=0 && m_outerIndex[i]==0)
         {
-          m_colPtrs[i] = m_data.size();
+          m_outerIndex[i] = m_data.size();
           --i;
         }
-        m_colPtrs[col+1] = m_colPtrs[col];
+        m_outerIndex[outer+1] = m_outerIndex[outer];
       }
-      assert(m_colPtrs[col+1] == m_data.size());
-      int id = m_colPtrs[col+1];
-      m_colPtrs[col+1]++;
+      assert(m_outerIndex[outer+1] == m_data.size());
+      int id = m_outerIndex[outer+1];
+      m_outerIndex[outer+1]++;
 
-      m_data.append(0, row);
+      m_data.append(0, inner);
       return m_data.value(id);
     }
 
     inline void endFill()
     {
       int size = m_data.size();
-      int i = m_cols;
+      int i = m_outerSize;
       // find the last filled column
-      while (i>=0 && m_colPtrs[i]==0)
+      while (i>=0 && m_outerIndex[i]==0)
         --i;
       i++;
-      while (i<=m_cols)
+      while (i<=m_outerSize)
       {
-        m_colPtrs[i] = size;
+        m_outerIndex[i] = size;
         ++i;
       }
     }
 
     void resize(int rows, int cols)
     {
-      if (m_cols != cols)
+      const int outerSize = RowMajor ? rows : cols;
+      m_innerSize = RowMajor ? cols : rows;
+      m_data.clear();
+      if (m_outerSize != outerSize)
       {
-        delete[] m_colPtrs;
-        m_colPtrs = new int [cols+1];
-        m_rows = rows;
-        m_cols = cols;
+        delete[] m_outerIndex;
+        m_outerIndex = new int [outerSize+1];
+        m_outerSize = outerSize;
       }
     }
 
     inline SparseMatrix(int rows, int cols)
-      : m_rows(0), m_cols(0), m_colPtrs(0)
+      : m_outerSize(0), m_innerSize(0), m_outerIndex(0)
     {
       resize(rows, cols);
     }
 
     template<typename OtherDerived>
     inline SparseMatrix(const MatrixBase<OtherDerived>& other)
-      : m_rows(0), m_cols(0), m_colPtrs(0)
+      : m_outerSize(0), m_innerSize(0), m_outerIndex(0)
     {
       *this = other.derived();
     }
 
-    inline void shallowCopy(const SparseMatrix& other)
+    inline void swap(SparseMatrix& other)
     {
-      EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: shallowCopy\n");
-      delete[] m_colPtrs;
-      m_colPtrs = 0;
-      m_rows = other.rows();
-      m_cols = other.cols();
-      m_colPtrs = other.m_colPtrs;
-      m_data.shallowCopy(other.m_data);
-      other.markAsCopied();
+      EIGEN_DBG_SPARSE(std::cout << "SparseMatrix:: swap\n");
+      std::swap(m_outerIndex, other.m_outerIndex);
+      std::swap(m_innerSize, other.m_innerSize);
+      std::swap(m_outerSize, other.m_outerSize);
+      m_data.swap(other.m_data);
     }
 
     inline SparseMatrix& operator=(const SparseMatrix& other)
     {
       if (other.isRValue())
       {
-        shallowCopy(other);
+        swap(other.const_cast_derived());
       }
       else
       {
         resize(other.rows(), other.cols());
-        for (int col=0; col<=cols(); ++col)
-          m_colPtrs[col] = other.m_colPtrs[col];
+        for (int j=0; j<=m_outerSize; ++j)
+          m_outerIndex[j] = other.m_outerIndex[j];
         m_data = other.m_data;
         return *this;
       }
@@ -216,7 +227,7 @@ class SparseMatrix : public SparseMatrixBase<SparseMatrix<_Scalar, _Flags> >
         s << "Column pointers:\n";
         for (uint i=0; i<m.cols(); ++i)
         {
-          s << m.m_colPtrs[i] << " ";
+          s << m.m_outerIndex[i] << " ";
         }
         s << std::endl;
         s << std::endl;
@@ -228,8 +239,7 @@ class SparseMatrix : public SparseMatrixBase<SparseMatrix<_Scalar, _Flags> >
     /** Destructor */
     inline ~SparseMatrix()
     {
-      if (this->isNotShared())
-        delete[] m_colPtrs;
+      delete[] m_outerIndex;
     }
 };
 
@@ -237,8 +247,8 @@ template<typename Scalar, int _Flags>
 class SparseMatrix<Scalar,_Flags>::InnerIterator
 {
   public:
-    InnerIterator(const SparseMatrix& mat, int col)
-      : m_matrix(mat), m_id(mat.m_colPtrs[col]), m_start(m_id), m_end(mat.m_colPtrs[col+1])
+    InnerIterator(const SparseMatrix& mat, int outer)
+      : m_matrix(mat), m_id(mat.m_outerIndex[outer]), m_start(m_id), m_end(mat.m_outerIndex[outer+1])
     {}
 
     InnerIterator& operator++() { m_id++; return *this; }
diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h
index c663a6902..cc43e7e07 100644
--- a/Eigen/src/Sparse/SparseMatrixBase.h
+++ b/Eigen/src/Sparse/SparseMatrixBase.h
@@ -43,19 +43,16 @@ class SparseMatrixBase : public MatrixBase<Derived>
     { return *static_cast<Derived*>(const_cast<SparseMatrixBase*>(this)); }
 
     SparseMatrixBase()
-      : m_isRValue(false), m_hasBeenCopied(false)
+      : m_isRValue(false)
     {}
 
     bool isRValue() const { return m_isRValue; }
-    Derived& temporary() { m_isRValue = true; return derived(); }
+    Derived& markAsRValue() { m_isRValue = true; return derived(); }
 
     inline Derived& operator=(const Derived& other)
     {
       if (other.isRValue())
-      {
-        m_hasBeenCopied = true;
-        derived().shallowCopy(other);
-      }
+        derived().swap(other.const_cast_derived());
       else
         this->operator=<Derived>(other);
       return derived();
@@ -82,7 +79,7 @@ class SparseMatrixBase : public MatrixBase<Derived>
       }
       temp.endFill();
 
-      derived() = temp.temporary();
+      derived() = temp.markAsRValue();
       return derived();
     }
 
@@ -119,7 +116,6 @@ class SparseMatrixBase : public MatrixBase<Derived>
 
     friend std::ostream & operator << (std::ostream & s, const SparseMatrixBase& m)
     {
-
       if (Flags&RowMajorBit)
       {
         for (int row=0; row<m.outerSize(); ++row)
@@ -130,6 +126,7 @@ class SparseMatrixBase : public MatrixBase<Derived>
             for ( ; col<it.index(); ++col)
               s << "0 ";
             std::cout << it.value() << " ";
+            ++col;
           }
           for ( ; col<m.cols(); ++col)
             s << "0 ";
@@ -144,15 +141,12 @@ class SparseMatrixBase : public MatrixBase<Derived>
       return s;
     }
 
-  protected:
-
-    bool isNotShared() { return !m_hasBeenCopied; }
-    void markAsCopied() const { m_hasBeenCopied = true; }
+    template<typename OtherDerived>
+    OtherDerived inverseProduct(const MatrixBase<OtherDerived>& other) const;
 
   protected:
 
     bool m_isRValue;
-    mutable bool m_hasBeenCopied;
 };
 
 #endif // EIGEN_SPARSEMATRIXBASE_H
diff --git a/Eigen/src/Sparse/SparseProduct.h b/Eigen/src/Sparse/SparseProduct.h
index e921e8775..6dff64f30 100644
--- a/Eigen/src/Sparse/SparseProduct.h
+++ b/Eigen/src/Sparse/SparseProduct.h
@@ -123,9 +123,6 @@ template<typename LhsNested, typename RhsNested> class Product<LhsNested,RhsNest
     const RhsNested m_rhs;
 };
 
-const int RowMajor = RowMajorBit;
-const int ColMajor = 0;
-
 template<typename Lhs, typename Rhs, typename ResultType,
   int LhsStorageOrder = ei_traits<Lhs>::Flags&RowMajorBit,
   int RhsStorageOrder = ei_traits<Rhs>::Flags&RowMajorBit,
diff --git a/Eigen/src/Sparse/TriangularSolver.h b/Eigen/src/Sparse/TriangularSolver.h
new file mode 100644
index 000000000..8634e114c
--- /dev/null
+++ b/Eigen/src/Sparse/TriangularSolver.h
@@ -0,0 +1,170 @@
+// 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_SPARSETRIANGULARSOLVER_H
+#define EIGEN_SPARSETRIANGULARSOLVER_H
+
+template<typename Lhs, typename Rhs,
+  int TriangularPart = (int(Lhs::Flags) & LowerTriangularBit)
+                     ? Lower
+                     : (int(Lhs::Flags) & UpperTriangularBit)
+                     ? Upper
+                     : -1,
+  int StorageOrder = int(Lhs::Flags) & RowMajorBit ? RowMajor : ColMajor
+  >
+struct ei_inverse_product_selector;
+
+// forward substitution, row-major
+template<typename Lhs, typename Rhs>
+struct ei_inverse_product_selector<Lhs,Rhs,Lower,RowMajor>
+{
+  typedef typename Rhs::Scalar Scalar;
+  static void run(const Lhs& lhs, const Rhs& rhs, Rhs& res)
+  {
+    for(int col=0 ; col<rhs.cols() ; ++col)
+    {
+      for(int i=0; i<lhs.rows(); ++i)
+      {
+        Scalar tmp = rhs.coeff(i,col);
+        Scalar lastVal = 0;
+        int lastIndex = 0;
+        for(typename Lhs::InnerIterator it(lhs, i); it; ++it)
+        {
+          lastVal = it.value();
+          lastIndex = it.index();
+          tmp -= lastVal * res.coeff(lastIndex,col);
+        }
+        if (Lhs::Flags & UnitDiagBit)
+          res.coeffRef(i,col) = tmp;
+        else
+        {
+          ei_assert(lastIndex==i);
+          res.coeffRef(i,col) = tmp/lastVal;
+        }
+      }
+    }
+  }
+};
+
+// backward substitution, row-major
+template<typename Lhs, typename Rhs>
+struct ei_inverse_product_selector<Lhs,Rhs,Upper,RowMajor>
+{
+  typedef typename Rhs::Scalar Scalar;
+  static void run(const Lhs& lhs, const Rhs& rhs, Rhs& res)
+  {
+    for(int col=0 ; col<rhs.cols() ; ++col)
+    {
+      for(int i=lhs.rows()-1 ; i>=0 ; --i)
+      {
+        Scalar tmp = rhs.coeff(i,col);
+        typename Lhs::InnerIterator it(lhs, i);
+        for(++it; it; ++it)
+        {
+          tmp -= it.value() * res.coeff(it.index(),col);
+        }
+
+        if (Lhs::Flags & UnitDiagBit)
+          res.coeffRef(i,col) = tmp;
+        else
+        {
+          typename Lhs::InnerIterator it(lhs, i);
+          ei_assert(it.index() == i);
+          res.coeffRef(i,col) = tmp/it.value();
+        }
+      }
+    }
+  }
+};
+
+// forward substitution, col-major
+template<typename Lhs, typename Rhs>
+struct ei_inverse_product_selector<Lhs,Rhs,Lower,ColMajor>
+{
+  typedef typename Rhs::Scalar Scalar;
+  static void run(const Lhs& lhs, const Rhs& rhs, Rhs& res)
+  {
+    // NOTE we could avoid this copy using an in-place API
+    res = rhs;
+    for(int col=0 ; col<rhs.cols() ; ++col)
+    {
+      for(int i=0; i<lhs.cols(); ++i)
+      {
+        typename Lhs::InnerIterator it(lhs, i);
+        if(!(Lhs::Flags & UnitDiagBit))
+        {
+          ei_assert(it.index()==i);
+          res.coeffRef(i,col) /= it.value();
+        }
+        Scalar tmp = res.coeffRef(i,col);
+        for(++it; it; ++it)
+          res.coeffRef(it.index(), col) -= tmp * it.value();
+      }
+    }
+  }
+};
+
+// backward substitution, col-major
+template<typename Lhs, typename Rhs>
+struct ei_inverse_product_selector<Lhs,Rhs,Upper,ColMajor>
+{
+  typedef typename Rhs::Scalar Scalar;
+  static void run(const Lhs& lhs, const Rhs& rhs, Rhs& res)
+  {
+    // NOTE we could avoid this copy using an in-place API
+    res = rhs;
+    for(int col=0 ; col<rhs.cols() ; ++col)
+    {
+      for(int i=lhs.cols()-1; i>=0; --i)
+      {
+        if(!(Lhs::Flags & UnitDiagBit))
+        {
+          // FIXME lhs.coeff(i,i) might not be always efficient while it must simply be the
+          // last element of the column !
+          res.coeffRef(i,col) /= lhs.coeff(i,i);
+        }
+        Scalar tmp = res.coeffRef(i,col);
+        typename Lhs::InnerIterator it(lhs, i);
+        for(; it && it.index()<i; ++it)
+          res.coeffRef(it.index(), col) -= tmp * it.value();
+      }
+    }
+  }
+};
+
+template<typename Derived>
+template<typename OtherDerived>
+OtherDerived
+SparseMatrixBase<Derived>::inverseProduct(const MatrixBase<OtherDerived>& other) const
+{
+  ei_assert(derived().cols() == other.rows());
+  ei_assert(!(Flags & ZeroDiagBit));
+  ei_assert(Flags & (UpperTriangularBit|LowerTriangularBit));
+
+  OtherDerived res(other.rows(), other.cols());
+  ei_inverse_product_selector<Derived, OtherDerived>::run(derived(), other.derived(), res);
+  return res;
+}
+
+#endif // EIGEN_SPARSETRIANGULARSOLVER_H
diff --git a/bench/BenchSparseUtil.h b/bench/BenchSparseUtil.h
index 97838f4b1..9d88148d0 100644
--- a/bench/BenchSparseUtil.h
+++ b/bench/BenchSparseUtil.h
@@ -64,7 +64,8 @@ void eiToGmm(const EigenSparseMatrix& src, GmmSparse& dst)
 
 #ifndef NOMTL
 #include <boost/numeric/mtl/mtl.hpp>
-typedef mtl::compressed2D<Scalar> MtlSparse;
+typedef mtl::compressed2D<Scalar, mtl::matrix::parameters<mtl::tag::col_major> > MtlSparse;
+typedef mtl::compressed2D<Scalar, mtl::matrix::parameters<mtl::tag::row_major> > MtlSparseRowMajor;
 void eiToMtl(const EigenSparseMatrix& src, MtlSparse& dst)
 {
   mtl::matrix::inserter<MtlSparse> ins(dst);
diff --git a/bench/sparse_trisolver.cpp b/bench/sparse_trisolver.cpp
new file mode 100644
index 000000000..27362c11b
--- /dev/null
+++ b/bench/sparse_trisolver.cpp
@@ -0,0 +1,206 @@
+
+//g++ -O3 -g0 -DNDEBUG  sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.005 -DSIZE=10000 && ./a.out
+//g++ -O3 -g0 -DNDEBUG  sparse_product.cpp -I.. -I/home/gael/Coding/LinearAlgebra/mtl4/ -DDENSITY=0.05 -DSIZE=2000 && ./a.out
+// -DNOGMM -DNOMTL
+
+#ifndef SIZE
+#define SIZE 10000
+#endif
+
+#ifndef DENSITY
+#define DENSITY 0.01
+#endif
+
+#ifndef REPEAT
+#define REPEAT 1
+#endif
+
+#include "BenchSparseUtil.h"
+
+#ifndef MINDENSITY
+#define MINDENSITY 0.0004
+#endif
+
+#ifndef NBTRIES
+#define NBTRIES 10
+#endif
+
+#define BENCH(X) \
+  timer.reset(); \
+  for (int _j=0; _j<NBTRIES; ++_j) { \
+    timer.start(); \
+    for (int _k=0; _k<REPEAT; ++_k) { \
+        X  \
+  } timer.stop(); }
+
+typedef SparseMatrix<Scalar,Upper> EigenSparseTriMatrix;
+typedef SparseMatrix<Scalar,RowMajorBit|Upper> EigenSparseTriMatrixRow;
+
+void fillMatrix(float density, int rows, int cols,  EigenSparseTriMatrix& dst)
+{
+  dst.startFill(rows*cols*density);
+  for(int j = 0; j < cols; j++)
+  {
+    for(int i = 0; i < j; i++)
+    {
+      Scalar v = (ei_random<float>(0,1) < density) ? ei_random<Scalar>() : 0;
+      if (v!=0)
+        dst.fill(i,j) = v;
+    }
+    dst.fill(j,j) = ei_random<Scalar>();
+  }
+  dst.endFill();
+}
+
+int main(int argc, char *argv[])
+{
+  int rows = SIZE;
+  int cols = SIZE;
+  float density = DENSITY;
+  BenchTimer timer;
+  #if 1
+  EigenSparseTriMatrix sm1(rows,cols);
+  VectorXf b = VectorXf::random(cols);
+  VectorXf x = VectorXf::random(cols);
+
+  bool densedone = false;
+
+  for (float density = DENSITY; density>=MINDENSITY; density*=0.5)
+  {
+    EigenSparseTriMatrix sm1(rows, cols);
+    fillMatrix(density, rows, cols, sm1);
+
+    // dense matrices
+    #ifdef DENSEMATRIX
+    if (!densedone)
+    {
+      densedone = true;
+      std::cout << "Eigen Dense\t" << density*100 << "%\n";
+      DenseMatrix m1(rows,cols);
+      Matrix<Scalar,Dynamic,Dynamic,Dynamic,Dynamic,RowMajorBit> m2(rows,cols);
+      eiToDense(sm1, m1);
+      m2 = m1;
+
+      BENCH(x = m1.marked<Upper>().inverseProduct(b);)
+      std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
+      std::cerr << x.transpose() << "\n";
+
+      BENCH(x = m2.marked<Upper>().inverseProduct(b);)
+      std::cout << "   rowmajor^-1 * b:\t" << timer.value() << endl;
+      std::cerr << x.transpose() << "\n";
+    }
+    #endif
+
+    // eigen sparse matrices
+    {
+      std::cout << "Eigen sparse\t" << density*100 << "%\n";
+      EigenSparseTriMatrixRow sm2 = sm1;
+
+      BENCH(x = sm1.inverseProduct(b);)
+      std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
+      std::cerr << x.transpose() << "\n";
+
+      BENCH(x = sm2.inverseProduct(b);)
+      std::cout << "   rowmajor^-1 * b:\t" << timer.value() << endl;
+      std::cerr << x.transpose() << "\n";
+
+//       x = b;
+//       BENCH(sm1.inverseProductInPlace(x);)
+//       std::cout << "   colmajor^-1 * b:\t" << timer.value() << " (inplace)" << endl;
+//       std::cerr << x.transpose() << "\n";
+//
+//       x = b;
+//       BENCH(sm2.inverseProductInPlace(x);)
+//       std::cout << "   rowmajor^-1 * b:\t" << timer.value() << " (inplace)" << endl;
+//       std::cerr << x.transpose() << "\n";
+    }
+
+    // GMM++
+    #ifndef NOGMM
+    {
+      std::cout << "GMM++ sparse\t" << density*100 << "%\n";
+      GmmSparse m1(rows,cols);
+      gmm::csr_matrix<Scalar> m2;
+      eiToGmm(sm1, m1);
+      gmm::copy(m1,m2);
+      std::vector<Scalar> gmmX(cols), gmmB(cols);
+      Map<Matrix<Scalar,Dynamic,1> >(&gmmX[0], cols) = x;
+      Map<Matrix<Scalar,Dynamic,1> >(&gmmB[0], cols) = b;
+
+      gmmX = gmmB;
+      BENCH(gmm::upper_tri_solve(m1, gmmX, false);)
+      std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
+      std::cerr << Map<Matrix<Scalar,Dynamic,1> >(&gmmX[0], cols).transpose() << "\n";
+
+      gmmX = gmmB;
+      BENCH(gmm::upper_tri_solve(m2, gmmX, false);)
+      timer.stop();
+      std::cout << "   rowmajor^-1 * b:\t" << timer.value() << endl;
+      std::cerr << Map<Matrix<Scalar,Dynamic,1> >(&gmmX[0], cols).transpose() << "\n";
+    }
+    #endif
+
+    // MTL4
+    #ifndef NOMTL
+    {
+      std::cout << "MTL4\t" << density*100 << "%\n";
+      MtlSparse m1(rows,cols);
+      MtlSparseRowMajor m2(rows,cols);
+      eiToMtl(sm1, m1);
+      m2 = m1;
+      mtl::dense_vector<Scalar> x(rows, 1.0);
+      mtl::dense_vector<Scalar> b(rows, 1.0);
+
+      BENCH(x = mtl::upper_trisolve(m1,b);)
+      std::cout << "   colmajor^-1 * b:\t" << timer.value() << endl;
+//       std::cerr << x << "\n";
+
+      BENCH(x = mtl::upper_trisolve(m2,b);)
+      std::cout << "   rowmajor^-1 * b:\t" << timer.value() << endl;
+//       std::cerr << x << "\n";
+    }
+    #endif
+
+
+
+  }
+  #endif
+
+  #if 0
+    // bench small matrices (in-place versus return bye value)
+    {
+      timer.reset();
+      for (int _j=0; _j<10; ++_j) {
+        Matrix4f m = Matrix4f::random();
+        Vector4f b = Vector4f::random();
+        Vector4f x = Vector4f::random();
+        timer.start();
+        for (int _k=0; _k<1000000; ++_k) {
+          b = m.inverseProduct(b);
+        }
+        timer.stop();
+      }
+      std::cout << "4x4 :\t" << timer.value() << endl;
+    }
+
+    {
+      timer.reset();
+      for (int _j=0; _j<10; ++_j) {
+        Matrix4f m = Matrix4f::random();
+        Vector4f b = Vector4f::random();
+        Vector4f x = Vector4f::random();
+        timer.start();
+        for (int _k=0; _k<1000000; ++_k) {
+          m.inverseProductInPlace(x);
+        }
+        timer.stop();
+      }
+      std::cout << "4x4 IP :\t" << timer.value() << endl;
+    }
+  #endif
+
+    std::cout << "\n\n";
+
+  return 0;
+}
+