Cómo optimizar "U [0] * v [0] + u [2] * v [2]" línea de código con SSE o GLSL

Question

tengo la siguiente función (del proyecto de código abierto "recast navigation"):

/// Derives the dot product of two vectors on the xz-plane. (@p u . @p v) 
/// @param[in]  u  A vector [(x, y, z)] 
/// @param[in]  v  A vector [(x, y, z)] 
/// @return The dot product on the xz-plane. 
/// 
/// The vectors are projected onto the xz-plane, so the y-values are ignored. 
inline float dtVdot2D(const float* u, const float* v) 
{ 
    return u[0]*v[0] + u[2]*v[2]; 
} 

Lo ejecuté a través de la prueba de rendimiento de la CPU VS2010 y me mostró que en toda la línea de código de refundición de código en esta función u[0]*v[0] + u[2]*v[2] es la CPU más caliente.

¿Cómo puedo optimizar la CPU (a través de SSE o GLSL like GLM (if it is easier or faster and appropriate in such case) por ejemplo) esta línea?

Editar: El contexto en el que las llamadas aparecen:

bool dtClosestHeightPointTriangle(const float* p, const float* a, const float* b, const float* c, float& h) { 
    float v0[3], v1[3], v2[3]; 
    dtVsub(v0, c,a); 
    dtVsub(v1, b,a); 
    dtVsub(v2, p,a); 

    const float dot00 = dtVdot2D(v0, v0); 
    const float dot01 = dtVdot2D(v0, v1); 
    const float dot02 = dtVdot2D(v0, v2); 
    const float dot11 = dtVdot2D(v1, v1); 
    const float dot12 = dtVdot2D(v1, v2); 

    // Compute barycentric coordinates 
    const float invDenom = 1.0f/(dot00 * dot11 - dot01 * dot01); 
    const float u = (dot11 * dot02 - dot01 * dot12) * invDenom; 
    const float v = (dot00 * dot12 - dot01 * dot02) * invDenom; 

Answer 1

Después de intentar algunas cosas sobre el papel, he encontrado algo que podría funcionar. Es una implementación adecuadamente paralelizada/vectorizada de la función en SSE.

Sin embargo, requiere reorganización de datos porque haremos un procesamiento paralelo para 4 triángulos a la vez.

Lo dividiré en pasos y usaré los nombres de las instrucciones aquí y allá, pero use los intrínsecos C (_mm_load_ps(), _mm_sub_ps() et al, están en xmmintrin.h en VC) - cuando hablar de registros esto solo significa __m128.

PASO 1.

No necesitamos la coordenada Y en absoluto, por lo que hemos creado punteros a pares de X y Z. de suministro al menos 4 pares (es decir, 4 triángulos total) por llamada. Llamaré a cada par X y Z un vértice.

PASO 2.

Uso MOVAPS (requiere los punteros a estar alineados a 16 bits) para cargar los dos primeros vértices apuntado por cada puntero en registros.

un registro cargado desde un se verá así: [a0.x, a0.z, a1.x, a1.z]

PASO 3.

Ahora, utilizando una sola instrucción de resta, se puede calcular los deltas (su v0, v1 , v2) durante 2 vértices a la vez.

Calcular v0, v1 y v2 no sólo para los 2 primeros triángulos, sino también para el último 2! Como dije, debe proporcionar un total de 4 vértices u 8 flotantes por entrada. Solo repita los pasos 2 y 3 para los datos.

Ahora tenemos 2 pares de vx registros, cada par que contiene el resultado para 2 triángulos. Me referiré a ellos como vx_0 (primer par) y vx_1 (segundo par) donde x es de 0 a 2.

PASO 4.

productos escalares. Para paralelizar el cálculo baricéntrico (más adelante), requerimos el resultado de cada producto de puntos para cada uno de los 4 triángulos, en 1 registro individual.

Entonces, donde calcularía dot01 por ejemplo, haremos lo mismo, pero para 4 triángulos a la vez. Cada v -registro contiene el resultado para 2 vectores, así que comenzamos por multiplicarlos.

Digamos que u y v - parámetros en función de su producto escalar - son ahora v0_0 y v1_0 (para calcular dot01):

Multiplicar u y v obtener: [(v0_0.x0) * (v1_0.x0), (v0_0.z0) * (v1_0.z0), (v0_0.x1) * (v1_0.x1), (v0_0.z1) * (v1_0.z1)]

Esto puede parecer confuso, debido a .x0/.z0 y .x1/.z1, pero mira lo que estaba cargado en el paso 2: a0, a1 .

Si por ahora esto todavía se siente difusa, recoger un pedazo de papel y escriba lo largo desde el principio.

A continuación, aún para el mismo producto escalar, hacer la multiplicación para v0_1 y v1_1 (es decir, el segundo par de triángulos).

Ahora tenemos 2 registros con 2 X & Z pares cada uno (4 vértices en total), multiplicados y listos para sumarse para formar 4 productos de puntos separados.SSE3 tiene una instrucción para hacer exactamente esto, y se llama HADDPS:

xmm0 = [A, B, C, D] XMM1 = [E, F, G, H]

HADDPS xmm0, XMM1 hace esto:

xmm0 = [A + B, C + D, E + F, G + H]

se llevará a las X & pares Z de nuestro primer registro, los de la segunda, sumarlos y almacenarlos en el primer, segundo, tercer y cuarto componente del registro de destino. Ergo: ¡en este punto tenemos este producto en particular para los 4 triángulos!

** Ahora repita este proceso para todos los productos dot: dot00 et cetera. **

PASO 5.

El último cálculo (por lo que pude ver por el código suministrado) es la materia barycentric. Este es un cálculo 100% escalar en su código. Pero sus entradas ahora no son resultados escalares de productos de puntos (es decir flotadores individuales), son vectores/registros SSE con un producto de puntos para cada uno de los 4 triángulos.

Así que si vas a vectorizar esto usando las operaciones SSE paralelas que operan en los 4 flotadores, eventualmente terminarás con 1 registro (o resultado) llevando 4 alturas, 1 para cada triángulo.

Como mi hora del almuerzo ya ha pasado, no voy a deletrear esto, pero teniendo en cuenta la configuración/idea que he dado, este es un último paso y no debería ser difícil de entender.

Sé que esta idea es un poco exagerada y requiere un poco de amor del código que se encuentra arriba y tal vez algún tiempo de calidad con lápiz y papel, pero será rápido (e incluso puede agregar OpenMP después si te gustaría).

Buena suerte :)

