performance - Rsqrt vectorizado rápido y recíproco con SSE/AVX dependiendo de la precisión

simd

Question

Supongamos que es necesario calcular la raíz cuadrada recíproca o recíproca para los datos de punto flotante empaquetados. Ambos se pueden hacer fácilmente por:

__m128 recip_float4_ieee(__m128 x) { return _mm_div_ps(_mm_set1_ps(1.0f), x); }
__m128 rsqrt_float4_ieee(__m128 x) { return _mm_div_ps(_mm_set1_ps(1.0f), _mm_sqrt_ps(x)); }

Esto funciona perfectamente pero lento: de acuerdo con la guía , toman 14 y 28 ciclos en Sandy Bridge (rendimiento). Las versiones correspondientes de AVX toman casi el mismo tiempo en Haswell.

Por otro lado, las siguientes versiones se pueden usar en su lugar:

__m128 recip_float4_half(__m128 x) { return _mm_rcp_ps(x); }
__m128 rsqrt_float4_half(__m128 x) { return _mm_rsqrt_ps(x); }

Toman solo uno o dos ciclos de tiempo (rendimiento), dando un impulso de rendimiento mayor. Sin embargo, son MUY aproximadas: producen un resultado con un error relativo menor a 1.5 * 2 ^ -12. Dado que épsilon de máquina de números de flotador de precisión simple es 2 ^ -24, podemos decir que esta aproximación tiene aproximadamente la mitad de precisión.

Parece que la iteración de Newton-Raphson se puede agregar para producir un resultado con precisión simple (tal vez no sea tan exacta como lo exige el estándar IEEE), ver GCC , ICC , discusiones en LLVM . Teóricamente, el mismo método se puede usar para valores de doble precisión, produciendo la mitad o una precisión simple o doble .

Me interesan las implementaciones de trabajo de este enfoque para los tipos de datos flotantes y dobles y para todas las precisiones (media, simple, doble). No es necesario manejar casos especiales (división por cero, sqrt (-1), inf / nan y similares). Además, no está claro para mí cuál de estas rutinas sería más rápida que las soluciones triviales IEEE-Compilant, y que sería más lenta.

Aquí hay algunas restricciones menores en las respuestas, por favor:

1. Use elementos intrínsecos en sus muestras de código. El ensamblado depende del compilador, por lo que es menos útil.
2. Use una convención de nomenclatura similar para las funciones.
3. Implemente rutinas teniendo como único registro SSE / AVX único que contenga valores flotantes / dobles densamente empaquetados. Si hay una mejora considerable en el rendimiento, también puede publicar rutinas teniendo varios registros como entrada (dos regs pueden ser viables).
4. No publique ambas versiones SSE / AVX si son absolutamente iguales hasta que cambie _mm por _mm256 y viceversa.

Cualquier estimación de rendimiento, medidas, discusiones son bienvenidas.

RESUMEN

Aquí están las versiones para números flotantes de precisión simple con una iteración NR:

__m128 recip_float4_single(__m128 x) {
  __m128 res = _mm_rcp_ps(x);
  __m128 muls = _mm_mul_ps(x, _mm_mul_ps(res, res));
  return res =  _mm_sub_ps(_mm_add_ps(res, res), muls);
}
__m128 rsqrt_float4_single(__m128 x) {
  __m128 three = _mm_set1_ps(3.0f), half = _mm_set1_ps(0.5f);
  __m128 res = _mm_rsqrt_ps(x); 
  __m128 muls = _mm_mul_ps(_mm_mul_ps(x, res), res); 
  return res = _mm_mul_ps(_mm_mul_ps(half, res), _mm_sub_ps(three, muls)); 
}

La respuesta dada por Peter Cordes explica cómo crear otras versiones, y contiene un análisis de rendimiento teórico completo.

Puede encontrar todas las soluciones implementadas con benchmark aquí: recip_rsqrt_benchmark .

Los resultados de rendimiento obtenidos en Ivy Bridge se presentan a continuación. Solo las implementaciones SSE de registro único han sido evaluadas. El tiempo dedicado se da en ciclos por llamada. El primer número es para la mitad de la precisión (sin NR), el segundo es para la precisión simple (1 iteración de NR), el tercero es para las iteraciones de 2 NR.

1. recip en flotador toma 1, 4 ciclos versus 7 ciclos.
2. rsqrt en flotación toma 1, 6 ciclos versus 14 ciclos.
3. recip en doble toma 3, 6, 9 ciclos versus 14 ciclos.
4. rsqrt en doble toma 3, 8, 13 ciclos contra 28 ciclos.

Advertencia: tuve que redondear los resultados crudos de forma creativa ...

Answer 1

Hay muchos ejemplos del algoritmo en la práctica. Por ejemplo:

