[C] 跨平台使用Intrinsic函数范例2——使用SSE2、AVX指令集 处理 双精度浮点数组求和

作者:zyl910

  本文面对对SSE等SIMD指令集有一定基础的读者,以双精度浮点数组求和为例演示了如何跨平台使用SSE2、AVX指令集。支持vc、gcc编译器,在Windows、Linux、Mac这三大平台上成功运行。


一、关键讲解

  前文(http://www.cnblogs.com/zyl910/archive/2012/10/22/simdsumfloat.html)演示了如何使用SSE、AVX指令集 处理 单精度浮点数组求和。现在对其进行改造,使用SSE2、AVX指令集 处理 双精度浮点数组求和。因程序基本上差不多,文本就不详细讲解了,只说关键变化。

1.1 指令集简介

  先来看看支持双精度浮点数的SIMD的指令集——
SSE指令集只支持单精度浮点运算,直到SSE2指令集才支持双精度浮点数运算。SSE2定义了128位紧缩双精度浮点类型,对应Intrinsic中的__m128d类型,它能一次能处理2个双精度浮点数。
AVX指令集从一开始就支持单精度与双精度浮点运算(但整数运算要等AVX2)。AVX定义了256位紧缩双精度浮点类型,对应Intrinsic中的__m256d类型,它能一次能处理4个双精度浮点数。


1.2 改造为 SSE2、AVX的双精度浮点代码

  在使用Intrinsic函数时,将 SSE、AVX的单精度浮点代码 改造为 SSE2、AVX的双精度浮点代码是很方便的。对比前文与本文的数组求和代码,变更的地方有——

float double 备注
指令 Intrinsic Asm 指令 Intrinsic Asm
SSE __m128 XMMWORD SSE2 __m128d XMMWORD 类型
_mm_setzero_ps XORPS _mm_setzero_pd XORPD 赋0
_mm_load_ps MOVAPS _mm_load_pd MOVAPD 加载
_mm_add_ps ADDPS _mm_add_pd ADDPD 加法
AVX __m256 YMMWORD AVX __m256d YMMWORD 类型
_mm256_setzero_ps VXORPS _mm256_setzero_pd VXORPD 赋0
_mm256_load_ps VMOVAPS _mm256_load_pd VMOVAPD 加载
_mm256_add_ps VADDPS _mm256_add_pd VADDPD 加法

  其次,还需要调整一下地址计算。对于C语言来说,将float指针改为double指针就能解决大多数地址计算问题了。然后再修正一下块宽计算、改写一下合并时的代码,便大功告成了。
  例如sumfloat_sse与sumdouble_sse——

// 单精度浮点数组求和_SSE版.
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;
}

// 双精度浮点数组求和_SSE版.
double sumdouble_sse(const double* pbuf, size_t cntbuf)
{
	double s = 0;	// 求和变量.
	size_t i;
	size_t nBlockWidth = 2;	// 块宽. SSE寄存器能一次处理2个double.
	size_t cntBlock = cntbuf / nBlockWidth;	// 块数.
	size_t cntRem = cntbuf % nBlockWidth;	// 剩余数量.
	__m128d xfdSum = _mm_setzero_pd();	// 求和变量。[SSE2] XORPD. 赋初值0.
	__m128d xfdLoad;	// 加载.
	const double* p = pbuf;	// SSE批量处理时所用的指针.
	const double* q;	// 将SSE变量上的多个数值合并时所用指针.

	// SSE批量处理.
	for(i=0; i<cntBlock; ++i)
	{
		xfdLoad = _mm_load_pd(p);	// [SSE2] MOVAPD. 加载.
		xfdSum = _mm_add_pd(xfdSum, xfdLoad);	// [SSE2] ADDPD. 双精浮点紧缩加法.
		p += nBlockWidth;
	}
	// 合并.
	q = (const double*)&xfdSum;
	s = q[0] + q[1];

	// 处理剩下的.
	for(i=0; i<cntRem; ++i)
	{
		s += p[i];
	}

	return s;
}



 

