What is a matrix?
What is a matrix?
While we have CRTP base class DenseMatrix
specifying the interface for matrices, this base is incomplete or wrong documented and even inconsistent. The specification is also not applied to all matrixlike classes neither in dunecommon nor in duneistl. So, I wanted to open a specific discussion on what is a matrix?
1. Matrix is a 2d container of values
On Wikipedia, a matrix is defined as "a rectangular array of numbers (or other mathematical objects)". Thus we have the following requirements:
 (a) A matrix stores elements/entries of a specific type. Seeing matrices as collections, this would mean we have a single elementtype. For collections this type is typically denoted by
value_type
.  "The size of a matrix is defined by the number of rows and columns that it contains". Actually, this definition specifies three size informations:
 (b) There is a type representing the sizes, typically denoted by
size_type
 (c) There is number of rows and number of columns. In Dune this is (for historical reasons) denoted by
N()
andM()
respectively.  (d) There us the total size of a matrix. Since a matrix shape is denoted by
n x m
, I would suggest to define the size as the productn*m
. This would allow to see singlerow and singlecolumn matrices as rowvector and columnvector, correspondingly. Compare, e.g.,numpy.matrix.size
,mtl::size(matrix)
, orEigen::Matrix::size
.
 (b) There is a type representing the sizes, typically denoted by
 (e) The elements of a matrix should be accessible by indices. Since the row and column indices are in the range
[0,rows)
and[0,cols)
, its index type should also besize_type
.  (f) Following the Carray syntax, elements are accessed by
matrix[i][j]
.
2. Matrix is an element of a vector space
Matrices over a field build an algebraic set with vectorspace properties:
 matrices for an additive abelian group (commutative group)
 matrices are scalable by scalars of the field type
 Multiplication with scalar is distributive with respect to addition.
These mathematical concepts result in the following requirements for the matrix type M
:
 There is a field type associate to the matrix, i.e.
Dune::FieldType<M>::field_type
exists.  There are arithmetic operations for matrices
A
andB
: (g)
A + B
,A  B
 (h)
A += B
,A = B
 (i)
A
,+A
(optional)
 (g)
 Multiplication with a scalar
s
of field type: (j)
A * s
,s * A
 (k)
A *= s
 (l)
A / s
,A /= s
 (j)
 (m) Additionally, we typically require compound operations:
A.axpy(s, B)
2.1. Normed space and innerproduct space
By defining norms over matrices and inner products, the vector spcae could be turned into a normed vector space or an innerproduct space

Typical norm for matrices include:
 (n)
frobenius_norm()
 (o)
forbenius_norm2()
(squared frobenius norm)  (p)
infinity_norm()
 (q)
infinity_norm_real()
(simplified infinity norm)
 (n)

Innerproducts:
 (r) frobenius innerproduct
A:B = tr(A^H A)
(not implemented by any matrix)
 (r) frobenius innerproduct
3. Matrix is a linear map
A linear transformation or linear map between vector spaces can be represented by a matrix. The operations of a linear map are typically of the form of matrixvector products:
 (s)
A.mv(x, y)
:y = A*x
 (t)
A.umv(x, y)
:y+= A*x
 (u)
A.mmv(x, y)
:y= A*x
 (v)
A.usmv(alpha, x, y)
:y+= alpha * A*x
If the matrix additionally provide transposed operations:
 (w)
A.mtv(x, y)
:y = A^T * x
 (x)
A.umtv(x, y)
:y+= A^T * x
 (y)
A.mmtv(x, y)
:y= A^T * x
 (z)
A.usmtv(alpha, x, y)
:y+= alpha * A^T * x
If the matrix additionally provide hermitian operations:
 (A)
A.mhv(x, y)
:y = A^H * x
(currently not implemented by any matrix)  (B)
A.umhv(x, y)
:y+= A^H * x
 (C)
A.mmhv(x, y)
:y= A^H * x
 (D)
A.usmhv(alpha, x, y)
:y+= alpha * A^H * x
4. Other operations with matrices
 (E) Matrix multiplication:
A in R^nxm
andB in R^mxk
thenA * B in R^nxk
 Matrices must be compatible
 square matrices might form a ring under matrix addition and matrix multiplication
this would allow
A *= B
 (F) The transposed of a matrix:
A in R^nxm
>A^T in R^mxn
. Either an operation
A.transpose()
or a free functiontranspose(A)
?  The transposed matrix is itself a matrix
 Either an operation
 (G) Trace of a matrix
tr(A) = sum_i A_ii
 Assumption: matrix is square
 Either an operation
A.trace()
or a free functiontrace(A)
?  Recursive definiton:
tr(A) = sum_i tr(A_ii)
 (H) Determinant of a matrix
 Assumption: matrix is square
 Easy to define for small matrices
 Generalization using, e.g., Laplace expansion
 Special operation: gramian determinant
 (I) Inversion of a matrix or solution of a linear system
 Assumption: matrix is square
 easy to implement for small matrices
 LU decomposition (or other forms of factorization)
5. Matrices can be traversed
Visiting all the (stored) elements of a matrix is a useful operation to generalize some algorithms
 Whats is a proper call for choosing the way/order of traversal?
 It exists
begin()
andend()
functions. Are not documented properly but typically mean traversal of rows.
 Additionally reverse traversal is given by
beforeEnd()
andbeforeBegin()
 Why not using
rbegin()
andrend()
as for standard containers?
 Why not using
 Missing functions in the interface:
 Traversal over rows first: e.g.
begin_rows()
,end_rows()
 Traversal over columns first: e.g.
begin_cols()
,end_cols()
 Traversal over rows first: e.g.
6. Special matrix types with more properties
 Sparse matrices:
 Check whether elements exist in the store is required:
A.exists(i,j)
 Number of stored elements (nonzeros):
A.nonzeroes()
 Iterators only iterate over stored elements/nonzeros.
 Often direct elementaccess not implemented (neither const nor mutable)
 Check whether elements exist in the store is required:
 Banded/Diagonal matrices:
 Special type of sparse matrix
 Additional access to diagonal(s):
A.diagonal()
(current interface does not allow return offdiagonals)  Direct access to diagonal entries:
A.diagonal(i)
 Triangular matrices
 Tridiagonal matrices
Examples in dunecommon
DenseMatrix
A DenseMatrix
is the CRTB base class for some matrix implementation, but it lacks some of the properties from above. Also
its documentation is partially wrong.
 (c) There are multiple functions returning number of rows/columns:
A.N()
andA.rows()
,A.M()
andA.cols()
, whererows() > mat_rows()
andcols() > mat_cols()
call the implementation functions.  (d) is defined differently. For some reason
A.size()
returns the number of rows.  (g) arithmetic operations not completely available for matrices,
+A
not implemented (not required)  (j) and (l) not implemented, except for
A /= s
 (r) There is no inner product for matrices implemented
 (A) Missing linmap operation
A.mhv(x,y)
 (E) binary matrixmatrix product not implemented. For square matrices the multassign is implemented from left and
right, i.e.,
A.leftmultiply(B)
andA.rightmultiply(B)
. The general matrixmatrix multiplication is commented out.  (G) Trace is not implemented.
FieldMatrix
Is an implementation of DenseMatrix
thus inherits its properties. Some additional comments
 (c) Here
rows
andcols
are enum constants hiding the interface functionsrows()
andcols()
 (g) FieldMatrix implements
A + B
andA  B
 (j) and (l) are completely implemented
 (E) Multiplication is implemented as
A * B
andleftmultiply[any]()
andrightmultiply[any]()
. Additional operations are implemented in the namespaceDune::FMatrixHelp
DynamicMatrix
Is an implementation of DenseMatrix
thus inherits its properties.
DiagonalMatrix
 (d) is defined differently. For some reason
A.size()
returns the number of rows.  (f) The element access is only implemented for diagonal entries, i.e.
matrix[i][j]
throws ifi != j
 (g), (j) Operations
