Я ищу несколько советов о том, как сделать параллельную сумму префикса с SSE. Я заинтересован в том, чтобы делать это на массиве целых чисел, чисел с плавающей точкой или двойных чисел.

Я придумал два решения. Особый случай и общий случай. В обоих случаях решение проходит через массив в два прохода параллельно с OpenMP. Для особого случая я использую SSE на обоих проходах. Для общего случая я использую его только на втором проходе.

Мой главный вопрос - как я могу использовать SSE на первом проходе в общем случае? Следующая ссылкаSIMD-префиксы сумма-на-интел процессорам показать улучшение для байтов, но не для 32-битных типов данных.

Причина, по которой специальный случай называется специальным, заключается в том, что он требует, чтобы массив был в специальном формате. Например, предположим, что в массиве было только 16 элементов.aплавает. Затем, если массив был перестроен следующим образом (массив структур в структуру массивов):

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]

Вертикальные суммы SSE могут использоваться на обоих проходах. Однако это будет эффективно только в том случае, если массивы уже были в специальном формате и выходные данные можно будет использовать в специальном формате. В противном случае пришлось бы выполнять дорогостоящую перестройку как на входе, так и на выходе, что делало бы это намного медленнее, чем в общем случае.

Может быть, я должен рассмотреть другой алгоритм для суммы префикса (например, двоичное дерево)?

Код для общего случая:

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