Element-wise multidimensional array iteration regardless of storage order

I have a n-dimensional array class and I would like to be able to iterate though the elements in a specific order, let's say rows first, and thus regardless of the storage order (row or column major).

The signature of the array class is template<typename T, StorageOrder T_storageOrder, std::size_t ... T_dimensions> class Array, I have a storageIndex() function that converts a list of indices into a 1-D index to retrieve the internal data, and there is an ElementWiseIterator class for the iteration logic.

Everything works well but for the end iterator in the column-major case. The idea behind the end iterator is to use the storageIndex() function to retrieve a 1-D index from a list of indices that is out of range of the array shape defined by the T_dimensions template parameter. In other words, if I have a 2 × 3 array, then I'm trying to retrieve an index that is located at {2, 0}, but this (rightfully) returns me the 3rd element when using a column-major order array, instead of returning the 7th (non existing) element as I'd need to represent the end iterator.

This leads me to believe that my approach is flawed but I have no other idea at the moment and am looking for a bit of inspiration to come up with a clean and generic approach. I had a look at numpy.nditer but couldn't understand anything of the implementation.

I was also wondering if the ElementWiseIterator couldn't be simplified a bit because it currently requires to hold quite some data.

Please find below a reduction of the code.

#include <array>
#include <cstddef>
#include <iostream>
#include <iterator>

// http://ift.tt/16jFSAc
#include <m3ta/product>


enum class StorageOrder {
    rowMajor,
    columnMajor
};

template<StorageOrder> struct StorageOrderTag {};
using RowMajorStorageOrderTag = StorageOrderTag<StorageOrder::rowMajor>;
using ColumnMajorStorageOrderTag = StorageOrderTag<StorageOrder::columnMajor>;


// - Converts a list of indices for a specific shape array into a 1-D index.

template<typename T_Shape>
std::size_t storageIndex(const T_Shape &indices, const T_Shape &shape, RowMajorStorageOrderTag)
{
    std::size_t i = 0;
    std::size_t out = indices[i];
    while (i++ < indices.size() - 1) {
        out = indices[i] + shape[i] * out;
    }

    return out;
}

template<typename T_Shape>
std::size_t storageIndex(const T_Shape &indices, const T_Shape &shape, ColumnMajorStorageOrderTag)
{
    std::size_t i = indices.size() - 1;
    std::size_t out = indices[i];
    while (i-- > 0) {
        out = indices[i] + shape[i] * out;
    }

    return out;
}


//- Element-wise iterator.

template<
    typename T,
    typename T_Data,
    StorageOrder T_storageOrder,
    std::size_t T_dimensionality
>
class ElementWiseIterator
    : public std::iterator<std::bidirectional_iterator_tag, T>
{
private:
    using Shape = std::array<std::size_t, T_dimensionality>;

public:
    T & operator*() const
    { return *_currentElement; }

    ElementWiseIterator & operator++()
    {
        std::size_t i = _shape.size();
        while (i-- > 0) {
            if (_currentIndices[i] < _shape[i] - 1 || i == 0) {
                ++_currentIndices[i];
                break;
            }
        }

        for (++i; i < _currentIndices.size(); ++i) {
            _currentIndices[i] = 0;
        }

        setCurrentElement();
        return *this;
    }

    friend bool operator==(const ElementWiseIterator &iterator1, const ElementWiseIterator &iterator2)
    { return iterator1._currentElement == iterator2._currentElement; }

    friend bool operator!=(const ElementWiseIterator &iterator1, const ElementWiseIterator &iterator2)
    { return !(iterator1 == iterator2); }

private:
    ElementWiseIterator(T_Data *data, const Shape &indices, const Shape &shape)
        : _currentElement(nullptr),
          _data(data),
          _currentIndices(indices),
          _shape(shape)
    {
        setCurrentElement();
    }

    void setCurrentElement()
    {
        std::size_t index = storageIndex(
            _currentIndices,
            _shape,
            StorageOrderTag<T_storageOrder>()
        );

        _currentElement = &(*_data)[index];
    }

    T *_currentElement;
    T_Data *_data;
    Shape _currentIndices;
    Shape _shape;

    template<typename, StorageOrder, std::size_t ...> friend class Array;
};


//- Array class.

template<typename T, StorageOrder T_storageOrder, std::size_t ... T_dimensions>
class Array
{
public:
    static constexpr std::size_t size()
    { return m3ta::product(T_dimensions ...); }

    using Shape = std::array<std::size_t, sizeof ... (T_dimensions)>;

    static constexpr Shape shape()
    { return {T_dimensions ...}; }

protected:
    using Storage = std::array<T, size()>;

public:
    using Iterator = typename Storage::iterator;
    using ElementWiseIterator = ElementWiseIterator<
        T,
        Storage,
        T_storageOrder,
        sizeof ... (T_dimensions)
    >;

    Iterator begin()
    { return _data.begin(); }

    Iterator end()
    { return _data.end(); }

    ElementWiseIterator elementWiseBegin()
    { return ElementWiseIterator(&_data, {0}, shape()); }

    ElementWiseIterator elementWiseEnd()
    {
        // Set the current iterator indices to the first out of range element.
        // Ie: for an a 2x3 array, that would be {2, 0}.
        Shape shape = this->shape();
        return ElementWiseIterator(&_data, {shape[0]}, shape);
    }

private:
    Storage _data;
};



int main(int argc, char **argv)
{
    int i;

    Array<int, StorageOrder::rowMajor, 2, 3> rowArray;
    i = 0;
    for (auto &value : rowArray) {
        value = i++;
    }

    Array<int, StorageOrder::columnMajor, 2, 3> columnArray;
    i = 0;
    for (auto &value : columnArray) {
        value = i++;
    }

    // Below returns 0 1 2 3 4 5 as expected.
    for (auto it = rowArray.elementWiseBegin(); it != rowArray.elementWiseEnd(); ++it) {
        std::cout << *it << " ";
    }
    std::cout << std::endl;

    // Below returns only 0.
    for (auto it = columnArray.elementWiseBegin(); it != columnArray.elementWiseEnd(); ++it) {
        std::cout << *it << " ";
    }
    std::cout << std::endl;

    // This is because the end iterator, pointing to the indices {2, 0} is
    // (fairly enough) converted into the index 2 by the `storageIndex` function,
    // instead of an index beyond the last element.
    std::cout << storageIndex({columnArray.shape()[0]}, columnArray.shape(), ColumnMajorStorageOrderTag()) << std::endl;

    return 0;
}