Intel cpu上的SIMD前缀和

Question

我需要实现前缀和算法，并且需要它尽可能快。
Ex：

[3, 1,  7,  0,  4,  1,  6,  3]

应给出：

[3, 4, 11, 11, 15, 16, 22, 25]

有没有办法使用SSE SIMD CPU指令做到这一点？

我的第一个想法是递归地并行求和每对，直到所有求和都像下面这样计算出来！

//in parallel do 
for (int i = 0; i < z.length; i++) {
    z[i] = x[i << 1] + x[(i << 1) + 1];
}

为了使算法更清楚一点，z不是最终输出，而是用于计算输出。

int[] w = computePrefixSum(z);
for (int i = 1; i < ouput.length; i++) {
    ouput[i] = (i % 2 == 0) ? (x[i] + ouput[i - 1]) :  w[(i - 1) >> 1];
}

Answer 1

前缀和可以并行计算，它实际上是GPU编程的基本算法之一。如果您在Intel处理器上使用SIMD扩展，我不确定并行执行是否会给您带来很大好处，但请看一看nvidia关于实现并行前缀和的这篇文章（只看算法，忽略CUDA）：CUDA的并行前缀和（扫描）。

Answer 2

对于较大的寄存器长度和较小的和，可以利用一些较小的并行性。例如，将1个字节的16个值相加（恰好适合一个sse寄存器）只需要log₂16个加法和相等的移位次数<不多，但比15次依赖的添加和额外的内存访问更快。

__m128i x = _mm_set_epi8(3,1,7,0,4,1,6,3,3,1,7,0,4,1,6,3);
x = _mm_add_epi8(x, _mm_srli_si128(x, 1));
x = _mm_add_epi8(x, _mm_srli_si128(x, 2));
x = _mm_add_epi8(x, _mm_srli_si128(x, 4));
x = _mm_add_epi8(x, _mm_srli_si128(x, 8));

// x == 3, 4, 11, 11, 15, 16, 22, 25, 28, 29, 36, 36, 40, 41, 47, 50

如果求和较长，则可以通过利用指令级并行性和利用指令重新排序来隐藏依赖关系。

编辑：类似

__m128i x0 = _mm_set_epi8(3,1,7,0,4,1,6,3,3,1,7,0,4,1,6,3);
__m128i x1 = _mm_set_epi8(3,1,7,0,4,1,6,3,3,1,7,0,4,1,6,3);
__m128i x2 = _mm_set_epi8(3,1,7,0,4,1,6,3,3,1,7,0,4,1,6,3);
__m128i x3 = _mm_set_epi8(3,1,7,0,4,1,6,3,3,1,7,0,4,1,6,3);

__m128i mask = _mm_set_epi8(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);

x0 = _mm_add_epi8(x0, _mm_srli_si128(x0, 1));
x1 = _mm_add_epi8(x1, _mm_srli_si128(x1, 1));
x2 = _mm_add_epi8(x2, _mm_srli_si128(x2, 1));
x3 = _mm_add_epi8(x3, _mm_srli_si128(x3, 1));

x0 = _mm_add_epi8(x0, _mm_srli_si128(x0, 2));
x1 = _mm_add_epi8(x1, _mm_srli_si128(x1, 2));
x2 = _mm_add_epi8(x2, _mm_srli_si128(x2, 2));
x3 = _mm_add_epi8(x3, _mm_srli_si128(x3, 2));

x0 = _mm_add_epi8(x0, _mm_srli_si128(x0, 4));
x1 = _mm_add_epi8(x1, _mm_srli_si128(x1, 4));
x2 = _mm_add_epi8(x2, _mm_srli_si128(x2, 4));
x3 = _mm_add_epi8(x3, _mm_srli_si128(x3, 4));

x0 = _mm_add_epi8(x0, _mm_srli_si128(x0, 8));
x1 = _mm_add_epi8(x1, _mm_srli_si128(x1, 8));
x2 = _mm_add_epi8(x2, _mm_srli_si128(x2, 8));
x3 = _mm_add_epi8(x3, _mm_srli_si128(x3, 8));

x1 = _mm_add_epi8(_mm_shuffle_epi8(x0, mask), x1);
x2 = _mm_add_epi8(_mm_shuffle_epi8(x1, mask), x2);
x3 = _mm_add_epi8(_mm_shuffle_epi8(x2, mask), x3);

Answer 3

