c++ SSE2的内在特性

C++ SSE2 intrinsics

本文关键字:SSE2 c++      更新时间:2023-10-16

我刚刚了解到有一种方法可以使用intrinsic实现一些并行化。我发现下面的代码,并想通过它,但我可以理解很多。我试图使运算达到单精度,但我怎么能做到呢?

#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 = %lfn", 4*res);
        return 0;
}

我正在考虑将一切从double改为float,并调用正确的内在函数,例如,而不是_mm_set_pd1 -> _mm_set_ps1。我不知道这是否会使程序从双精度变成单精度。

我尝试如下,但我得到一个分割错误

#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 = %lfn", 4*res);
        return 0;
}

需要更多的修复:

  • x1需要用4个元素声明
  • 第二个for循环需要增加4(这就是导致segfault的原因)。
  • 需要4次赋值给x1数组。

这些变化都是因为单精度将4个值装入16字节的矢量寄存器中,而双精度只装入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 = %lfn", 4*res);
   return 0;
}

滚筒滚……

$ ./foo
pi = 3.141597

关于malloc()的使用。我认为大多数实现将返回在SSE加载和存储所需的16字节边界上对齐的内存,但这可能无法保证,因为__m128不是C/c++类型(它保证为"正常"类型对齐)。使用memalign()posix_memalign()会更安全。