A + B
,A  B
, andA
not implemented.  (i) compound assignment operators implemented for scalars:
A += s
andA = s
meaning elementwise addition and subtraction of that scalar.  (m) Operation
axpy
not implemented.  (r) There is no inner product for matrices implemented
 (A) Missing linmap operation
A.mhv(x,y)
 (E) No matrixmatrix products implemented
 (G) Trace is not implemented.
TransposedMatrixWrapper
The wrapper is created by transpose(matrix)
, but does not produce a matrix in any sense as described above. It
has just one single operation: A * B^T
. Note, this wrapper only works for transposing a FieldMatrix
.
Test of the interface
There is a test in dunecommon: checkmatrixinterface.hh
but it does not test the complete matrix interface, but just the functions
that were already implemented. Thus, this test is probably written after the implementation of the DenseMatrix
class and
not before.
An alternative test, that is implemented in terms of separating the various aspects of a matrix, is implemented in !953 (closed) in terms of c++20 concepts. These tests can run in the ci jobs for toolchains supporting this new c++ standard and can report interface failures.
Major problems with the current interface
The existing matrix implementations in dunecommon (but also in duneistl) only partially fulfil the interface definition from above and are inconsistent. Major problems are:
 Element access is not a defined property of a matrix. It is only guaranteed to succeed in some matrix implementations. This is mainly for performance reasons. This argument is a bit weak, because there are other ways to ensure performant data access, e.g. using iterators or specialized access functions.
 Traversal of matrices is not properly specified and cannot easily be extended to columnmajor traversal.
begin()
andend()
have no knowledge about its traversal direction. There should be a better name for rowwise traversal and columnwise traversal. Proposal:begin_rows()
andbegin_cols()
. Then, the user has to explicitly use a rangegenerator for the iteration, e.g.,rows(matrix).begin()
andcolumns(matrix).begin()
 The iterators itself are not properly specified. Is the dereferenced iterator again a range? What does it traverse?
 Some matrices implement a single index access
matrix[r]
as a rowview onto the matrix, i.e., this operation returns a vectorlike container representing ther
th row of the matrix. This is not clearly specified. Is this a required property of a matrix? What are the costs of this operation? Is the row a propery vector or just a proxy providing minimal access to the row entries?  Especially in the implementation of geometries, one needs matrixmatrix multiplication. While this is implemented for
FieldMatrix
, all the other matrix types lack this operation. But, e.g.,YaspGrid
returns for its jacobian aDiagonalMatrix
. Other grid implementations maybe something else. We needA * B
, but alsoA^T * B
andA * B^T
for various combinations of matrices. This is partially implemented in some "Helper" namespaces in dunegeometry or forFieldMatrix
inFMatrixHelper
namespace. Multiplication seems not properly designed!  Sparse matrices are a special type of matrices. The interface should mirror this mathematical property. Actually, sparse
matrices are just matrices with a different storage type for its entries. Thus, sparse matrices should only extend the
matrix interface with methods providing fast access to the stored elements. Either by iterators or something else. Currently,
sparse matrices only implement a partial matrix interface and some methods even have different behaviour, e.g.
matrix[i][j]
is not guaranteed to succeed, not even for constaccess, although the matrix value is clearly specified.  Views over matrices, like transposedview or hermitianview are not properly implemented. They do not provide a matrix interface although this could easily be provided. One could discuss whether views are readonly or provide mutable access to the matrix elements as well.
 Currently, we can only test that a matrix supports all the different aspects as described above. Not the individual aspect. But, often an algorithm just needs some properties, e.g., in iterative krylov methods, one needs the linearmap properties for a matrix, for ludecomposition maybe fast element access. In geometry parametrizations we need matrixmatrix multiplication and trace/Gramian determinant computation. It is required to have individual concept/interface checks.
 Singlerow and Singlecolumn matrices are not vectors. This often leads to problems, e.g., there is no
two_norm()
implemented and alsodot()
andaxpy
with another vector does not work, even if the shape would allow this. This problem is that vectors are not matrices and also not compatible in its interface to matrices.