Draft: Matrix-Vector concepts

What is a matrix? What is a vector?

These questions are not so very easy to answer and there was already some discussions in some issues about what should a matrix provide as interface and also about some special types of matrices, e.g. sparse matrices. See the referenced issues. This MR is intended as a structured documentation in code of what a matrix (and vector) should provide. It is more the base for a discussion than the concrete realization. I'm using c++20 concepts here to describe the matrix and vector concepts but this is just a technical detail. It could be written in c++17 SFINAE-concepts as well.

The definition of these concepts is inspired by some MTL4 concepts, see, e.g., http://old.simunova.com/docs/mtl4/html/group___concepts.html, but it applied to the Dune interface. Not everything is properly specified yet.

Motivation

• For some algorithms in dune-common and dune-istl involving matrices and vectors, it might be necessary to distinguish different types / properties of these types. E.g. some algorithms that traverse a matrix must know whether the iterators begin() and end() provide iteration over the rows of the matrix or the columns - is it a row-major or column-major matrix. For other algorithms we need to know whether we get a vector or a scalar in terms of accessibility of elements.
• In dune-istl we have several adapters and wrappers taking a matrix as argument, but actually only very special types of matrices are allowed.
• Not all matrices are derived from DenseMatrix and not all vectors are derived from DenseVector.
• It gets more complicated if we allow multi-type containers, like MultiTypeBlockMatrix or TupleVector. Are these types vectors or collections? Do these types provide the necessary interface for an algorithm?
• There are several container-like data-structures in dune-common (ArrayList, SLList, TupleVector, ReservedVector, DiagonalMatrix, TransposedMatrixWrapper,...). Neither from the name nor from the documentation of these classes it can be deduced what interface is provided, which guarantees are given and in which way these data-structures can be used. So I end up with reading directly the code of these implementation to find out.
• Reading the doxygen documentation for vector and matrices often is a good start to understand the interface, but is also very complex, sometimes incorrect or misleading and partially undocumented. It must be repeated for each implementation what is error-prone and often is not up-to-date.

Open Questions

• What is a sparse matrix? I would say: 1. it is a matrix, 2. it provides an extended interface to get a notion of the sparsity, e.g. nonzeroes() and exists(i,j)
• Is there some concept for "sparse traversable"? We have DenseIterators. Should sparse vectors/matrices provide also this interface or is it different for sparse data-structures? The iterators at least should provide the notion of an index, i.e., iter.index()
• There are wrapper/view types, like TransposedMatrixWrapper or ScalarVectorView. Those should also provide some minimal interface, to be usable in algorithms requesting matrices/vectors. If not at least a minimal interface is provided, they should not be called Matrix or Vector and maybe should be moved in the Impl namespace (if not already there)
• [...]

Details

The implementation is still incomplete, but already shows some missing interface method in the dune classes.

The concepts are implemented in terms of orthogonal properties:

• A matrix and a vector both are "collections" of elements.
• They have algebraic properties, like size/shape
• Matrices are (following a definition on wikipedia) tables of elements. Thus all matrices have an interface to access these elements (const/mutable) by [i][j], whereas vectors are similar and have access to its elements by [i]. Note, by this defintion all matrices have these access. It is just an interface conept. Not necessary the optimal way of working with these data-structures.
• An orthogonal concept is to see the matrix as a linear operator. This defines the interface using mv and umv...
• Another orthogonal property is to see matrices/vectors as elements of a vector-space with all its vector-space operations, i.e. +, -, ...
• One more orthogonal concept is how to traverse the containers. Does the matrix provide iterators? Is this rowwise or columnwise iteration? Maybe is there an iterator to traverse all nonzeros of the matrix directly?

By separating orthogonal concepts, one could see a sparse matrix like DiagonalMatrix as a Matrix, as and element of a vector space and as a linear-operator. Additionally, it could fulfill a concept of a sparse-matrix or it can be iterated like a sparse matrix. So, an algorithm could specify the requirements for its arguments in terms of these orthogonal properties.

Applications

Krylov methods

For most iterative krylov methods we just require some operations for the operator and the vectors, it does not need to be a full matrix implementation, for example:

template <class L, Concept::HilbertSpace X, Concept::HilbertSpace Y>
  requires Concept::LinearMap<L,X,Y>
void cg_solver(const L& linOp, X& x, const Y& y)
{
  // implementation...
}

Other algorithms require more, e.g. the trasponsed operations:

template <class L, Concept::HilbertSpace X, Concept::HilbertSpace Y>
  requires Concept::TransposableLinearMap<L,X,Y>
void qmr_solver(const L& linOp, X& x, const Y& y)
{
  // implementation...
}

Matrix factorization

Depending on your algorithm, a factorization requires either a traversable collection or a matrix with index acces, e.g.

template <Concept::Matrix M>
void lu_factorization(M& matrix)
{
  // implementation...
}

C++20 Concepts

While the real concept requirements/constraints from c++20 can only be used with a recent compiler, one could provide concept-replacements using SFINAE with a similar syntax, i.e. as constexpr variable templates. Those could be used inside static_assert, e.g.

template <class L, class X, class Y>
void cg_solver(const L& linOp, X& x, const Y& y)
{
  static_assert(Concept::LinearMap<L,X,Y>);
  static_assert(Concept::HilbertSpace<X>);
  static_assert(Concept::HilbertSpace<Y>);

  // implementation...
}

Related issues

• #145 (Access to off-diagonal entries of a DiagonalMatrix)
• !552 (make access to off diagonal entries of DiagonalMatrix return zero)
• #74 (Future of IsMatrix traits / Providing a better alternative)
• dune-istl!424 (merged) (Cholmod: Add support for almost all matrix types )
• #131 (Complete binary operators for dense matrix/vector types)
• !816 (closed) (Add support for C++20 concepts)

Lot's of these issues have started a long discussion, mainly because it is not clear in Dune what we mean when speaking about (sparse) matrices and vectors.