SSE练习：单精度浮点数组求和

客官°小女子只卖身不卖艺 2022-07-16 04:49 165阅读 0赞

SSE（Streaming SIMD Extensions）指令是一种SIMD 指令， [ Intrinsics][Intrinsics]函数则是对SSE指令的函数封装，利用C语言形式来调用SIMD指令集，大大提高了易读性和可维护。[Intrinsics][]函数的使用可查看手册[Intel Intrinsics Guide][]。

关于本文实现了单精度浮点数组的求和，切实感受SSE带来的速度提升。本文代码主要来自[\[1\]][1].

首先是不使用任何加速手段的求和函数:

//普通版
    float sumfloat_base(const float *pbuf,unsigned int cntbuf)
    {
    	float s=0;
    	for (unsigned int i=0;i<cntbuf;i++)
    	{
    		s+=pbuf[i];
    	}
    	return s;
    }

在程序优化中有一种经常使用的方法：循环展开。循环展开可以降低循环开销，提高指令级并行性能,此处使用四路展开（测试表明，更多路展开难以带来更快的速度）:

//四路展开
    float sumfloat_base_4loop(const float*pbuf,unsigned int cntbuf)
    {
    	float s=0;
    	float fSum0=0,fSum1=0,fSum2=0,fSum3=0;
    	unsigned int i=0;
    	const float *p=pbuf;
    	for (;i<=cntbuf-4;i+=4)//cntbuf-4不会每次都计算，编译器实现会计算好循环次数！
    	{
    		fSum0+=p[i];
    		fSum1+=p[i+1];
    		fSum2+=p[i+2];
    		fSum3+=p[i+3];
    	}
    	s=fSum0+fSum1+fSum2+fSum3;
    	for (;i<cntbuf;i++)
    	{
    		s+=p[i];
    	}
    	return s;
    }

接着使用SSE 进行加速,由于SSE寄存器位宽128，因此一次能处理4个float类型的数据。SSE指令要求内存地址按16字节对齐，因此在声明缓冲区时使用了\_\_declspec(align(16))。对于动态申请的内存可使用\_aligned\_malloc。

//SSE版 
    float sumfloat_sse(const float *pbuf ,unsigned int cntbuf)
    {
    	float s=0;
    	int nBlockWidth=4;//SSE一次处理4个float
    	int cntBlock=cntbuf/nBlockWidth;
    	int cntRem=cntbuf%nBlockWidth;
    	__m128 fSum=_mm_setzero_ps();//求和变量，初值清零
    	__m128 fLoad;
    	const float*p=pbuf;
    	for (unsigned int i=0;i<cntBlock;i++)
    	{
    		fLoad=_mm_load_ps(p);//加载
    		fSum=_mm_add_ps(fSum,fLoad);//求和
    		p+=nBlockWidth;
    	}
    	const float *q=(const float*)&fSum;
    	s=q[0]+q[1]+q[2]+q[3];		//合并
    	for (int i=0;i<cntRem;i++)//处理尾部剩余数据
    	{
    		s+=p[i];
    	}
    	return s;
    }

本程序中使用的[Intrinsics][]函数为：

\_\_m128\_mm\_load\_ps(float const\*mem\_addr)：

从16字节对齐的内存mem\_addr中加载128位（4个单精度浮点数）到寄存器。对应指令 **movaps** xmm,m128

\_\_m128\_mm\_setzero\_ps(void)：返回一个\_m128类型的全零向量。对应指令：xorps xmm, xmm

\_\_m128\_mm\_add\_ps(\_\_m128a, \_\_m128b)：将4对32位浮点数同时进行相加操作。这4对32位浮点数来自两个128位的存储单元，再把计算结果（相加之和）赋给一个128位的存储单元。对应指令：**addps** xmm,xmm

void\_mm\_store\_ps(float\*mem\_addr, \_\_m128a)：将128位数据存入16字节对齐的内存中。对应指令：movaps m128, xmm

最后在SSE版本中再次使用循环展开:

