복합 곱셈 감소를 위한 휴대용 심드 코드 작성 방법

복합 곱셈 감소를 위한 휴대용 심드 코드 작성 방법

나는 복잡한 배열의 곱셈 감소를 계산하기 위해 빠른 심드 코드를 작성하고 싶습니다.표준 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를 대상으로 하고 싶다는 것을 알았다면 gcc에서 sse3 자동 벡터화를 활성화하는 방법을 팔로우하고 다음과 같이 쓸 수 있습니다.

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를 사용하여 빠른 벡터화된 어셈블리 코드를 생성합니다.비록 당신은 입력을 4의 배수로 조절해야 하지만요.조립은 다음과 같습니다.

.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에는 최적이 아닙니다.

sse, avx2 또는 avx512 중 하나에 대해 컴파일될 때 gcc가 최적의 코드를 생성할 C 또는 C++ 코드를 어떻게 작성할 수 있습니까?즉, SIMD 레지스터의 너비가 다를 때마다 항상 별도의 기능을 수기로 작성해야 합니까?

이것을 더 쉽게 해주는 오픈 소스 라이브러리가 있습니까?

Eigen 라이브러리를 사용한 예는 다음과 같습니다.

#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를 활성화(및 -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();
}

참고: 위의 기능은 해당 기능과 직접적으로 일치합니다.fOP 의한 바와 같이 않습니다. 곱셈을 하기 때문입니다.그러나 @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;
}

또는 더 일반적인(C++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

제게 명확하지 않은 이유로, 이것은 실제로 어떤 메모리 주위를 이동하는 실제 방법으로 불리는 방법에 숨겨져 있습니다. 왜 Igen/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;
}

여기서 다른 점은 클랜이 해당 외부 루프를 루프당 2회 곱하기로 해제한다는 것입니다.반면, gcc는 다음과 같이 컴파일될 때 fused-multiply-add 명령어를 사용할 것입니다.-mfma.

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++를 사용하느냐에 따라 assembler 출력이 상당히 달라집니다.전체적으로 g++는 입력 인수를 인라인 로드하는 데 실패하는 경향이 있고, clang은 FMA 연산을 사용하지 못합니다(YMMV...). 근본적으로 생성된 어셈블러 코드를 검사해야 합니다.그리고 더 중요한 것은 코드를 벤치마킹해야 한다는 것입니다(이 루틴이 전체 문제에 얼마나 영향을 미치는지 확실하지 않습니다).

또한 아이겐이 실제로 선형대수 라이브러리라는 점에 주목하고 싶었습니다.순수 휴대용 SIMD 코드 생성을 위해 이것을 이용하는 것은 실제로 의도된 것이 아닙니다.

휴대성이 주된 관심사라면, 독자적인 구문으로 SIMD 지침을 제공하는 라이브러리가 많이 있습니다.그들 대부분은 본질적인 것보다 더 단순하고 휴대하기 쉬운 명시적인 벡터화를 합니다.이 라이브러리(UME::SIMD)는 최근에 출시되었으며 성능이 매우 뛰어납니다.

본 논문(UME::SIMD)에서는 Vc 기반의 인터페이스를 구축하였으며, 이 인터페이스를 UME::SIMD라고 합니다. 이 인터페이스를 통해 프로그래머는 SIMD ISA에 대한 광범위한 지식 없이 SIMD 기능에 접근할 수 있습니다.UME::SIMD는 고유성에 비해 성능 손실 없이 명시적 벡터화를 위한 단순하고 유연하며 휴대용 추상화를 제공합니다.

이것에 대한 일반적인 해결책은 없다고 생각합니다."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개 어레이 요소 포함)을 모두 유지하여 플래그별로 전환할 수는 있지만, 이를 너무 지루하게 유지하여 큰 이점을 얻을 수도 있습니다.

언급URL : https://stackoverflow.com/questions/45298855/how-to-write-portable-simd-code-for-complex-multiplicative-reduction