Complex Mul and Div using sse Instructions

Question

¿Es beneficioso realizar multiplicaciones y divisiones complejas a través de las instrucciones de SSE? Sé que la suma y la resta funcionan mejor cuando se usa SSE. ¿Puede alguien decirme cómo puedo usar SSE para realizar una multiplicación compleja para obtener un mejor rendimiento?

Answer 1

multiplicación Bueno complejo se define como:

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

Así que sus 2 componentes en un número complejo serían

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

Así que asumiendo que usted está utilizando 8 flotadores para representar 4 números complejos definidos de la siguiente manera :

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

Y quiere hacer simultáneamente (c1 * c3) y (c2 * c4) su código SSE se vería "algo" como el siguiente:

(Nota: utilicé MSVC en Windows pero el principio SERÁ el mismo).

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

Lo que he hecho anteriormente es que he simplificado las matemáticas un poco. Suponiendo que la siguiente:

c1a c1b c2a c2b 
c3a c3b c4a c4b 

Reordenando termino con los siguientes vectores

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

I luego se multiplica 0 y 2 en conjunto para obtener:

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

Siguiente multiplico 3 y 1 juntos para obtener:

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

Finalmente me tapa los signos de un par de los flotadores en 3

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

Así que puedo sumarlos y obtener

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

que es lo que buscábamos :)

Answer 2

Sólo por integridad, el Manual de referencia de optimización de arquitecturas Intel® 64 e IA-32 que se puede descargar here contiene el ensamblaje para la multiplicación compleja (Ejemplo 6-9) y la división compleja (Ejemplo 6-10).

Así es, por ejemplo, el código se multiplican:

// 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 

El montaje se correlaciona directamente con gccs X86 intrinsics (solo predicado cada instrucción con __builtin_ia32_).

Answer 3

El algoritmo en la referencia de optimización Intel no maneja los desbordamientos y NaN s en la entrada correctamente.

Una sola NaN en la parte real o imaginaria del número se propagará incorrectamente a la otra parte.

Como varias operaciones con infinito (por ejemplo, infinito * 0) terminan en NaN, los desbordamientos pueden hacer que NaN s aparezcan en sus datos que por lo demás se comportan bien.

Si se desborda y NaN s son raros, una forma sencilla de evitar esto es simplemente comprobar si hay NaN en el resultado y volver a calcular con los compiladores IEEE aplicación compatible:

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