Come ottimizzare "u [0] * v [0] + u [2] * v [2]" linea di codice con SSE o GLSL

Question

Ho la seguente funzione (dal progetto opensource "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]; 
} 

L'ho eseguito attraverso il test delle prestazioni della CPU VS2010 e mi ha mostrato che in tutta la codeline di codebase di rifusione in questa funzione u[0]*v[0] + u[2]*v[2] è la CPU più calda.

Come posso ottimizzare la CPU (tramite SSE o GLSL like GLM (if it is easier or faster and appropriate in such case) per esempio) su questa linea?

Edit: Il contesto in cui le chiamate vengono visualizzate:

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

Dopo aver provato un po 'di cose sulla carta, ho trovato qualcosa che potrebbe funzionare per voi. È un'implementazione perfettamente parallelizzata/vettoriale della funzione in SSE.

Tuttavia richiede la riorganizzazione dei dati perché eseguiremo l'elaborazione parallela per 4 triangoli contemporaneamente.

Scomporlo in passaggi e utilizzare i nomi delle istruzioni qua e là, ma si prega di utilizzare i C intrinsechi (_mm_load_ps(), _mm_sub_ps() et al, sono in xmmintrin.h in VC) - quando io parlare di registri significa solo __m128.

FASE 1.

Non abbiamo bisogno della coordinata Y del tutto, così abbiamo istituito puntatori a coppie di X e Z. fornire almeno 4 paia (cioè 4 triangoli Totale) per chiamata. Chiamerò ogni coppia X e Z un vertice.

FASE 2.

Uso MOVAPS (richiede i puntatori da allineare a 16-bit) per caricare i primi due vertici puntati da ogni puntatore in registri.

Un registro caricato da un sarà simile a questa: [a0.x, a0.z, a1.x, a1.z]

FASE 3.

Ora, utilizzando una singola istruzione sottrarre, è possibile calcolare i delta (il vostro V0, v1, v2) per 2 vertici in una sola volta.

Calcolare V0, v1 e v2 non solo per i primi 2 triangoli, ma anche per questi ultimi 2! Come ho detto dovresti fornire un totale di 4 vertici, o 8 galleggianti, per input. Solo ripetere i passaggi 2 e 3 per i dati.

Ora abbiamo 2 coppie di registri vx, ciascuna coppia contenente il risultato per 2 triangoli. Mi riferirò a loro come vx_0 (prima coppia) e vx_1 (seconda coppia) in cui x va da 0 a 2.

punto 4.

Dot prodotti. Per parallelizzare il calcolo baricentrico (in seguito) è necessario il risultato di ciascun prodotto punto per ciascuno dei 4 triangoli, in 1 registro singolo.

Quindi, dove calcolereste dot01 ad esempio, faremo lo stesso, ma per 4 triangoli in una volta. Ogni v -registrazione contiene il risultato per 2 vettori, quindi iniziamo moltiplicandoli.

Diciamo che u e v - parametri nella vostra funzione prodotto scalare - ora sono v0_0 e v1_0 (da calcolare dot01):

Multiply u e v per ottenere: [(v0_0.x0) * (v1_0.x0), (v0_0.z0) * (v1_0.z0), (v0_0.x1) * (v1_0.x1), (v0_0.z1) * (v1_0.z1)]

Questo può sembrare confuso a causa della .x0/.z0 e .x1/.z1, ma guardare a ciò che è stato caricato al punto 2: a0, a1.

Se ormai questo si sente ancora confuso, prendi un pezzo di carta e scrivi dall'inizio.

successivo, sempre per lo stesso prodotto scalare, fare la moltiplicazione per v0_1 e v1_1 (cioè la seconda coppia di triangoli).

Ora abbiamo 2 registri con 2 coppie X & Z (4 vertici totali), moltiplicati e pronti per essere aggiunti insieme per formare 4 punti separati.SSE3 ha un'istruzione per fare esattamente questo, e si chiama HADDPS:

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

HADDPS xmm0, XMM1 fa questo:

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

Ci vorranno le X & coppie Z dei primo registro, quelle della seconda, aggiungere insieme e memorizzarli nel primo, secondo, terzo e quarto componente del registro di destinazione. Ergo: a questo punto abbiamo questo particolare prodotto dot per tutti e 4 i triangoli!

** Ora ripetere questo processo per tutti i prodotti dot: dot00 eccetera. **

STEP 5.

L'ultimo calcolo (per quanto ho potuto vedere dal codice in dotazione) è la roba baricentrica. Questo è un calcolo scalare al 100% nel codice. Ma i tuoi input ora non sono risultati di prodotti a punti scalari (cioè singoli float), sono vettori/registri SSE con un prodotto dot per ciascuno dei 4 triangoli.

Quindi, se si esegue il flatout di questa operazione utilizzando le operazioni SSE parallele che operano su tutti e 4 i float, alla fine si otterrà 1 registro (o risultato) che trasporta 4 altezze, 1 per ciascun triangolo.

Dal momento che la mia pausa pranzo è scaduta, non ho intenzione di scrivere questo, ma dato il setup/idea che ho dato questo è un ultimo passo e non dovrebbe essere difficile da capire.

So che questa idea è un po 'lunga e richiede un po' di amore dal codice che si trova sopra di esso e forse un po 'di tempo con carta e matita, ma sarà veloce (e puoi anche aggiungere OpenMP in seguito se ti piacerebbe).

Buona fortuna :)

