Trasposizione memoria veloce con SSE, AVX e OpenMP

Question

Ho bisogno di un algoritmo di trasposizione della memoria veloce per la mia funzione di convoluzione gaussiana in C/C++. Quello che faccio ora è

convolute_1D 
transpose 
convolute_1D 
transpose 

Risulta che con questo metodo la grandezza del filtro deve essere grande (o maggiore del previsto) o la trasposta richiede più la spira (ad esempio per una matrice 1920x1080 la convoluzione prende allo stesso tempo della trasposizione per una dimensione del filtro di 35). L'algoritmo di trasposizione corrente che sto usando utilizza il loop blocking/tiling insieme a SSE e OpenMP. Ho provato una versione usando AVX ma non è più veloce. Qualche suggerimento su come posso accelerarlo?

inline void transpose4x4_SSE(float *A, float *B, const int lda, const int ldb) { 
    __m128 row1 = _mm_load_ps(&A[0*lda]); 
    __m128 row2 = _mm_load_ps(&A[1*lda]); 
    __m128 row3 = _mm_load_ps(&A[2*lda]); 
    __m128 row4 = _mm_load_ps(&A[3*lda]); 
    _MM_TRANSPOSE4_PS(row1, row2, row3, row4); 
    _mm_store_ps(&B[0*ldb], row1); 
    _mm_store_ps(&B[1*ldb], row2); 
    _mm_store_ps(&B[2*ldb], row3); 
    _mm_store_ps(&B[3*ldb], row4); 
} 
//block_size = 16 works best 
inline void transpose_block_SSE4x4(float *A, float *B, const int n, const int m, const int lda, const int ldb ,const int block_size) { 
    #pragma omp parallel for 
    for(int i=0; i<n; i+=block_size) { 
     for(int j=0; j<m; j+=block_size) { 
      int max_i2 = i+block_size < n ? i + block_size : n; 
      int max_j2 = j+block_size < m ? j + block_size : m; 
      for(int i2=i; i2<max_i2; i2+=4) { 
       for(int j2=j; j2<max_j2; j2+=4) { 
        transpose4x4_SSE(&A[i2*lda +j2], &B[j2*ldb + i2], lda, ldb); 
       } 
      } 
     } 
    } 

}

Trasposizione 8x8 matrice float utilizzando AVX. Non è più veloce di quattro trasposizioni 4x4.

inline void transpose8_ps(__m256 &row0, __m256 &row1, __m256 &row2, __m256 &row3, __m256 &row4, __m256 &row5, __m256 &row6, __m256 &row7) { 
__m256 __t0, __t1, __t2, __t3, __t4, __t5, __t6, __t7; 
__m256 __tt0, __tt1, __tt2, __tt3, __tt4, __tt5, __tt6, __tt7; 
__t0 = _mm256_unpacklo_ps(row0, row1); 
__t1 = _mm256_unpackhi_ps(row0, row1); 
__t2 = _mm256_unpacklo_ps(row2, row3); 
__t3 = _mm256_unpackhi_ps(row2, row3); 
__t4 = _mm256_unpacklo_ps(row4, row5); 
__t5 = _mm256_unpackhi_ps(row4, row5); 
__t6 = _mm256_unpacklo_ps(row6, row7); 
__t7 = _mm256_unpackhi_ps(row6, row7); 
__tt0 = _mm256_shuffle_ps(__t0,__t2,_MM_SHUFFLE(1,0,1,0)); 
__tt1 = _mm256_shuffle_ps(__t0,__t2,_MM_SHUFFLE(3,2,3,2)); 
__tt2 = _mm256_shuffle_ps(__t1,__t3,_MM_SHUFFLE(1,0,1,0)); 
__tt3 = _mm256_shuffle_ps(__t1,__t3,_MM_SHUFFLE(3,2,3,2)); 
__tt4 = _mm256_shuffle_ps(__t4,__t6,_MM_SHUFFLE(1,0,1,0)); 
__tt5 = _mm256_shuffle_ps(__t4,__t6,_MM_SHUFFLE(3,2,3,2)); 
__tt6 = _mm256_shuffle_ps(__t5,__t7,_MM_SHUFFLE(1,0,1,0)); 
__tt7 = _mm256_shuffle_ps(__t5,__t7,_MM_SHUFFLE(3,2,3,2)); 
row0 = _mm256_permute2f128_ps(__tt0, __tt4, 0x20); 
row1 = _mm256_permute2f128_ps(__tt1, __tt5, 0x20); 
row2 = _mm256_permute2f128_ps(__tt2, __tt6, 0x20); 
row3 = _mm256_permute2f128_ps(__tt3, __tt7, 0x20); 
row4 = _mm256_permute2f128_ps(__tt0, __tt4, 0x31); 
row5 = _mm256_permute2f128_ps(__tt1, __tt5, 0x31); 
row6 = _mm256_permute2f128_ps(__tt2, __tt6, 0x31); 
row7 = _mm256_permute2f128_ps(__tt3, __tt7, 0x31); 
} 

