Vectorization with loops

[FrancoisFayard]

Hi,

I would like to to use xt::xarray as a multidimensional array in my code where some algorithm are written with loops. For instance, here is an algorithm that compute the diameter of a set of points in the plane. The problem I found with xt::xarray is that the data structure does not let the compiler know that on this code point(j, 0) and point(j + 1, 0) are next to each other in memory. It generates suboptimal code. For instance it does not let them compiler vectorize the code.

What is the recommended way to do that with xtensor?

François

const std::size_t n = 10000;
std::vector<std::size_t> shape = {n, 2};
std::vector<std::size_t> stride = {1, n};
xt::xarray<float, xt::layout_type::dynamic> point(shape, stride);
// Fill point

float ans = 0.0f;
for (std::size_t i = 0; i < n; ++i) {
  for (std::size_t j = i + 1; j < n; ++j) {
    const float dx = point(i, 0) - point(j, 0);
    const float dy = point(i, 1) - point(j, 1);
    const float distance = dx * dx + dy * dy;
    if (distance > ans) {
      ans = distance;
    }
  }
}
ans = std::sqrt(ans);

[JohanMabille]

Hi,

Sorry for the late response. xtensor optimization is done when operating on full arrays and expressions, it is more difficult to achieve when operating on single values in the array. I don't see what kind of vectorization you expect from the compiler in your example, even if point(j, 0) and point(j+1, 0) are contiguous in memory, they are involved in a same expression. If point(i, 0) and point(i, 1) were contiguous, we could expect the compiler to vectorize the computation of dx and dy.

Anyway we provide access to the one-dimensional container holding the data of the arrays:

using container_type = xt::array<float, xt::layout_type::dynamic>::container_type;
const container_type& dc = point.data();
for(std:size_t i = 0; i < dc.size(); ++i)

We also provide access to the raw buffer for low-level operations:

using value_type = xt::array<float, xt::layout_type::dynamic>::value_type;
value_type* vc = point.raw_data(); 

You can then access the elements of the container or the buffer with the traditional operator[]. Note that you'll have to deal with the strides and the shape to retrieve the data corresponding to your points.

[FrancoisFayard]

Hi Johan,

This code can be perfectly vectorized if the compiler knows at compiler time that time(j, 0) and time(j + 1, 0) are close together in memory. It will make the inner loop work by chunk of 8 (with AVX2) or 16 (with AVX512).

Thanks for the tip in how to get the raw pointer.

Best regards,
François

[wolfv]

I checked the following minimal example in godbolt: https://godbolt.org/g/jzGwqT
Will do the same for the xtensor solution to see if the loop code is very different.

Note that once we have merged xfixed we are able compile strides at compile time which will make these optimizations definitly feasible in xtensor as well.

[FrancoisFayard]

Hi. The godbolt code is kind of cheating because on the the file, the compiler can see the value of stride, so it can optimize everything. But if you compile some code that receives a Mat object without knowing anything about it, it won't be able to vectorize it.

It you want vectorization using elements access, you must have the "unit stride" information in the class itself. I don't see any workaround.

[wolfv]

It was not by accident that it was cheating. here is a side by side comparison of the inner loop from the example above (on the right) and on the left xtensor code + xtensor assembly.

Note that i used xtensor instead of xarray (with compile time known dimensions) and I used the xt::layout_type::column_major which your strides effectively are. The loop is almost identical to my eyes.

Here is the code (let me know if you find a mistake! :)

#include <xtensor/xtensor.hpp>

double test()
{
  const std::size_t n = 10000;
  std::array<std::size_t, 2> shape = {n, 2};
  xt::xtensor<float, 2, 
              xt::layout_type::column_major> 
              point(std::move(shape));
  // Fill point

  float ans = 0.0f;
  for (std::size_t i = 0; i < n; ++i) {
    for (std::size_t j = i + 1; j < n; ++j) {
      const float dx = point(i, 0) - point(j, 0);
      const float dy = point(i, 1) - point(j, 1);
      const float distance = dx * dx + dy * dy;
      if (distance > ans) {
        ans = distance;
      }
    }
  }
  ans = std::sqrt(ans);
  return ans;
}

[wolfv]

Also, if we find a workaround for the MSVC bug, we have the #586 (xfixed container) where the shape (and thus the strides) are part of the type and a lot of the accessor functions are constexpr. This should allow for very optimal code to be generated by the compiler.

An example would be that you can do something like: xfixed<double, xt::xshape<4, 3, 5>>( ... ).

[wolfv]

by the way, while it's not explicit and we don't do any optimizations based on that, if you select layout_type:row_major or layout_type::column_major the type will actually be constrained to have a unit stride (that's also indicated by the static variable "is_contigous" on the xarray and xtensor type.

[wolfv]

Hi @FrancoisFayard ,

I've just published a blog post with some performance tricks & tips around xtensor, which also introduces the xtensorf container (renamed from xfixed) which statically knows it's strides at compile time: https://medium.com/@wolfv/faster-xtensor-87fcbe7f5293

Maybe it's interesting to you?
We're really trying to make xtensor faster, so if you have any tricks or tips from your experience with HPC, please don't hesitate to let us know!

[FrancoisFayard]

Hi @wolfv. Sorry for being late as I was skiing and the next week will be very busy.

Having a static tensor does solve the problem for size that are known at compilation time but does not solve the general case. Moreover, size which are known at compilation time are often small and usually do not benefit that much from vectorization.

What is really needed, is that the compiler knowns that tensor(i, j) and tensor(i, j+ 1) are right next to each other in memory at compile time (if the tensor is C ordered). You end up having to define a tensor that is C-ordered and a tensor that is Fortran-ordered. I'll give you benchmark when I have time, maybe sometimes next week.

[SylvainCorlay]

@FrancoisFayard the "unit-stride" information is part of the type in the case of xt::layout_type::column_major and xt::layout_type::row_major layouts, and it is exploited in the xsimd-based vectorization.

[FrancoisFayard]

If I remember correctly, those layouts (such as xt::layout_type::row_major) impose the sizes of the tensor to be known at compile time, don't they?

[SylvainCorlay]

If I remember correctly, those layouts (such as xt::layout_type::row_major) impose the sizes of the tensor to be known at compile time, don't they?

Nope they don't require sizes or strides to be known at compile time, they only allow to make assumptions such as consequent elements being adjacent in memory.

We have three types of layout: row-major, column-major and dynamic. Default in xtensor containers are row-major.

[SylvainCorlay]

The new features added by Wolf are about compile-time dimensions, which enables even more optimization by the compiler, such as completely unrolling certain loops etc...