Newton Raphson con SSE2: alguien puede explicarme que estas 3 líneas tienen una respuesta que explica la iteración utilizada por uno de los ejemplos de Intel .

Para el análisis de rendimiento en, digamos, Haswell (que puede FP mul en dos puertos de ejecución, a diferencia de diseños anteriores), reproduciré el código aquí (con una operación por línea). Consulte http://agner.org/optimize/ para obtener información sobre el rendimiento y la latencia de las tablas, y sobre cómo comprender más antecedentes.

// execute (aka dispatch) on cycle 1, results ready on cycle 6
nr = _mm_rsqrt_ps( x );

// both of these execute on cycle 6, results ready on cycle 11
xnr = _mm_mul_ps( x, nr );         // dep on nr
half_nr = _mm_mul_ps( half, nr );  // dep on nr

// can execute on cycle 11, result ready on cycle 16
muls = _mm_mul_ps( xnr , nr );     // dep on xnr

// can execute on cycle 16, result ready on cycle 19
three_minus_muls = _mm_sub_ps( three, muls );  // dep on muls

// can execute on cycle 19, result ready on cycle 24
result = _mm_mul_ps( half_nr, three_minus_muls ); // dep on three_minus_muls
// result is an approximation of 1/sqrt(x), with ~22 to 23 bits of precision in the mantissa.

Hay MUCHO espacio para que otros cálculos se superpongan aquí, si no es parte de la cadena de dependencia. Sin embargo, si los datos para cada iteración de su código dependen de los datos de la anterior, es posible que esté mejor con la latencia de 11 ciclos de sqrtps . O incluso si cada iteración de bucle es lo suficientemente larga como para que la ejecución fuera de orden no pueda ocultarlo solapando iteraciones independientes.

Para obtener sqrt(x) lugar de 1/sqrt(x) , multiplica por x (y fixup si x puede ser cero, por ejemplo, enmascarando ( _mm_andn_ps ) con el resultado de un CMPPS contra 0.0). La forma óptima es reemplazar half_nr con half_xnr = _mm_mul_ps( half, xnr ); . Eso no alarga la cadena de dep, porque half_xnr puede comenzar en el ciclo 11 pero no es necesario hasta el final (ciclo 19). Lo mismo con FMA disponible: sin aumento en la latencia total.

Si 11 bits de precisión son suficientes (sin iteración de Newton), el manual de optimización de Intel (sección 11.12.3) sugiere usar rcpps (rsqrt (x)), que es peor que multiplicar por la x original, al menos con AVX. Puede guardar una instrucción movdqa con SSE de 128 bits, pero 256b rcpps es más lenta que 256b mul o fma. (Y le permite agregar el resultado sqrt a algo gratis con FMA en lugar de mul para el último paso).

La versión SSE de este ciclo, sin tener en cuenta ninguna instrucción de movimiento, es de 6 uops. Esto significa que debe tener un rendimiento en Haswell de uno por 3 ciclos (dado que dos puertos de ejecución pueden manejar FP mul, y rsqrt está en el puerto opuesto al FP add / sub). En SnB / IvB (y probablemente en Nehalem), debería tener un rendimiento de uno por 5 ciclos , ya que los mulps y rsqrtps compiten por el puerto 0. los subps están en el puerto 1, que no es el cuello de botella.

Para Haswell, podemos usar FMA para combinar el restar con el mul. Sin embargo, la unidad divisores / sqrt no tiene 256 b de ancho, por lo que a diferencia de todo lo demás, divps / sqrtps / rsqrtps / rcpps en ymm regs toma uops extra y ciclos extra.

// vrsqrtps ymm has higher latency
// execute on cycle 1, results ready on cycle 8
nr = _mm256_rsqrt_ps( x );

// both of can execute on cycle 8, results ready on cycle 13
xnr = _mm256_mul_ps( x, nr );         // dep on nr
half_nr = _mm256_mul_ps( half, nr );  // dep on nr

// can execute on cycle 13, result ready on cycle 18
three_minus_muls = _mm256_fnmadd_ps( xnr, nr, three );  // -(xnr*nr) + 3

// can execute on cycle 18, result ready on cycle 23
result = _mm256_mul_ps( half_nr, three_minus_muls ); // dep on three_minus_muls

Ahorramos 3 ciclos con FMA. Esto se compensa utilizando un rsqrt 256b de 2 ciclos más lento, para una ganancia neta de 1 ciclo menos de latencia (bastante bueno para el doble de ancho). SnB / IvB AVX sería la latencia de 24c para 128b, 26c de latencia para 256b.

