Skip to content

Interlacciare velocemente 2 array doppi in un array di strutture con 2 membri float e 1 int (invariante del ciclo), con conversione SIMD double->float?

Dopo aver ricercato con specialisti in materia, programmatori di diverse aree e professori, abbiamo trovato la risposta al problema e la condividiamo in questa pubblicazione.

Soluzione:

Ecco un tentativo con SSE4.1, senza AVX (è più difficile da fare e per ora mi vengono in mente ancora più rimescolamenti), e usando il formato 12byte/punto: (non testato)

void test3(MyStruct * _pPoints, double * pInputValues1, double * pInputValues2) {
        // struct MyStruct 
        // { 
        //    float O1;
        //    float O2;
        //    unsigned int Offset;
        // };
    __m128 offset = _mm_castsi128_ps(_mm_cvtsi32_si128(_uiDefaultOffset));
    int i;
    for (i = 0; i < _iNum - 2; i += 2)
    {
        // read inputs and convert to float
        __m128d inA = _mm_loadu_pd(&pInputValues1[i]);
        __m128d inB = _mm_loadu_pd(&pInputValues2[i]);
        __m128 inAf = _mm_cvtpd_ps(inA);    // 0 0 A1 A0
        __m128 inBf = _mm_cvtpd_ps(inB);    // 0 0 B1 B0
        // shuffle B0 from place 0 to place 1, merge with offset
        __m128 tempA = _mm_shuffle_ps(inBf, offset, _MM_SHUFFLE(1, 0, 0, 0)); // 0 OF B0 B0
        // shuffle A1 from place 1 to place 0, merge with offset
        __m128 tempB = _mm_shuffle_ps(inAf, offset, _MM_SHUFFLE(1, 0, 1, 1)); // 0 OF A1 A1
        // replace B0 at place 0 with A0
        __m128 outA = _mm_blend_ps(tempA, inAf, 1);  // 0 OF B0 A0
        // replace A1 at place 1 with B1
        __m128 outB = _mm_blend_ps(tempB, inBf, 2);  // 0 OF B1 A1
        // store results
        _mm_storeu_ps(&_pPoints[i].O1, outA);
        _mm_storeu_ps(&_pPoints[i + 1].O1, outB);
    }
    // remaining iteration if _iNum is not even
    for (; i < _iNum; i++)
    {
        _pPoints[i].O1 = static_cast(pInputValues1[i]);
        _pPoints[i].O2 = static_cast(pInputValues2[i]);
        _pPoints[i].Offset = _uiDefaultOffset;
    }
}

Questo utilizza la capacità di shufps di scegliere da due fonti diverse per fare la fusione dei dati dinamici e l'offset costante, gli stessi rimescolamenti spostano anche il galleggiante in ogni gruppo che deve essere spostato. Poi si usano le miscele per sostituire un singolo galleggiante con un altro galleggiante che era già al posto giusto. Questo richiede 2 shuffle e 2 blend, c'è anche un modo con 3 shuffle e zero blend, ma gli shuffle vanno tutti in p5 sugli attuali processori Intel mentre i blend possono andare in una porta diversa. Anche le conversioni usano già la p5, quindi si sta ingolfando; usare le miscele dovrebbe essere meglio. Sono comunque 4 µops p5 per iterazione, quindi ci vogliono almeno 2 cicli per ogni elemento processato, il che non è il massimo.

Il ciclo principale salta gli ultimi elementi in modo da non scrivere al di fuori dei limiti, ma esegue memorizzazioni leggermente sovrapposte di 16 byte che scrivono 4 byte oltre la fine della struttura. Questa parte viene sovrascritta con il risultato reale dalla memorizzazione successiva, ma potrebbe essere pericoloso farlo alla fine dell'array.

Questo problema non è molto simile a memcpy. Si tratta di ottimizzare l'interleaving con shuffle e/o memorizzazione scalare del membro intero invariante del ciclo. Questo rende difficile la SIMD.

