C++ SSE运算例子

参考：http://blog.163.com/chenqneu@126/blog/static/45738484200781493846369/

C++ SSE浮点运算实例(__m128)

实验环境VC++ 8.0
（1）
#include "stdafx.h"
#include "omp.h"
#include <xmmintrin.h>
#include <iostream>
using namespace std;
const int N=120;

int _tmain(int argc, _TCHAR* argv[])
{
    __m128 *p1,*p2,*p3;
    float f1[N],f2[N],f3[N];
    for(int i=0;i<N;i++){
        f1[i]=i+0.12;
        f2[i]=i+0.16;
    };
    p1=(__m128*)f1;
    p2=(__m128*)f2;
    p3=(__m128*)f3;
     for(int i=0;i<N;i+=4){
        *p3=_mm_add_ps(*p1,*p2);
        p3++;
        p2++;
        p1++;
    }
    for(int i=0;i<N;i++)
        cout<<f3[i]<<endl;
    return 0;
}
程序执行时错误，数据读取时异常，发生访问冲突。
（2）
#include "stdafx.h"
#include "omp.h"
#include <xmmintrin.h>
#include <iostream>
using namespace std;
const int N=120;

int _tmain(int argc, _TCHAR* argv[])
{
    __m128 *p1,*p2,*p3;
    __declspec(align(16)) float f1[N],f2[N],f3[N];
    for(int i=0;i<N;i++){
        f1[i]=i+0.12;
        f2[i]=i+0.16;
    };
    p1=(__m128*)f1;
    p2=(__m128*)f2;
    p3=(__m128*)f3;
    for(int i=0;i<N;i+=4){
        *p3=_mm_add_ps(*p1,*p2);
        p3++;
        p2++;
        p1++;
    }
    for(int i=0;i<N;i++)
        cout<<f3[i]<<endl;
    return 0;
}
程序正常运行，并且结果正确，原因是SSE浮点运算要求按16字（128位）对齐，所以 __declspec(align(16)) float f1[N],f2[N],f3[N];便可以了。
（3）
#include "stdafx.h"
#include "omp.h"
#include <xmmintrin.h>
#include <iostream>
using namespace std;
const int N=120;

int _tmain(int argc, _TCHAR* argv[])
{
    __m128 *p1,*p2,*p3;
    __declspec(align(16)) float *pf1,*pf2,*pf3;
    pf1=(__declspec(align(16))float*)_mm_malloc(sizeof(float)*N,16);
    pf2=(__declspec(align(16))float*)_mm_malloc(sizeof(float)*N,16);
    pf3=(__declspec(align(16))float*)_mm_malloc(sizeof(float)*N,16);
    for(int i=0;i<N;i++){
        pf1[i]=i+0.12;
        pf2[i]=i+0.16;
    };
    p1=(__m128*)pf1;
    p2=(__m128*)pf2;
    p3=(__m128*)pf3;
    for(int i=0;i<N;i+=4){
        *p3=_mm_add_ps(*p1,*p2);
        p3++;
        p2++;
        p1++;
    }
    for(int i=0;i<N;i++)
        cout<<pf3[i]<<endl;
    _mm_free(pf1);
    _mm_free(pf2);
    _mm_free(pf3);
    return 0;
}
这个程序也是正确的，通过动态分配内存来实现，_mm_malloc(sizeof,16)和_mm_free()。
注意：下面的动态内存申请方法是有问题的：
pf1=(__declspec(align(16))(float*))_mm_malloc(sizeof(float)*N,16);
提示：error C2440: “=”: 无法从“int (__cdecl *)(float *)”转换为“float *”
pf1=(__declspec(align(16)))(float*)_mm_malloc(sizeof(float)*N,16);
提示：error C2440: “=”: 无法从“int”转换为“float *”
（4）
#include "stdafx.h"
#include "omp.h"
#include <xmmintrin.h>
#include <iostream>
using namespace std;
const int N=120;

int _tmain(int argc, _TCHAR* argv[])
{
    __m128 m1,m2,m3;
    __declspec(align(16)) float *pf1,*pf2,*pf3;
    pf1=(__declspec(align(16))float*)_mm_malloc(sizeof(float)*N,16);
    pf2=(__declspec(align(16))float*)_mm_malloc(sizeof(float)*N,16);
    pf3=(__declspec(align(16))float*)_mm_malloc(sizeof(float)*N,16);
    for(int i=0;i<N;i++){
        pf1[i]=i+0.12;
        pf2[i]=i+0.16;
    };
    for(int i=0;i<N;i+=4){
        m1=_mm_load_ps(pf1+i);
        m2=_mm_load_ps(pf2+i);
        m3=_mm_add_ps(m1,m2);
        _mm_store_ps(pf3+i,m3);
    }
    for(int i=0;i<N;i++)
        cout<<pf3[i]<<endl;
    _mm_free(pf1);
    _mm_free(pf2);
    _mm_free(pf3);
    return 0;
}
使用_mm_load_ps(float *)和_mm_store_ps(float *,__m128)也可以用来直接加载和保存4组浮点数，替换__m128*来进行数据引用。当然，这里仍旧要求浮点数组是按16字对齐的。
（5）
#include "stdafx.h"
#include "omp.h"
#include <xmmintrin.h>
#include <iostream>
using namespace std;
const int N=120;