(y perdona mi explicación difusa, que puede azotar encima de la función si es necesario =))

ACTUALIZACIÓN

he escrito una implementación y no funcionó como esperaba, principalmente porque el componente Y se involucró más allá del fragmento de código que inicialmente pegó en su pregunta (lo busqué). Lo que he hecho aquí no es solo pedirte que reorganices todos los puntos en pares XZ y los alimentes por 4, sino que también alimentes 3 punteros (para los puntos A, B y C) con los valores Y de cada uno de los 4 triángulos. Desde una perspectiva local, esto es más rápido. Todavía puedo modificarlo para que requiera cambios menos intrusivos del final del llamado, pero háganme saber lo que es deseable.

A continuación, un descargo de responsabilidad: este código es directo como el infierno (algo que he encontrado que funciona bastante bien con los compiladores en términos de SSE ... pueden reorganizarse como corresponde y las CPU x86/x64 toman su parte también) También el nombre, no es mi estilo, no es de nadie, solo hazlo como mejor te parezca.

Espero que ayude y si no, con mucho gusto lo repasaremos de nuevo.Y si esto es un proyecto comercial también existe la opción de conseguir mí a bordo supongo;)

De todos modos, lo he puesto en Pastebin: http://pastebin.com/20u8fMEb

Answer 2

puede implementar su producto escalar individual con instrucciones SSE, pero el resultado no será significativamente más rápido (e incluso puede ser más lento) que el código como escrito ahora. Su reescritura puede vencer las optimizaciones del compilador que están ayudando a la versión actual.

Para aprovechar cualquier beneficio de reescribir eso con SSE o CUDA, tiene que optimizar el bucle que está llamando para ese producto de punto. Esto es particularmente cierto para CUDA, donde la sobrecarga de la ejecución de un producto de punto sería enorme. Solo puede ver una aceleración si envía miles de vectores a la GPU para calcular miles de productos de puntos. La misma idea vale para SSE en la CPU, pero es posible que pueda ver una mejora en un número menor de operaciones. Sin embargo, seguirá siendo más de un producto de dot.

Lo más fácil de probar es g++ -ftree-vectorize. GCC podrá alinear su pequeña función y luego tratar de optimizar el circuito por usted (de hecho, probablemente ya lo esté, pero sin instrucciones de SSE). El vectorizador de árbol intentará hacer automáticamente lo que propone hacer a mano. No siempre es exitoso.

Answer 3

que solicitó la versión de SSE de su algoritmo de modo que aquí es:

// Copied and modified from xnamathvector.inl 
XMFINLINE XMVECTOR XMVector2DotXZ 
(
    FXMVECTOR V1, 
    FXMVECTOR V2 
) 
{ 
#if defined(_XM_NO_INTRINSICS_) 

    XMVECTOR Result; 

    Result.vector4_f32[0] = 
    Result.vector4_f32[1] = 
    Result.vector4_f32[2] = 
    Result.vector4_f32[3] = V1.vector4_f32[0] * V2.vector4_f32[0] + V1.vector4_f32[2] * V2.vector4_f32[2]; 

    return Result; 

#elif defined(_XM_SSE_INTRINSICS_) 
    // Perform the dot product on x and z 
    XMVECTOR vLengthSq = _mm_mul_ps(V1,V2); 
    // vTemp has z splatted 
    XMVECTOR vTemp = _mm_shuffle_ps(vLengthSq,vLengthSq,_MM_SHUFFLE(2,2,2,2)); 
    // x+z 
    vLengthSq = _mm_add_ss(vLengthSq,vTemp); 
    vLengthSq = _mm_shuffle_ps(vLengthSq,vLengthSq,_MM_SHUFFLE(0,0,0,0)); 
    return vLengthSq; 
#else // _XM_VMX128_INTRINSICS_ 
#endif // _XM_VMX128_INTRINSICS_ 
} 

bool dtClosestHeightPointTriangle(FXMVECTOR p, FXMVECTOR a, FXMVECTOR b, FXMVECTOR c, float& h) 
{ 
    XMVECTOR v0 = XMVectorSubtract(c,a); 
    XMVECTOR v1 = XMVectorSubtract(b,a); 
    XMVECTOR v2 = XMVectorSubtract(p,a); 

    XMVECTOR dot00 = XMVector2DotXZ(v0, v0); 
    XMVECTOR dot01 = XMVector2DotXZ(v0, v1); 
    XMVECTOR dot02 = XMVector2DotXZ(v0, v2); 
    XMVECTOR dot11 = XMVector2DotXZ(v1, v1); 
    XMVECTOR dot12 = XMVector2DotXZ(v1, v2); 

    // Compute barycentric coordinates 
    XMVECTOR invDenom = XMVectorDivide(XMVectorReplicate(1.0f), XMVectorSubtract(XMVectorMultiply(dot00, dot11), XMVectorMultiply(dot01, dot01))); 

    XMVECTOR u = XMVectorMultiply(XMVectorSubtract(XMVectorMultiply(dot11, dot02), XMVectorMultiply(dot01, dot12)), invDenom); 
    XMVECTOR v = XMVectorMultiply(XMVectorSubtract(XMVectorMultiply(dot00, dot12), XMVectorMultiply(dot01, dot02)), invDenom); 
} 

El XMVector2Dot se toma de xnamathvector.inl, que renombraron y modificar para funcionar en las coordenadas X/Z.

XNAMath es una gran biblioteca de matemáticas vectorizada multiplataforma de Microsoft; También lo uso en OS X importando el encabezado sal.h (no estoy seguro acerca del problema de licencia allí, así que ten cuidado).
De hecho, cualquier plataforma que admita los intrínsecos SSE debería ser compatible con.

Un par de cosas a tener en cuenta:

• Es necesario para cargar sus flotadores en XMVECTORs utilizando el método XMLoadFloat3; esto cargará float3 no alineado en una estructura __m128.
• Probablemente no verá ninguna mejora en el rendimiento de este código SSE (perfil !!) ya que existe una penalización de rendimiento para cargar flotantes alineados en los registros SSE.
• Esta es una conversión de fuerza bruta del algoritmo a SSE, tendrá mejor suerte al ser más inteligente que yo y tratar de entender el algoritmo e implementar una versión compatible con SSE.
• Tendrá mejor suerte al convertir toda la aplicación para usar el código XNA Math/SSE en lugar de solo esa pequeña porción. Al menos imponer el uso de tipos de vectores alineados (XMFLOAT3A o struct __declspec (align (16)) myvectortype {};).
• Se desaconseja el ensamblaje SSE directo, especialmente en x64, a favor de los intrínsecos.

Answer 4

Las instrucciones SSE están diseñadas para optimizar alghoritms que procesan grandes bloques de datos representados como enteros o números de coma flotante. Los tamaños típicos son millones y miles de millones de números que deben procesarse de alguna manera. No tiene sentido optimizar la función que procesa solo cuatro (o veinte) escalares. Lo que ganes con SSE podrías perder con la sobrecarga de llamadas de funciones. La cantidad razonable de números procesados ​​por una llamada de función es de al menos mil. Es posible que pueda obtener una gran ganancia de rendimiento mediante el uso de intrínsecos SSE. Pero es difícil darle consejos específicos adaptados a sus necesidades en función de la información que proporcionó. Debe editar su pregunta y proporcionar una vista de alto nivel de su problema (funciones ubicadas más profundamente en su pila de llamadas). Por ejemplo, no es obvio cuántas veces se llama el método dtClosestHeightPointTriangle por segundo. Este número es crítico para juzgar objetivamente si la transición a SSE sería de valor práctico. La organización de datos también es muy importante. Lo ideal es que sus datos se almacenen en la menor cantidad posible de segmentos lineales de memoria para utilizar el subsistema de caché de la CPU de manera eficiente.