1.3 环境检查

  最后,别忘了检查环境——
INTRIN_SSE2、INTRIN_AVX 宏是 zintrin.h 提供的,可用来在编译时检测编译器是否支持SSE2、AVX指令集。
simd_sse_level、simd_avx_level函数是 ccpuid.h 提供的,可用来在运行时检测当前系统环境是否支持SSE2、AVX指令集。


二、全部代码

2.1 simdsumdouble.c

  全部代码——

#define __STDC_LIMIT_MACROS	1	// C99整数范围常量. [纯C程序可以不用, 而C++程序必须定义该宏.]

#include <stdlib.h>
#include <stdio.h>
#include <time.h>

#include "zintrin.h"
#include "ccpuid.h"


// Compiler name
#define MACTOSTR(x)	#x
#define MACROVALUESTR(x)	MACTOSTR(x)
#if defined(__ICL)	// Intel C++
#  if defined(__VERSION__)
#    define COMPILER_NAME	"Intel C++ " __VERSION__
#  elif defined(__INTEL_COMPILER_BUILD_DATE)
#    define COMPILER_NAME	"Intel C++ (" MACROVALUESTR(__INTEL_COMPILER_BUILD_DATE) ")"
#  else
#    define COMPILER_NAME	"Intel C++"
#  endif	// #  if defined(__VERSION__)
#elif defined(_MSC_VER)	// Microsoft VC++
#  if defined(_MSC_FULL_VER)
#    define COMPILER_NAME	"Microsoft VC++ (" MACROVALUESTR(_MSC_FULL_VER) ")"
#  elif defined(_MSC_VER)
#    define COMPILER_NAME	"Microsoft VC++ (" MACROVALUESTR(_MSC_VER) ")"
#  else
#    define COMPILER_NAME	"Microsoft VC++"
#  endif	// #  if defined(_MSC_FULL_VER)
#elif defined(__GNUC__)	// GCC
#  if defined(__CYGWIN__)
#    define COMPILER_NAME	"GCC(Cygmin) " __VERSION__
#  elif defined(__MINGW32__)
#    define COMPILER_NAME	"GCC(MinGW) " __VERSION__
#  else
#    define COMPILER_NAME	"GCC " __VERSION__
#  endif	// #  if defined(_MSC_FULL_VER)
#else
#  define COMPILER_NAME	"Unknown Compiler"
#endif	// #if defined(__ICL)	// Intel C++


//////////////////////////////////////////////////
// sumdouble: 双精度浮点数组求和的函数
//////////////////////////////////////////////////

// 双精度浮点数组求和_基本版.
//
// result: 返回数组求和结果.
// pbuf: 数组的首地址.
// cntbuf: 数组长度.
double sumdouble_base(const double* pbuf, size_t cntbuf)
{
	double s = 0;	// 求和变量.
	size_t i;
	for(i=0; i<cntbuf; ++i)
	{
		s += pbuf[i];
	}
	return s;
}

#ifdef INTRIN_SSE2
// 双精度浮点数组求和_SSE版.
double sumdouble_sse(const double* pbuf, size_t cntbuf)
{
	double s = 0;	// 求和变量.
	size_t i;
	size_t nBlockWidth = 2;	// 块宽. SSE寄存器能一次处理2个double.
	size_t cntBlock = cntbuf / nBlockWidth;	// 块数.
	size_t cntRem = cntbuf % nBlockWidth;	// 剩余数量.
	__m128d xfdSum = _mm_setzero_pd();	// 求和变量。[SSE2] XORPD. 赋初值0.
	__m128d xfdLoad;	// 加载.
	const double* p = pbuf;	// SSE批量处理时所用的指针.
	const double* q;	// 将SSE变量上的多个数值合并时所用指针.

	// SSE批量处理.
	for(i=0; i<cntBlock; ++i)
	{
		xfdLoad = _mm_load_pd(p);	// [SSE2] MOVAPD. 加载.
		xfdSum = _mm_add_pd(xfdSum, xfdLoad);	// [SSE2] ADDPD. 双精浮点紧缩加法.
		p += nBlockWidth;
	}
	// 合并.
	q = (const double*)&xfdSum;
	s = q[0] + q[1];

	// 处理剩下的.
	for(i=0; i<cntRem; ++i)
	{
		s += p[i];
	}

	return s;
}

