Эффективно вычислить абсолютные значения вектора std :: complex <float>с AVX2 - PullRequest
Новая
       9

Эффективно вычислить абсолютные значения вектора std :: complex <float>с AVX2

0 голосов
/ 03 декабря 2018

Для некоторых приложений DSP в реальном времени мне нужно вычислить абсолютные значения вектора со сложным значением.

Простая реализация будет выглядеть так: 

computeAbsolute (std::complex<float>* complexSourceVec,
                 float* realValuedDestinationVec,
                 int vecLength)
{
    for (int i = 0; i < vecLength; ++i)
        realValuedDestinationVec[i] = std::abs (complexSourceVec[i]);
}

Я хочу заменить этореализация с оптимизированной версией AVX2, основанной на инстинктах AVX2.Что было бы наиболее эффективным способом реализовать это таким образом?

Примечание. Исходные данные передаются мне через API, к которому у меня нет доступа, поэтому нет возможности изменить макет сложного ввода.вектор для повышения эффективности.

Ответы [ 2 ]

0 голосов
/ 03 декабря 2018

Вдохновленный ответом Дана М. Сначала я реализовал его версию с некоторыми изменениями:

Сначала изменил ее, чтобы использовать более широкие 256-битные регистры, затем пометил временные массивы re и im__attribute__((aligned (32))), чтобы можно было использовать выровненную нагрузку 

void computeAbsolute1 (const std::complex<float>* cplxIn, float* absOut, const int length)
{
    for (int i = 0; i < length; i += 8)
    {
        float re[8] __attribute__((aligned (32))) = {cplxIn[i].real(), cplxIn[i + 1].real(), cplxIn[i + 2].real(), cplxIn[i + 3].real(), cplxIn[i + 4].real(), cplxIn[i + 5].real(), cplxIn[i + 6].real(), cplxIn[i + 7].real()};
        float im[8] __attribute__((aligned (32))) = {cplxIn[i].imag(), cplxIn[i + 1].imag(), cplxIn[i + 2].imag(), cplxIn[i + 3].imag(), cplxIn[i + 4].imag(), cplxIn[i + 5].imag(), cplxIn[i + 6].imag(), cplxIn[i + 7].imag()};
        __m256 x4 = _mm256_load_ps (re);
        __m256 y4 = _mm256_load_ps (im);
        __m256 b4 = _mm256_sqrt_ps (_mm256_add_ps (_mm256_mul_ps (x4,x4), _mm256_mul_ps (y4,y4)));
        _mm256_storeu_ps (absOut + i, b4);
    }
}

Однако ручная перетасовка значений таким образом казалась задачей, которую можно как-то ускорить.Теперь это решение, которое я придумал, которое работает в 2–3 раза быстрее в быстром тесте, скомпилированном clang с полной оптимизацией: 

#include <complex>
#include <immintrin.h>

void computeAbsolute2 (const std::complex<float>* __restrict cplxIn, float* __restrict absOut, const int length)
{    
    for (int i = 0; i < length; i += 8)
    {
        // load 8 complex values (--> 16 floats overall) into two SIMD registers
        __m256 inLo = _mm256_loadu_ps (reinterpret_cast<const float*> (cplxIn + i    ));
        __m256 inHi = _mm256_loadu_ps (reinterpret_cast<const float*> (cplxIn + i + 4));

        // seperates the real and imaginary part, however values are in a wrong order
        __m256 re = _mm256_shuffle_ps (inLo, inHi, _MM_SHUFFLE (2, 0, 2, 0));
        __m256 im = _mm256_shuffle_ps (inLo, inHi, _MM_SHUFFLE (3, 1, 3, 1));

        // do the heavy work on the unordered vectors
        __m256 abs = _mm256_sqrt_ps (_mm256_add_ps (_mm256_mul_ps (re, re), _mm256_mul_ps (im, im)));

        // reorder values prior to storing
        __m256d ordered = _mm256_permute4x64_pd (_mm256_castps_pd(abs), _MM_SHUFFLE(3,1,2,0));
        _mm256_storeu_ps (absOut + i, _mm256_castpd_ps(ordered));
    }
}

Я думаю, что я пойду с этой реализацией, если никто не придетбыстрое решение

Эффективно компилируется с gcc и clang ( в проводнике компилятора Godbolt ).

0 голосов
/ 03 декабря 2018

Очень сложно (если возможно) написать «высокооптимизированную версию AVX2» комплексного пресса, поскольку способ определения комплексных чисел в стандарте предотвращает (особенно из-за всех угловых случаев inf / nan) большую оптимизацию.

Однако, если вас не волнует правильность, вы можете просто использовать -ffast-math, и некоторые компиляторы оптимизируют код для вас.Смотрите вывод gcc: https://godbolt.org/z/QbZlBI

Вы также можете взять этот вывод и создать свою собственную функцию abs со встроенной сборкой.Но да, как уже упоминалось, если вам действительно нужна производительность, вы, вероятно, захотите поменять std::complex на что-то другое.

Я смог получить приличный вывод для вашего конкретного случая со всеми необходимыми перемешиваниями:заполнение вручную небольших массивов re и im.См .: https://godbolt.org/z/sWAAXo Это может быть тривиально расширено для ymm регистров.

В любом случае, вот окончательное решение, адаптированное из этого ответа SO , в котором используются встроенные функции в сочетании с умнымоптимизация компилятора: 

#include <complex>
#include <cassert>
#include <immintrin.h>

static inline void cabs_soa4(const float *re, const float *im, float *b) {
    __m128 x4 = _mm_loadu_ps(re);
    __m128 y4 = _mm_loadu_ps(im);
    __m128 b4 = _mm_sqrt_ps(_mm_add_ps(_mm_mul_ps(x4,x4), _mm_mul_ps(y4,y4)));
    _mm_storeu_ps(b, b4);
}

void computeAbsolute (const std::complex<float>* src,
                 float* realValuedDestinationVec,
                 int vecLength)
{
    for (int i = 0; i < vecLength; i += 4) {
        float re[4] = {src[i].real(), src[i + 1].real(), src[i + 2].real(), src[i + 3].real()};
        float im[4] = {src[i].imag(), src[i + 1].imag(), src[i + 2].imag(), src[i + 3].imag()};
        cabs_soa4(re, im, realValuedDestinationVec);
    }
}

, которая компилируется в простые 

_Z15computeAbsolutePKSt7complexIfEPfi:
        test    edx, edx
        jle     .L5
        lea     eax, [rdx-1]
        shr     eax, 2
        sal     rax, 5
        lea     rax, [rdi+32+rax]
.L3:
        vmovups xmm0, XMMWORD PTR [rdi]
        vmovups xmm2, XMMWORD PTR [rdi+16]
        add     rdi, 32
        vshufps xmm1, xmm0, xmm2, 136
        vmulps  xmm1, xmm1, xmm1
        vshufps xmm0, xmm0, xmm2, 221
        vfmadd132ps     xmm0, xmm1, xmm0
        vsqrtps xmm0, xmm0
        vmovups XMMWORD PTR [rsi], xmm0
        cmp     rax, rdi
        jne     .L3
.L5:
        ret

https://godbolt.org/z/Yu64Wg

...