SIMD优化难题
SIMD优化难题
提问于 2010-10-29 05:29:14
我想使用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)*pi就可以用四条指令(pcmpgt，pxor，psub，pand)来完成。
回答 4
Stack Overflow用户
发布于 2010-11-01 21:18:35
这是我能做的最接近向量化的代码了。我真的不期望它会更快。我只是在尝试写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;
}Stack Overflow用户
发布于 2010-11-05 08:06:59
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值。诚然，如果你不小心，你会被零除。
Stack Overflow用户
发布于 2010-10-29 15:09:03
我建议您使用浮点SIMD操作-根据您的精度要求，可以是单精度还是双精度。使用SSE，从整型到浮点型或双精度型的转换相对较快。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/4047319
复制相关文章
相似问题
腾讯云开发者
