Fast element access for multi-dimensional representation of contiguous array

Question

I have a multi-dimensional array represented contiguously in memory. I want to keep this representation hidden and just let the user access the array elements as if it were a multi-dimensional one: e.g. my_array[0][3][5] or my_array(0,3,5) or something similar. The dimensions of the object are not determined until runtime, but the object is created with a type that specifies how many dimensions it has. This element look-up will need to be called billions of times, and so should hopefully involve minimal overhead for each call.

I have looked at similar questions but not really found a good solution. Using the [] operator requires the creation of N-1 dimensional objects, which is fine for multi-dimensional structures like vectors-of-vectors because the object already exists, but for a contiguous array it seems like it would get convoluted very quickly and require some kind of slicing through the original array.

I have also looked at overloading (), which seems more promising, but requires specifying the number of arguments, which will vary depending upon the number of dimensions of the array. I have thought about using list initialization or vectors, but wanted to avoid instantiating objects.

I am only a little familiar with templates and figure that there should be some way with C++'s majestic template powers to specify a unique overloading of () for arrays with different types (e.g. different numbers of dimensions). But I have only used templates in really basic generic cases like making a function use both float and double.

I am imagining something like this:

template<typename TDim>
class MultiArray {
public:
  MultiArray() {} //build some things
  ~MultiArray() {} //destroy some things

  // The number of arguments would be == to TDim for the instantiated class
  float& operator() (int dim1, int dim2, ...) {
    //convert to contiguous index and return ref to element
    // I believe the conversion equation is something like:
    // dim1 + Maxdim1 * ( dim2 + MaxDim2 * ( dim3 + MaxDim3 * (...)))
  }

private:
  vector<float> internal_array;
  vector<int> MaxDimX; // Each element says how large each corresponding dim is.
};

So if I initialize this class and attempted to access an element, it would look something like this:

my_array = MultiArray<4>();
element = my_array(2,5,4,1);

How might I go about doing this using templates? Is this even possible?

Answer 1

template<class T>
struct slice {
    T* data = 0;
    std::size_t const* stride = 0;
    slice operator[](std::size_t I)const {
        return{ data + I* *stride, stride + 1 };
    }
    operator T&()const {
      return *data;
    }
    T& operator=(typename std::remove_const<T>::type in)const {
      *data = std::move(in); return *data;
    }
};

store a vector<T> of data, and an std::vector<std::size_t> stride of strides, where stride[0] is the element-stride that the first index wants.

template<class T>
struct buffer {
  std::vector<T> data;
  std::vector<std::size_t> strides;

  buffer( std::vector<std::size_t> sizes, std::vector<T> d ):
    data(std::move(d)),
    strides(sizes)
  {
    std::size_t scale = 1;
    for (std::size_t i = 0; i<sizes.size(); ++i){
      auto next = scale*strides[sizes.size()-1-i];
      strides[sizes.size()-1-i] = scale;
      scale=next;
    }
  }
  slice<T> get(){ return {data.data(), strides.data()}; }
  slice<T const> get()const{ return {data.data(), strides.data()}; }
};

c++14. Live example.

If you use not enough []s it refers to the first element of the subarray in question. If you use too many it does UB. It does zero dimension checking, both in count of dimensions and in size.

Both can be added, but would cost performance.

The number of dimensions is dynamic. You can split buffer into two types, one that owns the buffer and the other that provides the dimensioned view of it.

Answer 2

It seems to me that you could use Boost.MultiArray, boost::multi_array_ref to be more specific. boost::multi_array_ref does exactly what you want: it wraps continuous data array into an object that may be treated as a multidimensional array. You may also use boost::multi_array_ref::array_view for slicing purposes.

I cannot provide you with any benchmark results, but from my experience, I can say that boost::multi_array_ref works pretty fast.

Answer 3

If you can use C++17, so variadic template folding, and row major order, I suppose you can write something like (caution: not tested)

template <template ... Args>
float & operator() (Args ... dims)
 {
   static_assert( sizeof...(Args) == TDim , "wrong number of indexes" );
   // or SFINAE enable instead of static_assert()

   std::size_t pos { 0U };
   std::size_t i   { 0U };

   ( pos *= MaxDimX[i++], pos += dims, ... );

   return internal_array[pos];
 }

OTPS (Off Topic Post Scriptum): your MaxDimX, if I understand correctly, is a vector of dimensions; so should be an unsigned integer, non a signed int; usually, for indexes, is used std::size_t [see Note 1].

OTPS 2: if you know compile time the number of dimensions (TDim, right?) instead of a std::vector, I suggest the use of a std::array; I mean

std::array<std::size_t, TDim>  MaxDimX;

-- EDIT --

If you can't use C++17, you can use the trick of the unused array initialization to obtain something similar.

I mean

template <template ... Args>
float & operator() (Args ... dims)
 {
   using unused = int[];

   static_assert( sizeof...(Args) == TDim , "wrong number of indexes" );
   // or SFINAE enable instead of static_assert()

   std::size_t pos { 0U };
   std::size_t i   { 0U };

   (void)unused { (pos *= MaxDimX[i++], pos += dims, 0) ... };

   return internal_array[pos];
 }

Note 1: as pointed by Julius, the use of a signed or an unsigned integer for indexes is controversial.

So I try to explain better why I suggest to use an unsigned value (std::size_t, by example) for they.

The point is that (as far I know) all Standard Template Library is designed to use unsigned integer for index values. You can see it by the value returned by the size() method and by the fact that access methods that receive an index, as at() or operator[], receive an unsigned value.