(e perdona la mia spiegazione fuzzy, posso montare la funzione se necessario =))

UPDATE

ho scritto un'implementazione ed non è andato come mi aspettavo, principalmente perché la componente Y è stata coinvolta oltre il pezzo di codice che hai inizialmente incollato nella tua domanda (l'ho cercato). Quello che ho fatto qui non è solo chiedervi di riorganizzare tutti i punti in coppie XZ e dar loro da mangiare per 4, ma anche di alimentare 3 puntatori (per i punti A, B e C) con i valori Y per ciascuno dei 4 triangoli. Da una prospettiva locale questo è il più veloce. Posso ancora modificarlo per richiedere modifiche meno invadenti dalla fine del callee, ma per favore fammi sapere cosa è desiderabile.

Quindi una dichiarazione di non responsabilità: questo codice è semplice come l'inferno (qualcosa che ho trovato che funziona molto bene con i compilatori in termini di SSE ... possono riorganizzarsi come si vede e le CPU x86/x64 prendono la loro parte anche in esso). Anche la denominazione, non è il mio stile, non è di nessuno, fallo con quello che ritieni opportuno.

speriamo vi sia utile e se non sarò contento di andare su di esso di nuovo.E se questo è un progetto commerciale c'è anche la possibilità di me di salire a bordo immagino;)

In ogni caso, ho messo su pastebin: http://pastebin.com/20u8fMEb

Answer 2

è possibile implementare il prodotto singolo punto con istruzioni SSE, ma il risultato non sarà molto più veloce (e può anche essere più lento) rispetto al codice come scritto ora La tua riscrittura potrebbe vanificare le ottimizzazioni del compilatore che stanno aiutando la versione corrente.

Per trarre beneficio dalla riscrittura che con SSE o CUDA è necessario ottimizzare il ciclo che richiede quel prodotto punto. Questo è particolarmente vero per CUDA, dove il sovraccarico di esecuzione di un prodotto punto sarebbe enorme. Puoi vedere un aumento della velocità solo se hai inviato migliaia di vettori nella GPU per calcolare migliaia di prodotti punto. La stessa idea vale per SSE sulla CPU, ma potresti essere in grado di vedere un miglioramento rispetto a un numero minore di operazioni. Sarà comunque superiore a un prodotto punto.

La cosa più semplice da provare potrebbe essere g++ -ftree-vectorize. GCC sarà in grado di allineare la tua piccola funzione e quindi provare a ottimizzare il ciclo per te (probabilmente lo è già, ma senza le istruzioni SSE). Il vectorizer dell'albero tenterà di eseguire automaticamente ciò che si propone di fare a mano. Non ha sempre successo.

Answer 3

hai chiesto una versione SSE del vostro algoritmo ecco è:

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

La XMVector2Dot è preso da xnamathvector.inl, ho rinominato e modificato per operare sulle coordinate X/Z.

XNAMath è una grande libreria matematica cross-platform vettoriale di Microsoft; Lo uso anche su OS X importando l'intestazione sal.h (non sono sicuro del problema delle licenze, quindi fai attenzione).
Infatti, qualsiasi piattaforma che supporti intrinseci SSE dovrebbe supportarlo.

Un paio di cose da guardare per:

• È necessario caricare i carri in XMVECTORs utilizzando il metodo XMLoadFloat3; questo caricherà float3 non allineato in una struttura __m128.
• Probabilmente non vedrete alcun miglioramento delle prestazioni da questo codice SSE (profilo !!) in quanto vi è una penalizzazione delle prestazioni per caricare i float non allineati nei registri SSE.
• Questa è una conversione a forza bruta dell'algoritmo in SSE, si avrà maggiore fortuna essendo più intelligenti di me e in realtà cercando di capire l'algoritmo e implementando una versione amichevole di SSE.
• Avrai maggiore fortuna convertendo l'intera applicazione in modo che utilizzi il codice XNA Math/SSE anziché solo quella piccola porzione. Implicare almeno l'uso di tipi di vettore allineati (XMFLOAT3A o struct __declspec (align (16)) myvectortype {};).
• L'assemblaggio SSE diritto è sconsigliato soprattutto in x64, a favore di elementi intrinseci.

Answer 4

Le istruzioni SSE sono pensate per ottimizzare gli algoritmi che elaborano grandi blocchi di dati rappresentati come numeri interi o in virgola mobile. Le dimensioni tipiche sono milioni e miliardi di numeri che devono essere elaborati in qualche modo. Non ha senso ottimizzare la funzione che elabora solo quattro (o venti) scalari. Quello che guadagni con SSE potresti perdere con la funzione chiamata overhead. Il numero ragionevole di numeri elaborati da una chiamata di funzione è almeno di mille. È possibile che si ottenga un enorme guadagno di prestazioni usando l'intrinseco SSE. Ma è difficile darti un consiglio specifico su misura per le tue esigenze in base alle informazioni che hai fornito. Dovresti modificare la tua domanda e fornire una visione più ad alto livello del tuo problema (le funzioni si trovano più in profondità nel tuo callstack). Ad esempio, non è ovvio quante volte il metodo dtClosestHeightPointTriangle viene chiamato al secondo? Questo numero è fondamentale per giudicare obiettivamente se la transizione verso SSE sarebbe di valore pratico. Anche l'organizzazione dei dati è molto importante. Idealmente, i tuoi dati dovrebbero essere archiviati nel minor numero possibile di segmenti lineari di memoria per utilizzare efficientemente il sottosistema di cache della CPU.