C++ Iterators for multi-dimensional C arrays

I have a large number of 3 to 6-dimensional C arrays I need to iterate through. More C++'y representation like boost::multi_array isn't an option as these arrays come via the C framework PETSc (using fortran ordering, hence the backward indexing). Straightforward loops end up looking something like this:

  for (int i=range.ibeg; i<=range.iend; ++i){
   for (int j=range.jbeg; j<=range.jend; ++j){
     for (int k=range.kbeg; k<=range.kend; ++k){
       (...)

or even worse:

  for (int i=range.ibeg-1; i<=range.iend+1; ++i){
    for (int j=range.jbeg-1; j<=range.jend+1; ++j){
      for (int k=range.kbeg-1; k<=range.kend+1; ++k){
       for (int ii=0; ii<Np1d; ++ii){
        for (int jj=0; jj<Np1d; ++jj){
         for (int kk=0; kk<Np1d; ++kk){
           data[k][j][i].member[kk][jj][ii] = 
            func(otherdata[k][j][i].member[kk][jj][ii],
                 otherdata[k][j][i].member[kk][jj][ii+1]);

There are many instances like this, with varying ranges on the loop indexes, and it all gets very ugly and potentially error prone. How should one construct iterators for multi-dimensional arrays like this?

A fully templated version was not so hard after all, so here it is in a separate answer, again with live example. If I'm not mistaken, this should have zero overhead on top of custom nested loops. You could measure and let me know. I intend to implement this for my own purposes anyway, that's why I put this effort here.

template<size_t N>
using size = std::integral_constant<size_t, N>;

template<typename T, size_t N>
class counter : std::array<T, N>
{
    using A = std::array<T, N>;
    A b, e;

    template<size_t I = 0>
    void inc(size<I> = size<I>())
    {
        if (++_<I>() != std::get<I>(e))
            return;

        _<I>() = std::get<I>(b);
        inc(size<I+1>());
    }

    void inc(size<N-1>) { ++_<N-1>(); }

public:
    counter(const A& b, const A& e) : A(b), b(b), e(e) { }

    counter& operator++() { return inc(), *this; }

    operator bool() const { return _<N-1>() != std::get<N-1>(e); }

    template<size_t I>
    T& _() { return std::get <I>(*this); }

    template<size_t I>
    constexpr const T& _() const { return std::get <I>(*this); }
};

Instead of operator[] I now have method _ (feel free to rename), which is just a shortcut for std::get, so usage is not so much more verbose than with operator[]:

    for (counter<int, N> c(begin, end); c; ++c)
        cout << c._<0>() << " " << c._<1>() << " " << c._<2>() << endl;

In fact, you may try the previous version

    for (counter<int, N> c(begin, end); c; ++c)
        cout << c[0] << " " << c[1] << " " << c[2] << endl;

and measure, because it may be equivalent. For this to work, switch std::array inheritance to public or declare using A::operator[]; in counter's public section.

What is definitely different is operator++, which is now based on recursive template function inc() and the problematic condition if (n < N - 1) is replaced by a specialization (actually, overload) that has no overhead.

If it turns out that there is overhead eventually, an ultimate attempt would be to replace std::array by std::tuple. In this case, std::get is the only way; there is no operator[] alternative. It will also be weird that type T is repeated N times. But I hope this won't be needed.

Further generalizations are possible, e.g. specifying a (compile-time) increment step per dimension or even specifying arbitrary indirect arrays per dimension, e.g. to simulate

a([3 5 0 -2 7], -4:2:20)

in Matlab-like syntax.

But this needs even more work, and I think you can take it on from here if you like the approach.