I have been lurking on the C++ side of stack-overflow only long enough to know there are a lot of beginners and intermediate programmers baffled by multi-dimensional arrays. I've seen a lot of monstrosities : three-star programmers, cache-unfriendly implementations, vectors of vectors of vectors ..., etc.

I thought to provide them with a basic multi-dimensional array:

  • open for extensions,
  • fit for production and experimentation,
  • working without allocations,
  • and cache-friendly.

Before I do so, I'd like some advice from the community.

  1. What can I improve in the provided features of this class?
  2. What new feature could one add, and why would it be useful?
#include <array>
#include <numeric>

namespace ysc
namespace _details
    template<class InputIt, class OutputIt>
    OutputIt partial_product(InputIt first, InputIt last, OutputIt output)
    { *output++ = 1; return partial_sum(first, last, output, std::multiplies<>{}); }

    // cache-friendly:
    // neighbor objects within the right-most coordinate are neighbors in memory
    template<class TDim, class TCoord>
    auto coordinates_to_index(TDim const& dimensions, TCoord const& coords)
        std::array<std::size_t, dimensions.size()> dimension_product;
        using std::crbegin, std::crend, std::prev;
        partial_product(crbegin(dimensions), prev(crend(dimensions)), begin(dimension_product));
        return std::inner_product(cbegin(dimension_product), cend(dimension_product), crbegin(coords), 0);

constexpr struct matrix_zero_t {} matrix_zero;

template<class T, std::size_t... Dimensions>
class matrix
template<class, std::size_t...> friend class matrix;

    static constexpr std::size_t order      = sizeof...(Dimensions);
    static constexpr std::array  dimensions = { Dimensions... };

    std::array<T, (Dimensions * ...)> _data;

    friend void swap(matrix& lhs, matrix& rhs)
        using std::swap;
        swap(lhs._data, rhs._data);

    matrix()                          = default;
    matrix(matrix&& other)            = default;
    matrix& operator=(matrix&& other) = default;

    matrix(matrix_zero_t) : _data({}) {}

    template<class U>
    matrix(matrix<U, Dimensions...> const& other) { std::copy(cbegin(other._data), cend(other._data), begin(_data)); }

    template<class U>
    matrix& operator=(matrix<U, Dimensions...> const& other)
        matrix o{other};
        swap(*this, o);
        return *this;

    template<class... Args>
    T const& operator()(Args... coordinates) const
    { return _data[_details::coordinates_to_index(dimensions, std::array{coordinates...})]; }

    template<class... Args>
    T& operator()(Args... coordinates)
    { return _data[_details::coordinates_to_index(dimensions, std::array{coordinates...})]; }
} // namespace ysc

Usage demo: http://coliru.stacked-crooked.com/a/1652451c78275436

This is a C++17 implementation; this itself is not set in stone.

  • \$\begingroup\$ How would this be used? \$\endgroup\$
    – Mast
    Commented Oct 26, 2018 at 8:45
  • \$\begingroup\$ @Mast I've included a small example in the usage demo I link under the code. This would be used as a class with value semantics: matrix<...> m = matrix_zero; m(x,y,z) = some_value; m(a,b,c) = m(x,y,z) + m(y,z,x);. \$\endgroup\$
    – YSC
    Commented Oct 26, 2018 at 8:50
  • 1
    \$\begingroup\$ While I like the idea behind this class, I don't think it's better than some nested std::arrays (for example like this rough sketch) with its current feature set. // In my experience, when working with multidimensional arrays you often want some way to refer to some "subslice" of the covered space. That would be a nice feature ;) \$\endgroup\$
    – hoffmale
    Commented Oct 26, 2018 at 20:46

1 Answer 1


Here are some suggestions:

  1. coordinates_to_index should return std::size_t instead of int because int may not be able to hold the required values.

    return std::inner_product(cbegin(dimension_product), cend(dimension_product),
                              crbegin(coords), std::size_t{0}); // instead of int
  2. You don't need to compute the strides and store them. You can use Horner's rule directly.

    std::size_t result = 0;
    for (std::size_t i = 0; i < N; ++i)
        result = result * dimensions[i] + coords[i];
    return result;
  3. Your calculation of the total size fails to consider the case where order == 0. Instead of

    (Dimensions * ...)                   // ill-formed if Dimensions is empty


    (Dimensions * ... * std::size_t{1})

    Moreover, I would suggest exposing the total size, along with the order and dimensions. Like this:

    static constexpr std::size_t size = (Dimensions * ... * std::size_t(1));

    And use it in place of ad-hoc computations further in code.

  4. Your matrix is not Copyable. You declare move operations but not copy operations (the template doesn't make a difference), and the compiler synthesizes deleted copy operations which take precedence over the copying template. You should either default the copy operations explicitly (and maybe destructor as well), or omit the redundant default declarations directly.

  5. Why use :_data({}) when :_data{} is more readable and less verbose?

  6. Consider adding out-of-range detection for operator() in some form.


  1. Currently, the converting copy constructor first default-initializes the elements and then copy assign from the initializer matrix, which is kinda surprising. It is more reasonable to initialize directly (with some template hacks):

    namespace detail {
        template <class T, std::size_t N, class U, std::size_t... Is>
        std::array<T, N> construct(const std::array<U, N>& arr, std::index_sequence<Is...>)
            return std::array<T, N>{static_cast<T>(arr[Is])...};
    template <class T, std::size_t... Dimensions>
    class matrix {
        // ...
        // ...
        template <class U>
        matrix(matrix<U, Dimensions...> const& other)
            : _data{detail::construct<T, N>(other._data, std::make_index_sequence<size>{})}
        // ...

    The above implementation requires access to relies on access to _data of another matrix specialization, which can be granted with a friend declaration or an access function. It also relies on the C++17 guaranteed copy elision.

    Now, the converting copy assignment operator cannot delegate to the converting copy constructor:

    template <class U>
    matrix& operator=(const matrix<U, Dimensions...>& other)
        std::copy_n(_data.begin(), size, other._data.begin());
        return *this;
  2. It makes sense to constrain the converting copy constructor and the converting copy assignment operator with SFINAE. Also, they should be made conditionally explicit and conditionally noexcept:

    template <class U, class = std::enable_if_t<
        std::is_constructible_v<T, const U&> && std::is_convertible_v<const U&, T>
    >> matrix(matrix<U, Dimensions...> const& other)
        noexcept(std::is_nothrow_constructible_v<T, const U&>);
    template <class U, class = std::enable_if_t<
        std::is_constructible_v<T, const U&> && !std::is_convertible_v<const U&, T>
    >> explicit matrix(matrix<U, Dimensions...> const& other)
        noexcept(std::is_nothrow_constructible_v<T, const U&>);

    Similarly for the converting copy assignment operator.

  3. Consider support converting move construction and converting move assignment.

  4. The class can be modified to support constexpr.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.