sparse module:

 - remove some useless stuff => let's focus on a single sparse matrix format
 - finalize the new RandomSetter

3 files changed, 621 insertions, 0 deletions

diff --git a/disabled/HashMatrix.h b/disabled/HashMatrix.h
new file mode 100644
index 000000000..1d6d33cce
--- /dev/null
+++ b/disabled/HashMatrix.h
@@ -0,0 +1,166 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_HASHMATRIX_H
+#define EIGEN_HASHMATRIX_H
+
+template<typename _Scalar, int _Flags>
+struct ei_traits<HashMatrix<_Scalar, _Flags> >
+{
+  typedef _Scalar Scalar;
+  enum {
+    RowsAtCompileTime = Dynamic,
+    ColsAtCompileTime = Dynamic,
+    MaxRowsAtCompileTime = Dynamic,
+    MaxColsAtCompileTime = Dynamic,
+    Flags = SparseBit | _Flags,
+    CoeffReadCost = NumTraits<Scalar>::ReadCost,
+    SupportedAccessPatterns = RandomAccessPattern
+  };
+};
+
+// TODO reimplement this class using custom linked lists
+template<typename _Scalar, int _Flags>
+class HashMatrix
+  : public SparseMatrixBase<HashMatrix<_Scalar, _Flags> >
+{
+  public:
+    EIGEN_GENERIC_PUBLIC_INTERFACE(HashMatrix)
+    class InnerIterator;
+  protected:
+
+    typedef typename std::map<int, Scalar>::iterator MapIterator;
+    typedef typename std::map<int, Scalar>::const_iterator ConstMapIterator;
+
+  public:
+    inline int rows() const { return m_innerSize; }
+    inline int cols() const { return m_data.size(); }
+
+    inline const Scalar& coeff(int row, int col) const
+    {
+      const MapIterator it = m_data[col].find(row);
+      if (it!=m_data[col].end())
+        return Scalar(0);
+      return it->second;
+    }
+
+    inline Scalar& coeffRef(int row, int col)
+    {
+      return m_data[col][row];
+    }
+
+  public:
+
+    inline void startFill(int /*reserveSize = 1000 --- currently unused, don't generate a warning*/) {}
+
+    inline Scalar& fill(int row, int col) { return coeffRef(row, col); }
+
+    inline void endFill() {}
+
+    ~HashMatrix()
+    {}
+
+    inline void shallowCopy(const HashMatrix& other)
+    {
+      EIGEN_DBG_SPARSE(std::cout << "HashMatrix:: shallowCopy\n");
+      // FIXME implement a true shallow copy !!
+      resize(other.rows(), other.cols());
+      for (int j=0; j<this->outerSize(); ++j)
+        m_data[j] = other.m_data[j];
+    }
+
+    void resize(int _rows, int _cols)
+    {
+      if (cols() != _cols)
+      {
+        m_data.resize(_cols);
+      }
+      m_innerSize = _rows;
+    }
+
+    inline HashMatrix(int rows, int cols)
+      : m_innerSize(0)
+    {
+      resize(rows, cols);
+    }
+
+    template<typename OtherDerived>
+    inline HashMatrix(const MatrixBase<OtherDerived>& other)
+      : m_innerSize(0)
+    {
+      *this = other.derived();
+    }
+
+    inline HashMatrix& operator=(const HashMatrix& other)
+    {
+      if (other.isRValue())
+      {
+        shallowCopy(other);
+      }
+      else
+      {
+        resize(other.rows(), other.cols());
+        for (int col=0; col<cols(); ++col)
+          m_data[col] = other.m_data[col];
+      }
+      return *this;
+    }
+
+    template<typename OtherDerived>
+    inline HashMatrix& operator=(const MatrixBase<OtherDerived>& other)
+    {
+      return SparseMatrixBase<HashMatrix>::operator=(other);
+    }
+
+  protected:
+
+    std::vector<std::map<int, Scalar> > m_data;
+    int m_innerSize;
+
+};
+
+template<typename Scalar, int _Flags>
+class HashMatrix<Scalar,_Flags>::InnerIterator
+{
+  public:
+
+    InnerIterator(const HashMatrix& mat, int col)
+      : m_matrix(mat), m_it(mat.m_data[col].begin()), m_end(mat.m_data[col].end())
+    {}
+
+    InnerIterator& operator++() { m_it++; return *this; }
+
+    Scalar value() { return m_it->second; }
+
+    int index() const { return m_it->first; }
+
+    operator bool() const { return m_it!=m_end; }
+
+  protected:
+    const HashMatrix& m_matrix;
+    typename HashMatrix::ConstMapIterator m_it;
+    typename HashMatrix::ConstMapIterator m_end;
+};
+
+#endif // EIGEN_HASHMATRIX_H
diff --git a/disabled/LinkedVectorMatrix.h b/disabled/LinkedVectorMatrix.h
new file mode 100644
index 000000000..4e5a4862e
--- /dev/null
+++ b/disabled/LinkedVectorMatrix.h
@@ -0,0 +1,317 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_LINKEDVECTORMATRIX_H
+#define EIGEN_LINKEDVECTORMATRIX_H
+
+template<typename _Scalar, int _Flags>
+struct ei_traits<LinkedVectorMatrix<_Scalar,_Flags> >
+{
+  typedef _Scalar Scalar;
+  enum {
+    RowsAtCompileTime = Dynamic,
+    ColsAtCompileTime = Dynamic,
+    MaxRowsAtCompileTime = Dynamic,
+    MaxColsAtCompileTime = Dynamic,
+    Flags = SparseBit | _Flags,
+    CoeffReadCost = NumTraits<Scalar>::ReadCost,
+    SupportedAccessPatterns = InnerCoherentAccessPattern
+  };
+};
+
+template<typename Element, int ChunkSize = 8>
+struct LinkedVectorChunk
+{
+  LinkedVectorChunk() : next(0), prev(0), size(0) {}
+  Element data[ChunkSize];
+  LinkedVectorChunk* next;
+  LinkedVectorChunk* prev;
+  int size;
+  bool isFull() const { return size==ChunkSize; }
+};
+
+template<typename _Scalar, int _Flags>
+class LinkedVectorMatrix
+  : public SparseMatrixBase<LinkedVectorMatrix<_Scalar,_Flags> >
+{
+  public:
+    EIGEN_GENERIC_PUBLIC_INTERFACE(LinkedVectorMatrix)
+    class InnerIterator;
+  protected:
+
+    enum {
+      RowMajor = Flags&RowMajorBit ? 1 : 0
+    };
+
+    struct ValueIndex
+    {
+      ValueIndex() : value(0), index(0) {}
+      ValueIndex(Scalar v, int i) : value(v), index(i) {}
+      Scalar value;
+      int index;
+    };
+    typedef LinkedVectorChunk<ValueIndex,8> VectorChunk;
+
+    inline int find(VectorChunk** _el, int id)
+    {
+      VectorChunk* el = *_el;
+      while (el && el->data[el->size-1].index<id)
+        el = el->next;
+      *_el = el;
+      if (el)
+      {
+        // binary search
+        int maxI = el->size-1;
+        int minI = 0;
+        int i = el->size/2;
+        const ValueIndex* data = el->data;
+        while (data[i].index!=id)
+        {
+          if (data[i].index<id)
+          {
+            minI = i+1;
+            i = (maxI + minI)+2;
+          }
+          else
+          {
+            maxI = i-1;
+            i = (maxI + minI)+2;
+          }
+          if (minI>=maxI)
+            return -1;
+        }
+        if (data[i].index==id)
+          return i;
+      }
+      return -1;
+    }
+
+  public:
+    inline int rows() const { return RowMajor ? m_data.size() : m_innerSize; }
+    inline int cols() const { return RowMajor ? m_innerSize : m_data.size(); }
+
+    inline const Scalar& coeff(int row, int col) const
+    {
+      const int outer = RowMajor ? row : col;
+      const int inner = RowMajor ? col : row;
+
+      VectorChunk* el = m_data[outer];
+      int id = find(&el, inner);
+      if (id<0)
+        return Scalar(0);
+      return el->data[id].value;
+    }
+
+    inline Scalar& coeffRef(int row, int col)
+    {
+      const int outer = RowMajor ? row : col;
+      const int inner = RowMajor ? col : row;
+
+      VectorChunk* el = m_data[outer];
+      int id = find(&el, inner);
+      ei_assert(id>=0);
+//       if (id<0)
+//         return Scalar(0);
+      return el->data[id].value;
+    }
+
+  public:
+
+    inline void startFill(int reserveSize = 1000)
+    {
+      clear();
+      for (unsigned int i=0; i<m_data.size(); ++i)
+        m_ends[i] = m_data[i] = 0;
+    }
+
+    inline Scalar& fill(int row, int col)
+    {
+      const int outer = RowMajor ? row : col;
+      const int inner = RowMajor ? col : row;
+//       std::cout << " ll fill " << outer << "," << inner << "\n";
+      if (m_ends[outer]==0)
+      {
+        m_data[outer] = m_ends[outer] = new VectorChunk();
+      }
+      else
+      {
+        ei_assert(m_ends[outer]->data[m_ends[outer]->size-1].index < inner);
+        if (m_ends[outer]->isFull())
+        {
+
+          VectorChunk* el = new VectorChunk();
+          m_ends[outer]->next = el;
+          el->prev = m_ends[outer];
+          m_ends[outer] = el;
+        }
+      }
+      m_ends[outer]->data[m_ends[outer]->size].index = inner;
+      return m_ends[outer]->data[m_ends[outer]->size++].value;
+    }
+
+    inline void endFill() { }
+
+    void printDbg()
+    {
+      for (int j=0; j<m_data.size(); ++j)
+      {
+        VectorChunk* el = m_data[j];
+        while (el)
+        {
+          for (int i=0; i<el->size; ++i)
+            std::cout << j << "," << el->data[i].index << " = " << el->data[i].value << "\n";
+          el = el->next;
+        }
+      }
+      for (int j=0; j<m_data.size(); ++j)
+      {
+        InnerIterator it(*this,j);
+        while (it)
+        {
+          std::cout << j << "," << it.index() << " = " << it.value() << "\n";
+          ++it;
+        }
+      }
+    }
+
+    ~LinkedVectorMatrix()
+    {
+      clear();
+    }
+
+    void clear()
+    {
+      for (unsigned int i=0; i<m_data.size(); ++i)
+      {
+        VectorChunk* el = m_data[i];
+        while (el)
+        {
+          VectorChunk* tmp = el;
+          el = el->next;
+          delete tmp;
+        }
+      }
+    }
+
+    void resize(int rows, int cols)
+    {
+      const int outers = RowMajor ? rows : cols;
+      const int inners = RowMajor ? cols : rows;
+
+      if (this->outerSize() != outers)
+      {
+        clear();
+        m_data.resize(outers);
+        m_ends.resize(outers);
+        for (unsigned int i=0; i<m_data.size(); ++i)
+          m_ends[i] = m_data[i] = 0;
+      }
+      m_innerSize = inners;
+    }
+
+    inline LinkedVectorMatrix(int rows, int cols)
+      : m_innerSize(0)
+    {
+      resize(rows, cols);
+    }
+
+    template<typename OtherDerived>
+    inline LinkedVectorMatrix(const MatrixBase<OtherDerived>& other)
+      : m_innerSize(0)
+    {
+      *this = other.derived();
+    }
+
+    inline void swap(LinkedVectorMatrix& other)
+    {
+      EIGEN_DBG_SPARSE(std::cout << "LinkedVectorMatrix:: swap\n");
+      resize(other.rows(), other.cols());
+      m_data.swap(other.m_data);
+      m_ends.swap(other.m_ends);
+    }
+
+    inline LinkedVectorMatrix& operator=(const LinkedVectorMatrix& other)
+    {
+      if (other.isRValue())
+      {
+        swap(other.const_cast_derived());
+      }
+      else
+      {
+        // TODO implement a specialized deep copy here
+        return operator=<LinkedVectorMatrix>(other);
+      }
+      return *this;
+    }
+
+    template<typename OtherDerived>
+    inline LinkedVectorMatrix& operator=(const MatrixBase<OtherDerived>& other)
+    {
+      return SparseMatrixBase<LinkedVectorMatrix>::operator=(other.derived());
+    }
+
+  protected:
+
+    // 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;
+
+};
+
+
+template<typename Scalar, int _Flags>
+class LinkedVectorMatrix<Scalar,_Flags>::InnerIterator
+{
+  public:
+
+    InnerIterator(const LinkedVectorMatrix& mat, int col)
+      : m_matrix(mat), m_el(mat.m_data[col]), m_it(0)
+    {}
+
+    InnerIterator& operator++()
+    {
+      m_it++;
+      if (m_it>=m_el->size)
+      {
+        m_el = m_el->next;
+        m_it = 0;
+      }
+      return *this;
+    }
+
+    Scalar value() { return m_el->data[m_it].value; }
+
+    int index() const { return m_el->data[m_it].index; }
+
+    operator bool() const { return m_el && (m_el->next || m_it<m_el->size); }
+
+  protected:
+    const LinkedVectorMatrix& m_matrix;
+    VectorChunk* m_el;
+    int m_it;
+};
+
+#endif // EIGEN_LINKEDVECTORMATRIX_H
diff --git a/disabled/SparseSetter.h b/disabled/SparseSetter.h
new file mode 100644
index 000000000..20569d3ba
--- /dev/null
+++ b/disabled/SparseSetter.h
@@ -0,0 +1,138 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_SPARSESETTER_H
+#define EIGEN_SPARSESETTER_H
+
+template<typename MatrixType, int AccessPattern,
+  int IsSupported = ei_support_access_pattern<MatrixType,AccessPattern>::ret>
+struct ei_sparse_setter_selector;
+
+/** \class SparseSetter
+  *
+  * Goal: provides a unified API to fill/update a dense or sparse matrix.
+  *
+  * Usage:
+  * \code
+  * {
+  *   SparseSetter<MatrixType, RandomAccessPattern> w(m);
+  *   for (...) w->coeffRef(rand(),rand()) = rand();
+  * }
+  * \endcode
+  *
+  * In the above example we want to fill a matrix m (could be a SparseMatrix or whatever other matrix type)
+  * in a random fashion (whence the RandomAccessPattern). Internally, if \a MatrixType supports random writes
+  * then \c w behaves as a pointer to m, and m is filled directly. Otherwise, a temporary matrix supporting
+  * random writes is created and \c w behaves as a pointer to this temporary object. When the object \c w
+  * is deleted (at the end of the block), then the temporary object is assigned to the matrix m.
+  * 
+  * So far we can distinghished 4 types of access pattern:
+  * - FullyCoherentAccessPattern (if col major, i+j*rows must increase)
+  * - InnerCoherentAccessPattern (if col major, i must increase for each column j)
+  * - OuterCoherentAccessPattern (if col major, the column j is set in a random order, but j must increase)
+  * - RandomAccessPattern
+  * 
+  * See the wiki for more details.
+  * 
+  * The template class ei_support_access_pattern is used to determine the type of the temporary object (which
+  * can be a reference to \a MatrixType if \a MatrixType support \a AccessPattern)
+  *
+  * Currently only the RandomAccessPattern seems to work as expected.
+  *
+  * \todo define the API for each kind of access pattern
+  * \todo allows both update and set modes (set start a new matrix)
+  * \todo implement the OuterCoherentAccessPattern
+  *
+  */
+template<typename MatrixType,
+         int AccessPattern,
+         typename WrapperType = typename ei_sparse_setter_selector<MatrixType,AccessPattern>::type>
+class SparseSetter
+{
+    typedef typename ei_unref<WrapperType>::type _WrapperType;
+  public:
+
+    inline SparseSetter(MatrixType& matrix) : m_wrapper(matrix), mp_matrix(&matrix) {}
+
+    ~SparseSetter()
+    { *mp_matrix = m_wrapper; }
+
+    inline _WrapperType* operator->() { return &m_wrapper; }
+
+    inline _WrapperType& operator*() { return m_wrapper; }
+
+  protected:
+
+    WrapperType m_wrapper;
+    MatrixType* mp_matrix;
+};
+
+template<typename MatrixType, int AccessPattern>
+struct ei_sparse_setter_selector<MatrixType, AccessPattern, AccessPatternSupported>
+{
+  typedef MatrixType& type;
+};
+
+// forward each derived of SparseMatrixBase to the generic SparseMatrixBase specializations
+template<typename Scalar, int Flags, int AccessPattern>
+struct ei_sparse_setter_selector<SparseMatrix<Scalar,Flags>, AccessPattern, AccessPatternNotSupported>
+: public ei_sparse_setter_selector<SparseMatrixBase<SparseMatrix<Scalar,Flags> >,AccessPattern, AccessPatternNotSupported>
+{};
+
+template<typename Scalar, int Flags, int AccessPattern>
+struct ei_sparse_setter_selector<LinkedVectorMatrix<Scalar,Flags>, AccessPattern, AccessPatternNotSupported>
+: public ei_sparse_setter_selector<LinkedVectorMatrix<SparseMatrix<Scalar,Flags> >,AccessPattern, AccessPatternNotSupported>
+{};
+
+template<typename Scalar, int Flags, int AccessPattern>
+struct ei_sparse_setter_selector<HashMatrix<Scalar,Flags>, AccessPattern, AccessPatternNotSupported>
+: public ei_sparse_setter_selector<HashMatrix<SparseMatrix<Scalar,Flags> >,AccessPattern, AccessPatternNotSupported>
+{};
+
+// generic SparseMatrixBase specializations
+template<typename Derived>
+struct ei_sparse_setter_selector<SparseMatrixBase<Derived>, RandomAccessPattern, AccessPatternNotSupported>
+{
+  typedef HashMatrix<typename Derived::Scalar, Derived::Flags> type;
+};
+
+template<typename Derived>
+struct ei_sparse_setter_selector<SparseMatrixBase<Derived>, OuterCoherentAccessPattern, AccessPatternNotSupported>
+{
+  typedef HashMatrix<typename Derived::Scalar, Derived::Flags> type;
+};
+
+template<typename Derived>
+struct ei_sparse_setter_selector<SparseMatrixBase<Derived>, InnerCoherentAccessPattern, AccessPatternNotSupported>
+{
+  typedef LinkedVectorMatrix<typename Derived::Scalar, Derived::Flags> type;
+};
+
+template<typename Derived>
+struct ei_sparse_setter_selector<SparseMatrixBase<Derived>, FullyCoherentAccessPattern, AccessPatternNotSupported>
+{
+  typedef SparseMatrix<typename Derived::Scalar, Derived::Flags> type;
+};
+
+#endif // EIGEN_SPARSESETTER_H