Comment écrire un code simd portable pour la réduction multiplicative complexe ?

Question

Je veux écrire un code simd rapide pour calculer la réduction multiplicative d'un tableau complexe. En C standard c'est :

#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 sera au maximum de 50.

Gcc ne peut pas vectoriser automatiquement les multiplications complexes mais, comme je suis heureux d'assumer le compilateur gcc et si je savais que je voulais cibler sse3 je pourrais suivre Comment activer l'autovectorisation sse3 dans gcc et écrire :

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

Cela produit en effet du code assembleur vectorisé rapide en utilisant gcc. Bien que vous ayez toujours besoin de rembourrer votre entrée à un multiple de 4. L'assemblage que vous obtenez est :

.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

Cependant, il est conçu pour le jeu d'instructions exact de simd et n'est pas optimal pour avx2 ou avx512 par exemple pour lesquels vous devez changer le code.

Comment pouvez-vous écrire du code C ou C++ pour lequel gcc produira du code optimal lorsqu'il est compilé pour l'un des systèmes sse, avx2 ou avx512 ? En d'autres termes, devez-vous toujours écrire des fonctions séparées à la main pour chaque largeur différente de registre SIMD ?

Existe-t-il des bibliothèques open source qui facilitent cette tâche ?

Answer 1

Voici un exemple utilisant le Bibliothèque Eigen :

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

Si vous compilez ceci avec clang ou g++ et sse ou avx activé (et -O2), vous devriez obtenir un code machine assez décent. Cela fonctionne également pour d'autres architectures comme Altivec ou NEON. Si vous savez que la première entrée de x est aligné, vous pouvez utiliser MapAligned au lieu de Map .

Vous obtiendrez un code encore meilleur, si vous connaissez la taille de votre vecteur au moment de la compilation :

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

Note : Les fonctions ci-dessus correspondent directement à la fonction f de l'OP. Cependant, comme @PeterCordes l'a souligné, il est généralement mauvais de stocker des nombres complexes entrelacés, car cela nécessitera beaucoup de brassage pour la multiplication. Au lieu de cela, on devrait stocker les parties réelles et imaginaires de manière à ce qu'elles puissent être directement chargées d'un paquet à la fois.

Édition/Addendum : Pour implémenter une structure de tableaux comme la multiplication complexe, vous pouvez en fait écrire quelque chose comme :

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

Ou plus générique (en utilisant 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;
}

Lorsqu'il est compilé avec -mavx -O2 la compilation donne quelque chose comme ceci (avec 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

Pour des raisons qui ne me sont pas évidentes, ceci est en fait caché dans une méthode qui est appelée par la méthode actuelle, qui ne fait que déplacer de la mémoire -- je ne sais pas pourquoi Eigen/gcc ne suppose pas que les arguments sont déjà correctement alignés. Si je compile la même chose avec clang 3.8.0 (et les mêmes arguments), elle est compilée en juste :

    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

Encore une fois, le déplacement de la mémoire au début est bizarre, mais au moins il est vectorisé. Pour gcc et clang, cela est optimisé lorsqu'il est appelé dans une boucle, cependant :

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 différence ici est que clang va dérouler cette boucle extérieure à 2 multiplications par boucle. D'un autre côté, gcc utilisera des instructions de type fused-multiply-add lorsqu'il est compilé avec le paramètre -mfma .

Le site f8 peut bien sûr aussi être généralisée à des dimensions arbitraires :

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

Et pour réduire la complexX<N> à un seul std::complex la fonction suivante peut être utilisée :

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

Cependant, selon que j'utilise clang ou g++, j'obtiens une sortie assembleur très différente. En général, g++ a tendance à ne pas charger en ligne les arguments d'entrée, et clang n'utilise pas les opérations FMA (YMMV ...). Essentiellement, vous devez inspecter le code assembleur généré de toute façon. Et plus important encore, vous devriez évaluer le code (je ne suis pas sûr de l'impact de cette routine sur votre problème global).

Aussi, je voulais noter que Eigen est en fait une bibliothèque d'algèbre linéaire. L'exploiter pour la génération de code SIMD purement portable n'est pas vraiment ce pour quoi elle a été conçue.

Answer 2

Si la portabilité est votre principale préoccupation, il existe de nombreuses bibliothèques ici qui fournissent des instructions SIMD dans leur propre syntaxe. La plupart d'entre eux font la vectorisation explicite plus simple et portable que les intrinsèques. Cette bibliothèque (UME::SIMD) a été publié récemment et est très performant

Sur ce document(UME::SIMD) une interface basée sur Vc a été établi qui s'appelle UME::SIMD. Il permet au programmeur d'accéder aux capacités SIMD sans avoir besoin d'une connaissance approfondie des ISA SIMD. UME::SIMD fournit une abstraction simple, flexible et portable pour la vectorisation explicite. vectorisation explicite sans perte de performance par rapport aux intrinsèques

Answer 3

Je ne pense pas que vous ayez une solution générale pour cela. Vous pouvez augmenter votre "vector_size" à 32 :

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

Augmentez également tous les tableaux pour qu'ils aient 8 éléments :

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

Cela rendra le compilateur capable de générer du code AVX512 (n'oubliez pas d'ajouter -mavx512f ), mais rendra votre code légèrement plus mauvais en SSE en rendant les transferts de mémoire sous-optimaux. Cependant, il ne désactivera certainement pas la vectorisation SSE.

Vous pourriez conserver les deux versions (avec 4 et 8 éléments de tableau), en passant de l'une à l'autre par un drapeau, mais cela pourrait être trop fastidieux pour un bénéfice limité.