[C] 跨平台使用Intrinsic函数范例2——使用SSE2、AVX指令集 处理 双精度浮点数组求和
作者：zyl910。
本文面对对SSE等SIMD指令集有一定基础的读者，以双精度浮点数组求和为例演示了如何跨平台使用SSE2、AVX指令集。支持vc、gcc编译器，在Windows、Linux、Mac这三大平台上成功运行。
一、关键讲解
前文（）演示了如何使用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
其次，还需要调整一下地址计算。对于C语言来说，将float指针改为double指针就能解决大多数地址计算问题了。然后再修正一下块宽计算、改写一下合并时的代码，便大功告成了。 例如sumfloat_sse与sumdouble_sse——
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; * p = pbuf; * q; // 将SSE变量上的多个数值合并时所用指针.(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];(i=0; i<cntRem; ++i){s += p[i];}return s;}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; * p = pbuf; * q; // 将SSE变量上的多个数值合并时所用指针.(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];(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 #include MACTOSTR(x) #x#define MACROVALUESTR(x) MACTOSTR(x)#if defined(__ICL) // Intel C++# if defined(__VERSION__)# define COMPILER_NAME __VERSION__# elif defined(__INTEL_COMPILER_BUILD_DATE)# define COMPILER_NAME MACROVALUESTR(__INTEL_COMPILER_BUILD_DATE) # else# define COMPILER_NAME # endif defined(_MSC_VER) // Microsoft VC++# if defined(_MSC_FULL_VER)# define COMPILER_NAME MACROVALUESTR(_MSC_FULL_VER) # elif defined(_MSC_VER)# define COMPILER_NAME MACROVALUESTR(_MSC_VER) # else# define COMPILER_NAME # endif defined(__GNUC__) // GCC# if defined(__CYGWIN__)# define COMPILER_NAME __VERSION__# elif defined(__MINGW32__)# define COMPILER_NAME __VERSION__# else# define COMPILER_NAME __VERSION__# endif # define COMPILER_NAME 双精度浮点数组求和_基本版.//// result: 返回数组求和结果.// pbuf: 数组的首地址.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_SSE2sumdouble_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; * p = pbuf; * q; // 将SSE变量上的多个数值合并时所用指针.(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];(i=0; i<cntRem; ++i){s += p[i];}return s;}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;* q; // 将SSE变量上的多个数值合并时所用指针.(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];(i=0; i<cntRem; ++i){s += p[i];}return s;}#ifdef INTRIN_AVXsumdouble_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; * p = pbuf; * q; // 将AVX变量上的多个数值合并时所用指针.(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];(i=0; i<cntRem; ++i){s += p[i];}return s;}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;* q; // 将AVX变量上的多个数值合并时所用指针.(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];(i=0; i<cntRem; ++i){s += p[i];}return s;}#ifndef ATTR_ALIGN# if defined(__GNUC__) // GCC# define ATTR_ALIGN(n) __attribute__((aligned(n)))# # define ATTR_ALIGN(n) __declspec(align(n))# endifBUFSIZE 2048 // = 32KB{L1 Cache} / (2 * sizeof(double))ATTR_ALIGN(32) double buf[BUFSIZE];// 测试时的函数类型typedef double (*TESTPROC)(const double* pbuf, size_t cntbuf);runTest(const char* szname, TESTPROC proc){clock_t TIMEOUT = CLOCKS_PER_SEC/ i,j,k;clock_t tm0, dt; mps; mps_good = n=(i=1; i<=3; ++i) // 多次测试. {tm0 = clock();// maink=0;do{for(j=1; j<=testloop; ++j) // 重复运算几次延长时间，避免计时开销带来的影响.{n = proc(buf, BUFSIZE); // 避免内循环被编译优化消掉.}++k;dt = clock() – tm0;}while(dt<TIMEOUT);// showmps = ((mps_good<mps) mps_good=mps; // 选取最佳值. }printf(, szname, mps_good, n);}int main(int argc, char* argv[]){char szBuf[64];int i;printf(, INTRIN_WORDSIZE);printf(, COMPILER_NAME);cpu_getbrand(szBuf);printf(, szBuf);printf();// init buf srand( (unsigned)time( NULL ) );for (i = 0; i < BUFSIZE; i++) buf[i] = (double)(rand() & 0x7fff); // 使用&0x7fff是为了让求和后的数值在一定范围内，便于观察结果是否正确.// testrunTest(, sumdouble_base); // 双精度浮点数组求和_基本版.#ifdef INTRIN_SSE2if (simd_sse_level(NULL) >= SIMD_SSE_2){runTest(, sumdouble_sse); // 双精度浮点数组求和_SSE版.runTest(, sumdouble_sse_4loop); // 双精度浮点数组求和_SSE四路循环展开版. }#ifdef INTRIN_AVXif (simd_avx_level(NULL) >= SIMD_AVX_1){runTest(, sumdouble_avx); // 双精度浮点数组求和_SSE版.runTest(, sumdouble_avx_4loop); // 双精度浮点数组求和_SSE四路循环展开版. };}
2.2 makefile
全部代码——
好想从现在开始抱着你，紧紧地抱着你，一直走到上帝面前。