// 双精度浮点数组求和_SSE四路循环展开版.
double sumdouble_sse_4loop(const double* pbuf, size_t cntbuf)
{
	double s = 0;	// 返回值.
	size_t i;
	size_t nBlockWidth = 2*4;	// 块宽. SSE寄存器能一次处理2个double,然后循环展开4次.
	size_t cntBlock = cntbuf / nBlockWidth;	// 块数.
	size_t cntRem = cntbuf % nBlockWidth;	// 剩余数量.
	__m128d xfdSum = _mm_setzero_pd();	// 求和变量。[SSE2] XORPD. 赋初值0.
	__m128d xfdSum1 = _mm_setzero_pd();
	__m128d xfdSum2 = _mm_setzero_pd();
	__m128d xfdSum3 = _mm_setzero_pd();
	__m128d xfdLoad;	// 加载.
	__m128d xfdLoad1;
	__m128d xfdLoad2;
	__m128d xfdLoad3;
	const double* p = pbuf;	// SSE批量处理时所用的指针.
	const double* q;	// 将SSE变量上的多个数值合并时所用指针.

	// SSE批量处理.
	for(i=0; i<cntBlock; ++i)
	{
		xfdLoad = _mm_load_pd(p);	// [SSE2] MOVAPD. 加载.
		xfdLoad1 = _mm_load_pd(p+2);
		xfdLoad2 = _mm_load_pd(p+4);
		xfdLoad3 = _mm_load_pd(p+6);
		xfdSum = _mm_add_pd(xfdSum, xfdLoad);	// [SSE2] ADDPD. 双精浮点紧缩加法.
		xfdSum1 = _mm_add_pd(xfdSum1, xfdLoad1);
		xfdSum2 = _mm_add_pd(xfdSum2, xfdLoad2);
		xfdSum3 = _mm_add_pd(xfdSum3, xfdLoad3);
		p += nBlockWidth;
	}
	// 合并.
	xfdSum = _mm_add_pd(xfdSum, xfdSum1);	// 两两合并(0~1).
	xfdSum2 = _mm_add_pd(xfdSum2, xfdSum3);	// 两两合并(2~3).
	xfdSum = _mm_add_pd(xfdSum, xfdSum2);	// 两两合并(0~3).
	q = (const double*)&xfdSum;
	s = q[0] + q[1];

	// 处理剩下的.
	for(i=0; i<cntRem; ++i)
	{
		s += p[i];
	}

	return s;
}
#endif	// #ifdef INTRIN_SSE2


#ifdef INTRIN_AVX
// 双精度浮点数组求和_AVX版.
double sumdouble_avx(const double* pbuf, size_t cntbuf)
{
	double s = 0;	// 求和变量.
	size_t i;
	size_t nBlockWidth = 4;	// 块宽. AVX寄存器能一次处理4个double.
	size_t cntBlock = cntbuf / nBlockWidth;	// 块数.
	size_t cntRem = cntbuf % nBlockWidth;	// 剩余数量.
	__m256d yfdSum = _mm256_setzero_pd();	// 求和变量。[AVX] VXORPD. 赋初值0.
	__m256d yfdLoad;	// 加载.
	const double* p = pbuf;	// AVX批量处理时所用的指针.
	const double* q;	// 将AVX变量上的多个数值合并时所用指针.

	// AVX批量处理.
	for(i=0; i<cntBlock; ++i)
	{
		yfdLoad = _mm256_load_pd(p);	// [AVX] VMOVAPD. 加载.
		yfdSum = _mm256_add_pd(yfdSum, yfdLoad);	// [AVX] VADDPD. 双精浮点紧缩加法.
		p += nBlockWidth;
	}
	// 合并.
	q = (const double*)&yfdSum;
	s = q[0] + q[1] + q[2] + q[3];

	// 处理剩下的.
	for(i=0; i<cntRem; ++i)
	{
		s += p[i];
	}

	return s;
}