Lo fate anche voi bisogno di di avere questo formato di memorizzazione con il int interleaved con il float membri? L'interleaving dei float è già abbastanza grave. Suppongo che un codice successivo modificherà il membro intin diverse strutture, altrimenti non ha senso duplicarlo per ogni elemento.

Si potrebbe lavorare a gruppi di 4 elementi, come struct { float a[4], b[4]; int i[4]; }; in modo da poter caricare+convertire 4x elementi contigui double in 4x float e fare una memorizzazione SIMD a 128 bit? Si avrebbe comunque una certa località spaziale quando si accede a tutti e 3 i membri di una singola "struct" di output-array.


Comunque, supponendo che il formato di output debba essere completamente interleaved, non abbiamo bisogno di imbottirlo a 16 byte. Le CPU x86 possono gestire in modo efficiente la sovrapposizione di store a 16 byte per lavorare con strutture a 12 byte, come mostra la risposta di @harold. Le suddivisioni delle linee di cache probabilmente costano meno della larghezza di banda extra della memoria necessaria per memorizzare il padding.

Oppure un'altra strategia potrebbe essere quella di usare magazzini separati per i float rispetto al intin modo da non avere bisogno di sovrapposizioni. Probabilmente possiamo ottimizzarla fino al punto in cui dovrebbe avere un collo di bottiglia con 1 memorizzazione per ciclo di clock per 1 struct ogni 2 cicli. (O leggermente inferiore, perché i negozi di cache-split devono riprodurre l'operazione di memorizzazione, almeno sulle CPU Intel). Potremmo anche srotolare per 4*12 = 3*16 byte e risparmiare 2 memorizzazioni di interi utilizzando memorizzazioni SIMD che vengono sovrapposte ai dati float. 48 byte = xyIx|yIxy|IxyI ha quattro I come parte di quattro strutture, ma sono abbastanza vicini da poterli memorizzare tutti e quattro con due _mm_storeu_si128( set1(offset) ) intrinseche. Quindi memorizzare gli elementi xy che si sovrappongono. I limiti di 16 byte sono contrassegnati da |. Se le suddivisioni delle linee di cache sono un problema, si potrebbero fare 2 scalari e un SIMD per l'ultimo vettore che è allineato (se l'array di output è allineato a 16 byte). Oppure, sulle CPU Intel Haswell e successive, un archivio allineato a 32 byte potrebbe andare bene.


Se non si fa attenzione, è molto facile che si verifichi un collo di bottiglia nel throughput dello shuffle sulle CPU Intel, in particolare quelle della famiglia Sandybridge (da SnB a Skylake/Coffee Lake), dove gli shuffle FP possono essere eseguiti solo sulla porta 5. Questo è il motivo per cui sto considerando non di mischiare tutto insieme per avere un solo archivio per struttura.

La conversione SIMD double->float costa 2 uops: shuffle + FP-math, perché il float è largo la metà e l'istruzione impacchetta i float nella parte inferiore del registro vettoriale.

AVX è utile per convertire 4 doublein un vettore SIMD di 4 floats.

A parte questo, sono d'accordo con @harold che i vettori a 128 bit sono probabilmente una buona scommessa. Anche l'AVX2 non ha un buon sistema di rimescolamento a 2 ingressi, e l'AVX1 è molto limitato. Quindi possiamo usare la conversione 256-bit -> 128-bit double->float per alimentare una strategia di interleave basata su __m128.

vmovhps [mem], xmm non costa un uop di rimescolamento sulle CPU Intel, ma solo un puro store, quindi rimescola insieme 2 vettori e ottiene [ B1 A1 B0 A0 ] in un singolo vettore ci permette di ottenere due memorizzazioni a 64 bit delle metà bassa e alta senza ulteriori rimescolamenti.

