Подтвердить что ты не робот

Как написать переносный simd-код для комплексного мультипликативного сокращения

Я хочу написать быстрый simd-код для вычисления мультипликативной редукции сложного массива. В стандарте C это:

#include <complex.h>
complex float f(complex float x[], int n ) {
   complex float p = 1.0;
   for (int i = 0; i < n; i++)
      p *= x[i];
   return p;
}

n будет не более 50.

Gcc не может автоматически векторизовать сложное умножение, но, поскольку я с удовольствием принимаю компилятор gcc, и если бы я знал, что хочу нацелить sse3, я мог бы следовать Как включить autovectorization sse3 в gcc и напишите:

typedef float v4sf __attribute__ ((vector_size (16)));
typedef union {
  v4sf v;
  float e[4];
} float4
typedef struct {
  float4 x;
  float4 y;
} complex4;
static complex4 complex4_mul(complex4 a, complex4 b) {
  return (complex4){a.x.v*b.x.v -a.y.v*b.y.v, a.y.v*b.x.v + a.x.v*b.y.v};
}
complex4 f4(complex4 x[], int n) {
  v4sf one = {1,1,1,1};
  complex4 p = {one,one};
  for (int i = 0; i < n; i++) p = complex4_mul(p, x[i]);
  return p;
}

Это действительно дает быстрый векторный код сборки с использованием gcc. Несмотря на то, что вам по-прежнему нужно вводить свой ввод в несколько раз. Сборка, которую вы получаете, это:

.L3:
    vmovaps xmm0, XMMWORD PTR 16[rsi]
    add     rsi, 32
    vmulps  xmm1, xmm0, xmm2
    vmulps  xmm0, xmm0, xmm3
    vfmsubps        xmm1, xmm3, XMMWORD PTR -32[rsi], xmm1
    vmovaps xmm3, xmm1
    vfmaddps        xmm2, xmm2, XMMWORD PTR -32[rsi], xmm0
    cmp     rdx, rsi
    jne     .L3

Однако он предназначен для точного набора команд simd и не является оптимальным для avx2 или avx512, например, для которого вам нужно изменить код.

Как вы можете написать код C или С++, для которого gcc будет создавать оптимальные код при компиляции для любого из sse, avx2 или avx512? То есть, вам всегда приходится писать отдельные функции вручную для каждой разной ширины регистра SIMD?

Есть ли библиотеки с открытым исходным кодом, которые облегчают это?

4b9b3361

Ответ 1

Вот пример с помощью Eigen library:

#include <Eigen/Core>
std::complex<float> f(const std::complex<float> *x, int n)
{
    return Eigen::VectorXcf::Map(x, n).prod();
}

Если вы скомпилируете это с помощью clang или g++ и sse или avx enabled (и -O2), вы должны получить довольно приличный машинный код. Он также работает для некоторых других архитектур, таких как Altivec или NEON. Если вы знаете, что первая запись x выровнена, вы можете использовать MapAligned вместо Map.

Вы получаете еще лучший код, если вы случайно знаете размер вашего вектора при компиляции:

template<int n>
std::complex<float> f(const std::complex<float> *x)
{
    return Eigen::Matrix<std::complex<float>, n, 1> >::MapAligned(x).prod();
}

Примечание. Вышеуказанные функции непосредственно соответствуют функции f OP. Однако, как указывал @PeterCordes, обычно сложно хранить сложные числа, чередующиеся, так как это потребует много перетасовки для умножения. Вместо этого нужно хранить реальные и мнимые части так, чтобы они могли сразу загружать один пакет за один раз.

Изменить/добавить. Чтобы реализовать структуру массивов, например комплексное умножение, вы можете написать что-то вроде:

typedef Eigen::Array<float, 8, 1> v8sf; // Eigen::Array allows element-wise standard operations
typedef std::complex<v8sf> complex8;
complex8 prod(const complex8& a, const complex8& b)
{
    return a*b;
}

Или более общий (с использованием С++ 11):

template<int size, typename Scalar = float> using complexX = std::complex<Eigen::Array<Scalar, size, 1> >;

template<int size>
complexX<size> prod(const complexX<size>& a, const complexX<size>& b)
{
    return a*b;
}

При компиляции с помощью -mavx -O2 это скомпилируется примерно так (с помощью g++ - 5.4):

    vmovaps 32(%rsi), %ymm1
    movq    %rdi, %rax
    vmovaps (%rsi), %ymm0
    vmovaps 32(%rdi), %ymm3
    vmovaps (%rdi), %ymm4
    vmulps  %ymm0, %ymm3, %ymm2
    vmulps  %ymm4, %ymm1, %ymm5
    vmulps  %ymm4, %ymm0, %ymm0
    vmulps  %ymm3, %ymm1, %ymm1
    vaddps  %ymm5, %ymm2, %ymm2
    vsubps  %ymm1, %ymm0, %ymm0
    vmovaps %ymm2, 32(%rdi)
    vmovaps %ymm0, (%rdi)
    vzeroupper
    ret