// 双精度浮点数组求和_AVX四路循环展开版.
double sumdouble_avx_4loop(const double* pbuf, size_t cntbuf)
{
	double s = 0;	// 求和变量.
	size_t i;
	size_t nBlockWidth = 4*4;	// 块宽. AVX寄存器能一次处理8个double,然后循环展开4次.
	size_t cntBlock = cntbuf / nBlockWidth;	// 块数.
	size_t cntRem = cntbuf % nBlockWidth;	// 剩余数量.
	__m256d yfdSum = _mm256_setzero_pd();	// 求和变量。[AVX] VXORPD. 赋初值0.
	__m256d yfdSum1 = _mm256_setzero_pd();
	__m256d yfdSum2 = _mm256_setzero_pd();
	__m256d yfdSum3 = _mm256_setzero_pd();
	__m256d yfdLoad;	// 加载.
	__m256d yfdLoad1;
	__m256d yfdLoad2;
	__m256d yfdLoad3;
	const double* p = pbuf;	// AVX批量处理时所用的指针.
	const double* q;	// 将AVX变量上的多个数值合并时所用指针.

	// AVX批量处理.
	for(i=0; i<cntBlock; ++i)
	{
		yfdLoad = _mm256_load_pd(p);	// [AVX] VMOVAPD. 加载.
		yfdLoad1 = _mm256_load_pd(p+4);
		yfdLoad2 = _mm256_load_pd(p+8);
		yfdLoad3 = _mm256_load_pd(p+12);
		yfdSum = _mm256_add_pd(yfdSum, yfdLoad);	// [AVX] VADDPD. 双精浮点紧缩加法.
		yfdSum1 = _mm256_add_pd(yfdSum1, yfdLoad1);
		yfdSum2 = _mm256_add_pd(yfdSum2, yfdLoad2);
		yfdSum3 = _mm256_add_pd(yfdSum3, yfdLoad3);
		p += nBlockWidth;
	}
	// 合并.
	yfdSum = _mm256_add_pd(yfdSum, yfdSum1);	// 两两合并(0~1).
	yfdSum2 = _mm256_add_pd(yfdSum2, yfdSum3);	// 两两合并(2~3).
	yfdSum = _mm256_add_pd(yfdSum, yfdSum2);	// 两两合并(0~3).
	q = (const double*)&yfdSum;
	s = q[0] + q[1] + q[2] + q[3];

	// 处理剩下的.
	for(i=0; i<cntRem; ++i)
	{
		s += p[i];
	}

	return s;
}

#endif	// #ifdef INTRIN_AVX



//////////////////////////////////////////////////
// main
//////////////////////////////////////////////////

// 变量对齐.
#ifndef ATTR_ALIGN
#  if defined(__GNUC__)	// GCC
#    define ATTR_ALIGN(n)	__attribute__((aligned(n)))
#  else	// 否则使用VC格式.
#    define ATTR_ALIGN(n)	__declspec(align(n))
#  endif
#endif	// #ifndef ATTR_ALIGN


#define BUFSIZE	2048	// = 32KB{L1 Cache} / (2 * sizeof(double))
ATTR_ALIGN(32) double buf[BUFSIZE];

// 测试时的函数类型
typedef double (*TESTPROC)(const double* pbuf, size_t cntbuf);

