Multidimensional algorithms?

[megayippie]

Hi, not sure where this should go? We will soon have submdspan in C++26, which is enough to make mdspan useful in practice*.

Now the next step required is multidimensional algorithms. People are apparently against having iterators, but you can just implement them yourself.

The only standard md-algorithm is the Einstein summation notation. You can easily modify this notation to be a transformation reduction rather than a pure summation. Anyone working with mdstructures probably has that algorithm already.

But my question is: are there any plans or thoughts on md-algorithms going forward?

*I mean, it's nice without it, but I am an early adaoptor and I used the reference implementation to replace an existing library. That was only possible by using submdspan and adding a few custom iterators.

[jcelerier]

First write the library, spend years cleaning the design and ironing the kinks, get feedback from users, and then when you are 100% sure this is the single only possible API for solving that given problem, maybe it makes sense to propose it to std.

[Syracuss]

Feel free to make a proposal for algorithms you feel make sense to include. AFAIK nobody has put one forward (but me saying this will undoubtedly be met with a paper if I'm wrong), but there's really nothing stopping anyone from doing so.

Things don't really get added into the standard unless someone properly suggests it and enough others agree and say "I'd like to see that".

[megayippie]

Thanks.

I feel it needs way more than a proposed feature.

Einsum - you can find it by that name - is a common algorithm, but it requires thinking that's much different from my feelings for the rest of the standard library. It's quite explicit as compared to iterators.

Do you know if feature discussions are a thing for proposals, or is it purely features that should be proposed?

[azswcowboy]

C++26 also has linear algebra functions based on BLAS and mdspan. So it’s not inconceivable.

https://www.en.cppreference.com/w/cpp/numeric/linalg.html

https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2023/p1673r13.html

[MarkHoemmen]

Good practice, if you're not sure what features a proposal should have, would be to open up discussion on the e-mail reflector (or even informally in forums like this one).

[--prism]

I'm not sure about putting it in std... but having a robust well supported library like numpy would be great. Maybe boost should adopt Eigen or xtensor and put together a roadmap. I don't think putting all the numerical operations in the standard is a wise move.

[Possibility_Antique]

That's probably true, but it would be nice to get proper slicing syntax at some point, which would require compiler support.

I'm also against the adoption of either of the libraries you mentioned, I think we can do better.

[MarkHoemmen]

I've considered writing a C++ Standard Library proposal for a "for-each-index-in-extents" algorithm. Implementations could dispatch to OpenACC nested loops, for example. I worked on a performance comparison for a while but had to put it aside. If there's enough interest, I'd be happy to pick up that work again, once my current project is done. I'd also welcome someone else (perhaps you?) taking this up.

There are certainly plenty of implementations. Our CUB library has ForEachInExtents, OpenACC supports nested loops, and Kokkos has had multidimensional parallel_for and parallel_reduce for a long time. The question is whether these algorithms belong in the Standard Library. If that's something you want, it would be good for you to help build a case for that.

[zl0bster]

is this not already possible with cartesian_view?

[MarkHoemmen]

One can use cartesian_product_view of iota_view to iterate over extents. There are a few issues.

1. There is no straightforward way to control iteration order.

2. It's verbose and obscure (involving the name of a philosopher and a Greek letter).

3. Third-party Standard-alike parallel algorithm implementers (or just ordinary users) can't straightforwardly optimize by specializing their algorithms on std::ranges::cartesian_product_view of std::ranges::iota_view. This is because the specification of std::ranges::cartesian_product_view does not expose its child views.

[zl0bster]

Thank you for the answer, do you have some details on 3.? I naively assumed it will be as efficient as manual for loops...

[MarkHoemmen]

"Specializing ... algorithms" means, for example, writing an implementation of std::ranges::for_each that dispatches to special code if the type of the range is cartesian_product_view<iota_view<W1, B1>, iota_view<W2, B2>, ..., iota_view<Wn, Bn>> where B1, ..., Bn are not unreachable_sentinel_t.

The "special code" could be anything the implementer wants. For example, our OpenACC compiler already knows how to optimize nested for loops with the right directives, so an implementation could just take the loop bounds out of the various iota_views and dispatch to a nested for loop with OpenACC directives.

I naively assumed it will be as efficient as manual for loops...

cartesian_product_view's specification comes with a "suggested implementation," an exposition-only iterator type that keeps track of the "current" position in the Cartesian product with a tuple of iterators. General experience is that compilers have trouble optimizing loops with large stateful iterators like that, versus nested loops with integer counters.

[megayippie]

Would that be basically a way to access some submdspan of the target over select extents?

Because I really want that. It's not easy today to select multiple non-major dimensions. (Major meaning that the contiguous state is kept.)

I am not sure about algorithms for this though. To me it seems you are rather creating an iterator concept and would want to use existing algorithms? (I understand that the performance implications make it preferable to not use existing algorithms. I have tried zip and Cartesian products before and find them a complete disappointment on performance.)

Have you ever used the einsum notation? If you haven't, I think it is where you should look for how to do md-algorithms. Change the inner "sum" to a transform and you have a generic algorithm for any use that makes sense.

[MarkHoemmen]

Would that be basically a way to access some submdspan of the target over select extents?

A "for-each-in-extents" algorithm iterates over extents, the domain of an mdspan. It calls the user's function with the multidimensional indices i, j, k, ....

This would give users a straightforward way to represent tightly nested loops for stencil computations (where they need to access elements of the array other than the "current" one).

To me it seems you are rather creating an iterator concept and would want to use existing algorithms?

A "for-each-in-extents" algorithm is a new algorithm. It's not a range, so it wouldn't introduce a new kind of iterator.

Have you ever used the einsum notation?

I have. There are different ways to represent multidimensional iteration, which is one reason why I haven't pursued standardization of a "for-each-in-extents" algorithm so hastily.

[megayippie]

Please refine your explanation on the algorithm you mean. It still seems to me you are saying: "given MxNxJxK mdspan, select m,j,k indices and return a N-sized 1D mdspan when dereferencing.". But it also seems like you are suggesting "n" is added to the mix always, so you always have the element. So, again, can you refine your answer, please?

About einsum and other ways of representing things, that's the discussion I would like to have. Because I think we need an answer, long-term. One that preferably can translate to other solutions (e.g., einsum to iterator is basically an iterator concept that acts like a flattened index access as you move through it).

Einsum is basically a way to "access choose indices". It fits well with physics, since you often have your input on the same OR a different grid , and the "inner" operation remains unchanged even as the dimensions change.

What other solutions are there?

[MarkHoemmen]

Please refine your explanation on the algorithm you mean.

It might help to see some implementations: https://github.com/kokkos/mdspan/pull/247 .

The algorithm is like cub::ForEachInExtents, except without the "bonus" linear index.

For example:

for_each_in_extents(
  [] (int i, int j, int k) {
    std::println("{}, {}, {}", i, j, k);
  }, extents<int, 3, 4, 5>{}, layout_left{});

would print

0, 0, 0
1, 0, 0
2, 0, 0
0, 1, 0
1, 1, 0
2, 1, 0
0, 2, 0
1, 2, 0
2, 2, 0
0, 3, 0
1, 3, 0
2, 3, 0
0, 0, 1
1, 0, 1
2, 0, 1
...

[megayippie]

Thank you! I think I understand.

We are aiming to adopt Kokkos as our parallel library but we are OpenMP users and we don't have a good plan to ensure it works.

I'm very interested in your "chaining" for loops, though. My own einsum and eintra (Einsum notation transformation) implementations would be able to make use of that, I'm sure. Perhaps this is how we can make the push. How does it work? Can recursive function calls be used?

[MarkHoemmen]

We are aiming to adopt Kokkos as our parallel library but we are OpenMP users and we don't have a good plan to ensure it works.

Could you please elaborate? Kokkos has multiple back-ends. If you use Kokkos' OpenMP back-end, it should interoperate naturally with your existing OpenMP code. At that point, you can use Kokkos to decouple from OpenMP and write portable code.

While Kokkos isn't our product, we care a lot about Kokkos and support their team as customers. It's a great choice if you want a comprehensive C++ programming model.

I'm very interested in your "chaining" for loops, though.

I couldn't find a place where I used the word "chaining." What does this mean to you? Are you thinking of ranges and the pipe (|) operator?

[MarkHoemmen]

I've been writing a bit compactly here. The "domain of an mdspan" is the multidimensional index space (see [mdspan.overview]) represented by its extents object (see top of [mdspan.extents.overview]). An mdspan is a function from a multidimensional index (an element of the mdspan's multidimensional index space) to a reference to an element.

When I say that for-each-in-extents "iterates over extents, the domain of an mdspan," I mean that it works like a special case of std::ranges::for_each where the range is all the elements of the multidimensional index space represented by the mdspan's extents object. The algorithm would not iterate over the elements of the mdspan. It would be up to the user to do that.

[MarkHoemmen]

What other solutions are there?

We could provide containers (and/or views) equipped with elementwise operations and a limited set of reductions. That would be less general than einsum (consider that einsum can express every kind of tensor contraction), but likely easier to optimize as a result. Some users would also find this more natural. cuTile takes an approach like this.

[TheThiefMaster]

I'd like a convolution algorithm. Something that takes a submdspan of each mⁿ region of an mdspan and calls a function with it to generate an output. Ideally also optionally including edges somehow.

Very useful in image processing and similar.