Return to Question

#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);
            for (size_t i = 0; i < ret.rows * ret.cols; i++)
                ret.pv[i] = -ret.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;
        }
    };

    // 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);
        if (lhs.rows * rhs.cols < 200)
        {
            for (size_t row = 0; row < lhs.rows; row++)
                for (size_t col = 0; col < rhs.cols; col++)
                {
                    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];
                }
        }
        else
        {
            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;
    }


    // subtract
    template <typename T>
    Matrix2D<T> operator-(const Matrix2D<T>& lhs, const Matrix2D<T>& rhs)
    {
        if (lhs.cols != rhs.cols && lhs.rows != rhs.rows)
            throw std::range_error("matrixes must have same dims");
        return lhs + (-rhs);
    }

    // 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;
    }
}

#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;
    }
}

Edit: I've added code to address G. Sliepen's excellent performance comments. I also investigated performance to determine the crossover where first transposing results in slower performance. This was done by creating 100,000+ random matrixes then walking through pairs of them and accumulating execution time. Turns out the point where the transpose starts to outperform the standard multiply algo is where the matrix size is only 4x4. At 8x8 the transpose is over twice as fast. This makes sense given the transpose is O(n^2) while the multiply is O(n^3). In both cases the largest performance hit is from the new dynamic allocation. Given this, I've removed the non-transpose portion of the multiply.

CompilerExplorer, gcc, clang, msvc

#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;
    }
}

CompilerExplorer, gcc, clang, msvc

#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);
            for (size_t i = 0; i < ret.rows * ret.cols; i++)
                ret.pv[i] = -ret.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;
        }
    };

    // 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);
        if (lhs.rows * rhs.cols < 200)
        {
            for (size_t row = 0; row < lhs.rows; row++)
                for (size_t col = 0; col < rhs.cols; col++)
                {
                    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];
                }
        }
        else
        {
            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;
    }


    // subtract
    template <typename T>
    Matrix2D<T> operator-(const Matrix2D<T>& lhs, const Matrix2D<T>& rhs)
    {
        if (lhs.cols != rhs.cols && lhs.rows != rhs.rows)
            throw std::range_error("matrixes must have same dims");
        return lhs + (-rhs);
    }

    // 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;
    }
}

CompilerExplorer, gcc, clang, msvc

#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);
            for (size_t i = 0; i < ret.rows * ret.cols; i++)
                ret.pv[i] = -ret.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;
        }
    };

    // 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);
        if (lhs.rows * rhs.cols < 200)
        {
            for (size_t row = 0; row < lhs.rows; row++)
                for (size_t col = 0; col < rhs.cols; col++)
                {
                    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];
                }
        }
        else
        {
            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;
    }


    // subtract
    template <typename T>
    Matrix2D<T> operator-(const Matrix2D<T>& lhs, const Matrix2D<T>& rhs)
    {
        if (lhs.cols != rhs.cols && lhs.rows != rhs.rows)
            throw std::range_error("matrixes must have same dims");
        return lhs + (-rhs);
    }

    // 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;
    }
}

Edit: I've added code to address G. Sliepen's excellent performance comments. I also investigated performance to determine the crossover where first transposing results in slower performance. This was done by creating 100,000+ random matrixes then walking through pairs of them and accumulating execution time. Turns out the point where the transpose starts to outperform the standard multiply algo is where the matrix size is only 4x4. At 8x8 the transpose is over twice as fast. This makes sense given the transpose is O(n^2) while the multiply is O(n^3). In both cases the largest performance hit is from the new dynamic allocation. Given this, I've removed the non-transpose portion of the multiply.

CompilerExplorer, gcc, clang, msvc

#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;
    }
}