// 进行测试
void runTest(const char* szname, TESTPROC proc)
{
	const int testloop = 4000;	// 重复运算几次延长时间,避免计时精度问题.
	const clock_t TIMEOUT = CLOCKS_PER_SEC/2;	// 最短测试时间.
	int i,j,k;
	clock_t	tm0, dt;	// 存储时间.
	double mps;	// M/s.
	double mps_good = 0;	// 最佳M/s. 因线程切换会导致的数值波动, 于是选取最佳值.
	volatile double n=0;	// 避免内循环被优化.
	for(i=1; i<=3; ++i)	// 多次测试.
	{
		tm0 = clock();
		// main
		k=0;
		do
		{
			for(j=1; j<=testloop; ++j)	// 重复运算几次延长时间,避免计时开销带来的影响.
			{
				n = proc(buf, BUFSIZE);	// 避免内循环被编译优化消掉.
			}
			++k;
			dt = clock() - tm0;
		}while(dt<TIMEOUT);
		// show
		mps = (double)k*testloop*BUFSIZE*CLOCKS_PER_SEC/(1024.0*1024.0*dt);	// k*testloop*BUFSIZE/(1024.0*1024.0) 将数据规模换算为M,然后再乘以 CLOCKS_PER_SEC/dt 换算为M/s .
		if (mps_good<mps)	mps_good=mps;	// 选取最佳值.
		//printf("%s:\t%.0f M/s\t//%f\n", szname, mps, n);
	}
	printf("%s:\t%.0f M/s\t//%f\n", szname, mps_good, n);
}

int main(int argc, char* argv[])
{
	char szBuf[64];
	int i;

	printf("simdsumdouble v1.00 (%dbit)\n", INTRIN_WORDSIZE);
	printf("Compiler: %s\n", COMPILER_NAME);
	cpu_getbrand(szBuf);
	printf("CPU:\t%s\n", szBuf);
	printf("\n");

	// init buf
	srand( (unsigned)time( NULL ) );
	for (i = 0; i < BUFSIZE; i++) buf[i] = (double)(rand() & 0x7fff);	// 使用&0x7fff是为了让求和后的数值在一定范围内,便于观察结果是否正确.

	// test
	runTest("sumdouble_base", sumdouble_base);	// 双精度浮点数组求和_基本版.
#ifdef INTRIN_SSE2
	if (simd_sse_level(NULL) >= SIMD_SSE_2)
	{
		runTest("sumdouble_sse", sumdouble_sse);	// 双精度浮点数组求和_SSE版.
		runTest("sumdouble_sse_4loop", sumdouble_sse_4loop);	// 双精度浮点数组求和_SSE四路循环展开版.
	}
#endif	// #ifdef INTRIN_SSE2
#ifdef INTRIN_AVX
	if (simd_avx_level(NULL) >= SIMD_AVX_1)
	{
		runTest("sumdouble_avx", sumdouble_avx);	// 双精度浮点数组求和_SSE版.
		runTest("sumdouble_avx_4loop", sumdouble_avx_4loop);	// 双精度浮点数组求和_SSE四路循环展开版.
	}
#endif	// #ifdef INTRIN_AVX

	return 0;
}



2.2 makefile

  全部代码——

# flags
CC = g++
CFS = -Wall -msse2

# args
RELEASE =0
BITS =
CFLAGS =

# [args] 生成模式. 0代表debug模式, 1代表release模式. make RELEASE=1.
ifeq ($(RELEASE),0)
	# debug
	CFS += -g
else
	# release
	CFS += -O3 -DNDEBUG
	//CFS += -O3 -g -DNDEBUG
endif

# [args] 程序位数. 32代表32位程序, 64代表64位程序, 其他默认. make BITS=32.
ifeq ($(BITS),32)
	CFS += -m32
else
	ifeq ($(BITS),64)
		CFS += -m64
	else
	endif
endif

# [args] 使用 CFLAGS 添加新的参数. make CFLAGS="-mavx".
CFS += $(CFLAGS)


.PHONY : all clean

# files
TARGETS = simdsumdouble
OBJS = simdsumdouble.o

all : $(TARGETS)

simdsumdouble : $(OBJS)
	$(CC) $(CFS) -o $@ $^


simdsumdouble.o : simdsumdouble.c zintrin.h ccpuid.h
	$(CC) $(CFS) -c $<


clean :
	rm -f $(OBJS) $(TARGETS) $(addsuffix .exe,$(TARGETS))



三、编译测试

3.1 编译

  在以下编译器中成功编译——