//SSE+四路展开
    float sumfloat_sse_4loop(const float *pbuf,unsigned int cntbuf)
    {
    	float s=0;
    	unsigned int nBlockWidth=4*4;
    	unsigned int cntBlock=cntbuf/nBlockWidth;
    	unsigned int cntRem=cntbuf%nBlockWidth;
    	__m128 fSum0=_mm_setzero_ps();//求和变量，初值清零
    	__m128 fSum1=_mm_setzero_ps();
    	__m128 fSum2=_mm_setzero_ps();
    	__m128 fSum3=_mm_setzero_ps();
    	__m128 fLoad0,fLoad1,fLoad2,fLoad3;
    	const float *p=pbuf;
    	for (unsigned int i=0;i<cntBlock;i++)
    	{
    		fLoad0=_mm_load_ps(p);//加载
    		fLoad1=_mm_load_ps(p+4);
    		fLoad2=_mm_load_ps(p+8);
    		fLoad3=_mm_load_ps(p+12);
    		fSum0=_mm_add_ps(fSum0,fLoad0);//求和
    		fSum1=_mm_add_ps(fSum1,fLoad1);
    		fSum2=_mm_add_ps(fSum2,fLoad2);
    		fSum3=_mm_add_ps(fSum3,fLoad3);
    		p+=nBlockWidth;
    	}
    	fSum0=_mm_add_ps(fSum0,fSum1);
    	fSum2=_mm_add_ps(fSum2,fSum3);
    	fSum0=_mm_add_ps(fSum0,fSum2);
    	const float*q=(const float*)&fSum0;
    	s=q[0]+q[1]+q[2]+q[3];				//合并
    	for (unsigned int i=0;i<cntRem;i++)//处理尾部剩余数据
    	{
    		s+=p[i];
    	}
    	return s;
    }

**完整程序：**

Timing.h

#include <windows.h>
    static _LARGE_INTEGER time_start, time_over;
    static double dqFreq;
    static inline void startTiming()
    {
    	_LARGE_INTEGER f;
    	QueryPerformanceFrequency(&f);
    	dqFreq=(double)f.QuadPart;
    
    	QueryPerformanceCounter(&time_start);
    }
    static inline double stopTiming()
    {
    	QueryPerformanceCounter(&time_over);
    	return ((double)(time_over.QuadPart-time_start.QuadPart)/dqFreq*1000);
    }

