I've expanded, and enhanced my matrix class which uses, aside from dynamic memory for the elements, only 3, const literals. These are rows, cols, and a ptr. The methods are all constexpr and the example code includes compile time verification of tests which compare compile time hashes against hashes generated at runtime. Since compile time execution is much better at detecting UB this seemed desirable as using consts in classes is rather new and not really viable pre c++20. This code no longer uses construct_at
on sub-objects but only on full objects which is more clearly legal in c++20.
this was the previous version I have incorporated many of the suggestions other commenters have made.
I also made a transpose method using a single loop with pointers expecting it to be faster than the version looping over both rows and cols. To my surprise it wasn't and came in 25% slower! That code is in the transpose2 method.
matrix2d.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;
using type = T;
constexpr Matrix2D() noexcept : pv(nullptr), cols(0), rows(0) {}
// uninitialized matrix. All values must be set priot to reading/using any
explicit constexpr Matrix2D(size_t a_rows, size_t a_cols) :
pv(new T[a_rows * a_cols]), cols(a_cols), rows(a_rows) {}
// value initialized matrix
explicit constexpr Matrix2D(size_t a_rows, size_t a_cols, T init) :
pv(new T[a_rows * a_cols]), cols(a_cols), rows(a_rows) {
std::fill(pv, pv + rows * cols, init);
}
// object only raw ctor for use with std::construct_at
explicit constexpr Matrix2D(T* pv, size_t cols, size_t rows) : pv(pv), cols(cols), rows(rows) {}
explicit constexpr 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);
}
constexpr Matrix2D& operator=(const Matrix2D& rhs) {
if (this != &rhs)
{
if (cols * rows != rhs.cols * rhs.rows) // skip delete/new if mem size unchanged
{
T* newp = new T[rhs.rows * rhs.cols]; // if new fails, *this is unchanged
delete[] pv;
std::construct_at(this, newp, rhs.cols, rhs.rows);
}
std::copy(rhs.pv, rhs.pv + rows * cols, pv);
}
return *this;
}
constexpr Matrix2D(Matrix2D&& rhs) noexcept : pv(rhs.pv), cols(rhs.cols), rows(rhs.rows) {
std::construct_at(&rhs, nullptr, 0, 0); // set rhs to a valid, null matrix
}
constexpr Matrix2D& operator=(Matrix2D&& rhs) noexcept {
delete[] pv;
std::construct_at(this, rhs.pv, rhs.cols, rhs.rows); // grab rhs's guts
std::construct_at(&rhs, nullptr, 0, 0); // set rhs to a valid, null matrix
return *this;
}
constexpr ~Matrix2D() {
delete[] pv;
}
// additional ctor to provide list initialization
constexpr 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(this, nullptr, 0, 0);
}
// bracket access v[r][c]
constexpr T* operator[](size_t row) noexcept { return pv + row * cols; }
constexpr const T* operator[](size_t row) const noexcept { return pv + row * cols; }
// paren access v(r,c);
constexpr T& operator()(size_t row, size_t col) noexcept { return pv[row * cols + col]; }
constexpr 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
constexpr Matrix2D insert_matrix(const Matrix2D& mat, size_t row_start, size_t col_start) const
{
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
constexpr Matrix2D get_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
constexpr Matrix2D operator-() const
{
Matrix2D ret(rows, cols);
for (size_t i = 0; i < rows * cols; i++)
ret.pv[i] = -pv[i];
return ret;
}
// add
constexpr 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 +=
constexpr 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
constexpr 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 -=
constexpr 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;
}
// This optimizes into just one loop yet is slower than
// the simple version using ret(c, r) = rhs(r, c);
constexpr Matrix2D transpose2() const
{
Matrix2D ret(cols, rows); // shape is flipped
T* p = pv;
T* pend = pv + rows * cols;
T* tp = ret.pv;
T* tpstart = tp;
T* tpend = ret.pv + rows * cols;
while (p < pend)
{
*tp = *p++;
tp += rows;
if (tp >= tpend)
tp = ++tpstart;
}
return ret;
}
// transpose
constexpr Matrix2D transpose() const
{
Matrix2D<T> ret(cols, rows, 0);
for (size_t r = 0; r < rows; r++)
for (size_t c = 0; c < cols; c++)
{
auto r1 = r;
auto c1 = c;
ret(c, r) = operator()(r, c);
}
return ret;
}
// append new columns
constexpr Matrix2D append_columns(const Matrix2D& new_cols) const
{
if (rows != new_cols.rows)
throw std::range_error("appended matrix must have same number of rows");
Matrix2D<T> ret(rows, cols + new_cols.cols);
for (size_t r = 0; r < ret.rows; r++)
for (size_t c = 0; c < ret.cols; c++)
if (c < cols)
ret(r, c) = operator()(r, c);
else
ret(r, c) = new_cols(r, c - cols);
return ret;
}
// append new rows
constexpr Matrix2D append_rows(const Matrix2D& new_rows) const
{
if (cols != new_rows.cols)
throw std::range_error("appended matrix must have same number of columns");
Matrix2D<T> ret(rows + new_rows.rows, cols);
for (size_t r = 0; r < ret.rows; r++)
for (size_t c = 0; c < ret.cols; c++)
if (r < rows)
ret(r, c) = operator()(r, c);
else
ret(r, c) = new_rows(r - rows, c);
return ret;
}
constexpr Matrix2D invert() const
{
if (rows != cols)
throw std::range_error("cols and rows must be the same");
auto ident = Matrix2D(rows, rows, 0);
for (size_t i = 0; i < rows; i++)
ident(i, i) = 1;
Matrix2D m = this->append_columns(ident);
for (size_t i = m.rows - 1; i > 0; i--)
//Swap rows to put largest initial element first
if (m[i - 1][0] < m[i][0])
for (size_t j = 0; j < m.cols; j++)
std::swap(m[i][j], m[i - 1][j]);
for (size_t i = 0; i < m.rows; i++)
{
for (size_t j = 0; j < m.rows; j++)
{
if (j != i)
{
if (m[i][i] == 0)
throw std::domain_error("singular");
T temp = m[j][i] / m[i][i];
for (size_t k = 0; k < 2 * m.rows; k++)
m[j][k] -= m[i][k] * temp;
}
}
}
for (size_t i = 0; i < m.rows; i++)
{
T temp = m[i][i];
for (size_t j = 0; j < 2 * m.rows; j++)
m[i][j] = m[i][j] / temp;
}
return m.get_sub_matrix(0, m.rows, m.rows, m.rows);
}
};
// multiply matrixes
template <typename T>
constexpr 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(); // 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;
}
// equality test
template <typename T>
constexpr bool operator==(const Matrix2D<T>& lhs, const Matrix2D<T>& rhs)
{
if (lhs.cols != rhs.cols || lhs.rows != rhs.rows)
return false;
for (size_t row = 0; row < lhs.rows; row++)
for (size_t col = 0; col < lhs.cols; col++)
if (lhs(row,col) != rhs(row, col))
return false;
return true;
}
template <typename T>
constexpr bool operator!=(const Matrix2D<T>& lhs, const Matrix2D<T>& rhs)
{
return !(lhs==rhs);
}
}
main.cpp
#include <iostream>
#include <array>
#include <string>
#include <vector>
#include "matrix2d.h"
using namespace matrix_const;
struct HashList {
unsigned int hash[100];
};
// set #if 0 to force printout of all matrixes
#if 1
constexpr HashList hashlist{
4231457547u,559445916u,3889815258u,559445916u,559445916u,559445916u,1196435762u,146730798u,2933268210u,2648446219u,
1101633479u,2648446219u,1934638098u,2735629034u,3480429190u,2106445168u,559445916u,559445916u,3480429190u,2106445168u,
559445916u,2648446219u,1101633479u,2648446219u,3682454706u,3247604110u,191436758u,2648446219u,1101633479u,2648446219u,
2526570054u, };
#else
constexpr HashList hashlist{};
#endif
template <typename T>
void print(const Matrix2D<T>& mat, const std::string& message = "")
{
std::cout << message << "\n";
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 error_count{};
template <typename T>
constexpr unsigned int matrix_hash(const Matrix2D<T>& mat, size_t &index)
{
auto crc_cycle = [](unsigned int v, auto x) {return (std::bit_cast<unsigned int>(x) + v) * 1664525u + 1013904223u; };
unsigned int hash{};
hash = crc_cycle(hash, static_cast<unsigned int>(mat.rows));
hash = crc_cycle(hash, static_cast<unsigned int>(mat.cols));
for (size_t r = 0; r < mat.rows; r++)
for (size_t c = 0; c < mat.cols; c++)
hash = crc_cycle(hash, mat(r, c));
if (!std::is_constant_evaluated() && hash != hashlist.hash[index])
{
print(mat, std::to_string(index) + " Matrix Difference Detected");
error_count++;
}
index++;
return hash;
}
// Generate hashes of matrixes after various tests at both runtime and compile time
// so they can be checked against each other
constexpr std::array<unsigned int, 100> do_tests()
{
std::vector<unsigned int> hashes;
std::array<unsigned int, 100> hashes_array{};
size_t index{ 0 };
using Mati = Matrix2D<int>;
// index:0 Test initialize with fill value
{
Mati x(1, 2, 10);
hashes.push_back(matrix_hash(x, index));
}
// index:1-3 Test list ctor, access with paren and [], move assignment
{
Mati x{ {1, 2}, { 3, 4 } };
hashes.push_back(matrix_hash(x, index));
x[0][1]++;
x(1, 0)++;
hashes.push_back(matrix_hash(x, index));
x = Mati({ {1, 2}, { 3, 4 } });
hashes.push_back(matrix_hash(x, index));
}
// index:4-6 Test of Move assigment operator and verification guts moved
{
Mati x{ {1, 2}, { 3, 4 } };
Mati y;
y = x;
hashes.push_back(matrix_hash(y, index));
y = std::move(x);
hashes.push_back(matrix_hash(y, index));
hashes.push_back(matrix_hash(x, index)); // compiler warns but this is to check gutted object remains valid
}
// index 7:8 Test insert and extract submatrixes
{
Mati x{ {1,2,3},{4,5,6},{7,8,9} };
Mati y{ {9,10},{11,12} };
auto z = x.get_sub_matrix(0, 1, 2, 2);
hashes.push_back(matrix_hash(z, index));
auto zz = x.insert_matrix(y, 0, 1); // insert at row 0, col 1
hashes.push_back(matrix_hash(zz, index));
}
// Index 9:11 Test transpose
{
Mati x = Mati{ { 2, 3, 4 }, { 5, 6, 7 }, {8, 9, 10},{11,12,13} };
hashes.push_back(matrix_hash(x, index));
auto x1 = x.transpose();
hashes.push_back(matrix_hash(x1, index));
auto x2 = x1.transpose();
hashes.push_back(matrix_hash(x2, index));
}
// Index 12:13 Test append cols, append rows
{
Mati x{ {1, 2}, { 3, 4 } };
Mati y{ {5}, {6} };
Mati z;
z = x.append_columns(y);
hashes.push_back(matrix_hash(z, index));
z = x.append_rows(y.transpose());
hashes.push_back(matrix_hash(z, index));
}
// Index 14:17 Test +-assignments
{
Mati x{ {1, 2}, { 3, 4 } }, y1, y2;
y1 = x; y2 = x;
y2 += y1 += x; // note prec is left to right
hashes.push_back(matrix_hash(y1, index));
hashes.push_back(matrix_hash(y2, index));
y2 -= y1 -= x; // note prec is left to right
y2 -= x; // y2 == x
hashes.push_back(matrix_hash(y1, index));
hashes.push_back(matrix_hash(y2, index));
}
{
// Index 18:20 Test + and - operators
Mati x{ {1, 2}, { 3, 4 } }, y1, y2;
y1 = x + x;
hashes.push_back(matrix_hash(y1, index));
y1 = y1 + x;
hashes.push_back(matrix_hash(y1, index));
y1 = y1 - x - x; // y1 == x
hashes.push_back(matrix_hash(y1, index));
}
{
// Index 21:23 Test transpose
Mati x{ { 2, 3, 4 }, { 5, 6, 7 }, {8, 9, 10},{11,12,13} };
hashes.push_back(matrix_hash(x, index));
auto x1 = x.transpose();
hashes.push_back(matrix_hash(x1, index));
auto x2 = x1.transpose();
hashes.push_back(matrix_hash(x2, index));
}
{
// Index 24:26, float matrix inversion test
// inversion produces small errors so any difference between compile time
// and runtime float calcs will print associated tests
// This may occur with different fp settings
Matrix2D<float> x{ {1, 2}, { 3, 4 } };
hashes.push_back(matrix_hash(x, index));
auto x1 = x.invert();
hashes.push_back(matrix_hash(x1, index));
auto x2 = x1.invert();
hashes.push_back(matrix_hash(x2, index));
}
{
// Index 27:30 Test Transpose and multiply
Mati x{ {1} };
x = Mati{ { 2, 3, 4 }, { 5, 6, 7 }, {8, 9, 10},{11,12,13} };
hashes.push_back(matrix_hash(x, index));
auto x1 = x.transpose();
hashes.push_back(matrix_hash(x1, index));
auto x2 = x1.transpose();
hashes.push_back(matrix_hash(x2, index));
auto x3 = x * x1;
hashes.push_back(matrix_hash(x3, index));
}
// only print hashes if called at runtime and hashlist not initialized with hashes
if (!std::is_constant_evaluated() && hashlist.hash[0] == 0)
{
int count{};
for (auto x : hashes)
{
std::cout << x << "u,";
if (count++ % 10 == 9)
std::cout << "\n";
}
}
for (size_t i = 0; i < hashes.size(); i++)
hashes_array[i] = hashes[i];
return hashes_array;
}
int main()
{
// check compiler requirements for this code
static_assert(sizeof(unsigned int) == 4 && sizeof(float) == 4);
using namespace matrix_const;
{
constexpr std::array<unsigned int, 100> hash_static = do_tests();
std::array<unsigned int, 100> hash_dynamic = do_tests();
if (hash_dynamic != hash_static)
throw "dynamic/static hashes don't match";
if (hashlist.hash[0] != 0)
{
if (error_count)
std::cout << error_count << " hash error(s) detected\n";
else
std::cout << "passed\n";
}
}
}