A matrix class using const members

Question

I've long used a simple matrix class with a vector and rows/cols for the shape. But I've always disliked getters even when necessary to prevent exposing invariants. Then came c++20 which provided the ability to modify consts in a class with the change in Basic.Life. So then I modified my matrix class to get rid of the rows and cols getters. However, I was unable to use a const vector as the move operations could not be done hence pessimizing performance.

So I just use new with a const ptr instead of a vector. Now the entire matrix class, except for the dynamic memory (contents of the matrix), are declared const so safer programming in class code, friend functions, and no longer need getters for safe user code. And the sizeof each matrix object is 40% smaller than the vector version.

Note that declaring the values in the matrix class const is consistent with them being unchanged from start of lifetime to end of lifetime. I believe the code provides a string exception guarantee but haven't tested this. Also, the matrix multiply, has a transpose mode that aligns the matrixes such that multi-threading (not shown) them is fairly trivial and it also keeps the cache happy. Would like any comments, especially of flaws that might be left.

CompilerExplorer, gcc, clang, msvc

header file: matrix2db.h

#pragma once

#include <memory>
#include <numeric>
#include <initializer_list>
#include <exception>
#include <stdexcept>

namespace matrix_const
{

    template <typename T>
    class Matrix2D {
        T*  const pv;
    public:
        const size_t cols;
        const size_t rows;
        Matrix2D() noexcept : pv(nullptr), cols(0), rows(0) {}
        Matrix2D(size_t a_rows, size_t a_cols) :
            pv(new T[a_rows * a_cols]), cols(a_cols), rows(a_rows) {}

        Matrix2D(const Matrix2D& rhs) : pv(new T[rhs.rows * rhs.cols]), cols(rhs.cols), rows(rhs.rows) {
            std::copy(rhs.pv, rhs.pv + rows * cols, pv);
        }
        Matrix2D& operator=(const Matrix2D& rhs) {
            T* newp = new T[rhs.rows * rhs.cols];    // if new fails, *this is unchanged
            delete[] pv;
            std::construct_at(&rows, rhs.rows);
            std::construct_at(&cols, rhs.cols);
            std::construct_at(&pv, newp);
            std::copy(rhs.pv, rhs.pv + rows * cols, pv);
            return *this;
        }
        Matrix2D(Matrix2D&& rhs) noexcept : pv(rhs.pv), cols(rhs.cols), rows(rhs.rows) {
            std::construct_at(&rhs.rows, 0);    // maintain invariants in remnant
            std::construct_at(&rhs.pv, nullptr);
        }
        Matrix2D& operator=(Matrix2D&& rhs) noexcept {
            delete[] pv;
            std::construct_at(&rows, rhs.rows);
            std::construct_at(&cols, rhs.cols);
            std::construct_at(&rhs.rows, 0);    // maintain invariants
            std::construct_at(&pv, rhs.pv);
            std::construct_at(&rhs.pv, nullptr);
            return *this;
        }
        ~Matrix2D() {
            delete[] pv;
        }
    
        // additional ctor to provide list initialization
        Matrix2D(std::initializer_list<std::initializer_list<T>> list) :
            pv(new T[list.size() * (list.begin())->size()]),
            cols((list.begin())->size()),
            rows(list.size())
        {
            if (rows > 0 && cols > 0)
                for (size_t row = 0; row < list.size(); row++)
                {
                    if (list.begin()->size() != (list.begin() + row)->size())
                        throw std::runtime_error("number of columns in each row must be the same");
                    for (size_t col = 0; col < cols; col++)
                        pv[cols * row + col] = *((list.begin() + row)->begin() + col);
                }
            else
                std::construct_at(&pv, nullptr);
        }

        // bracket access v[r][c]
        T* operator[](size_t row) noexcept { return pv + row * cols; }
        const T* operator[](size_t row) const noexcept { return pv + row * cols; }

        // paren access v(r,c);
        T& operator()(size_t row, size_t col) noexcept { return pv[row * cols + col]; }
        const T& operator()(size_t row, size_t col) const noexcept { return pv[row * cols + col]; }

        // insert a sub-matrix by overlaying a selected portion, returns the whole matrix
        Matrix2D insert_matrix(const Matrix2D& mat, size_t row_start, size_t col_start)
        {
            if (row_start + mat.rows > rows || col_start + mat.cols > cols)
                throw std::range_error("Requested extents exceed bounds");
            Matrix2D ret(*this);
            for (size_t row = row_start; row < row_start + mat.rows; row++)
                for (size_t col = col_start; col < col_start + mat.cols; col++)
                    ret(row, col) = mat(row - row_start, col - col_start);
            return ret;
        }