inline void transpose8x8_avx(float *A, float *B, const int lda, const int ldb) { 
    __m256 row0 = _mm256_load_ps(&A[0*lda]); 
    __m256 row1 = _mm256_load_ps(&A[1*lda]); 
    __m256 row2 = _mm256_load_ps(&A[2*lda]); 
    __m256 row3 = _mm256_load_ps(&A[3*lda]); 
    __m256 row4 = _mm256_load_ps(&A[4*lda]); 
    __m256 row5 = _mm256_load_ps(&A[5*lda]); 
    __m256 row6 = _mm256_load_ps(&A[6*lda]); 
    __m256 row7 = _mm256_load_ps(&A[7*lda]);  
    transpose8_ps(row0, row1, row2, row3, row4, row5, row6, row7); 
    _mm256_store_ps(&B[0*ldb], row0); 
    _mm256_store_ps(&B[1*ldb], row1); 
    _mm256_store_ps(&B[2*ldb], row2); 
    _mm256_store_ps(&B[3*ldb], row3); 
    _mm256_store_ps(&B[4*ldb], row4); 
    _mm256_store_ps(&B[5*ldb], row5); 
    _mm256_store_ps(&B[6*ldb], row6); 
    _mm256_store_ps(&B[7*ldb], row7); 

} 

Answer 1

Direi che la soluzione migliore sarebbe quella di cercare di combinare la convoluzione e la trasposta - cioè scrivere i risultati del convolve trasposto, come si va. Quasi certamente la larghezza di banda della memoria è limitata sulla trasposizione, quindi ridurre il numero di istruzioni utilizzate per la trasposizione non è di grande aiuto (da qui la mancanza di miglioramenti nell'uso di AVX). Riducendo il numero di passaggi sui tuoi dati otterrai i migliori miglioramenti delle prestazioni.

Answer 2

FWIW, su una CPU per laptop Core i7 M di 3 anni, questa trasposizione 4x4 ingenua era appena più lenta della versione SSE, mentre quasi il 40% più veloce su una nuova CPU desktop Intel Xeon E5-2630 v2 @ 2,60 GHz.

inline void transpose4x4_naive(float *A, float *B, const int lda, const int ldb) { 
    const float r0[] = { A[0], A[1], A[2], A[3] }; // memcpy instead? 
    A += lda; 
    const float r1[] = { A[0], A[1], A[2], A[3] }; 
    A += lda; 
    const float r2[] = { A[0], A[1], A[2], A[3] }; 
    A += lda; 
    const float r3[] = { A[0], A[1], A[2], A[3] }; 

    B[0] = r0[0]; 
    B[1] = r1[0]; 
    B[2] = r2[0]; 
    B[3] = r3[0]; 
    B += ldb; 
    B[0] = r0[1]; 
    B[1] = r1[1]; 
    B[2] = r2[1]; 
    B[3] = r3[1]; 
    B += ldb; 
    B[0] = r0[2]; 
    B[1] = r1[2]; 
    B[2] = r2[2]; 
    B[3] = r3[2]; 
    B += ldb; 
    B[0] = r0[3]; 
    B[1] = r1[3]; 
    B[2] = r2[3]; 
    B[3] = r3[3]; 
} 

Stranamente, la CPU vecchio computer portatile è più veloce del desktop v2 dual E5-2630 con il doppio della base, ma questa è un'altra storia :)

In caso contrario, si potrebbe anche essere interessato a http://research.colfaxinternational.com/file.axd?file=2013%2F8%2FColfax_Transposition-7110P.pdf http://colfaxresearch.com/multithreaded-transposition-of-square-matrices-with-common-code-for-intel-xeon-processors-and-intel-xeon-phi-coprocessors/ (richiede login ora ...)

Answer 3

Considerate questa trasposizione 4x4.

struct MATRIX { 
    union { 
     float f[4][4]; 
     __m128 m[4]; 
     __m256 n[2]; 
    }; 
}; 
MATRIX myTranspose(MATRIX in) { 

    // This takes 15 assembler instructions (compile not inline), 
    // and is faster than XMTranspose 
    // Comes in like this 1 2 3 4 5 6 7 8 
    //      9 10 11 12 13 14 15 16 
    // 
    // Want the result  1 5 9 13 2 6 10 14 
    //      3 7 11 15 4 8 12 16 

    __m256 t0, t1, t2, t3, t4, t5, n0, n1; 
    MATRIX result; 

    n0 = in.n[0];            // n0 = 1, 2, 3, 4, 5, 6, 7, 8 
    n1 = in.n[1];            // n1 = 9, 10, 11, 12, 13, 14, 15, 16 
    t0 = _mm256_unpacklo_ps(n0, n1);       // t0 = 1, 9, 2, 10, 5, 13, 6, 14 
    t1 = _mm256_unpackhi_ps(n0, n1);       // t1 = 3, 11, 4, 12, 7, 15, 8, 16 

    t2 = _mm256_permute2f128_ps(t0, t1, 0x20);     // t2 = 1, 9, 2, 10, 3, 11, 4, 12 
    t3 = _mm256_permute2f128_ps(t0, t1, 0x31);     // t3 = 5, 13, 6, 14, 7, 15, 8, 16 

    t4 = _mm256_unpacklo_ps(t2, t3);       // t2 = 1, 5, 9, 13, 3, 7, 11, 15 
    t5 = _mm256_unpackhi_ps(t2, t3);       // t3 = 2, 6, 10, 14, 4, 8, 12, 16 

    result.n[0] = _mm256_permute2f128_ps(t4, t5, 0x20);   // t6 = 1, 5, 9, 13, 2, 6, 10, 14 
    result.n[1] = _mm256_permute2f128_ps(t4, t5, 0x31);   // t7 = 3, 7, 11, 15, 4, 8, 12, 16 
    return result; 
}