Как использовать std::experimental::simd?

Question

Я попытался выполнить пример, приведенный на github std::simd, но моя векторизованная версия оказалась в 2-3 раза медленнее. Как правильно пользоваться?

В документальном фильме отсутствует серьезная документация и другие примеры использования. Никаких конструкторов в списке нет и т. Д. Я уверен, что использую его неправильно, но с ограниченной документацией я не знаю, как действовать.

g++ -o test test.cpp --std= C++2a -O0

#include <array>
#include <chrono>
#include <cstdlib>
#include <experimental/simd>
#include <iostream>
#include <random>

using std::experimental::native_simd;
using Vec3D_v = std::array<native_simd<float>, 3>;
native_simd<float> scalar_product(const Vec3D_v& a, const Vec3D_v& b) {
  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
using Vec3D = std::array<float, 3>;
float scalar_product(const std::array<float, 3>& a, const std::array<float, 3>& b) {
  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

int main(){
  constexpr std::size_t VECREG_SIZE = native_simd<float>::size();
  std::array<Vec3D, VECREG_SIZE * 1000> arr;
  std::array<Vec3D_v, VECREG_SIZE * 1000> arr_v;
  std::random_device rd;
  std::mt19937 generator(rd());
  std::uniform_real_distribution<float> distribution(0.f, 1.f);
  for( std::size_t i = 0; i < arr.size(); ++i ){
    arr[i] = {distribution(generator), distribution(generator), distribution(generator)};
    arr_v[i] = {distribution(generator), distribution(generator), distribution(generator)};
  }
  float result = 0.f;
  auto start = std::chrono::high_resolution_clock::now();
  for( std::size_t i = 1; i < arr.size(); ++i ){
    result += scalar_product(arr_v[i-1], arr_v[i])[0];
  }
  auto end = std::chrono::high_resolution_clock::now();
  auto elapsed = end - start;
  std::cout << "VC: " << elapsed.count() << '\n' << std::endl;

  result = 0;
  start = std::chrono::high_resolution_clock::now();
  for( std::size_t i = 1; i < arr.size(); ++i ){
    result += scalar_product(arr[i-1], arr[i]);
  }
  end = std::chrono::high_resolution_clock::now();
  elapsed = end - start;
  std::cout << "notVC: " << elapsed.count() << '\n';
  return EXIT_SUCCESS;
}

Answer 1

Проблема 1. При использовании инструкций SIMD взимается начальная стоимость. Возьмите свой код и зациклите его три раза (я компилирую с -O3и распечатать result в противном случае большая часть кода удаляется):

$ ./test
VC: 37240 (result: 5986.1)
notVC: 18668 (result: 5983.29)
VC: 26177 (result: 5986.1)
notVC: 18516 (result: 5983.29)
VC: 25895 (result: 5986.1)
notVC: 18083 (result: 5983.29)

Сборка основной петли для _v версия теперь гласит:

    1840:       c5 fc 28 d5             vmovaps %ymm5,%ymm2
    1844:       c5 fc 28 28             vmovaps (%rax),%ymm5
    1848:       c5 fc 28 cc             vmovaps %ymm4,%ymm1
    184c:       c5 fc 28 c3             vmovaps %ymm3,%ymm0
    1850:       c5 fc 28 60 20          vmovaps 0x20(%rax),%ymm4
    1855:       c5 fc 28 58 40          vmovaps 0x40(%rax),%ymm3
    185a:       48 83 c0 60             add    $0x60,%rax
    185e:       c5 d4 59 d2             vmulps %ymm2,%ymm5,%ymm2
    1862:       c4 e2 6d 98 cc          vfmadd132ps %ymm4,%ymm2,%ymm1
    1867:       c4 e2 75 98 c3          vfmadd132ps %ymm3,%ymm1,%ymm0
    186c:       c5 ca 58 f0             vaddss %xmm0,%xmm6,%xmm6
    1870:       48 39 d8                cmp    %rbx,%rax
    1873:       75 cb                   jne    1840 <main+0x6f0>

Проблема 2: при каждом повороте петли вы переводите native_simd<float> привести к float используя [0]оператор. Это может иметь ужасные последствия, но компилятор достаточно умен, чтобы этого не делать, как показывает приведенная выше сборка.

Проблема 3. Как видим, nativeпросто инструктирует компилятор поместить значения в регистры SIMD. В этом нет особой выгоды: где здесь аспект множественных данных? Что вы хотите сделать, так это упаковать свой 3D-вектор в один регистр SIMD и переписать цикл, чтобы собрать каждое измерение скалярного произведения в одном компоненте. Наконец, вы должны взять сумму всех компонентов:

using std::experimental::fixed_size_simd;
using Vec3D_v = fixed_size_simd<float, 3>;

и

  for( std::size_t i = 1; i < arr.size(); ++i ){
    result_v += arr_v[i-1] * arr_v[i];
  }
  float result = std::experimental::reduce (result_v);

Запустив это, у нас есть:

$ ./test
VC: 14958 (result: 2274.7)
notVC: 5279 (result: 2274.7)
VC: 4718 (result: 2274.7)
notVC: 5177 (result: 2274.7)
VC: 4720 (result: 2274.7)
notVC: 5132 (result: 2274.7)

А сборка основной петли - вот такая красивая деталь:

    1588:       c5 f8 28 d0             vmovaps %xmm0,%xmm2
    158c:       c5 f8 28 00             vmovaps (%rax),%xmm0
    1590:       48 83 c0 10             add    $0x10,%rax
    1594:       c4 e2 79 b8 ca          vfmadd231ps %xmm2,%xmm0,%xmm1
    1599:       48 39 c3                cmp    %rax,%rbx
    159c:       75 ea                   jne    1588 <main+0x438>

Здесь каждый %xmmрегистр хранит сразу 3 значения с плавающей запятой. Кроме того, компилятор сильно оптимизирует второй цикл для использования инструкций AVX, поэтому выигрыш не так уж и важен (но все еще существует!).

---

Полный код:

#include <array>
#include <chrono>
#include <cstdlib>
#include <experimental/simd>
#include <iostream>
#include <random>

using std::experimental::fixed_size_simd;
using Vec3D_v = fixed_size_simd<float, 3>;

using Vec3D = std::array<float, 3>;
float scalar_product (const std::array<float, 3> &a, const std::array<float, 3> &b) {
  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

int main () {
  constexpr std::size_t VECREG_SIZE = fixed_size_simd<float, 3>::size ();
  std::array<Vec3D, VECREG_SIZE * 1000> arr;
  std::array<Vec3D_v, VECREG_SIZE * 1000> arr_v;
  std::random_device rd;
  std::mt19937 generator (rd ());
  std::uniform_real_distribution<float> distribution (0.f, 1.f);

  for (std::size_t i = 0; i < arr.size (); ++i) {
    arr[i] = {distribution (generator), distribution (generator), distribution (generator) };

    for (int j = 0; j < 3; ++j)
      arr_v[i][j] = arr[i][j];
  }

  Vec3D_v result_v;

  for (int iter = 0; iter < 3; ++iter) {

    for (int j = 0; j < 3; ++j)
      result_v[j] = 0.f;

    auto start = std::chrono::high_resolution_clock::now ();

    for (std::size_t i = 1; i < arr.size (); ++i) {
      result_v += arr_v[i - 1] * arr_v[i];
    }

    float result = std::experimental::reduce (result_v);
    auto end = std::chrono::high_resolution_clock::now ();
    auto elapsed = end - start;
    std::cout << "VC: " << elapsed.count () << " (result: " << result << ")" << std::endl;

    result = 0;
    start = std::chrono::high_resolution_clock::now ();

    for (std::size_t i = 1; i < arr.size (); ++i) {
      result += scalar_product (arr[i - 1], arr[i]);
    }

    end = std::chrono::high_resolution_clock::now ();
    elapsed = end - start;
    std::cout << "notVC: " << elapsed.count () << " (result: " << result << ")" << std::endl;
  }

  return EXIT_SUCCESS;
}