#include <stdio.h>
    #include <intrin.h>
    #include <stdlib.h>
    #include <time.h>
    #include "Timing.h"
    #define BUFSIZE	4096 // = 32KB{L1 Cache} / (2 * sizeof(float))
    __declspec(align(16))float buf[BUFSIZE];//内存对齐
    typedef float (*TESTPROC)(const float* pbuf, unsigned int cntbuf);//函数指针（用于测试时统一表示）
    void RunTest(const char *szname,TESTPROC proc);
    float sumfloat_base(const float *pbuf,unsigned int cntbuf);
    float sumfloat_base_4loop(const float*pbuf,unsigned int cntbuf);
    float sumfloat_sse(const float *pbuf ,unsigned int cntbuf);
    float sumfloat_sse_4loop(const float *pbuf,unsigned int cntbuf);
    
    
    int main()
    {
    	srand( (unsigned)time( NULL ) );
    	for (int i = 0; i < BUFSIZE; i++) 
    		buf[i] = (float)(rand() & 0x3f);// 使用&0x3f是为了让求和后的数值不会超过float类型的有效位数，便于观察结果是否正确.
    	RunTest("sumfloat_base",sumfloat_base);
    	RunTest("sumfloat_base_4loop",sumfloat_base_4loop);
    	RunTest("sumfloat_sse",sumfloat_sse);
    	RunTest("sumfloat_sse_4loop",sumfloat_sse_4loop);
    	return 0;
    }
    
    
    //测试函数
    void RunTest(const char *szname,TESTPROC proc)
    {
    	unsigned int testloop=4000;//循环次数
    
    	volatile float result;//volatile类型放止编译器优化使得循环内部不执行！
    	double mpsgood=0;
    	double mps;
    	for (int k=0;k<=3;k++)//循环多次，选取最好情况
    	{
    		startTiming();
    		for(unsigned int i=0;i<testloop;i++)
    		{
    			result=proc(buf,BUFSIZE);
    		}
    		double interval=stopTiming();
    		mps=testloop*BUFSIZE*4*1000.0/(interval*1024*1024);//单位MB/s
    		if (mpsgood<mps)mpsgood=mps;
    	}
    	
    	printf("%s:\t%f,\t%.0lfMB/s\n",szname,result,mpsgood);//测速单位MB/s
    }
    //普通版
    float sumfloat_base(const float *pbuf,unsigned int cntbuf)
    {
    	float s=0;
    	for (unsigned int i=0;i<cntbuf;i++)
    	{
    		s+=pbuf[i];
    	}
    	return s;
    }
    
    //四路展开
    float sumfloat_base_4loop(const float*pbuf,unsigned int cntbuf)
    {
    	float s=0;
    	float fSum0=0,fSum1=0,fSum2=0,fSum3=0;
    	unsigned int i=0;
    	const float *p=pbuf;
    	for (;i<=cntbuf-4;i+=4)//cntbuf-4不会每次都计算，编译器实现会计算好循环次数！
    	{
    		fSum0+=p[i];
    		fSum1+=p[i+1];
    		fSum2+=p[i+2];
    		fSum3+=p[i+3];
    	}
    	s=fSum0+fSum1+fSum2+fSum3;
    	for (;i<cntbuf;i++)
    	{
    		s+=p[i];
    	}
    	return s;
    }
    //SSE版 
    float sumfloat_sse(const float *pbuf ,unsigned int cntbuf)
    {
    	float s=0;
    	int nBlockWidth=4;//SSE一次处理4个float
    	int cntBlock=cntbuf/nBlockWidth;
    	int cntRem=cntbuf%nBlockWidth;
    	__m128 fSum=_mm_setzero_ps();//求和变量，初值清零
    	__m128 fLoad;
    	const float*p=pbuf;
    	for (unsigned int i=0;i<cntBlock;i++)
    	{
    		fLoad=_mm_load_ps(p);//加载
    		fSum=_mm_add_ps(fSum,fLoad);//求和
    		p+=nBlockWidth;
    	}
    	const float *q=(const float*)&fSum;
    	s=q[0]+q[1]+q[2]+q[3];		//合并
    	for (int i=0;i<cntRem;i++)//处理尾部剩余数据
    	{
    		s+=p[i];
    	}
    	return s;
    }
    
    //SSE+四路展开
    float sumfloat_sse_4loop(const float *pbuf,unsigned int cntbuf)
    {
    	float s=0;
    	unsigned int nBlockWidth=4*4;
    	unsigned int cntBlock=cntbuf/nBlockWidth;
    	unsigned int cntRem=cntbuf%nBlockWidth;
    	__m128 fSum0=_mm_setzero_ps();//求和变量，初值清零
    	__m128 fSum1=_mm_setzero_ps();
    	__m128 fSum2=_mm_setzero_ps();
    	__m128 fSum3=_mm_setzero_ps();
    	__m128 fLoad0,fLoad1,fLoad2,fLoad3;
    	const float *p=pbuf;
    	for (unsigned int i=0;i<cntBlock;i++)
    	{
    		fLoad0=_mm_load_ps(p);//加载
    		fLoad1=_mm_load_ps(p+4);
    		fLoad2=_mm_load_ps(p+8);
    		fLoad3=_mm_load_ps(p+12);
    		fSum0=_mm_add_ps(fSum0,fLoad0);//求和
    		fSum1=_mm_add_ps(fSum1,fLoad1);
    		fSum2=_mm_add_ps(fSum2,fLoad2);
    		fSum3=_mm_add_ps(fSum3,fLoad3);
    		p+=nBlockWidth;
    	}
    	fSum0=_mm_add_ps(fSum0,fSum1);
    	fSum2=_mm_add_ps(fSum2,fSum3);
    	fSum0=_mm_add_ps(fSum0,fSum2);
    	const float*q=(const float*)&fSum0;
    	s=q[0]+q[1]+q[2]+q[3];				//合并
    	for (unsigned int i=0;i<cntRem;i++)//处理尾部剩余数据
    	{
    		s+=p[i];
    	}
    	return s;
    }

**测试结果：**

![20161006144411100][]

\[1\][http://www.cnblogs.com/zyl910/archive/2012/10/22/simdsumfloat.html\#undefined][1]

[Intrinsics]: http://blog.csdn.net/lee353086/article/details/38437469
[Intel Intrinsics Guide]: https://software.intel.com/sites/landingpage/IntrinsicsGuide/
[1]: http://www.cnblogs.com/zyl910/archive/2012/10/22/simdsumfloat.html#undefined
[20161006144411100]: /images/20220716/e7b0ef17c85748448c9968ccddc6d301.png