Right or wrong, the language itself is designed to return a std::size_t from the old sizeof() and from more recent variadic sizeof...(). The same class std::index_sequence is an alias for std::integer_sequence with a fixed unsigned, again std::size_t, type.

In a world designed to use unsigned integers for indexes, the use of a signed integer for they it's possible but, IMHO, dangerous (because error prone).

Answer 4

I've used this pattern several times when creating a class templates of a matrix class with variable dimensions.

Matrix.h

#ifndef MATRIX_H

template<typename Type, size_t... Dims>
class Matrix {
public:
    static const size_t numDims_ = sizeof...(Dims);

private:
    size_t numElements_;

    std::vector<Type> elements_;
    std::vector<size_t> strides_; // Technically this vector contains the size of each dimension... (its shape)
    // actual strides would be the width in memory of each element to that dimension of the container.
    // A better name for this container would be dimensionSizes_ or shape_ 

public:
    Matrix() noexcept;

    template<typename... Arg>
    Matrix( Arg&&... as  ) noexcept;

    const Type& operator[]( size_t idx ) const;

    size_t numElements() const {
        return elements_.size();
    }

    const std::vector<size_t>& strides() const {
        return strides_;
    }

    const std::vector<Type>& elements() const {
        return elements_;
    }

}; // matrix

#include "Matrix.inl"

#endif // MATRIX_H

Matrix.inl

template<typename Type, size_t... Dims>
Matrix<Type, Dims...>::Matrix() noexcept :
strides_( { Dims... } ) {
    using std::begin;
    using std::end;
    auto mult = std::accumulate( begin( strides_ ), end( strides_ ), 1, std::multiplies<>() );
    numElements_ = mult;
    elements_.resize( numElements_ );
} // Matrix


template<typename Type, size_t... Dims>
template<typename... Arg>
Matrix<Type, Dims...>::Matrix( Arg&&... as ) noexcept   : 
elements_( { as... } ),
strides_( { Dims... } ){
    numElements_ = elements_.size();
} // Matrix

template<typename T, size_t... d>
const T& Matrix<T,d...>::operator[]( size_t idx ) const {
    return elements_[idx];
} // Operator[]

Matrix.cpp

#include "Matrix.h"
#include <vector>
#include <numeric>
#include <functional>
#include <algorithm>

main.cpp

#include <vector>
#include <iostream>
#include "matrix.h"

int main() {

    Matrix<int, 3, 3> mat3x3( 1, 2, 3, 4, 5, 6, 7, 8, 9 );

    for ( size_t idx = 0; idx < mat3x3.numElements(); idx++ ) {
        std::cout << mat3x3.elements()[idx] << " ";
    }

    std::cout << "\n\nnow using array index operator\n\n";

    for ( size_t idx = 0; idx < mat3x3.numElements(); idx++ ) {
        std::cout << mat3x3[idx] << " ";
    }

    std::cout << "\n\ncheck the strides\n\n";
    for ( size_t idx = 0; idx < mat3x3.numDims_; idx++ ) {
        std::cout << mat3x3.strides()[idx] << " ";
    }
    std::cout << "\n\n";

    std::cout << "=================================\n\n";

    Matrix<float, 5, 2, 9, 7> mf5x2x9x7;
    // Check Strides
    // Num Elements
    // Total Size
    std::cout << "The total number of dimensions are: " << mf5x2x9x7.numDims_ << "\n";
    std::cout << "The total number of elements are: " << mf5x2x9x7.numElements() << "\n";
    std::cout << "These are the strides: \n";
    for ( size_t n = 0; n < mf5x2x9x7.numDims_; n++ ) {
        std::cout << mf5x2x9x7.strides()[n] << " ";
    }
    std::cout << "\n";

    std::cout << "The elements are: ";
    for ( size_t n = 0; n < mf5x2x9x7.numElements(); n++ ) {
        std::cout << mf5x2x9x7[n] << " ";
    }
    std::cout << "\n";      

    std::cout << "\nPress any key and enter to quit." << std::endl;
    char c;
    std::cin >> c;

    return 0;

} // main

This is a simple variable multidimensional matrix class of the Same Type <T>

You can create a matrix of floats, ints, chars etc of varying sizes such as a 2x2, 2x3, 5x3x7, 4x9x8x12x2x19. This is a very simple but versatile class.

It is using std::vector<> so the search time is linear. The larger the multi - dimensional matrix grows in dimensions the larger the internal container will grow depending on the size of each dimension; this can "explode" fairly quick if each individual dimensions are of a large dimensional size for example: a 9x9x9 is only a 3 dimensional volumetric matrix that has many more elements than a 2x2x2x2x2 which is a 5 dimensional volumetric matrix. The first matrix has 729 elements where the second matrix has only 32 elements.

I did not include a default constructor, copy constructor, move constructor, nor any overloaded constructors that would accept either a std::container<T> or another Matrix<T,...>. This can be done as an exercise for the OP.

I also did not include any simple functions that would give the size of total elements from the main container, nor the number of total dimensions which would be the size of the strides container size. The OP should be able to implement these very simply.

As for the strides and for indexing with multiple dimensional coordinates the OP would need to use the stride values to compute the appropriate indexes again I leave this as the main exercise.

EDIT - I went ahead and added a default constructor, moved some members to private section of the class, and added a few access functions. I did this because I just wanted to demonstrate in the main function the power of this class even when creating an empty container of its type.

Even more you can take user Yakk's answer with his "stride & slice" algorithm and should easily be able to plug it right into this class giving you the full functionality of what you are looking for.