Tuttavia, la versione di @harold potrebbe essere migliore. 4 shuffle per 2 struct potrebbero essere migliori di 4 store per 2 struct, dato che gli store a volte devono essere riprodotti per la divisione delle linee di cache, mentre gli shuffle non lo fanno. Ma con il trucco della sovrapposizione dei negozi, 3,5 o 3 negozi per 2 strutture sembrano fattibili.


Oppure, ecco un'altra idea che utilizza alcuni dei metodi sopra descritti, ma che prevede una miscelazione per salvare i negozi

In pratica mi è venuto in mente mentre modificavo il codice di @harold per implementare l'idea di cui ho scritto nel testo precedente. L'uso di un blend è un buon modo per ridurre la pressione sulle porte store e shuffle.

Alcune delle idee precedenti valgono ancora la pena di essere esplorate, soprattutto per quanto riguarda la memorizzazione di un'ampia porzione di set1(offset) e poi sovrapporla a quella a 64 bit di vmovlps a 64 bit. (Dopo aver srotolato per 3x2 = 6 o 3x4 = 12 le strutture di output, per renderle un multiplo dei 4 doppi che convertiamo in una volta sola). 12 * 12 = 144 byte, che è un multiplo di 16 ma non di 32 o 64, in modo da poter sapere in ogni momento dove ci troviamo rispetto a un confine di 16 byte, ma non rispetto alle linee di cache, a meno di non srotolare ancora di più. (potenzialmente lasciando altro lavoro da ripulire e gonfiando la dimensione del codice).

#include 
#include 
#include 

struct f2u { 
  float O1, O2;
  unsigned int Offset;
};

// names with a leading _ at file scope are reserved for the implementation.
// fixed that portability problem for you.
static const unsigned uiDefaultOffset = 123;

// only requires AVX1
// ideally pA and pB should be 32-byte aligned.
// probably also dst 16-byte aligned is good.
void cvt_interleave_avx(f2u *__restrict dst, double *__restrict pA, double *__restrict pB, ptrdiff_t len)
{
    __m128 voffset = _mm_castsi128_ps(_mm_set1_epi32(uiDefaultOffset));

    // 48 bytes per iteration: 3x16 = 4x12
    ptrdiff_t i;
    for (i = 0; i < len - 3; i += 4)
    {
        // read inputs and convert to float
        __m256d inA = _mm256_loadu_pd(&pA[i]);
        __m256d inB = _mm256_loadu_pd(&pB[i]);
        __m128 inAf = _mm256_cvtpd_ps(inA);    // A3 A2 A1 A0
        __m128 inBf = _mm256_cvtpd_ps(inB);    // B3 B2 B1 B0

        // interleave to get XY pairs
        __m128 lo = _mm_unpacklo_ps(inAf, inBf); // B1 A1 B0 A0
        __m128 hi = _mm_unpackhi_ps(inAf, inBf); // B3 A3 B2 A2

        // blend integer into place
        __m128 out0 = _mm_blend_ps(lo, voffset, 1<<2);  // x OF B0 A0
        __m128 out2 = _mm_blend_ps(hi, voffset, 1<<2);  // x OF B2 A2

        // TODO: _mm_alignr_epi8 to create OF OF B1 A1 spending 1 more shuffle to save a store.

        // store results
        _mm_storeu_ps(&dst[i + 0].O1, out0);  // 16 bytes with blended integer
        _mm_storeh_pi((__m64*)&dst[i + 1].O1, lo);    // 8 bytes from top half of reg, partial overlap
        dst[i + 1].Offset = uiDefaultOffset;

        _mm_storeu_ps(&dst[i + 2].O1, out2);  // 16 bytes with blended integer
        _mm_storeh_pi((__m64*)&dst[i + 3].O1, hi);    // 8 bytes from top half of reg, partial overlap
        dst[i + 3].Offset = uiDefaultOffset;
    }

    // scalar cleanup for  if _iNum is not even
    for (; i < len; i++)
    {
        dst[i].O1 = static_cast(pA[i]);
        dst[i].O2 = static_cast(pB[i]);
        dst[i].Offset = uiDefaultOffset;
    }
}

