当前位置： 文江博客 话题详情

SIMD 优化难题

发布于 2024-09-29 10:45:46 字数 384 浏览 0 评论 0 原文

我想使用 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) 和大小~sqrt(N)。我们还可以假设p只能取-1、0和1的值；所以我们不需要真正的乘法，如果我们能以某种方式计算 N/i，则 (N/i)*p[i] 可以用四个指令（pcmpgt、pxor、psub、pand）来完成。

原文

I Want to optimize the following function using SIMD (SSE2 & such):

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

This seems like an eminently vectorizable task, except that the needed instructions just aren't there ...

We can assume that N is very large (10^12 to 10^18) and size~sqrt(N). We can also assume that p can only take values of -1, 0, and 1; so we don't need a real multiplication, the (N/i)*p[i] can be done with four instructions (pcmpgt, pxor, psub, pand), if we could just somehow compute N/i.

如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论， 你可以免费 注册 一个本站的账号。

评论（4）

作妖 2024-10-06 10:45:46

这是我能对代码进行矢量化的最接近的结果。我真的不希望它会更快。我只是尝试编写 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;
}

This is as close as I could get to vectorizing that code. I don't really expect it to be faster. I was just trying my hand at writting SIMD code.

#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;
}
回复 原文
筱果果 2024-10-06 10:45:46

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

对于已知的 N/x 值（我们称该值 r），我们可以确定 x 的下一个值（我们称该值x' 使得 N/x'：由于

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 有很多次相同的值。当然，如果你不小心，你就会被零除。

The derivative of 1/x is -1/x^2, which means as x gets bigger, N/x==N/(x + 1).

For a known value of N/x (let's call that value r), we can determine the next value of x (let's call that value x' such that N/x'<r:

x'= N/(r - 1)

And since we are dealing with integers:

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

So, the loop becomes something like this:

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;

For sufficiently large N, you will have many many runs of identical values for N/i. Granted, you will hit a divide by zero if you aren't careful.

回复 原文
情徒 2024-10-06 10:45:46

我建议您使用浮点 SIMD 运算来执行此操作 - 根据您的精度要求，可以是单精度也可以是双精度。使用 SSE 从 int 到 float 或 double 的转换相对较快。

I suggest you do this with floating point SIMD operations - either single or double precision depending on your accuracy requirements. Conversion from int to float or double is relatively fast using SSE.

回复 原文
写下不归期 2024-10-06 10:45:46

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

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

减少总除法数看起来是一种更有前途的方法。人们可能会注意到，对于正整数x和y，则floor(x/(2y)) = Floor(floor(x/y)/2)< /em>.在 C 术语中，一旦计算出 N/i（截断除法），只需将其右移一位即可获得 N/(2*i)。使用得当，这使得一半的分区几乎是免费的（“正确”还包括以不会对缓存造成严重破坏的方式访问数十亿个 p[i] 值，因此它不会看起来很容易）。

The cost is concentrated in computing the divisions. There is no opcode in SSE2 for integral divisions, so you would have to implement a division algorithm yourself, bit by bit. I do not think it would be worth the effort: SSE2 allow you to perform two instances in parallel (you use 64-bit numbers, and SSE2 registers are 128-bit) but I find it likely that a handmade division algorithm would be at least twice as slow as the CPU idiv opcode.

(By the way, do you compile in 32-bit or 64-bit mode ? The latter will be more comfortable with 64-bit integers.)

Reducing the overall number of divisions looks like a more promising way. One may note that for positive integers x and y, then floor(x/(2y)) = floor(floor(x/y)/2). In C terminology, once you have computed N/i (truncated division) then you just have to shift it right by one bit to obtain N/(2*i). Used properly, this makes half of your divisions almost free (that "properly" also includes accessing the billions of p[i] values in a way which does not wreak havoc with the caches, so it does not seem very easy).

回复 原文
~没有更多了~

关于作者

温柔戏命师

暂无简介

0 文章
0 评论
22 人气
发私信

相关话题

更多

热门标签

操作系统 程序设计 IT运维 Linux系统管理 JavaScript 服务器应用 solaris C/C++ PHP Shell BSD Vue.js aix Oracle Python HTML 系统管理 HTML5 CSS 前端
更多

推荐作者

qq_FjTq5B

文章 0 评论 0

18273202778

文章 0 评论 0

WordPress小学生

文章 0 评论 0

〃温暖了心ぐ

文章 0 评论 0

迷乱花海

文章 0 评论 0

niuniu

文章 0 评论 0

更多

友情链接

文江博客
    我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站，即表示您同意使用 Cookies 和您的相关数据。
    原文