SIMD 优化难题

Posted 2023-02-16

【中文标题】SIMD 优化难题

【英文标题】：SIMD optimization puzzle

【发布时间】：2010-10-28 21:29:13

【问题描述】：

我想使用 SIMD（SSE2 等）优化以下功能：

int64_t fun(int64_t N, int size, int* p)

    int64_t sum = 0;
    for(int i=1; i<size; i++)
       sum += (N/i)*p[i];
    return sum;

这似乎是一个非常可向量化的任务，只是所需的指令不在那里......

我们可以假设 N 非常大（10^12 到 10^18）并且 size~sqrt(N)。我们还可以假设 p 只能取 -1、0 和 1 的值；所以我们不需要真正的乘法，如果我们能以某种方式计算 N/i，则 (N/i)*p[i] 可以用四个指令（pcmpgt、pxor、psub、pand）来完成。

【问题讨论】：

@user434507 我想如果你听起来不那么傲慢，你会得到更多的回应。

@srean 如果您不知道如何解决，请直说。 :)

我不太了解 SIMD 约束，但我了解一些数学。如果是i/3，你能做到吗？也就是说，如果您需要计算许多不同的数字，每个数字除以一个常数，而不是相反，那是否可以按照您想要的方式并行化？

整数除法感觉像是近似值。如果 N/i 的计算是估计的，因此偶尔会偏离一两个单位，是否可以？还是必须在现场？

@user434507 哈哈哈 随便你的船夫

【参考方案1】：

这是尽可能接近矢量化该代码的方法。我真的不指望它会更快。我只是在尝试编写 SIMD 代码。

#include <stdint.h>

int64_t fun(int64_t N, int size, const int* p)

    int64_t sum = 0;
    int i;
    for(i=1; i<size; i++) 
        sum += (N/i)*p[i];
    
    return sum;


typedef int64_t v2sl __attribute__ ((vector_size (2*sizeof(int64_t))));

int64_t fun_simd(int64_t N, int size, const int* p)

    int64_t sum = 0;
    int i;
    v2sl v_2 =  2, 2 ;
    v2sl v_N =  N, N ;
    v2sl v_i =  1, 2 ;
    union  v2sl v; int64_t a[2];  v_sum;

    v_sum.a[0] = 0;
    v_sum.a[1] = 0;
    for(i=1; i<size-1; i+=2) 
        v2sl v_p =  p[i], p[i+1] ;
        v_sum.v += (v_N / v_i) * v_p;
        v_i += v_2;
    
    sum = v_sum.a[0] + v_sum.a[1];
    for(; i<size; i++) 
        sum += (N/i)*p[i];
    
    return sum;


typedef double v2df __attribute__ ((vector_size (2*sizeof(double))));

int64_t fun_simd_double(int64_t N, int size, const int* p)

    int64_t sum = 0;
    int i;
    v2df v_2 =  2, 2 ;
    v2df v_N =  N, N ;
    v2df v_i =  1, 2 ;
    union  v2df v; double a[2];  v_sum;

    v_sum.a[0] = 0;
    v_sum.a[1] = 0;
    for(i=1; i<size-1; i+=2) 
        v2df v_p =  p[i], p[i+1] ;
        v_sum.v += (v_N / v_i) * v_p;
        v_i += v_2;
    
    sum = v_sum.a[0] + v_sum.a[1];
    for(; i<size; i++) 
        sum += (N/i)*p[i];
    
    return sum;


#include <stdio.h>

static const int test_array[] = 
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0,
 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0, 1, 0, -1, 0
;
#define test_array_len (sizeof(test_array)/sizeof(int))

#define big_N (1024 * 1024 * 1024)

int main(int argc, char *argv[]) 
    int64_t res1;
    int64_t res2;
    int64_t res3;
    v2sl a =  123, 456 ;
    v2sl b =  100, 200 ;
    union  v2sl v; int64_t a[2];  tmp;

    a = a + b;
    tmp.v = a;
    printf("a =  %ld, %ld \n", tmp.a[0], tmp.a[1]);

    printf("test_array size = %zd\n", test_array_len);

    res1 = fun(big_N, test_array_len, test_array);
    printf("fun() = %ld\n", res1);

    res2 = fun_simd(big_N, test_array_len, test_array);
    printf("fun_simd() = %ld\n", res2);

    res3 = fun_simd_double(big_N, test_array_len, test_array);
    printf("fun_simd_double() = %ld\n", res3);

    return 0;

【参考方案2】：

1/x 的导数是-1/x^2，这意味着随着x 变大，N/x==N/(x + 1)。

对于N/x 的已知值（我们称其为r），我们可以确定x 的下一个值（我们称该值为x'，这样N/x'<r：

x'= N/(r - 1)

因为我们处理的是整数：

x'= ceiling(N/(r - 1))

所以，循环变成了这样：

int64_t sum = 0;   
int i=1; 
int r= N;
while (i<size)

    int s= (N + r - 1 - 1)/(r - 1);

    while (i<s && i<size)
    
        sum += (r)*p[i];
        ++i;
    

    r= N/s;

return sum;   

对于足够大的 N，N/i 将有许多相同值的运行。当然，如果你不小心，你会被零除。

【讨论】：

这是一个非常敏锐的观察，不幸的是，它对我没有帮助，因为 size ~ sqrt(N) 并且在那之后你不会得到任何重复。

嘿，你说得对。我从不费心去弄清楚“足够大的 N”。

【参考方案3】：

我建议您使用浮点 SIMD 操作来执行此操作 - 单精度或双精度，具体取决于您的精度要求。使用 SSE 从 int 到 float 或 double 的转换相对较快。

【讨论】：

使用浮点是有问题的，因为浮点 SIMD 操作有 52 位分辨率，而我的 N 可能比那个大...

@user434507：这取决于您的精度要求 - 您实际上并没有说明您需要多少精度。鉴于 N/i 是一个截断整数除法，那么您似乎并没有过度关注结果中的绝对数值精度？

@user434507：在这种情况下，我建议使用标量代码计算前 M 项，直到使用 SIMD 浮点数可以管理 N/i。然后，您可以处理点块并在每个块之后定期更新标量 64 位整数和。

【参考方案4】：

成本集中在计算除法上。 SSE2 中没有用于整数除法的操作码，因此您必须自己实现除法算法，一点一点。我认为这不值得付出努力：SSE2 允许您并行执行两个实例（您使用 64 位数字，而 SSE2 寄存器是 128 位）但我发现手工除法算法可能至少是比 CPU 慢两倍 idiv 操作码。

（顺便问一下，你是在 32 位还是 64 位模式下编译？后者更适合 64 位整数。）

减少部门总数似乎是一种更有希望的方法。有人可能会注意到，对于正整数 x 和 y，则 floor(x/(2y)) = floor(floor(x/y)/2)。在 C 术语中，一旦您计算了 N/i（截断除法），那么您只需将其右移一位即可获得 N/(2*i)。如果使用得当，这会使您的一半部门几乎免费（“正确”还包括以不会对缓存造成严重破坏的方式访问数十亿 p[i] 值，因此看起来并不容易）。