réduction avec OpenMP avec SSE/AVX

Question

Je veux faire une réduction sur un tableau en utilisant OpenMP et SIMD. J'ai lu qu'une réduction en OpenMP est équivalente à :

inline float sum_scalar_openmp2(const float a[], const size_t N) {
    float sum = 0.0f;
    #pragma omp parallel
    {
        float sum_private = 0.0f;
        #pragma omp parallel for nowait
        for(int i=0; i<N; i++) {
            sum_private += a[i];
        }
        #pragma omp atomic
        sum += sum_private;
    }
    return sum;
}

J'ai eu cette idée grâce au lien suivant : http://bisqwit.iki.fi/story/howto/openmp/#ReductionClause Mais atomic ne supporte pas non plus les opérateurs complexes. Ce que j'ai fait, c'est remplacer atomique par critique et implémenter la réduction avec OpenMP et SSE comme ceci :

#define ROUND_DOWN(x, s) ((x) & ~((s)-1))
inline float sum_vector4_openmp(const float a[], const size_t N) {
    __m128 sum4 = _mm_set1_ps(0.0f);
    #pragma omp parallel 
    {
        __m128 sum4_private = _mm_set1_ps(0.0f);
        #pragma omp for nowait
        for(int i=0; i < ROUND_DOWN(N, 4); i+=4) {
            __m128 a4 = _mm_load_ps(a + i);
            sum4_private = _mm_add_ps(a4, sum4_private);
        }
        #pragma omp critical
        sum4 = _mm_add_ps(sum4_private, sum4);
    }
    __m128 t1 = _mm_hadd_ps(sum4,sum4);
    __m128 t2 = _mm_hadd_ps(t1,t1);
    float sum = _mm_cvtss_f32(t2);  
    for(int i = ROUND_DOWN(N, 4); i < N; i++) {
        sum += a[i];
    }
    return sum;
} 

Cependant, cette fonction ne fonctionne pas aussi bien que je l'espère. J'utilise Visual Studio 2012 Express. Je sais que je peux améliorer un peu les performances en déroulant plusieurs fois le load/add SSE, mais c'est toujours moins que ce que j'attends.

J'obtiens de bien meilleures performances en fonctionnant sur des tranches de tableaux égales au nombre de threads :

inline float sum_slice(const float a[], const size_t N) {
    int nthreads = 4;
    const int offset = ROUND_DOWN(N/nthreads, nthreads);
    float suma[8] = {0};
    #pragma omp parallel for num_threads(nthreads) 
    for(int i=0; i<nthreads; i++) {
        suma[i] = sum_vector4(&a[i*offset], offset);
    }
    float sum = 0.0f;
    for(int i=0; i<nthreads; i++) {
        sum += suma[i]; 
    }
    for(int i=nthreads*offset; i < N; i++) {
        sum += a[i];
    }
    return sum;    
}

inline float sum_vector4(const float a[], const size_t N) {
    __m128 sum4 = _mm_set1_ps(0.0f);
    int i = 0;
    for(; i < ROUND_DOWN(N, 4); i+=4) {
        __m128 a4 = _mm_load_ps(a + i);
        sum4 = _mm_add_ps(sum4, a4);
    }
    __m128 t1 = _mm_hadd_ps(sum4,sum4);
    __m128 t2 = _mm_hadd_ps(t1,t1);
    float sum = _mm_cvtss_f32(t2);
    for(; i < N; i++) {
        sum += a[i];
    }
    return sum;

}

Quelqu'un sait-il s'il existe une meilleure façon de faire des réductions avec des opérateurs plus compliqués en OpenMP ?

Answer 1

Je pense que la réponse à votre question est Non. Je ne pense pas qu'il y ait une meilleure façon de faire de la réduction avec des opérateurs plus compliqués dans OpenMP.

En supposant que le tableau est aligné sur 16 bits, que le nombre de threads OpenMP est de 4, on peut s'attendre à un gain de performance de 12x - 16x avec OpenMP + SIMD. En réalité, il se peut que le gain de performance ne soit pas suffisant car

1. Il y a un surcoût dans la création des threads openmp.
2. Le code effectue une opération d'addition pour une opération de chargement. Par conséquent, le CPU ne fait pas assez de calculs. Il semble donc que le CPU passe la plupart du temps à charger les données, ce qui limite la bande passante de la mémoire.