* 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

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;
 };
 

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;
 
 };
 

Eigen/src/Sparse/SparseMatrix.h

@@ -58,41 +58,50 @@ class SparseMatrix : public SparseMatrixBase<SparseMatrix<_Scalar, _Flags> >
     EIGEN_GENERIC_PUBLIC_INTERFACE(SparseMatrix)
 
   protected:
+  public:
 
-    int* m_colPtrs;
+    typedef SparseMatrixBase<SparseMatrix> SparseBase;
+    enum {
+      RowMajor = SparseBase::RowMajor
+    };
+
+    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; }

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

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,

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