int _tmain(int argc, _TCHAR* argv[])
{
    __m128 m1,m2,m3;
    __declspec(align(16)) float *pf1,*pf2,*pf3;
    pf1=(float*)malloc(sizeof(float)*N);
    pf2=(float*)malloc(sizeof(float)*N);
    pf3=(float*)malloc(sizeof(float)*N);
    for(int i=0;i<N;i++){
        pf1[i]=i+0.12;
        pf2[i]=i+0.16;
    };
    for(int i=0;i<N;i+=4){
        m1=_mm_load_ps(pf1+i);
        m2=_mm_load_ps(pf2+i);
        m3=_mm_add_ps(m1,m2);
        _mm_store_ps(pf3+i,m3);
    }
    for(int i=0;i<N;i++)
        cout<<pf3[i]<<endl;
    free(pf1);
    free(pf2);
    free(pf3);
    return 0;
}
这里的程序就会出现运行时内存访问错误，原因就是浮点数组没有16字对齐时，不能用_mm_load_ps()和_mm_store_ps()，而是需要用 _mm_loadu_ps，_mm_storeu_ps()来处理，不过因为没对齐，所以处理速度会慢一些。
（6）
#include "stdafx.h"
#include "omp.h"
#include <xmmintrin.h>
#include <iostream>
using namespace std;
const int N=120;

int _tmain(int argc, _TCHAR* argv[])
{
    __m128 m1,m2,m3;
    __declspec(align(16)) float *pf1,*pf2,*pf3;
    pf1=(float*)malloc(sizeof(float)*N);
    pf2=(float*)malloc(sizeof(float)*N);
    pf3=(float*)malloc(sizeof(float)*N);
    for(int i=0;i<N;i++){
        pf1[i]=i+0.12;
        pf2[i]=i+0.16;
    };
    for(int i=0;i<N;i+=4){
        m1=_mm_loadu_ps(pf1+i);
        m2=_mm_loadu_ps(pf2+i);
        m3=_mm_add_ps(m1,m2);
        _mm_storeu_ps(pf3+i,m3);
    }
    for(int i=0;i<N;i++)
        cout<<dec<<pf3[i]<<endl;
    free(pf1);
    free(pf2);
    free(pf3);
    return 0;


}




下面是几个简单的例子： 

float sumfloat_sse(const float* pbuf, size_t cntbuf)
{
    float s = 0;    // 求和变量.
    size_t i;
    size_t nBlockWidth = 4; // 块宽. SSE寄存器能一次处理4个float.
    size_t cntBlock = cntbuf / nBlockWidth; // 块数.
    size_t cntRem = cntbuf % nBlockWidth;   // 剩余数量.
    __m128 xfsSum = _mm_setzero_ps();   // 求和变量。[SSE] 赋初值0
    __m128 xfsLoad; // 加载.
    const float* p = pbuf;  // SSE批量处理时所用的指针.
    const float* q; // 将SSE变量上的多个数值合并时所用指针.
    // SSE批量处理.
    for(i=0; i<cntBlock; ++i)
    {
        xfsLoad = _mm_load_ps(p);   // [SSE] 加载
        xfsSum = _mm_add_ps(xfsSum, xfsLoad);   // [SSE] 单精浮点紧缩加法
        p += nBlockWidth;
    }
    // 合并.
    q = (const float*)&xfsSum;
    s = q[0] + q[1] + q[2] + q[3];
    // 处理剩下的.
    for(i=0; i<cntRem; ++i)
    {
        s += p[i];
    }
    return s;
}
int main(){
float norm=0.5f,S=0.25,M=0.8f;
float *a;
  *a=2.0f;
float* res=new float[4];
__m128 _res;
const __m128 _s=sse::SET(S);
const __m128 _m=sse::SET(M);
//乘法
_res= sse::MUL(*a,_m);
//    float *res2=(float *)_In;
 memcpy(res,&_res,16);
 printf("res:%.9f,%.9f,%.9f\n",res[0],res[1],res[2]);
 //开平方的导数
_res=sse::RCPSQRT(_s);
memcpy(res,&_res,16);
printf("RCPSQRT:%.9f,%.9f,%.9f\n",res[0],res[1],res[2]);
//除法
_res=sse::RCP(_s);//开平方的导数
memcpy(res,&_res,16);
printf("RCP:%.9f,%.9f,%.9f\n",res[0],res[1],res[2]);
        float *f;
        int  n=10; // 数组的位数
        f = new float[n];
        for(int i=0;i<n;i++)
        f[i] = i;
        
        //    for(int i=0;i<n;i++)  cout << f[i] << "  ";
        
        float  aaa=  sumfloat_sse(f,n);
        printf("result:%f",aaa);
        delete [] f;