Merge branch 'feature/improve-matrixindexset' into 'master'

Improve performance of MatrixIndexSet

See merge request core/dune-istl!527

CHANGELOG.md

@@ -5,6 +5,11 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
 
 # Master (will become release 2.10)
 
+- The internal storage in `MatrixIndexSet` was changed from `std::set` to a sorted `std::vector`
+  (with a `std::set` fallback for very dense rows) to improve performance. The stored index
+  type was changed from `std::size_t` to `uint32_t` to reduce memory consumption and improve
+  performance. Hence, `MatrixIndexSet` can no longer be used for very large matrices with more
+  than 2^32 columns.
 
 ## Deprecations and removals
 

dune/istl/matrixindexset.hh

@@ -5,55 +5,149 @@
 #ifndef DUNE_ISTL_MATRIXINDEXSET_HH
 #define DUNE_ISTL_MATRIXINDEXSET_HH
 
-#include <vector>
+#include <algorithm>
+#include <cstddef>
+#include <cstdint>
 #include <set>
+#include <variant>
+#include <vector>
+
+#include <dune/common/overloadset.hh>
 
 namespace Dune {
 
 
-  /** \brief Stores the nonzero entries in a sparse matrix */
+  /**
+   * \brief Stores the nonzero entries for creating a sparse matrix
+   *
+   * This stores std::set-like container for the column index
+   * of each row. A sorted std::vector is used for this container
+   * up to a customizable maxVectorSize. If this size is exceeded,
+   * storage of the corresponding row is switched to std::set.
+   *
+   * The default value for maxVectorSize works well and ensures
+   * that the slow std::set fallback is only used for very
+   * dense rows.
+   *
+   * This class is thread safe in the sense that concurrent calls
+   * to all const methods and furthermore to add(row,col) with different
+   * rows in each thread are safe.
+   */
   class MatrixIndexSet
   {
+    using Index = std::uint_least32_t;
+
+    // A vector that partly mimics a std::set by staying
+    // sorted on insert() and having unique values.
+    class FlatSet : public std::vector<Index>
+    {
+      using Base = std::vector<Index>;
+    public:
+      using Base::Base;
+      using Base::begin;
+      using Base::end;
+      void insert(const Index& value) {
+        auto it = std::lower_bound(begin(), end(), value);
+        if ((it == end() or (*it != value)))
+          Base::insert(it, value);
+      }
+      bool contains(const Index& value) const {
+        return std::binary_search(begin(), end(), value);
+      }
+    };
+
+    using RowIndexSet = std::variant<FlatSet, std::set<Index>>;
 
   public:
-    typedef std::size_t size_type;
 
-    /** \brief Default constructor */
-    MatrixIndexSet() : rows_(0), cols_(0)
+    using size_type = Index;
+
+    /**
+     * \brief Default value for maxVectorSize
+     *
+     * This was selected after benchmarking for the worst case
+     * of reverse insertion of column indices. In many applications
+     * this works well. There's no need to use a different value
+     * unless you have many dense rows with more than defaultMaxVectorSize
+     * nonzero entries. Even in this case defaultMaxVectorSize may work
+     * well and a finding a better value may require careful
+     * benchmarking.
+     */
+    static constexpr size_type defaultMaxVectorSize = 2048;
+
+    /**
+     * \brief Constructor with custom maxVectorSize
+     *
+     * \param maxVectorSize Maximal size for stored row vector (default is defaultMaxVectorSize).
+     */
+    MatrixIndexSet(size_type maxVectorSize=defaultMaxVectorSize) noexcept : rows_(0), cols_(0), maxVectorSize_(maxVectorSize)
     {}
 
-    /** \brief Constructor setting the matrix size */
-    MatrixIndexSet(size_type rows, size_type cols) : rows_(rows), cols_(cols) {
-      indices_.resize(rows_);
+    /**
+     * \brief Constructor setting the matrix size
+     *
+     * \param rows Number of matrix rows
+     * \param cols Number of matrix columns
+     * \param maxVectorSize Maximal size for stored row vector (default is defaultMaxVectorSize).
+     */
+    MatrixIndexSet(size_type rows, size_type cols, size_type maxVectorSize=defaultMaxVectorSize) : rows_(rows), cols_(cols), maxVectorSize_(maxVectorSize)
+    {
+      indices_.resize(rows_, FlatSet());
     }
 
     /** \brief Reset the size of an index set */
     void resize(size_type rows, size_type cols) {
       rows_ = rows;
       cols_ = cols;
-      indices_.resize(rows_);
+      indices_.resize(rows_, FlatSet());
     }
 
-    /** \brief Add an index to the index set */
-    void add(size_type i, size_type j) {
-      indices_[i].insert(j);
+    /**
+     * \brief Add an index to the index set
+     *
+     * It is safe to call add(row, col) for different rows in concurrent threads,
+     * but it is not safe to do concurrent calls for the same row, even for different
+     * columns.
+     */
+    void add(size_type row, size_type col) {
+      return std::visit(Dune::overload(
+        // If row is stored as set, call insert directly
+        [&](std::set<size_type>& set) {
+          set.insert(col);
+        },
+        // If row is stored as vector only insert directly
+        // if maxVectorSize_ is not reached. Otherwise switch
+        // to set storage first.
+        [&](FlatSet& sortedVector) {
+          if (sortedVector.size() < maxVectorSize_)
+            sortedVector.insert(col);
+          else if (not sortedVector.contains(col))
+          {
+            std::set<size_type> set(sortedVector.begin(), sortedVector.end());
+            set.insert(col);
+            indices_[row] = std::move(set);
+          }
+        }
+      ), indices_[row]);
     }
 
     /** \brief Return the number of entries */
     size_type size() const {
       size_type entries = 0;
       for (size_type i=0; i<rows_; i++)
-        entries += indices_[i].size();
-
+        entries += rowsize(i);
       return entries;
     }
 
     /** \brief Return the number of rows */
     size_type rows() const {return rows_;}
 
-
     /** \brief Return the number of entries in a given row */
-    size_type rowsize(size_type row) const {return indices_[row].size();}
+    size_type rowsize(size_type row) const {
+      return std::visit([&](const auto& rowIndices) {
+        return rowIndices.size();
+      }, indices_[row]);
+    }
 
     /** \brief Import all nonzero entries of a sparse matrix into the index set
         \tparam MatrixType Needs to be BCRSMatrix<...>
@@ -92,17 +186,16 @@ namespace Dune {
       matrix.setSize(rows_, cols_);
       matrix.setBuildMode(MatrixType::random);
 
-      for (size_type i=0; i<rows_; i++)
-        matrix.setrowsize(i, indices_[i].size());
+      for (size_type row=0; row<rows_; row++)
+        matrix.setrowsize(row, rowsize(row));
 
       matrix.endrowsizes();
 
-      for (size_type i=0; i<rows_; i++) {
-
-        typename std::set<size_type>::iterator it = indices_[i].begin();
-        for (; it!=indices_[i].end(); ++it)
-          matrix.addindex(i, *it);
-
+      for (size_type row=0; row<rows_; row++) {
+        std::visit([&](const auto& rowIndices) {
+          for(size_type col : rowIndices)
+             matrix.addindex(row, col);
+        }, indices_[row]);
       }
 
       matrix.endindices();
@@ -111,9 +204,10 @@ namespace Dune {
 
   private:
 
-    std::vector<std::set<size_type> > indices_;
+    std::vector<RowIndexSet> indices_;
 
     size_type rows_, cols_;
+    size_type maxVectorSize_;
 
   };
 

dune/python/istl/matrixindexset.hh

@@ -3,6 +3,7 @@
 #ifndef DUNE_PYTHON_ISTL_BCRSMATRIX_HH
 #define DUNE_PYTHON_ISTL_BCRSMATRIX_HH
 
+#include <cstdint>
 #include <memory>
 #include <stdexcept>
 #include <string>
@@ -29,7 +30,7 @@ namespace Dune
     void registerMatrixIndexSet(pybind11::handle scope,
                             pybind11::class_<MatrixIndexSet, options...> cls)
     {
-      typedef std::size_t size_type;
+      using size_type = Dune::MatrixIndexSet::size_type;
 
       // two different possible constructors
       cls.def( pybind11::init( [] () { return new MatrixIndexSet(); } ) );