我刚刚了解到有一种方法可以使用内在函数实现某种程度的并行化。我找到了以下代码并想通过它,但我能理解很多。我正在尝试使操作采用单精度,但我该怎么做?
#include <stdio.h>
#include <stdlib.h>
#include <xmmintrin.h>
inline double pi_4 (int n){
int i;
__m128d mypart2,x2, b, c, one;
double *x = (double *)malloc(n*sizeof(double));
double *mypart = (double *)malloc(n*sizeof(double));
double sum = 0.0;
double dx = 1.0/n;
double x1[2] __attribute__((aligned(16)));
one = _mm_set_pd1(1.0); // set one to (1,1)
for (i = 0; i < n; i++){
x[i] = dx/2 + dx*i;
}
for (i = 0; i < n; i+=2){
x1[0]=x[i]; x1[1]=x[i+1];
x2 = _mm_load_pd(x1);
b = _mm_mul_pd(x2,x2);
c = _mm_add_pd(b,one);
mypart2 = _mm_div_pd(one,c);
_mm_store_pd(&mypart[i], mypart2);
}
for (i = 0; i < n; i++)
sum += mypart[i];
return sum*dx;
}
int main(){
double res;
res=pi_4(128);
printf("pi = %lf\n", 4*res);
return 0;
}
我正在考虑将所有内容从 double 更改为 float 并调用正确的内部函数,例如,而不是 _mm_set_pd1 -> _mm_set_ps1。我不知道这是否会使程序从 double 变为单精度。
更新
我试过如下但我遇到了段错误
#include <stdio.h>
#include <stdlib.h>
#include <xmmintrin.h>
inline float pi_4 (int n){
int i;
__m128 mypart2,x2, b, c, one;
float *x = (float *)malloc(n*sizeof(float));
float *mypart = (float*)malloc(n*sizeof(float));
float sum = 0.0;
float dx = 1.0/n;
float x1[2] __attribute__((aligned(16)));
one = _mm_set_ps1(1.0); // set one to (1,1)
for (i = 0; i < n; i++){
x[i] = dx/2 + dx*i;
}
for (i = 0; i < n; i+=2){
x1[0]=x[i]; x1[1]=x[i+1];
x2 = _mm_load_ps(x1);
b = _mm_mul_ps(x2,x2);
c = _mm_add_ps(b,one);
mypart2 = _mm_div_ps(one,c);
_mm_store_ps(&mypart[i], mypart2);
}
for (i = 0; i < n; i++)
sum += mypart[i];
return sum*dx;
}
int main(){
float res;
res=pi_4(128);
printf("pi = %lf\n", 4*res);
return 0;
}
最佳答案
还需要一些修复:
x1
需要用 4 个元素声明。- 第二个 for 循环需要递增 4(这是导致段错误的原因)。
- 需要对
x1
数组进行 4 次赋值。
这些变化都是因为单精度将 4 个值打包到一个 16 字节的 vector 寄存器中,而 double 仅将 2 个值打包。我想就是这样:
#include <stdio.h>
#include <stdlib.h>
#include <xmmintrin.h>
inline float pi_4 (int n){
int i;
__m128 mypart2,x2, b, c, one;
float *x = (float *)malloc(n*sizeof(float));
float *mypart = (float*)malloc(n*sizeof(float));
float sum = 0.0;
float dx = 1.0/n;
float x1[4] __attribute__((aligned(16)));
one = _mm_set_ps1(1.0); // set one to (1,1,1,1)
for (i = 0; i < n; i++){
x[i] = dx/2 + dx*i;
}
for (i = 0; i < n; i+=4){
x1[0]=x[i]; x1[1]=x[i+1];
x1[2]=x[i+2]; x1[3]=x[i+3];
x2 = _mm_load_ps(x1);
b = _mm_mul_ps(x2,x2);
c = _mm_add_ps(b,one);
mypart2 = _mm_div_ps(one,c);
_mm_store_ps(&mypart[i], mypart2);
}
for (i = 0; i < n; i++)
sum += mypart[i];
return sum*dx;
}
int main(){
float res;
res=pi_4(128);
printf("pi = %lf\n", 4*res);
return 0;
}
鼓声...
$ ./foo
pi = 3.141597
关于 malloc()
的使用的一句话。我认为大多数实现将根据 SSE 加载和存储的要求返回在 16 字节边界上对齐的内存,但这可能无法保证,因为 __m128 不是 C/C++ 类型(它保证与“正常”类型对齐) .使用 memalign()
或 posix_memalign()
会更安全。
关于C++ SSE2 内在函数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/15260623/