SSE的并行前缀（累加）总和

Question

我正在寻找有关如何使用SSE进行并行前缀求和的建议。我有兴趣在整数，浮点数或双精度数组上执行此操作。

我提出了两种解决方案。特殊情况和一般情况。在这两种情况下，解决方案都需要与OpenMP并行进行两次遍历。对于特殊情况，我在两次通过时都使用SSE。对于一般情况，我仅在第二遍使用它。

我的主要问题是，在一般情况下，如何在第一次通过时使用SSE？ 以下链接simd-prefix-sum-intel-cpu显示了对字节的改进，但对32位数据类型却没有。

特殊情况被称为特殊情况的原因是，它要求数组采用特殊格式。例如，假设a浮点数组只有16个元素。然后，如果数组是这样重新排列的（从结构数组到数组结构）：

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;

}