C++ Iterators for multi-dimensional C arrays

Question

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?

Answer 1

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.

Answer 2

A full-blown n-dimensional iterator is not necessary in your simple case of nested for loops. Since a single traversal is only needed, a simple counter is enough, which is easily custom-made like this:

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

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

    counter& operator++()
    {
        for (size_t n = 0; n < N; ++n)
        {
            if (++i[n] == e[n])
            {
                if (n < N - 1)
                    i[n] = b[n];
            }
            else
                break;
        }

        return *this;
    }

    operator bool() { return i[N - 1] != e[N - 1]; }

    T&       operator[](size_t n)       { return i[n]; }
    const T& operator[](size_t n) const { return i[n]; }
};

It is then very easy to use this counter like this:

int main()
{
    constexpr size_t N = 3;
    using A = std::array<int, N>;

    A begin = {{0, -1,  0}};
    A end   = {{3,  1,  4}};

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

assuming there's an operator << for counter. See live example for full code.

The innermost condition if (n < N - 1) accounts for being able to check for termination and is not so efficient to always check. It was not so apparent to me how to factor it out, but anyhow it only takes place when we advance to the next "digit" of the counter, not at every increment operation.

Instead of using c[0], c[1], c[2] etc., it is more efficient to use std::get if counter derives std::array instead of having member i (while b,e remain members). This idea can be extended towards a compile-time recursive implementation of operator++ (operator bool as well) that would eliminate the for loop inside it, along with the problematic check discussed above. operator[] would be discarded in this case. But all this would make counter code more obscure and I just wanted to highlight the idea. It would also make usage of counter a bit more verbose, but that's a price you'd need to pay for efficiency.

Of course, a full-blown n-dimensional iterator can be built by extending counter with more methods and traits. But to make it generic enough may be a huge undertaking.