Complexe Mul et Div utilisant sse Instructions

Question

La multiplication et la division complexes sont-elles bénéfiques grâce aux instructions SSE? Je sais que l'addition et la soustraction fonctionnent mieux lorsque vous utilisez SSE. Quelqu'un peut-il me dire comment je peux utiliser SSE pour effectuer une multiplication complexe pour obtenir de meilleures performances?

Answer 1

bien multiplication complexe est défini comme:

((c1a * c2a) - (c1b * c2b)) + ((c1b * c2a) + (c1a * c2b))i 

Ainsi, vos 2 composants dans un nombre complexe seraient

((c1a * c2a) - (c1b * c2b)) and ((c1b * c2a) + (c1a * c2b))i 

Donc, en supposant que vous utilisez 8 flotteurs pour représenter 4 nombres complexes définis comme suit :

c1a, c1b, c2a, c2b 
c3a, c3b, c4a, c4b 

Et vous voulez faire simultanément (c1 * c3) et (c2 * c4) votre code SSE ressemblerait à "quelque chose" comme le suivant:

(Note J'ai utilisé MSVC sous windows mais le principe sera le même).

__declspec(align(16)) float c1c2[]  = { 1.0f, 2.0f, 3.0f, 4.0f }; 
__declspec(align(16)) float c3c4[]   = { 4.0f, 3.0f, 2.0f, 1.0f }; 
__declspec(align(16)) float mulfactors[] = { -1.0f, 1.0f, -1.0f, 1.0f }; 
__declspec(align(16)) float res[]   = { 0.0f, 0.0f, 0.0f, 0.0f }; 

__asm 
{ 
    movaps xmm0, xmmword ptr [c1c2]   // Load c1 and c2 into xmm0. 
    movaps xmm1, xmmword ptr [c3c4]   // Load c3 and c4 into xmm1. 
    movaps xmm4, xmmword ptr [mulfactors] // load multiplication factors into xmm4 

    movaps xmm2, xmm1      
    movaps xmm3, xmm0      
    shufps xmm2, xmm1, 0xA0     // Change order to c3a c3a c4a c4a and store in xmm2 
    shufps xmm1, xmm1, 0xF5     // Change order to c3b c3b c4b c4b and store in xmm1 
    shufps xmm3, xmm0, 0xB1     // change order to c1b c1a c2b c2a abd store in xmm3 

    mulps xmm0, xmm2       
    mulps xmm3, xmm1      
    mulps xmm3, xmm4      // Flip the signs of the 'a's so the add works correctly. 

    addps xmm0, xmm3      // Add together 

    movaps xmmword ptr [res], xmm0   // Store back out 
}; 

float res1a = (c1c2[0] * c3c4[0]) - (c1c2[1] * c3c4[1]); 
float res1b = (c1c2[1] * c3c4[0]) + (c1c2[0] * c3c4[1]); 

float res2a = (c1c2[2] * c3c4[2]) - (c1c2[3] * c3c4[3]); 
float res2b = (c1c2[3] * c3c4[2]) + (c1c2[2] * c3c4[3]); 

if (res1a != res[0] || 
    res1b != res[1] || 
    res2a != res[2] || 
    res2b != res[3]) 
{ 
    _exit(1); 
} 

Ce que j'ai fait ci-dessus, c'est que j'ai simplifié un peu les calculs. En supposant que les éléments suivants:

c1a c1b c2a c2b 
c3a c3b c4a c4b 

En réorganisant je finis avec les vecteurs suivants

0 => c1a c1b c2a c2b 
1 => c3b c3b c4b c4b 
2 => c3a c3a c4a c4a 
3 => c1b c1a c2b c2a 

Je puis multiplier 0 et 2 ensemble pour obtenir:

0 => c1a * c3a, c1b * c3a, c2a * c4a, c2b * c4a 

Ensuite, je multiplierai 3 et 1 ensemble pour obtenir:

Enfin je feuillette les signes d'un couple des flotteurs dans 3

3 => -(c1b * c3b), c1a * c3b, -(c2b * c4b), c2a * c4b 

Je peux les ajouter ensemble et obtenir

(c1a * c3a) - (c1b * c3b), (c1b * c3a) + (c1a * c3b), (c2a * c4a) - (c2b * c4b), (c2b * c4a) + (c2a * c4b) 

Ce qui est ce que nous recherchions :)

Answer 2

Juste pour l'exhaustivité, le Manuel de référence de l'optimisation des architectures Intel® 64 et IA-32 téléchargeable here contient un assemblage pour la division complexe (exemple 6-9) et complexe (exemple 6-10).

est ici par exemple le code se multiplient:

// Multiplication of (ak + i bk) * (ck + i dk) 
// a + i b can be stored as a data structure 
movsldup xmm0, src1; load real parts into the destination, a1, a1, a0, a0 
movaps xmm1, src2; load the 2nd pair of complex values, i.e. d1, c1, d0, c0 
mulps xmm0, xmm1; temporary results, a1d1, a1c1, a0d0, a0c0 
shufps xmm1, xmm1, b1; reorder the real and imaginary parts, c1, d1, c0, d0 
movshdup xmm2, src1; load imaginary parts into the destination, b1, b1, b0, b0 
mulps xmm2, xmm1; temporary results, b1c1, b1d1, b0c0, b0d0 
addsubps xmm0, xmm2; b1c1+a1d1, a1c1 -b1d1, b0c0+a0d0, ; a0c0-b0d0 

L'ensemble des cartes directement à gccs X86 intrinsics (juste prédicat chaque instruction avec __builtin_ia32_).

Answer 3

L'algorithme de la référence d'optimisation Intel ne gère pas correctement les débordements et NaN dans l'entrée.

Un seul NaN dans la partie réelle ou imaginaire du numéro se propagera incorrectement à l'autre partie.

Étant donné que plusieurs opérations avec l'infini (par exemple, l'infini * 0) se terminent par NaN, les débordements peuvent provoquer l'apparition de NaN s dans vos données par ailleurs saines.

Si déversoirs et NaN s sont rares, d'une manière simple d'éviter cela est de vérifier juste pour NaN dans le résultat et recalcule avec les compilateurs mise en œuvre conforme IEEE:

float complex a[2], b[2]; 
__m128 res = simd_fast_multiply(a, b); 

/* store unconditionally, can be executed in parallel with the check 
* making it almost free if there is no NaN in data */ 
_mm_store_ps(dest, res); 

/* check for NaN */ 
__m128 n = _mm_cmpneq_ps(res, res); 
int have_nan = _mm_movemask_ps(n); 
if (have_nan != 0) { 
    /* do it again unvectorized */ 
    dest[0] = a[0] * b[0]; 
    dest[1] = a[1] * b[1]; 
}