diff --git a/dune/istl/matrixindexset.hh b/dune/istl/matrixindexset.hh
index 10baccf35c6e56b11da8b1ada0c5a8130d8dce63..acc4466b4a093ff0923042c36ef3fe08a01fa2d1 100644
--- a/dune/istl/matrixindexset.hh
+++ b/dune/istl/matrixindexset.hh
@@ -5,55 +5,139 @@
#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.
+ */
class MatrixIndexSet
{
+ using Index = std::size_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);
+ 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 +176,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 +194,10 @@ namespace Dune {
private:
- std::vector<std::set<size_type> > indices_;
+ std::vector<RowIndexSet> indices_;
size_type rows_, cols_;
+ size_type maxVectorSize_;
};