gcc9.1 -O3 -march=skylake su Godbolt compila il ciclo principale con 19 uops a dominio fuso per il front-end. (Né vcvtpd2ps possono essere microfuse, perché GCC non ha fatto nulla di intelligente come l'indirizzamento di pB rispetto a pA per evitare una modalità di indirizzamento indicizzata per una di esse. Quindi sono 3 uops ciascuno: load + convert + shuffle)

Ma il collo di bottiglia si verifica comunque con i depositi nel back-end, anche se ci vogliono ben 5 cicli per iterazione per emettere dal front-end a 4 larghezze.

Con 6 negozi (per 4 strutture) per iterazione, il collo di bottiglia si riduce al massimo a 1 iterazione per 6 cicli, con un collo di bottiglia sulla porta dei negozi/dati/unità di esecuzione. (Fino a Ice Lake, che può eseguire 2 memorizzazioni per clock). Quindi questo ottiene 1 struct ogni 1,5 cicli nel caso teorico migliore, lo stesso che stimavo prima per l'idea dell'overlapping-store.

(Sappiamo già che i negozi divisi in linee di cache devono essere rigiocati, il che costa in termini di throughput, quindi sappiamo che questo non riuscirà a raggiungere 1,5 cicli per struttura anche senza missioni di cache. Ma probabilmente è ancora meglio del collo di bottiglia di Harold, 4 cicli per 2 strutture = 2 cicli per struttura. Questa velocità dovrebbe comunque essere raggiungibile, perché il collo di bottiglia si verifica con gli shuffle, che non hanno bisogno di essere rigiocati in caso di split di cache-line).

Mi aspetto che il throughput su Ryzen sia simile, con un collo di bottiglia sul throughput dei negozi. Utilizziamo principalmente vettori a 128 bit e Ryzen ha un throughput di shuffle migliore di Intel. Su SnB-family, ci sono 4 uops di shuffle nel loop.

Se potessi mischiare in modo diverso in modo da poter ottenere due strutture contigue come metà alta della coppia di vettori, si aprirebbe la possibilità di combinare le due assegnazioni scalari in una sola _mm_storeu_si128 che sovrappongo a due _mm_storeh_pi (movhps) a 64 bit. (Facendo ancora due miscele per le altre due strutture di output). Questo ridurrebbe il tutto a 5 memorizzazioni totali.

Ma shufps ha delle restrizioni su dove prende i dati di partenza, quindi non è possibile usarlo per emulare unpcklps o interlacciare in modo diverso.

Probabilmente varrebbe la pena di usare palignr per la struttura B1 A1, spendendo un uop in più di shuffle per risparmiare una memorizzazione.

Non ho fatto un benchmark né ho calcolato la frequenza con cui i depositi non allineati attraversano il confine di una linea della cache (e quindi il costo del throughput).


AVX512

Se avessimo AVX512, avremmo uno shuffle a 2 ingressi che ci permetterebbe di costruire vettori di dati float+int in modo più efficiente, con un minor numero di istruzioni di shuffle e store per struttura. (Potremmo usare vpermt2ps con la maschera di fusione in set1(integer) per interlacciare 2 vettori di risultati di conversione insieme a numeri interi nei posti giusti).

Ispirata all'esempio di trasposizione 4x3 di Intel e basata sulla soluzione di @PeterCordes, ecco una soluzione AVX1, che dovrebbe ottenere un throughput di 8 strutture in 8 cicli (il collo di bottiglia è ancora p5):

#include 
#include 

struct f2u { 
  float O1, O2;
  unsigned int Offset;
};
static const unsigned uiDefaultOffset = 123;

void cvt_interleave_avx(f2u *__restrict dst, double *__restrict pA, double *__restrict pB, ptrdiff_t len)
{
    __m256 voffset = _mm256_castsi256_ps(_mm256_set1_epi32(uiDefaultOffset));

    // 8 structs per iteration
    ptrdiff_t i=0;
    for(; i(dst + i);

        // 4*vcvtpd2ps    --->  4*(p1,p5,p23)
        __m128 inA3210 = _mm256_cvtpd_ps(_mm256_loadu_pd(&pA[i]));
        __m128 inB3210 = _mm256_cvtpd_ps(_mm256_loadu_pd(&pB[i]));
        __m128 inA7654 = _mm256_cvtpd_ps(_mm256_loadu_pd(&pA[i+4]));
        __m128 inB7654 = _mm256_cvtpd_ps(_mm256_loadu_pd(&pB[i+4]));

        // 2*vinsertf128  --->  2*p5
        __m256 A76543210 = _mm256_set_m128(inA7654,inA3210);
        __m256 B76543210 = _mm256_set_m128(inB7654,inB3210);

        // 2*vpermilps    --->  2*p5
        __m256 A56741230 = _mm256_shuffle_ps(A76543210,A76543210,_MM_SHUFFLE(1,2,3,0));
        __m256 B67452301 = _mm256_shuffle_ps(B76543210,B76543210,_MM_SHUFFLE(2,3,0,1));

        // 6*vblendps     ---> 6*p015 (does not need to use p5)
        __m256 outA1__B0A0 = _mm256_blend_ps(A56741230,B67452301,2+16*2);
        __m256 outA1ccB0A0 = _mm256_blend_ps(outA1__B0A0,voffset,4+16*4);

        __m256 outB2A2__B1 = _mm256_blend_ps(B67452301,A56741230,4+16*4);
        __m256 outB2A2ccB1 = _mm256_blend_ps(outB2A2__B1,voffset,2+16*2);

        __m256 outccB3__cc = _mm256_blend_ps(voffset,B67452301,4+16*4);
        __m256 outccB3A3cc = _mm256_blend_ps(outccB3__cc,A56741230,2+16*2);

        // 3* vmovups     ---> 3*(p237,p4)
        _mm_storeu_ps(dst_f+ 0,_mm256_castps256_ps128(outA1ccB0A0));
        _mm_storeu_ps(dst_f+ 4,_mm256_castps256_ps128(outB2A2ccB1));
        _mm_storeu_ps(dst_f+ 8,_mm256_castps256_ps128(outccB3A3cc));
        // 3*vextractf128 ---> 3*(p23,p4)
        _mm_storeu_ps(dst_f+12,_mm256_extractf128_ps(outA1ccB0A0,1));
        _mm_storeu_ps(dst_f+16,_mm256_extractf128_ps(outB2A2ccB1,1));
        _mm_storeu_ps(dst_f+20,_mm256_extractf128_ps(outccB3A3cc,1));
    }

    // scalar cleanup for  if _iNum is not even
    for (; i < len; i++)
    {
        dst[i].O1 = static_cast(pA[i]);
        dst[i].O2 = static_cast(pB[i]);
        dst[i].Offset = uiDefaultOffset;
    }
}

Collegamento a Godbolt, con un codice di prova minimo alla fine: https://godbolt.org/z/0kTO2b

Per qualche ragione, a gcc non piace generare il file vcvtpd2ps che convertono direttamente dalla memoria a un registro. Questo potrebbe funzionare meglio con carichi allineati (avere input e output allineati è probabilmente vantaggioso in ogni caso). E clang apparentemente vuole superarmi in astuzia con uno dei comandi vextractf128 alla fine.

Valutazioni e recensioni

Se puoi, puoi lasciare un tutorial su cosa aggiungeresti a questa affermazione.



Utilizzate il nostro motore di ricerca

Ricerca
Generic filters

Lascia un commento

Il tuo indirizzo email non sarà pubblicato.