License : Creative Commons Attribution 4.0 International (CC BY-NC-SA 4.0)
Copyright : Frédéric Pennerath, CentraleSupelec
Last modified : April 9, 2024 11:25
Link to the source : matrix.md

Write your own templates:
The example of matrix template

The goal of this programming exercise is for you to learn how to write the main forms of templates and to use some of their features, like partial specialization or conditional branching (if constexpr).

To this end, we will develop a class template for mathematical matrices, that is very similar in its principle to the class template Vector introduced during the lectures. In particular we consider matrices whose size \(n \times m\) (\(n\) is the number of row, \(m\) the number of column) is fixed at compile time and can therefore be stored in a fixed-size array. Matrices can then be created very quickly as local variables directly on the call stack without storing components in the heap, using dynamic allocation (new). This feature can then be used intensively to proceed intermediary computations very efficiently. Integrating matrix dimension into the template parameters also makes these dimensions as an inherent part of the matrix type so that matrices with different dimensions are considered to be of different types. This is particularly interesting to detect dimension inconsistencies at compile time.

As a summary, the purpose of the template class Matrix is to provide a generic, while at the same time, efficient and strongly typed implementation of matrices, where genericity refers to the freedom for the users of your template, to choose

• the type Scalar of components on which the matrix is defined (like double, float, or any other algebraic field like std::complex)
• the dimensions n and m of the matrix of type int.

Class template declaration and definition

In a file Matrix.hpp implement from scratch a class template Matrix parameterized by n, m and Scalar (by default equal to double) that defines

• An alias type scalar_type to “export” type Scalar.
• A constructor that sets its matrix to a scalar constant. If omitted, the scalar constant will be initialized by calling its default constructor Scalar{}. One will assume that this default constructor will return the zero for the considered field, i.e. the identity with respect to addition (this is actually true for standard numerical types like int, float, etc).
• A copy constructor using the default definition.
• A read/write accessor operator[](std::pair<int,int> coords) that returns a reference of type Scalar to the component of the matrix of coordinates coords (i.e the pair \((i,j)\) where \(i\) is the row index and \(j\) the column index, both starting from \(0\)).

Few important preliminary comments:

• Do not forget that templates must be declared and defined in the header file .hpp (i.e do not use .cpp file to define template functions and methods).
• For code clarity, methods and internal operators of the class will be defined outside the class declaration, right after it.
• Scalar components of the matrix will be stored in a one-dimensional std::array, which is just a class wrapping a native array whose sole purpose is to provide utility methods like begin, end and size so that loops like for(auto& x : array) { ... } can be used. Components will be stored in this array, from top to bottom, left to right.
• Operator operator[] is only provided to the user in order to initialize a given component in a matrix. For sake of efficiency, you should not call it from any method or operator of Matrix you will develop. Write your loops by incrementing iterators instead (which are actually native C pointer for std::array).
• To avoid any naming conflict, one will declare all matrix classes and functions in a dedicated namespace matrix.
• One will test every new method/operator with some unit tests. A test file test-Matrix.cpp is provided but feel free to modify it or to write your own tests. If using the provided test file, you just need to update macro #define QUESTION to define which tests you want to compile (as their compilation will depend on your progress).
• A solution is provided here Matrix.hpp in case you get stuck.

Function / external operator template

Add an external operator template << that serializes a matrix into an output stream std::ostream.

Notes:

Do not forget to add the right friend statement (see the similar example of class Vector in the lectures). Test (see test-Matrix.cpp).

Method / internal operator templates

Implement internal operators computing the sum and product of two matrices.

Notes:

• One recalls that given two matrices \(A\) and \(B\) of respective size \((n_A,m_A)\) and \((n_B,m_B)\) where \(n\) is the number of rows and \(m\) is the number of columns, then

  • Sum \(A + B\) requires that \(n_A = n_B\) and \(m_A = m_B\). Then \(A + B\) is a matrix of size \((n_A,m_A)\) such that its coefficient \((i,j)\) is \[(A + B)_{i,j} = A_{i,j} + B_{i,j}\]

  • Product \(A \times B\) requires that \(m_A = n_B\). Then \(A \times B\) is a matrix of size \((n_A,m_B)\) such that \[(A \times B)_{i,j} = \sum_{k=1}^{m_A} A_{i,k} \times B_{k,j}\]

• Test what happens when one attempts to multiply two matrices with inconsistent dimensions?

Code branching with if constexpr

When computing a product \(A \times B\) of two matrices such that \(A\) is a column vector (i.e. \(m_A = 1\)) and thus \(B\) is a row vector (i.e. \(n_B = 1\)), the standard implementation is particularly unefficient. Can you see why?

Using if constexpr as introduced in C++17, specialize the code of the product operator when \(m_A = 1\) so that the code is more efficient.

Notes:

• You can print debug messages to check the right code branches are called.

Generation of custom compilation error using static_assert

In order to produce interpretable error message, implement an internal operator operator* that accept as right operand any matrix, including those with sizes incompatible with the right operand (i.e the current object *this). In that operator, introduce a static_assert statement that produces a compiler error if dimensions are inconsistent.

Partial specialization of class template

Matrices, as arrays (or std::array), store their components directly inside the object Matrix. This is only possible because the matrix has a fixed size known at compile time. This property allows the matrix components to be stored directly on the stack and thus to gain in performance (no dynamic allocation, no memory indirection). However, storing on the stack can quickly become a drawback when the matrix contains many coefficients to the point of causing stack overflow.

So we would like a class template Matrix to store its components in the heap instead of in the stack as soon as the number of its components exceeds a certain threshold (for example 10000). To do this, one simple solution is to replace std::array by std::vector when the size exeeds a threshold. We could specialize directly the class matrix but this would require to duplicate all methods and operators. Instead we can introduce a template class Storage whose sole purpose is to store components either in an array or in a vector. Class Matrix then inherits from Storage.

Modify your implementation accordingly, by defining class template Storage parameterized by Scalar, the dimension of the matrix and a fourth boolean argument stack that tells whether Storage must store components in an array (if stack is true) or in a vector (if stack is false). Use partial specialization of class template (see example in the lectures) to implement both versions (with array and with vector). Then make class template Matrix inherits from Storage by setting the stack argument to the right value. To to this, one will introduce a constexpr function constexpr bool best_on_stack(size_t n_rows, size_t n_columns) that will returns true if and only if the components are better stored on the stack.

Notes:

• Do not forget to inialize the vector to the right size using method resize.
• Observe that if constexpr is pointless here as it cannot help to specialize attributes. Only template specialization can do it.
• You can add a static boolean flag static const bool on_stack in each specialization, that tells whether the components are stored in an array or not.

Concepts

We now wants to formalize and document syntactic constraints that Scalar must fullfil to be an algebraic field.

Introduce a concept Field that lists syntactic constraints of an algebraic field (See the example in lectures on generic programming). Then constraint parameter Scalar of template class matrix to satisfy concept Field. Test by instantiating an absurd matrix (e.g a matrix defined on std::string). See error messages.

Notes:

• Concepts have been introduced in C++20. Be sure to compile with the option --std=c++20 and to use a recent compiler.

Metaprogramming

When computing the product of two matrices \(A\) and \(B\), the termination test of the internal loop is particularly costly (why?).

Update your class template to do loop unrolling in the internal loop.

• Be careful: the internal loop might have many iterations so that loop unrolling should not unroll all iterations. We only want to group k successive iterations using loop unrolling with k typically equal to 4 or 8.
• In a first version we can assume \(m_A\) is a multiple of k. Then manage properly rounding when k does not divide anymore \(m_A\).