prefixo paralelo (cumulativo) soma com SSE
Estou procurando alguns conselhos sobre como fazer uma soma de prefixo paralela com o SSE. Estou interessado em fazer isso em uma matriz de ints, floats ou doubles.
Eu criei duas soluções. Um caso especial e um caso geral. Em ambos os casos, a solução é executada sobre o array em duas passagens em paralelo com o OpenMP. Para o caso especial eu uso o SSE em ambos os passes. Para o caso geral eu uso apenas na segunda passagem.
Minha principal questão é como eu posso usar o SSE na primeira passagem no caso geral? O seguinte linksimd-prefixo-soma-on-intel-cpu mostra uma melhoria para bytes, mas não para tipos de dados de 32 bits.
A razão pela qual o caso especial é chamado de especial é que ele requer que o array esteja em um formato especial. Por exemplo, vamos supor que havia apenas 16 elementos de uma matriza
de flutuadores. Então, se o array foi reorganizado assim (array de structs para struct de arrays):
a[0] a[1] ...a[15] -> a[0] a[4] a[8] a[12] a[1] a[5] a[9] a[13]...a[3] a[7] a[11] a[15]
As somas verticais SSE podem ser usadas em ambos os passes. No entanto, isso só seria eficiente se as matrizes já estivessem no formato especial e a saída pudesse ser usada no formato especial. Caso contrário, dispendioso rearranjo teria que ser feito em entrada e saída, o que tornaria muito mais lento do que o caso geral.
Talvez eu deva considerar um algoritmo diferente para a soma do prefixo (por exemplo, uma árvore binária)?
Código para o caso geral:
void prefix_sum_omp_sse(double a[], double s[], int n) {
double *suma;
#pragma omp parallel
{
const int ithread = omp_get_thread_num();
const int nthreads = omp_get_num_threads();
#pragma omp single
{
suma = new double[nthreads + 1];
suma[0] = 0;
}
double sum = 0;
#pragma omp for schedule(static) nowait //first parallel pass
for (int i = 0; i<n; i++) {
sum += a[i];
s[i] = sum;
}
suma[ithread + 1] = sum;
#pragma omp barrier
#pragma omp single
{
double tmp = 0;
for (int i = 0; i<(nthreads + 1); i++) {
tmp += suma[i];
suma[i] = tmp;
}
}
__m128d offset = _mm_set1_pd(suma[ithread]);
#pragma omp for schedule(static) //second parallel pass with SSE as well
for (int i = 0; i<n/4; i++) {
__m128d tmp1 = _mm_load_pd(&s[4*i]);
tmp1 = _mm_add_pd(tmp1, offset);
__m128d tmp2 = _mm_load_pd(&s[4*i+2]);
tmp2 = _mm_add_pd(tmp2, offset);
_mm_store_pd(&s[4*i], tmp1);
_mm_store_pd(&s[4*i+2], tmp2);
}
}
delete[] suma;
}