c++ - Cómo escribir código simd portátil para la reducción multiplicativa compleja

gcc avx

Question

Quiero escribir código SIMD rápido para calcular la reducción multiplicativa de una matriz compleja. En la norma C esto es:

#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 será a lo sumo 50.

Gcc no puede auto-vectorizar la multiplicación compleja pero, como estoy feliz de asumir el compilador gcc y si supiera que quería apuntar a sse3 podría seguir Cómo habilitar la autovectorización de sse3 en gcc y escribir:

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;
}

De hecho, esto produce un código de ensamblaje vectorizado rápido usando gcc. Aunque todavía necesita rellenar su entrada a un múltiplo de 4. El ensamblaje que obtiene es:

.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

Sin embargo, está diseñado para el conjunto de instrucciones simd exacto y no es óptimo para avx2 o avx512, por ejemplo, para lo que necesita cambiar el código.

¿Cómo puede escribir código C o C ++ para el cual gcc producirá un código óptimo cuando se compile para cualquiera de sse, avx2 o avx512? Es decir, ¿siempre tiene que escribir funciones separadas a mano para cada ancho diferente del registro SIMD?

¿Hay bibliotecas de código abierto que hacen esto más fácil?

Answer 1

Aquí hay un ejemplo usando la librería Eigen :

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

Si compilas esto con clang o g ++ y sse o avx habilitados (y -O2), deberías obtener un código de máquina bastante decente. También funciona para algunas otras arquitecturas como Altivec o NEON. Si sabe que la primera entrada de x está alineada, puede usar MapAligned lugar de Map .

Obtendrá un código aún mejor, si conoce el tamaño de su vector en tiempo de compilación utilizando esto:

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

Nota: Las funciones anteriores corresponden directamente a la función f del OP. Sin embargo, como lo señaló @PeterCordes, generalmente es malo almacenar números complejos intercalados, ya que esto requerirá mucha reorganización para la multiplicación. En su lugar, uno debe almacenar partes reales e imaginarias de manera que puedan cargarse directamente un paquete a la vez.

Edición / Addendum : para implementar una estructura de matrices como la multiplicación compleja, puede escribir algo como:

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;
}

O más genérico (usando 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;
}

Cuando se compila con -mavx -O2 , esto se compila en algo como esto (usando 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

Por razones que no son obvias para mí, esto está realmente oculto en un método llamado por el método real, que simplemente se mueve alrededor de alguna memoria. No sé por qué Eigen / gcc no asume que los argumentos ya están alineados correctamente. Si compilo lo mismo con clang 3.8.0 (y los mismos argumentos), se compila para:

    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

Una vez más, el movimiento de la memoria al principio es extraño, pero al menos eso está vectorizado. Sin embargo, para gcc y clang, esto se optimiza cuando se llama en un bucle, sin embargo:

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;
}

La diferencia aquí es que Clang desenrollará ese bucle externo a 2 multiplicaciones por bucle. Por otro lado, gcc utilizará instrucciones fusionadas, multiplicadas y -mfma cuando se compile con -mfma .

La función f8 puede, por supuesto, también generalizarse a dimensiones arbitrarias:

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;
}

Y para reducir el complexX<N> a un solo std::complex , se puede utilizar la siguiente función:

// 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]);
}

Sin embargo, dependiendo de si utilizo clang o g ++, obtengo una salida de ensamblador bastante diferente. En general, g ++ tiene una tendencia a fallar en la carga en línea de los argumentos de entrada, y Clang no puede utilizar las operaciones FMA (YMMV ...) Esencialmente, debe inspeccionar el código del ensamblador generado de todos modos. Y, lo que es más importante, debe evaluar el código (no está seguro, cuánto impacto tiene esta rutina en su problema general).

Además, quería señalar que Eigen en realidad es una biblioteca de álgebra lineal. Explotarlo para la generación de código SIMD portátil puro no es realmente para lo que está diseñado.

Answer 2

No creo que tengas una solución completamente general para esto. Puede aumentar su "vector_size" a 32:

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

También aumenta todas las matrices para tener 8 elementos:

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};
}

Esto hará que el compilador sea capaz de generar código AVX512 (no olvide agregar -mavx512f ), pero hará que su código sea un poco peor en SSE al hacer que las transferencias de memoria no sean óptimas. Sin embargo, no deshabilitará la vectorización SSE.

Puedes mantener ambas versiones (con 4 y con 8 elementos de matriz), cambiando entre ellas por alguna bandera, pero puede ser demasiado tedioso para poco beneficio.

Answer 3

Si la portabilidad es su principal preocupación, aquí hay muchas bibliotecas que proporcionan instrucciones SIMD en su propia sintaxis. La mayoría de ellos hacen la vectorización explícita más simple y portátil que los intrínsecos. Esta biblioteca (UME :: SIMD) se publicó recientemente y tiene un gran rendimiento

En este documento (UME :: SIMD) se ha establecido una interfaz basada en Vc que se denomina UME :: SIMD. Permite al programador acceder a las capacidades de SIMD sin la necesidad de un amplio conocimiento de las ISA de SIMD. UME :: SIMD proporciona una abstracción simple, flexible y portátil para la vectorización explícita sin pérdidas de rendimiento en comparación con los intrínsecos