La versión 256b FMA usa 7 uops en total. ( VRSQRTPS es 3 uops, 2 para p0, y uno para p1 / 5). 256b mulps y fma son instrucciones single-uop, y ambos pueden ejecutarse en el puerto 0 o en el puerto 1. (p0 solo en pre-Haswell). Por lo tanto, el rendimiento debería ser: uno por 3c , si el motor OOO despacha uops a los puertos de ejecución óptimos. (es decir, el shuffle uop de rsqrt siempre va a p5, nunca a p1, donde ocuparía el ancho de banda mul / fma). En cuanto a la superposición con otros cálculos (no solo la ejecución independiente), de nuevo es bastante liviano. Solo 7 uops con una cadena de depósito de 23 ciclos dejan mucho espacio para que sucedan otras cosas mientras esos uops se sientan en el búfer de pedido.

Si este es un paso en una cadena de depósito gigante sin nada más (ni siquiera una próxima iteración independiente), entonces 256b vsqrtps es 19 latencia de ciclo, con un rendimiento de uno por 14 ciclos. (Haswell). Si todavía necesita el recíproco, 256b vdivps también tiene una latencia de 18-21c, con una por cada 14c de rendimiento. Entonces, para sqrt normal, la instrucción tiene una latencia menor. Para sqrt rec, una iteración de la aproximación es menos latencia. (Y MUCHO más rendimiento, si se superpone consigo mismo. Si se superpone con otras cosas que no hacen la unidad de división, sqrtps no es un problema).

sqrtps puede ser una iteración de rendimiento versus iteración de rsqrt + Newton, si es parte de un cuerpo de bucle con la suficiente cantidad de trabajo no dividido y no cuadrado que la unidad de división no esté saturada.

Esto es especialmente cierto si necesita sqrt(x) , no 1/sqrt(x) . por ejemplo, en Haswell con AVX2, un bucle copia + arcsinh sobre una matriz de flotantes que se adapta a caché L3, implementado como fastlog(v + sqrt(v*v + 1)) , se ejecuta aproximadamente al mismo rendimiento con VSQRTPS reales o con VRSQRTPS + una iteración Newton-Raphson. (Incluso con una aproximación muy rápida para log () , por lo que el cuerpo del bucle total es de aproximadamente 9 operaciones FMA / add / mul / convertir, y 2 booleanos, más el VSQRTPS ymm. Hay una aceleración al usar simplemente fastlog(v2_plus_1 * rsqrt(v2_plus_1) + v2_plus_1) , por lo que no tiene un cuello de botella en el ancho de banda de la memoria, pero puede ser un cuello de botella en la latencia (por lo que la ejecución fuera de orden no puede explotar todo el paralelismo de las iteraciones independientes).

Otras precisiones

Para la precisión media , no hay instrucciones para hacer cálculos en medio flotadores. Debe convertir sobre la marcha mientras carga / almacena, usando las instrucciones de conversión.

Para precisión doble , no hay rsqrtpd . Presumiblemente, el pensamiento es que si no necesitas precisión total, deberías estar usando flotador en primer lugar. Entonces puedes convertir para flotar y volver, luego hacer el mismo algoritmo exacto, pero con pd lugar de ps intrinsics. O bien, podría mantener sus datos flotando por un tiempo. por ejemplo, convertir dos registros ymm de dobles en un registro ymm de singles.

Desafortunadamente, no hay una sola instrucción que tome dos registros de ymm de dobles y salidas de ymm de singles. Tienes que ir ymm-> xmm dos veces, luego _mm256_insertf128_ps uno xmm al alto 128 del otro. Pero luego puedes alimentar ese vector 256b ymm al mismo algoritmo.

Si va a convertir de nuevo a doble inmediatamente después, puede tener sentido realizar la iteración rsqrt + Newton-Raphson en los dos registros 128b de singles por separado. Los 2 uops extra para insertar / extraer, y los 2 uops adicionales para rsqrt 256b, comienzan a sumar, sin mencionar la latencia de 3 ciclos de vinsertf128 / vextractf128 . La computación se superpondrá y ambos resultados estarán listos con un par de ciclos de diferencia. (O 1 ciclo aparte, si tiene una versión especial de su código para operaciones de entrelazado en 2 entradas a la vez).

Recuerde que la precisión simple tiene un rango menor de exponentes que el doble, por lo que la conversión puede desbordarse a + Inf o subdesbordamiento a 0.0. Si no está seguro, definitivamente solo use el _mm_sqrt_pd normal.

Con AVX512F, hay _mm512_rsqrt14_pd( __m512d a) . Con AVX512ER (KNL pero no SKX o Cannonlake) , hay _mm512_rsqrt28_pd . _ps versiones _ps por supuesto también existen. 14 bits de precisión de mantisa pueden ser suficientes para usar sin una iteración de Newton en más casos.

Una iteración de Newton todavía no le dará un flotador de precisión simple redondeado como lo hará el sqrt regular.