VC6:x86版。
VC2003:x86版。
VC2005:x86版。
VC2010:x86版、x64版。
GCC 4.7.0(Fedora 17 x64):x86版、x64版。
GCC 4.6.2(MinGW(20120426)):x86版。
GCC 4.7.1(TDM-GCC(MinGW-w64)):x86版、x64版。
llvm-gcc-4.2(Mac OS X Lion 10.7.4, Xcode 4.4.1):x86版、x64版。



3.2 测试

  因虚拟机上的有效率损失,于是仅在真实系统上进行测试。

  系统环境——
CPU:Intel(R) Core(TM) i3-2310M CPU @ 2.10GHz
操作系统:Windows 7 SP1 x64版

  然后分别运行VC与GCC编译的Release版可执行文件,即以下4个程序——
exe\simdsumdouble_vc32.exe:VC2010 SP1 编译的32位程序,/O2 /arch:SSE2。
exe\simdsumdouble_vc64.exe:VC2010 SP1 编译的64位程序,/O2 /arch:AVX。
exe\simdsumdouble_gcc32.exe:GCC 4.7.1(TDM-GCC(MinGW-w64)) 编译的32位程序,-O3 -mavx。
exe\simdsumdouble_gcc64.exe:GCC 4.7.1(TDM-GCC(MinGW-w64)) 编译的64位程序,-O3 -mavx。

  测试结果(使用cmdarg_ui)——

 

参考文献——
《Intel® 64 and IA-32 Architectures Software Developer’s Manual Combined Volumes:1, 2A, 2B, 2C, 3A, 3B, and 3C》044US. August 2012.http://www.intel.com/content/www/us/en/processors/architectures-software-developer-manuals.html
《Intel® Architecture Instruction Set Extensions Programming Reference》014. AUGUST 2012.http://software.intel.com/en-us/avx/
《AMD64 Architecture Programmer’s Manual Volume 4: 128-Bit and 256-Bit Media Instructions》. December 2011.http://developer.amd.com/documentation/guides/Pages/default.aspx#manuals
《[C] 让VC、BCB支持C99的整数类型(stdint.h、inttypes.h)(兼容GCC)》. http://www.cnblogs.com/zyl910/archive/2012/08/08/c99int.html
《[C] zintrin.h: 智能引入intrinsic函数 V1.01版。改进对Mac OS X的支持,增加INTRIN_WORDSIZE宏》. http://www.cnblogs.com/zyl910/archive/2012/10/01/zintrin_v101.html
《[C/C++] ccpuid:CPUID信息模块 V1.03版,改进mmx/sse指令可用性检查(使用signal、setjmp,支持纯C)、修正AVX检查Bug》.http://www.cnblogs.com/zyl910/archive/2012/10/13/ccpuid_v103.html
《[x86]SIMD指令集发展历程表(MMX、SSE、AVX等)》. http://www.cnblogs.com/zyl910/archive/2012/02/26/x86_simd_table.html
《SIMD(MMX/SSE/AVX)变量命名规范心得》. http://www.cnblogs.com/zyl910/archive/2012/04/23/simd_var_name.html
《GCC 64位程序的makefile条件编译心得——32位版与64位版、debug版与release版(兼容MinGW、TDM-GCC)》. http://www.cnblogs.com/zyl910/archive/2012/08/14/gcc64_make.html
《[C#] cmdarg_ui:“简单参数命令行程序”的通用图形界面》.  http://www.cnblogs.com/zyl910/archive/2012/06/19/cmdarg_ui.html
《[C] 跨平台使用Intrinsic函数范例1——使用SSE、AVX指令集 处理 单精度浮点数组求和(支持vc、gcc,兼容Windows、Linux、Mac)》. http://www.cnblogs.com/zyl910/archive/2012/10/22/simdsumfloat.html


源码下载——
http://files.cnblogs.com/zyl910/simdsumdouble.rar

阅读更多
个人分类: #All - 所有 SIMD
想对作者说点什么? 我来说一句

没有更多推荐了,返回首页

关闭
关闭
关闭