Draft: Introduce a new BCRSMatrix

Description

The goal of this MR introduces is to introduce a BCRSMatrix-like matrix, called IBCRSMatrix, that fixes several of the issues noticed with BCRSMatrix over the years:

  • Column indices can be chosen at compile time, with a default of uint32 which saves about 25% of memory usage on flat matrices.
  • Row offsets are now part of the shared pattern, and are now stored more compactly.
  • Sparse patterns are now an external class with a slim interface for general patterns, allowing the pattern to be shared with other classes, e.g. for mixed-precision algorithms, and to simplify implementation of dense, banded, and diagonal matrices.
  • There are separate builder classes for sparse pattern creation, making the matrix itself much smaller, which is important for sparse-sparse matrices.
  • The sparse pattern builders have proper buffer growth and compression, removing the need of the overflow map or scattered memory allocations.
  • Full support of move semantics with much better exception guarantees.
  • Full support of stateful allocators, even on move and swap operations.
  • Follows the AllocatorAwareContainer named requirement from the standard library, thus propagating allocators to its stored objects.

Trade-offs

It is designed to be an almost drop-in replacement of BCRSMatrix with a few caveats that are usually easy to modify:

  • Rows are proxy objects that are generated on-the-fly, this means that they cannot be bound to non-const references.
  • Row and const row iterators are different, respecting const correctness. This may cause some issues in code that used the relaxed const correctness guarantees of the BCRMatrix iterators.
  • Implicit build mode does not store values during pattern formation and only default initialized blocks are accepted.
  • The row-wise build mode has different begin and end iterator types.
  • The setSize() method should be used after selecting the build mode and before the pattern construction is finalized.
  • The N() and M() methods are not available during the pattern construction.
  • The build mode becomes BuildMode::unknown after the pattern is built, meaning that it needs to be reset before constructing a new pattern.
  • Setting a build mode cleans up the matrix.
  • There is no BuildStage and their related functions.
  • There is no one-before iterators as they would be undefined behavior or would make the implementation much more difficult (see dune-common#430).

Design

To separate concerns and make the implementation and maintenance simpler, I divided this into several key classes, each one with its own tasks and slim interface to interact with the others. This is a first implementation so whether the name is fine or whether they belong to the Dune or Dune::Imp namespace is open for discussion.

  • IndexedArrayVew: This is a class that unifies base_array_unmanaged and compressed_base_array_unmanaged to have a view over a stream of values indexed by another (maybe sparse) range. A key change to allow the unified object is that they operate on data iterators and a borrowed index range views.
  • VectorView: This is a class that unifies block_vector_unmanaged and CompressedBlockVectorWindow. This simply extends IndexedArrayVew to have vector operations. Importantly, since the base class unifies sparse and dense cases, they are interoperable. This can later be used to implement BlockVector in terms of VectorView and allow interoperability between matrix rows and block vectors.
  • MatrixView: This class takes a data iterator and a pattern view with nnz indices, which then exposes the data with matrix-like operations. Importantly, it doesn't own the data and the pattern is abstracted away with a very slim interface that exposes each row indices and the offset to the data. Since it uses IndexedArrayVew as a base class to create rows on-the-fly, patterns that create indices on-the-fly do not require to store indices at all, allowing to use this for dense matrices, banded matrices, or any other exotic pattern (e.g. compile time-known) that fulfills the pattern interface.
  • SparseIndexRanges and CompressedSparseIndexRanges: These are instances of patterns that expose each row indices, which are internally stored contiguously. The compressed version stores the indices in the smallest possible index type that fits all indices and convert them on-the-fly via a variant.
  • SequencedSparseIndexRangeBuilder and UnsequencedSparseIndexRangeBuilder: They take the responsibility to construct patterns and provide an instance of either SparseIndexRanges or CompressedSparseIndexRanges. The idea is to separate construction of patterns from the actual data used by other classes, and to provide an builder interface that is less error-prone to memory errors. A key part of their implementation is the growth of buffers when inserting indices. The sequenced buffer simply grows the back buffer and copies the data over, while the unsequenced builder grows all index sets based on the index set with the maximum number of entries. This is designed to be amortized when inserting random distributions of equal sized index sets and for sequential insertion. Note that the random and implicit matrix build modes are then unified in the unsequenced builders. This is because both need to insert indices in a pre-allocated buffer and need to grow once they are exhausted. Setting the index set sizes in the random mode before adding indices "simply" re-distributes the buffers to those sizes. Thus, not giving a size, like in implicit matrix build mode, just distributes the buffer evenly.
  • IMatrix: This is an matrix which owns the data. It basically extend MatrixView by properly allocating the data, and delegating everything else to the MatrixView. The name Matrix is already taken, but I would hope that gets replaced by this class which has much better allocation behaviors. The prefix I comes from "indexed" because it has a pattern as argument, which as said above, can be used to index matrix entries in a generic way. Other names are welcome!
  • IBCRSMatrix: This is a class based on IMatrix which explicitly uses SparseIndexRanges as a pattern. In particular, it holds an internal sparse index range builder that allows it to be an almost drop-in replacement of BCRSMatrix.

Many design choices have been made to make this possible, so discussions are welcome and encouraged. If a particular class needs a more detailed discussions, we can separate this into several MRs that address these issues individually. For instance, note that, aside some assumptions, the matrix and vector interfaces are relatively independent of the sparse index ranges implementation and grow mechanisms.

Requires

TODOs

  • Find work arounds for older compilers/standard libraries
  • Find a work around to implement setIndices and setIndicesNoSort.
  • Add CHANGELOG entry
Edited by Santiago Ospina de los Ríos

Merge request reports

Loading