我所知道的最快的并行前缀和算法是在两个过程中并行运行和，在第二个过程中也使用SSE。

在第一步中，并行计算部分和，并存储每个部分和的总和。在第二步中，将前一部分和的总和添加到下一部分和。您可以使用多个线程（例如使用OpenMP）并行运行这两个过程。第二步还可以使用SIMD，因为每个部分和都要加一个常数值。

假设阵列中有n个元素、m个核心，SIMD宽度为w，则时间成本应为

n/m + n/(m*w) = (n/m)*(1+1/w)

由于首次传递不使用SIMD，时间成本将始终大于n/m

例如，对于SIMD\U宽度为4的四个核（使用SSE的四个32位浮点），成本将为5n/16。或比时间成本为n的顺序码快3.2倍左右。使用超线程，速度会更快。

在特殊情况下，也可以在第一次通过时使用SIMD。那么时间成本就是

2*n/(m*w)

我发布了通用案例的代码，该代码使用OpenMP进行线程处理，并在以下链接中讨论了有关特殊案例的详细信息并行前缀累积总和与sse

编辑：我设法找到了第一遍的SIMD版本，速度大约是顺序代码的两倍。现在，我的四芯常春藤桥系统总共提升了7个左右。

编辑：对于较大的数组，一个问题是在第一次传递后，大多数值都已从缓存中逐出。我提出了一个解决方案，它在一个块中并行运行，但在每个块中串行运行。chunk\u大小是一个需要调整的值。例如，我将其设置为1MB=256K浮点。现在，当值仍在二级缓存中时，完成第二次传递。这样做可以大大改进大型阵列。

这是SSE的代码。AVX代码的速度差不多，所以我没有在这里发布。进行前缀和的函数是scan\u omp\u SSEp2\u SSEp1\u chunk。向其传递一个浮点数组，并用累积和填充数组。

__m128 scan_SSE(__m128 x) {
    x = _mm_add_ps(x, _mm_castsi128_ps(_mm_slli_si128(_mm_castps_si128(x), 4))); 
    x = _mm_add_ps(x, _mm_shuffle_ps(_mm_setzero_ps(), x, 0x40)); 
    return x;
}

float pass1_SSE(float *a, float *s, const int n) {
    __m128 offset = _mm_setzero_ps();
    #pragma omp for schedule(static) nowait
    for (int i = 0; i < n / 4; i++) {
        __m128 x = _mm_load_ps(&a[4 * i]);
        __m128 out = scan_SSE(x);
        out = _mm_add_ps(out, offset);
        _mm_store_ps(&s[4 * i], out);
        offset = _mm_shuffle_ps(out, out, _MM_SHUFFLE(3, 3, 3, 3));
    }
    float tmp[4];
    _mm_store_ps(tmp, offset);
    return tmp[3];
}

void pass2_SSE(float *s, __m128 offset, const int n) {
    #pragma omp for schedule(static)
    for (int i = 0; i<n/4; i++) {
        __m128 tmp1 = _mm_load_ps(&s[4 * i]);
        tmp1 = _mm_add_ps(tmp1, offset);
        _mm_store_ps(&s[4 * i], tmp1);
    }
}

void scan_omp_SSEp2_SSEp1_chunk(float a[], float s[], int n) {
    float *suma;
    const int chunk_size = 1<<18;
    const int nchunks = n%chunk_size == 0 ? n / chunk_size : n / chunk_size + 1;
    //printf("nchunks %d\n", nchunks);
    #pragma omp parallel
    {
        const int ithread = omp_get_thread_num();
        const int nthreads = omp_get_num_threads();

        #pragma omp single
        {
            suma = new float[nthreads + 1];
            suma[0] = 0;
        }

        float offset2 = 0.0f;
        for (int c = 0; c < nchunks; c++++) {
            const int start = c*chunk_size;
            const int chunk = (c + 1)*chunk_size < n ? chunk_size : n - c*chunk_size;
            suma[ithread + 1] = pass1_SSE(&a[start], &s[start], chunk);
            #pragma omp barrier
            #pragma omp single
            {
                float tmp = 0;
                for (int i = 0; i < (nthreads + 1); i++) {
                    tmp += suma[i];
                    suma[i] = tmp;
                }
            }
            __m128 offset = _mm_set1_ps(suma[ithread]+offset2);
            pass2_SSE(&s[start], offset, chunk);
            #pragma omp barrier
            offset2 = s[start + chunk-1];
        }
    }
    delete[] suma;
}