如何在 C 中为两个 for 循环编写 SSE 指令？

Question

使用库“immintrin.h”，我可以为简单的 for 循环和操作编写 SSE 指令。但是，如何为显示的语句编写 SSE 指令？

for (int i =0; i<n; i++){
   for (int j=0; j<n; j++) {
     x[i] += a[i] + a[j];
}}

x 和 a 是使用 _mm_malloc() 初始化的 float*。内存访问模式可用作 __m128 和 4 字节的展开策略。

对不起，如果我不太清楚，但就像

for (int i = 0; i < vecsize; i+=4) {
        __m128 a  = _mm_load_ps(a+i);
        __m128 x  = _mm_add_ps(x,a);
        _mm_store_ps(x+i, x);
    }

（仅适用于 1 个循环），我想要上面显示的循环类似的东西。

编辑：我（EricPostpischil）正在从评论中注入此文本，因为它对问题陈述重要。作者 NeilDA 应该对此进行扩展：

...在我的程序中，'a' 总是在变化，因此我希望 'x' 随之变化。

我已经做到了！我提交了答案..

Answer 1

这只是部分答案，但太长/太详细，无法评论。

我怀疑你的问题写得是否正确。如图所示，对于每个x[i]，它添加a[i]n次，每个a[j]添加一次，对于0≤jn。所以相当于：

sum = 0;
for (j = 0; j < n; ++j)
    sum += a[j];
for (i = 0; i < n; ++i)
    x[i] += n*a[i] + sum;

这将使用比其他可能的数组操作更简单的 SSE 代码来实现。当然，像上面那样简单地重写它会产生比原始公式更快的代码。

Answer 2

__m128 *mx = (__m128*)x;
__m128 *ma = (__m128*)a;
__m128 temp_a;

for (int i = 0; i < (n>>2); ++i) {
    for (int j = 0; j < (n>>2); ++j) {
        temp_a = _mm_add_ps(*(ma+i), *(ma+j));
        *mx = _mm_add_ps(*mx, temp_a);
    }
    mx++;
}

当然，您必须确保n 是 4 的倍数，并确保 x 初始化为 0，以便累加正确。

Answer 3

__m128 x_v, a_v, aj_v;
for (int i = 0; i < (vec_size); i+=4) {
            x_v = _mm_load_ps(x + i);
            a_v = _mm_load_ps(a+i);
        for (int j = 0; j < N; j++) {
            aj_v = _mm_set1_ps(a[j]);
            x_v = _mm_add_ps(x_v, _mm_add_ps(aj_v, a_v));
        }
            _mm_store_ps(x + i, x_v);
    }

我不知道它是否可以进一步改进，但这将我的时间从 0.17 秒减少到 0.04 秒：D 任何改进它的评论或方法都会很棒！

Answer 4

变换：

    for (i = 0; i < n; i++)
        for (j = 0; j < n; j++)
            x[i] += a[i] + a[j];

进入：

    for (i = 0; i < n; i++)
        x[i] += n*a[i] + sum(a));

请参阅下面代码中的f_sse()：

#include <stdio.h>
#include <string.h>
#include <immintrin.h>

enum {
    N = 4,
};

float x[N], a[N] = { .1, .2, .3, .4 }, y[N];

void f(float *x, float *a, int n)
{
    int i, j;
    for (i = 0; i < n; i++)
        for (j = 0; j < n; j++)
            x[i] += a[i] + a[j];
}

float array_sum(float *a, int n)
{
    /* could be vectorized as well */
    int i;
    float s;
    for (s = 0, i = 0; i < n; i++)
        s += a[i];
    return s;
}

void f_sse(float *x, float *a, int n)
{
    int i, l;
    float t;
    __m128 sum_a, n_vec;
    t = array_sum(a, n);
    sum_a = _mm_set1_ps(t);
    n_vec = _mm_set1_ps(n);
    l = n / 4;
    for (i = 0; i < l; i += 4) {
        __m128 ai, xi;
        ai = _mm_load_ps(a + i);
        xi = _mm_load_ps(x + i);
        ai = _mm_mul_ps(ai, n_vec);
        ai = _mm_add_ps(ai, sum_a);
        xi = _mm_add_ps(xi, ai);
        _mm_store_ps(x + i, xi);
    }
}

int main()
{
    int i, r;
    f(x, a, N);
    f_sse(y, a, N);
    r = memcmp(x, y, N);
    if (r == 0)
        return 0;
    printf("x: { ");
    for (i = 0; i < N; i++)
        printf("%f, ", x[i]);
    printf("}\n");
    printf("y: { ");
    for (i = 0; i < N; i++)
        printf("%f, ", y[i]);
    printf("}\n");
    return 3;
}

由于您声明a 正在同时更新，因此您无法进行上述循环转换。您需要有意识地决定何时从a 获取更新。它永远不会匹配原始的非矢量化版本，因为我们一次加载 4 个浮点数。

在内循环中重新加载a[j]：

void f_sse(float *x, float *a, int n)
{
    int i, j, l;
    l = n / 4;
    for (i = 0; i < l; i += 4) {
        __m128 ai, xi;
        ai = _mm_load_ps(a + i);
        xi = _mm_load_ps(x + i);
        for (j = 0; j < n; j++) {
            /* re-load a[j] to get updates */
            xi = _mm_add_ps(xi, _mm_add_ps(ai, _mm_set1_ps(((volatile float *)a)[j])));
        }
        _mm_store_ps(x + i, xi);
    }

}

在内循环中重新加载a[j] 和(a + i)：

void f_sse(float *x, float *a, int n)
{
    int i, j, l;
    l = n / 4;
    for (i = 0; i < l; i += 4) {
        __m128 ai, xi;
        xi = _mm_load_ps(x + i);
        for (j = 0; j < n; j++) {
            /* re-load (a + i) to get updates */
            ai = _mm_load_ps(a + i);
            /* re-load a[j] to get updates */
            xi = _mm_add_ps(xi, _mm_add_ps(ai, _mm_set1_ps(((volatile float *)a)[j])));
        }
        _mm_store_ps(x + i, xi);
    }
}