C++20 : N-dimensional minimal Matrix class

Question

I wrote a N-dim matrix (tensor) class, based on my previous 2D matrix implementation(2D Matrix in C++20 and Strassen's algorithm), accepting many helpful reviews from here.

MatrixBase.h

#ifndef FROZENCA_MATRIXBASE_H
#define FROZENCA_MATRIXBASE_H

#include <algorithm>
#include <array>
#include <cassert>
#include <cstddef>
#include <concepts>
#include <functional>
#include <initializer_list>
#include <iostream>
#include <numeric>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include "MatrixUtils.h"
#include "MatrixInitializer.h"

namespace frozenca {

template <typename T, std::regular R, std::size_t N>
class MatrixBase {
public:
    static constexpr std::size_t ndim = N;

private:
    std::size_t size_;
    std::array<std::size_t, N> dims_;
    std::array<std::size_t, N> strides_;

protected:
    MatrixBase() = delete;
    virtual ~MatrixBase() = default;

    template <typename... Dims>
    explicit MatrixBase(Dims... dims);

    MatrixBase(std::size_t size,
               const std::array<std::size_t, N> dims,
               const std::array<std::size_t, N> strides) : size_ {size}, dims_ {dims}, strides_ {strides} {}

public:
    using value_type = R;
    using reference = R&;
    using const_reference = const R&;
    using pointer = R*;

    friend void swap(MatrixBase& a, MatrixBase& b) noexcept {
        std::swap(a.size_, b.size_);
        std::swap(a.dims_, b.dims_);
        std::swap(a.strides_, b.strides_);
    }

    auto begin() { return static_cast<T&>(*this).begin(); }
    auto cbegin() const { return static_cast<const T&>(*this).cbegin(); }
    auto end() { return static_cast<T&>(*this).end(); }
    auto cend() const { return static_cast<const T&>(*this).cend(); }
    auto rbegin() { return static_cast<T&>(*this).rbegin(); }
    auto crbegin() const { return static_cast<const T&>(*this).crbegin(); }
    auto rend() { return static_cast<T&>(*this).rend(); }
    auto crend() const { return static_cast<const T&>(*this).crend(); }

    template <RequestingElement... Args>
    reference operator()(Args... args);

    template <RequestingElement... Args>
    const_reference operator()(Args... args) const;

    reference operator[](const std::array<std::size_t, N>& pos);
    const_reference operator[](const std::array<std::size_t, N>& pos) const;

    template <std::regular U>
    MatrixBase(std::initializer_list<U>) = delete;

    template <std::regular U>
    MatrixBase& operator=(std::initializer_list<U>) = delete;

    MatrixBase(typename MatrixInitializer<R, N>::type init);

    [[nodiscard]] std::size_t size() const { return size_;}

    [[nodiscard]] const std::array<std::size_t, N>& dims() const {
        return dims_;
    }

    [[nodiscard]] std::size_t dims(std::size_t n) const {
        if (n >= N) {
            throw std::out_of_range("Out of range in dims");
        }
        return dims_[n];
    }

    [[nodiscard]] const std::array<std::size_t, N>& strides() const {
        return strides_;
    }

    [[nodiscard]] std::size_t strides(std::size_t n) const {
        if (n >= N) {
            throw std::out_of_range("Out of range in strides");
        }
        return strides_[n];
    }

};

template <typename T, std::regular R, std::size_t N>
template <typename... Dims>
MatrixBase<T, R, N>::MatrixBase(Dims... dims) : size_{(... * static_cast<std::size_t>(dims))},
    dims_{static_cast<std::size_t>(dims)...} {
    static_assert(sizeof...(Dims) == N);
    if (std::ranges::find(dims_, 0lu) != std::end(dims_)) {
        throw std::invalid_argument("Zero dimension not allowed");
    }
    strides_ = computeStrides(dims_);
}

template <typename T, std::regular R, std::size_t N>
MatrixBase<T, R, N>::MatrixBase(typename MatrixInitializer<R, N>::type init) {
    dims_ = deriveDims<N>(init);
    strides_ = computeStrides(dims_);
    size_ = strides_[0] * dims_[0];
}

template <typename T, std::regular R, std::size_t N>
template <RequestingElement... Args>
typename MatrixBase<T, R, N>::reference MatrixBase<T, R, N>::operator()(Args... args) {
    return const_cast<typename MatrixBase<T, R, N>::reference>(std::as_const(*this).operator()(args...));
}

template <typename T, std::regular R, std::size_t N>
template <RequestingElement... Args>
typename MatrixBase<T, R, N>::const_reference MatrixBase<T, R, N>::operator()(Args... args) const {
    static_assert(sizeof...(args) == N);
    std::array<std::size_t, N> pos {std::size_t(args)...};
    return this->operator[](pos);
}

template <typename T, std::regular R, std::size_t N>
typename MatrixBase<T, R, N>::reference MatrixBase<T, R, N>::operator[](const std::array<std::size_t, N>& pos) {
    return const_cast<typename MatrixBase<T, R, N>::reference>(std::as_const(*this).operator[](pos));
}

template <typename T, std::regular R, std::size_t N>
typename MatrixBase<T, R, N>::const_reference MatrixBase<T, R, N>::operator[](const std::array<std::size_t, N>& pos) const {
    if (!std::equal(std::cbegin(pos), std::cend(pos), std::cbegin(dims_), std::less<>{})) {
        throw std::out_of_range("Out of range in element access");
    }
    return *(cbegin() + std::inner_product(std::cbegin(pos), std::cend(pos), std::cbegin(strides_), 0lu));
}

} // namespace frozenca

#endif //FROZENCA_MATRIXBASE_H

Matrix.h

#ifndef FROZENCA_MATRIX_H
#define FROZENCA_MATRIX_H

#include <array>
#include <iterator>
#include <memory>
#include "MatrixBase.h"
#include "MatrixView.h"
#include "MatrixInitializer.h"

namespace frozenca {

template <std::regular T, std::size_t N>
class Matrix final : public MatrixBase<Matrix<T, N>, T, N> {
private:
    std::unique_ptr<T[]> data_;

public:
    using Base = MatrixBase<Matrix<T, N>, T, N>;
    using Base::size;
    using Base::dims;
    using Base::strides;

    template <typename... Dims>
    explicit Matrix(Dims... dims);

    ~Matrix() override = default;

    friend void swap(Matrix& a, Matrix& b) noexcept {
        swap(static_cast<Base&>(a), static_cast<Base&>(b));
        std::swap(a.data_, b.data_);
    }

    using iterator = T*;
    using const_iterator = const T*;
    using reverse_iterator = std::reverse_iterator<iterator>;
    using const_reverse_iterator = std::reverse_iterator<const_iterator>;

    iterator begin() { return &data_[0];}
    const_iterator cbegin() const { return &data_[0]; }
    iterator end() { return &data_[size()];}
    const_iterator cend() const { return &data_[size()];}
    reverse_iterator rbegin() { return std::make_reverse_iterator(end());}
    const_reverse_iterator crbegin() const { return std::make_reverse_iterator(cend());}
    reverse_iterator rend() { return std::make_reverse_iterator(begin());}
    const_reverse_iterator crend() const { return std::make_reverse_iterator(cbegin());}

    template <typename M, std::regular U>
    Matrix(const MatrixBase<M, U, N>&);

    template <typename M, std::regular U>
    Matrix& operator=(const MatrixBase<M, U, N>&);

    Matrix(typename MatrixInitializer<T, N>::type init);
    Matrix& operator=(typename MatrixInitializer<T, N>::type init);

    MatrixView<T, N> submatrix(const std::array<std::size_t, N>& pos_begin, const std::array<std::size_t, N>& pos_end = dims()) {
        if (!std::equal(std::cbegin(pos_begin), std::cend(pos_begin), std::cbegin(pos_end), std::less<>{})) {
            throw std::out_of_range("submatrix begin/end position error");
        }
        std::array<std::size_t, N> view_dims;
        std::transform(std::cbegin(pos_end), std::cend(pos_end), std::cbegin(pos_begin), std::begin(view_dims), std::minus<>{});
        std::size_t view_size = std::accumulate(std::cbegin(view_dims), std::cend(view_dims), 1lu, std::multiplies<>{});
        MatrixView<T, N> view (view_size, view_dims, strides(), &this->operator[](pos_begin));
        return view;
    }
};

template <std::regular T, std::size_t N>
template <typename... Dims>
Matrix<T, N>::Matrix(Dims... dims) : Base(dims...), data_(std::make_unique<T[]>(size())) {
}

template <std::regular T, std::size_t N>
Matrix<T, N>::Matrix(typename MatrixInitializer<T, N>::type init) : Base(init),
        data_(std::make_unique<T[]>(size())) {
    insertFlat(data_, init);
}

template <std::regular T, std::size_t N>
Matrix<T, N>& Matrix<T, N>::operator=(typename MatrixInitializer<T, N>::type init) {
    Matrix<T, N> mat(init);
    swap(*this, mat);
    return *this;
}

} // namespace frozenca

#endif //FROZENCA_MATRIX_H

MatrixView.h

#ifndef FROZENCA_MATRIXVIEW_H
#define FROZENCA_MATRIXVIEW_H

#include <compare>
#include "MatrixBase.h"

namespace frozenca {

template <std::regular T, std::size_t N>
class MatrixView final : public MatrixBase<MatrixView<T, N>, T, N> {
private:
    T* data_view_;
    std::array<std::size_t, N> view_strides_;

public:
    using Base = MatrixBase<MatrixView<T, N>, T, N>;
    using Base::size;
    using Base::dims;
    using Base::strides;

    explicit MatrixView(std::size_t size,
                        const std::array<std::size_t, N>& dims,
                        const std::array<std::size_t, N>& strides,
                        T* data_view) : Base(size, dims, strides), data_view_ {data_view} {
        view_strides_ = computeStrides(dims);
    }

    ~MatrixView() override = default;

    friend void swap(MatrixView& a, MatrixView& b) noexcept {
        swap(static_cast<Base&>(a), static_cast<Base&>(b));
        std::swap(a.data_view_, b.data_view_);
        std::swap(a.view_strides_, b.view_strides_);
    }

    template <typename T_>
    struct MVIterator {
        MatrixView* ptr_ = nullptr;
        std::array<std::size_t, N> pos_ = {0};
        std::size_t offset_ = 0;
        std::size_t index_ = 0;

        using difference_type = std::ptrdiff_t;
        using value_type = T_;
        using pointer = T_*;
        using reference = T_&;
        using iterator_category = std::random_access_iterator_tag;

        MVIterator(MatrixView* ptr, std::array<std::size_t, N> pos = {0}) : ptr_ {ptr}, pos_ {pos} {
            ValidateOffset();
        }

        reference operator*() const {
            return ptr_->data_view_[offset_];
        }

        pointer operator->() const {
            return ptr_->data_view_ + offset_;
        }

        void ValidateOffset() {
            offset_ = std::inner_product(std::cbegin(pos_), std::cend(pos_), std::cbegin(ptr_->strides()), 0lu);
            index_ = std::inner_product(std::cbegin(pos_), std::cend(pos_), std::cbegin(ptr_->view_strides_), 0lu);
            if (index_ > ptr_->size()) {
                throw std::out_of_range("MatrixView iterator out of range");
            }
        }

        void Increment() {
            for (std::size_t i = N - 1; i < N; --i) {
                ++pos_[i];
                if (pos_[i] != ptr_->dims(i) || i == 0) {
                    break;
                } else {
                    pos_[i] = 0;
                }
            }
            ValidateOffset();
        }

        void Increment(std::ptrdiff_t n) {
            if (n < 0) {
                Decrement(-n);
                return;
            }
            auto carry = static_cast<std::size_t>(n);
            for (std::size_t i = N - 1; i < N; --i) {
                std::size_t curr_dim = ptr_->dims(i);
                pos_[i] += carry;
                if (pos_[i] < curr_dim || i == 0) {
                    break;
                } else {
                    carry = pos_[i] / curr_dim;
                    pos_[i] %= curr_dim;
                }
            }
            ValidateOffset();
        }

        void Decrement() {
            for (std::size_t i = N - 1; i < N; --i) {
                --pos_[i];
                if (pos_[i] != static_cast<std::size_t>(-1) || i == 0) {
                    break;
                } else {
                    pos_[i] = ptr_->dims(i) - 1;
                }
            }
            ValidateOffset();
        }

        void Decrement(std::ptrdiff_t n) {
            if (n < 0) {
                Increment(-n);
                return;
            }
            auto carry = static_cast<std::size_t>(n);
            for (std::size_t i = N - 1; i < N; --i) {
                std::size_t curr_dim = ptr_->dims(i);
                pos_[i] -= carry;
                if (pos_[i] < curr_dim || i == 0) {
                    break;
                } else {
                    carry = static_cast<std::size_t>(-quot(static_cast<long>(pos_[i]), static_cast<long>(curr_dim)));
                    pos_[i] = mod(static_cast<long>(pos_[i]), static_cast<long>(curr_dim));
                }
            }
            ValidateOffset();
        }

        MVIterator& operator++() {
            Increment();
            return *this;
        }

        MVIterator operator++(int) {
            MVIterator temp = *this;
            Increment();
            return temp;
        }

        MVIterator& operator--() {
            Decrement();
            return *this;
        }

        MVIterator operator--(int) {
            MVIterator temp = *this;
            Decrement();
            return temp;
        }

        MVIterator operator+(difference_type n) const {
            MVIterator temp = *this;
            temp.Increment(n);
            return temp;
        }

        MVIterator& operator+=(difference_type n) {
            Increment(n);
            return *this;
        }

        MVIterator operator-(difference_type n) const {
            MVIterator temp = *this;
            temp.Decrement(n);
            return temp;
        }

        MVIterator& operator-=(difference_type n) {
            Decrement(n);
            return *this;
        }

        reference operator[](difference_type n) const {
            return *(this + n);
        }

        template <typename T2>
        std::enable_if_t<std::is_same_v<std::remove_cv_t<T_>, std::remove_cv_t<T2>>, difference_type>
                operator-(const MVIterator<T2>& other) const {
            return offset_ - other.offset_;
        }

        // oh no.. *why* defining operator<=> doesn't work to automatically define these in gcc?
        template <typename T1, typename T2>
        friend std::enable_if_t<std::is_same_v<std::remove_cv_t<T1>, std::remove_cv_t<T2>>,
                bool>
        operator==(const MVIterator<T1>& it1, const MVIterator<T2>& it2) {
            return it1.offset_ == it2.offset_;
        }

        template <typename T1, typename T2>
        friend std::enable_if_t<std::is_same_v<std::remove_cv_t<T1>, std::remove_cv_t<T2>>,
                bool>
        operator!=(const MVIterator<T1>& it1, const MVIterator<T2>& it2) {
            return !(it1 == it2);
        }

        template <typename T1, typename T2>
        friend std::enable_if_t<std::is_same_v<std::remove_cv_t<T1>, std::remove_cv_t<T2>>,
        std::compare_three_way_result_t<std::size_t, std::size_t>>
        operator<=>(const MVIterator<T1>& it1, const MVIterator<T2>& it2) {
            return it1.offset_ <=> it2.offset_;
        }
    };

    using iterator = MVIterator<T>;
    using const_iterator = MVIterator<const T>;
    using reverse_iterator = std::reverse_iterator<iterator>;
    using const_reverse_iterator = std::reverse_iterator<const_iterator>;

    iterator begin() { return iterator(this);}
    const_iterator cbegin() const { return const_iterator(this); }
    iterator end() { return iterator(this, {dims(0), });}
    const_iterator cend() const { return iterator(this, {dims(0), });}
    reverse_iterator rbegin() { return std::make_reverse_iterator(end());}
    const_reverse_iterator crbegin() const { return std::make_reverse_iterator(cend());}
    reverse_iterator rend() { return std::make_reverse_iterator(begin());}
    const_reverse_iterator crend() const { return std::make_reverse_iterator(cbegin());}

};

} // namespace frozenca

#endif //FROZENCA_MATRIXVIEW_H

MatrixInitializer.h

#ifndef FROZENCA_MATRIXINITIALIZER_H
#define FROZENCA_MATRIXINITIALIZER_H

#include <concepts>
#include <initializer_list>

namespace frozenca {

template <std::regular T, std::size_t N>
struct MatrixInitializer {
    using type = std::initializer_list<typename MatrixInitializer<T, N - 1>::type>;
};

template <std::regular T>
struct MatrixInitializer<T, 1> {
    using type = std::initializer_list<T>;
};

template <std::regular T>
struct MatrixInitializer<T, 0>;

} // namespace frozenca

#endif //FROZENCA_MATRIXINITIALIZER_H

MatrixUtils.h

#ifndef FROZENCA_MATRIXUTILS_H
#define FROZENCA_MATRIXUTILS_H

#include <array>
#include <cassert>
#include <cstddef>
#include <iostream>
#include <iterator>
#include <memory>
#include <type_traits>

namespace frozenca {

template <typename... Args>
inline constexpr bool All(Args... args) { return (... && args); };

template <typename... Args>
inline constexpr bool Some(Args... args) { return (... || args); };

template <typename... Args>
concept RequestingElement = All(std::is_convertible_v<Args, std::size_t>...);

template <std::size_t N>
std::array<std::size_t, N> computeStrides(const std::array<std::size_t, N>& dims) {
    std::array<std::size_t, N> strides;
    std::size_t str = 1;
    for (std::size_t i = N - 1; i < N; --i) {
        strides[i] = str;
        str *= dims[i];
    }
    return strides;
}

template <std::size_t N, typename Initializer>
bool checkNonJagged(const Initializer& init) {
    auto i = std::cbegin(init);
    for (auto j = std::next(i); j != std::cend(init); ++j) {
        if (i->size() != j->size()) {
            return false;
        }
    }
    return true;
}

template <std::size_t N, typename Iter, typename Initializer>
void addDims(Iter& first, const Initializer& init) {
    if constexpr (N > 1) {
        assert(checkNonJagged<N>(init));
    }
    *first = std::size(init);
    ++first;
    if constexpr (N > 1) {
        addDims<N - 1>(first, *std::begin(init));
    }
}

template <std::size_t N, typename Initializer>
std::array<std::size_t, N> deriveDims(const Initializer& init) {
    std::array<std::size_t, N> dims;
    auto f = std::begin(dims);
    addDims<N>(f, init);
    return dims;
}

template <std::regular T>
void addList(std::unique_ptr<T[]>& data,
             const T* first, const T* last,
             std::size_t& index) {
    for (; first != last; ++first) {
        data[index] = *first;
        ++index;
    }
}

template <std::regular T, typename I>
void addList(std::unique_ptr<T[]>& data,
             const std::initializer_list<I>* first, const std::initializer_list<I>* last,
             std::size_t& index) {
    for (; first != last; ++first) {
        addList(data, first->begin(), first->end(), index);
    }
}

template <std::regular T, typename I>
void insertFlat(std::unique_ptr<T[]>& data, std::initializer_list<I> list) {
    std::size_t index = 0;
    addList(data, std::begin(list), std::end(list), index);
}

inline long quot(long a, long b) {
    return (a % b) >= 0 ? (a / b) : (a / b) - 1;
}

inline long mod(long a, long b) {
    return (a % b + b) % b;
}

} // namespace frozenca

#endif //FROZENCA_MATRIXUTILS_H

Test code:

#include "Matrix.h"
#include <iostream>
#include <numeric>

int main() {
    frozenca::Matrix<float, 2> m {{1, 2, 3}, {4, 5, 6}};
    m = {{4, 6}, {1, 3}};
    std::cout << m(1, 1) << '\n';
    std::cout << m[{1, 1}] << '\n';

    // 3 x 4 x 5
    frozenca::Matrix<int, 3> m2 (3, 4, 5);
    std::iota(std::begin(m2), std::end(m2), 0);

    // should output 0 ~ 59
    for (auto n : m2) {
        std::cout << n << ' ';
    }
    std::cout << '\n';

    auto sub = m2.submatrix({1, 1, 1}, {2, 3, 4});

    // 26 27 28
    // 31 32 33
    for (auto n : sub) {
        std::cout << n << ' ';
    }
    std::cout << '\n';

}

Works as intended.

Yet to come: Custom allocators, slices, numerical operations (ex. matrix multiplication), reshape(), 1-D template specialization, etc.

Feel free to comment anything as always!

Answer 1

This is looking very good. I especially like the improvement of submatrix(). Still, there are a few things I mentioned in my review of the 2D class that I'll repeat here, as well as add a few other remarks.

Avoid restricting template types too much

I see you are using std::regular now instead of concept Scalar. This still requires that the type is equality comparable, but nothing in your classes depend on that. You could perhaps use std::semiregular instead.

Avoid code duplication

Always look for opportunities to reduce code duplication, even if it looks like small stuff. In MatrixBase, begin() and related functions all do a static_cast<T&>(*this). You can make this code look nicer by writing:

template <typename T, std::semiregular R, std::size_t N>
class MatrixBase {
    T &self() { return static_cast<T&>(*this); }
    ...
public:
    auto begin() { return self().begin(); }
    ...
}

quot() in MatrixUtils.h could be rewritten as:

inline long quot(long a, long b) {
    return (a / b) - (a % b < 0);
}

Removing duplication reduces the chance of bugs.

Issues with concept RequestingElement

This concept is used in variadic templates, but it should not itself be a variadic template. I would rewrite it as:

template <typename Arg>
concept RequestingElement = std::is_convertible_v<Arg, std::size_t>;

With this, the concepts All and Some are not necessary. However, there are more issues with this concept. For example, this just checks that something is implicitly convertible to size_t, but a float would also implicitly convert to a size_t. It would be better to restrict this to integral types.

Finally, the name is is not very clear to me, I think it would be better to rename this to IndexType.

Make size_, dims_ and strides_ const

It looks like the dimensions of a matrix are fixed at construction time. You should therefore make the member variables that represent the size and shape of the matrix const.

Note that you can, and will need to, call functions inside the member initializer list, like so for example:

template <typename T, std::regular R, std::size_t N>
template <typename... Dims>
MatrixBase<T, R, N>::MatrixBase(Dims... dims) :
    size_{(... * static_cast<std::size_t>(dims))},
    dims_{static_cast<std::size_t>(dims)...},
    strides_{computeStrides(dims_)},
{
    static_assert(sizeof...(Dims) == N);
}

Allow zero length dimensions

It looks like your code handles matrices with one of the dimensions having length zero just fine. I would remove the check for that in the constructor of MatrixBase.

Note that allocating an array of length zero using new, and by extension also when using std::make_unique<T[]>(size), is legal in C++.

Answer 2

Since C++11 you can count on / and % behaving in a prescribed manner: Rather than being implementation-defined (whatever the underlying machine instruction does), it is defined to truncate toward zero. So, you don't need your quot wrapper, and mod can likewise assume that division worked that way already.

And before C++11, there was std::div.

---

 auto f = std::begin(dims);

You need to invoke begin and a few others by using the "std two-step", not by calling the qualified version directly. If the type of dims is not a built-in, and requires argument-dependant lookup, then this will fail.

using std::begin;
auto f = begin(dims);

or, use the newer versions in the ranges namespace, or your own wrappers.

But you know the type is std::array, so just like you did when you had an std::initializer_list, just use the member function form dims.begin() so it doesn't raise red flags to the reader.

---

BTW, #pragma once; works on all the major compilers.