Linear Algebra support in C++26

It’s likely that you’re already familiar with std::string_view and std::span<T> in Modern C++, each of which provides a lightweight “view” onto either string data or other contiguous typed data (such as made available by std::vector<T>::data()). A key difference is that std::string_view is immutable, even if initialized with non-const data. (You could instead use std::span<char> to get around this, if necessary.) A limitation of std::span<T> is that it only supports a single-dimensional view, so is not guaranteed to work with a construct like the following:

void f() {
  int array2d[2][3];
  std::span<int> arr_span{ reinterpret_cast<int*>(array2d), 6 };  
  // ...
}

Of course, you should try to avoid casting if at all possible, so you’ll probably be glad to know that Modern C++ will soon have another solution in the form of std::mdspan (“Multi-Dimensional Span”). It has been possible to “roll your own” matrix (or similar classes) using template syntax since C++98, for example (some newer syntax is present):

template<size_t rows, size_t columns, typename T>
class Matrix {
  T *data = nullptr;
public:
  Matrix() : data{ new T[rows * columns] } {}
  ~Matrix() { delete[] data; }
  const T& operator()(size_t row, size_t column) const {
    return data[row * columns + column];
  }
  T& operator()(size_t row, size_t column) {
    return data[row * columns + column];
  }
  // other functionality below ...
};

void g() {
  Matrix<2,3,int> matrix;
  matrix(1, 2) = 10;  // calls non-const operator()
  std::cout << matrix(1, 2) << '\n';  // calls const operator()
}

It is a design-time decision as to whether to use compile-time (template non-type parameters) or runtime (member variables) logic to allow for element access. Also the choice of row-major (as above) or column-major layout is available to the Matrix class designer, as is indexing from zero or one.

The good news is that this is all standardized and improved upon by std::mdspan, provided by header <mdspan> new with C++23. The above code written with this class would be:

#include <mdspan>
#include <iostream>

void h() {
  int array[] = { 1, 2, 3, 4, 5, 6 };
  std::mdspan matrix(array, 2, 3);
  matrix[1, 2] = 10;
  std::cout << matrix[1, 2];
}

A few things to take onboard from this:

1. std::mdspan is non-owning so you must allocate the memory it uses and assign to the data elements (including zeroing them) beforehand
2. A list of any number of dimensions (“extents”) follows the data pointer
3. (Not shown here but “dynamic”, runtime determined extents are possible)
4. Square brackets are used instead of parentheses (this implies that new language-level support is needed)
5. Indexing from zero is mandated

So you’ve read this far and are maybe wondering: “What has this to do with Linear Algebra?” The answer lies in another new header, <linalg> coming with C++26, with transformation and BLAS functions as additions to the (new) std::linalg namespace. The reason these are being added added to the C++ Standard Library is because element and stride accesses are standardized due to use of std::mdspan.

There are quite a few BLAS 1, 2 and 3 functions slated for inclusion, including dot and matrix_product (see links below for a complete list). These will doubtless be important to many students and scientists using Modern C++ for their work, and will assist in promoting use of the language as a whole.

That’s all for this article. We’ve briefly covered the <mdspan> and <linalg> headers, and the use of std::mdspan instead of a home-made multi-dimensional array class.

References:

1. https://en.cppreference.com/w/cpp/container/mdspan
2. https://en.cppreference.com/w/cpp/numeric/linalg