Compl Mul e Div utilizzando le istruzioni sse

Question

Sta eseguendo la moltiplicazione e la divisione complesse con le istruzioni SSE? So che addizioni e sottrazioni funzionano meglio quando si usa SSE. Qualcuno può dirmi come posso usare SSE per eseguire una moltiplicazione complessa per ottenere prestazioni migliori?

Answer 1

Bene complesso moltiplicazione è definita come:

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

In modo che le 2 componenti in un numero complesso sarebbero

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

Quindi, supponendo che si sta utilizzando 8 carri per rappresentare 4 numeri complessi definiti come segue :

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

E si desidera eseguire simultaneamente (c1 * c3) e (c2 * c4) il tuo codice SSE sembrerebbe "qualcosa" come il seguente:

(Nota ho usato MSVC sotto Windows ma il principio sarà lo stesso).

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

Quello che ho fatto sopra è che ho semplificato un po 'la matematica. Assumendo il seguente:

c1a c1b c2a c2b 
c3a c3b c4a c4b 

Riordinando alla fine con i seguenti vettori

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

ho quindi moltiplicare 0 e 2 insieme per ottenere:

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

Successivo Moltiplico 3 e 1 insieme per ottenere:

3 => c1b * c3b, c1a * c3b, c2b * c4b, c2a * c4b 

Infine ho capovolgere i segni di un paio di carri in 3

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

modo che io possa aggiungere insieme e

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

che è quello che cercavamo :)

Answer 2

Solo per completezza, il Manuale di riferimento sull'ottimizzazione delle architetture Intel® 64 e IA-32 che può essere scaricato here contiene l'assembly per il multiplo complesso (Esempio 6-9) e la divisione complessa (Esempio 6-10).

Ecco per esempio il codice si moltiplicano:

// 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'assemblea mappe direttamente a gccs X86 intrinsics (solo predicato ogni istruzione con __builtin_ia32_).

Answer 3

L'algoritmo nel riferimento di ottimizzazione Intel non gestisce gli overflow e NaN s nell'input correttamente.

Un singolo NaN nella parte reale o immaginaria del numero non verrà distribuito correttamente all'altra parte.

Poiché diverse operazioni con infinito (ad esempio infinito * 0) terminano con NaN, gli overflow possono causare l'apparire di NaN nei dati altrimenti ben funzionanti.

Se overflow e NaN s sono rari, un modo semplice per evitare questo è quello di controllare solo per NaN nel risultato e ricalcolare con i compilatori IEEE compatibile attuazione:

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