        // return a matrix subset
        Matrix2D sub_matrix(size_t row_start, size_t col_start, size_t row_count, size_t col_count) const
        {
            Matrix2D ret(row_count, col_count);
            if (row_start + row_count > rows || col_start + col_count > cols)
                throw std::range_error("Requested extents exceed bounds");
            for (size_t row = row_start; row < row_start + row_count; row++)
                for (size_t col = col_start; col < col_start + col_count; col++)
                    ret[row - row_start][col - col_start] = (*this)(row, col);
            return ret;
        }

        // negate
        friend Matrix2D operator-(const Matrix2D& rhs)
        {
            Matrix2D ret(rhs.rows, rhs.cols);
            for (size_t i = 0; i < ret.rows * ret.cols; i++)
                ret.pv[i] = -rhs.pv[i];
            return ret;
        }

        // add
        friend Matrix2D operator+(const Matrix2D& lhs, const Matrix2D& rhs)
        {
            if (lhs.cols != rhs.cols && lhs.rows != rhs.rows)
                throw std::range_error("matrixes must have same dims");
            Matrix2D ret(lhs);
            for (size_t i = 0; i < ret.rows * ret.cols; i++)
                ret.pv[i] += rhs.pv[i];
            return ret;
        }
        // add +=
        Matrix2D& operator+=(const Matrix2D& lhs)
        {
            if (lhs.cols != cols && lhs.rows != rows)
                throw std::range_error("matrixes must have same dims");
            for (size_t i = 0; i < rows * cols; i++)
                pv[i] += lhs.pv[i];
            return *this;
        }
        // subtract
        friend Matrix2D operator-(const Matrix2D& lhs, const Matrix2D& rhs)
        {
            if (lhs.cols != rhs.cols && lhs.rows != rhs.rows)
                throw std::range_error("matrixes must have same dims");
            Matrix2D ret(lhs);
            for (size_t i = 0; i < ret.rows * ret.cols; i++)
                ret.pv[i] -= rhs.pv[i];
            return ret;
        }
        // sub -=
        Matrix2D& operator-=(const Matrix2D& lhs)
        {
            if (lhs.cols != cols && lhs.rows != rows)
                throw std::range_error("matrixes must have same dims");
            for (size_t i = 0; i < rows * cols; i++)
                pv[i] -= lhs.pv[i];
            return *this;
        }
    };

    // multiply matrixes
    template <typename T>
    Matrix2D<T> operator*(const Matrix2D<T>& lhs, const Matrix2D<T>& rhs)
    {
        if (lhs.cols != rhs.rows)
            throw std::range_error("cols of first must == rows of second");
        Matrix2D<T> ret(lhs.rows, rhs.cols);
        Matrix2D<T> mat_tmp = ~rhs;     // transpose to improve cache locality
        for (size_t row = 0; row < lhs.rows; row++)
            for (size_t col = 0; col < rhs.cols; col++)
                ret[row][col] = std::inner_product(lhs[row], lhs[row] + lhs.cols, mat_tmp[col], T(0));
        return ret;
    }

    // transpose
    template <typename T>
    Matrix2D<T> operator ~(const Matrix2D<T>& rhs)
    {
        Matrix2D<T> ret(rhs.cols, rhs.rows);
        for (size_t r = 0; r < rhs.rows; r++)
            for (size_t c = 0; c < rhs.cols; c++)
                ret(c, r) = rhs(r, c);
        return ret;
    }
}

And here's some short test code that exercises the special member functions:

#include "matrix2db.h"
#include <iostream>

void print(const auto& mat)
{
    for (size_t r = 0; r < mat.rows; r++)
    {
        for (size_t c = 0; c < mat.cols; c++)
        {
            std::cout << mat[r][c] << " ";
        }
        std::cout << "\n";
    }
    std::cout << "\n";
}