По причинам, не очевидным для меня, это фактически скрыто в методе, который вызывается фактическим методом, который просто перемещается вокруг некоторой памяти - я не знаю, почему Eigen/gcc не предполагает, что аргументы уже правильно выровнен. Если я скомпилирую то же самое с clang 3.8.0 (и теми же аргументами), он скомпилируется только:

    vmovaps (%rsi), %ymm0
    vmovaps %ymm0, (%rdi)
    vmovaps 32(%rsi), %ymm0
    vmovaps %ymm0, 32(%rdi)
    vmovaps (%rdi), %ymm1
    vmovaps (%rdx), %ymm2
    vmovaps 32(%rdx), %ymm3
    vmulps  %ymm2, %ymm1, %ymm4
    vmulps  %ymm3, %ymm0, %ymm5
    vsubps  %ymm5, %ymm4, %ymm4
    vmulps  %ymm3, %ymm1, %ymm1
    vmulps  %ymm0, %ymm2, %ymm0
    vaddps  %ymm1, %ymm0, %ymm0
    vmovaps %ymm0, 32(%rdi)
    vmovaps %ymm4, (%rdi)
    movq    %rdi, %rax
    vzeroupper
    retq

Опять же, движение памяти в начале является странным, но, по крайней мере, оно векторизовано. Однако для gcc и clang это оптимизируется при вызове в цикле:

complex8 f8(complex8 x[], int n) {
    if(n==0)
        return complex8(v8sf::Ones(),v8sf::Zero()); // I guess you want p = 1 + 0*i at the beginning?

    complex8 p = x[0];
    for (int i = 1; i < n; i++) p = prod(p, x[i]);
    return p;
}

Различие заключается в том, что clang будет разворачивать этот внешний цикл на 2 умножения на цикл. С другой стороны, при компиляции с -mfma gcc будет использовать команды с плавным перемножением-добавлением.

Конечно, функция f8 также может быть обобщена на произвольные размеры:

template<int size>
complexX<size> fX(complexX<size> x[], int n) {
    using S= typename complexX<size>::value_type;
    if(n==0)
        return complexX<size>(S::Ones(),S::Zero());

    complexX<size> p = x[0];
    for (int i = 1; i < n; i++) p *=x[i];
    return p;
}

И для уменьшения complexX<N> до одного std::complex можно использовать следующую функцию:

// only works for powers of two
template<int size> EIGEN_ALWAYS_INLINE
std::complex<float> redux(const complexX<size>& var) {
    complexX<size/2> a(var.real().template head<size/2>(), var.imag().template head<size/2>());
    complexX<size/2> b(var.real().template tail<size/2>(), var.imag().template tail<size/2>());
    return redux(a*b);
}
template<> EIGEN_ALWAYS_INLINE
std::complex<float> redux(const complexX<1>& var) {
    return std::complex<float>(var.real()[0], var.imag()[0]);
}

Однако, в зависимости от того, использую ли я clang или g++, я получаю совсем другой выход ассемблера. В целом, g++ имеет тенденцию к сбою встроенной загрузки входных аргументов, а clang не использует операции FMA (YMMV...) По сути, вам все равно нужно проверить сгенерированный код ассемблера. И что еще более важно, вам следует сравнить код (не уверен, насколько сильно эта процедура работает в вашей общей проблеме).

Кроме того, я хотел бы отметить, что Eigen фактически является библиотекой линейной алгебры. Использование его для создания чистого портативного кода SIMD на самом деле не предназначено для.

Ответ 2

Если важна Портативность, существует множество библиотек здесь, которые предоставляют инструкции SIMD в их собственном синтаксисе. Большинство из них делают явную векторность более простой и переносимой, чем внутренняя. Эта библиотека (UME:: SIMD) недавно опубликована и имеет отличную производительность

В этот документ (UME:: SIMD) интерфейс, основанный на Vc, который называется UME:: SIMD. Это позволяет программисту получить доступ к SIMD возможности без необходимости в обширных знаниях SIMD ISA. UME:: SIMD обеспечивает простую, гибкую и портативную абстракцию для явная векторизация без потерь производительности по сравнению с встроенные функции

Ответ 3

Я не думаю, что у вас есть полное общее решение. Вы можете увеличить свой "vector_size" до 32:

typedef float v4sf __attribute__ ((vector_size (32)));

Также увеличьте все массивы, чтобы иметь 8 элементов:

typedef float v8sf __attribute__ ((vector_size (32)));

typedef union {
  v8sf v;
  float e[8];
} float8;
typedef struct {
  float8 x;
  float8 y;
} complex8;
static complex8 complex8_mul(complex8 a, complex8 b) {
  return (complex8){a.x.v*b.x.v -a.y.v*b.y.v, a.y.v*b.x.v + a.x.v*b.y.v};
}

Это сделает компилятор способным генерировать код AVX512 (не забудьте добавить -mavx512f), но сделает ваш код немного хуже в SSE, сделав невозможным передачу памяти. Тем не менее, это, конечно же, не приведет к отключению векторизации SSE.

Вы можете сохранить обе версии (с 4 и 8 элементами массива), переключаясь между ними с помощью некоторого флага, но это может быть слишком утомительно для небольшой пользы.