Hi Rob, Thanks for the response. I also think this is a pretty important use case for this great library. Let me justify why I think an Iterator based approach is the best: * For maximum efficiency (I will be using this in numerical simulations), we will want have an incrementing pointer within the data structure pointing to the actual object, and not just the iterator. * Also, if we have a single forward iterator, then we should be able to reuse a bunch of standard algorithms in the cases where we do not need the actual index values or the order of execution doesn't matter. * In the iterator algorithm, I would love to check the ordering of the data from the template parameters. If it is fortran ordering then we would do the loops in the opposite order for example, and this should be based on querying the data structure in the iterator code. * As for your range based use case, I don't know enough about how iterators and ranges interact to see how that would work. * As to implementation of this, thanks for the pointer to Fusion. This indeed seems to be the best approach since the loops can be generated at compile time, and compilers could then do some very fancy optimization on them. Sadly, I am not a template metaprogrammer, or a very good C++ programmer in general, and am terrified to learn it and lose sight of my non-CS research. Now for the semantics of this, it seems to me that the way the iterator concepts are designed, the distance() function is intended to pass back the related index information. Now my question is whether it is OK for the semantics of distance to pass back a boost::array<NUMDIMS> instead? Seems reasonable to me that it would not be based on passing back a scalar value in this case because that is arbitrarily dependent on the column vs. row ordering. And we can easily generate a boost::array from the distance function by looking at the distance between the .begin() for all of our current positions in the sub-iterators? If the iterator framework was there and the distance implemented as vectors, I can imagine the code looking like: int n1 = 3; int n2 = 4; int n3 = 5; boost::multi_array<double, 3> A(boost::extents[n1][n2][n3]); //NOTE SURE OF THE RIGHT SEMANTICS FOR GETTING AN ITERATOR OUTSIDE OF ONE GENERATED BY THE CONTAINER. multi_array_iterator<double, 3> begin = mult_sequential_iterator_begin(A); multi_array_iterator<double, 3> end = mult_sequential_iterator_end(A); for(mult_sequential_iterator<double, 3> iter = begin; iter != end; ++iter) *iter = f(distance(begin, iter); //POSSIBLE IMPLEMENTATION OF DISTANCE template<int NUMDIMS> boost::array<NUMDIMS> distance(const mult_sequential_iterator& iter1, const mult_sequential_iterator& iter1) const { static_assert(iter1::dimensionality == iter2::dimensionality) boost::array<NUMDIMS> dist; foreach dimension //PSEUDOCODE, FROM FUSION???? array<dimension> = distance(iter1.get_sub_iter<dimension>(), iter2.get_position<dimension>()); //here get_sub_iter returns the sub-iterator stored for that current dimension in the iterator. } template<int NUMDIMS> //HERE THE FUNCTION WOULD BE BASED ON A BOOST::ARRAY. double f(const boost::array<NUMDIMS>& indices) { ...} //AND LETS SAY WE WANTED TO NEGATE THE ENTIRE DATA STRUCTURE, NOW CAN USE STD: int n1 = 3; int n2 = 4; int n3 = 5; boost::multi_array<double, 3> A(boost::extents[n1][n2][n3]); boost::multi_array<double, 3> B(boost::extents[n1][n2][n3]); std::transform(mult_sequential_iterator_begin(A), mult_sequential_iterator_end(A), mult_sequential_iterator_begin(B), negate<double>()); //AND IF WE ADDED THIS ITERATOR TO THE CODE FOR THE MULTDIM: A.begin_sequential() could return the sequential iterator...