int main()
{
    using namespace matrix_const;
    {
        // Check move assignemnt operator
        Matrix2D<int> x { {1, 2}, { 3, 4 } };
        print(x);
        x = Matrix2D<int>{ { 2, 3 }, { 4, 5 }, {6, 7 } };
        print(x);
    }
    {
        // multiply a row vector by a column vector to produce a 1x1 result
        Matrix2D<int> x{ {1, 2, 3 } };
        Matrix2D<int> y{ {1}, {2}, {3}  };
        print(x*y);
    }
    {
        // check sub matrix operations
        Matrix2D<int> z{ {1,2,3},{4,5,6},{7,8,9} };
        Matrix2D<int> zi{ {9,10},{11,12} };
        print(z.sub_matrix(0, 1, 2, 2));     // sub matrix starting at row 0, col 1 with 2 rows 2 cols
        auto zz = z.insert_matrix(zi, 0, 1); // insert at row 0, col 1
        print(zz);
    }
    {
        // Check matrix of matrix operations
        Matrix2D<Matrix2D<int>> zz = Matrix2D<Matrix2D<int>>{ { Matrix2D<int>{ {1}}, Matrix2D<int>{ {3}} } };
        print(zz[0][0] + zz[0][1]); // add the two internal matrixes
    }
    {
        // Check matrix add and multiply
        Matrix2D<int> z1 = Matrix2D<int>{ {1,2},{3,4},{5,6},{7,8} };
        Matrix2D<int> z2{ {1,2,3},{4,5,6} };
        print(z1);
        print(z2);
        print(z2 + z2 + z2);
        Matrix2D<int> z3 = z1 * z2;     // z3: 4 rows, 3 cols
        print(z3);

        // product from separate math program product of z1 and z2
        Matrix2D<int> ref{ {9,12,15},{19,26,33},{29,40,51},{39,54,69} };
        print(z3 - ref);    // s/b all zeros
    }
}

Answer 1

Performance

Some operations look very inefficient. Consider subtracting two matrices, this involves:

1. Copying rhs into a new Matrix2D
2. Negating each element of that matrix
3. Copying lhs to a new Matrix2D
4. Adding the elements of the new matrix from step 1 to those of the matrix from step 3
5. Deleting the matrix from step 1

Then there is matrix multiplication, where you do a straightforward multiplication for matrices with less than 200 elements, but in the name of cache locality you first transpose the rhs. Sure that will make the multiplication more efficient, but the transposing operation itself might not be cache friendly (consider that memory stores might end up in the cache as well), and then you just wasted memory bandwidth making and reading mat_tmp.

Furthermore, where does the number 200 come from? Consider that there is a huge variation in the design of CPU caches, and that the size of T can vary a lot (consider an int16_t versus a std::complex<double>).

If you really want to optimize cache locality of your matrix multiplication function, consider storing the matrix using a Z-order curve.

For large matrices, the fact that you do not implement operator+=() and related functions might mean you are leaving some performance on the table.

Answer 2

G. Sliepen already pointed out a bunch of inefficiencies. I won't repeat what's already been said.

The inner loop of your matrix multiplication,

ret[row][col] = T{};
for (size_t row_col = 0; row_col < lhs.cols; row_col++)
   ret[row][col] += lhs[row][row_col] * rhs[row_col][col];

seems rather inefficient. (Note that lhs[row] calls a function that computes the pointer to the first element of the given row.) I hope an optimizing compiler can fix a lot of the redundant indexing you do here, but I would make those optimizations explicit. For example.

T res{};
T* lhs_it = lhs[row];
T* lhs_end = lhs_it + lhs.cols;
T* rhs_it = res[0] + col;
for (; lhs_it < lhs_end; ++lhs_it, rhs_it += rhs.cols) {
   res += *lhs_it * *rhs_it;
}
ret[row][col] = res;

---

Overloading operator~() as the transpose is unexpected. Yes, it is possible. No, it does not lead to clear code. Whenever you give a different meaning to an operator, you are making client code harder to understand. I would prefer to see mat.transposed() in code.

---

There's already been a discussion surrounding the use of std::construct_at() to change const variables. That discussion focuses on whether the compiler will put the variables in read-only memory or not. I think that is a red herring. What const says to the compiler is that the value will never change. The compiler can thus take shortcuts with this assumption. For example, if defining mat.cols = 3 at the top of a function, all throughout this function the compiler can replace the variable mat.cols with the constant 3. The program can avoid reading this memory location, since the value is known. Changing the value of mat.cols somewhere in the function hence leads to UB: the compiler can choose to ignore changes to the value!

Note that std::construct_at() is meant for use in compile-time memory allocation and construction (i.e. we can now use heap memory in a constexpr function).