• 并行编程之跨平台使用SSE、AVX指令集心得——以单精度浮点数组求和为例(支持vc、gcc,兼容Windows、Linux、Mac)


    作者:zyl910

      本文面对对SSE等SIMD指令集有一定基础的读者,以单精度浮点数组求和为例演示了如何跨平台使用SSE、AVX指令集。因使用了stdint、zintrin、ccpuid这三个模块,可以完全避免手工编写汇编代码,具有很高可移植性。支持vc、gcc编译器,在Windows、Linux、Mac这三大平台上成功运行。


    一、问题背景

      最初,我们只能使用汇编语言来编写SIMD代码。不仅写起来很麻烦,而且易读性、可维护性、移植性都较差。
      不久,VC、GCC等编译器相继支持了Intrinsic函数,使我们可以摆脱汇编,利用C语言来调用SIMD指令集,大大提高了易读性和可维护。而且移植性也有提高,能在同一编译器上实现32位与64位的平滑过渡。
      但当代码在另一种编译器编译时,会遇到一些问题而无法编译。甚至在使用同一种编译器的不同版本时,也会遇到无法编译问题。

      首先是整数类型问题——
      传统C语言的short、int、long等整数类型是与平台相关的,不同平台上的位长是不同的(例如Windows是LLP64模型,Linux、Mac等Unix系统多采用LP64模型)。而使用SSE等SIMD指令集时需要精确计算数据的位数,不同位长的数据必须使用不同的指令来处理。
      有一个解决办法,就是使用C99标准中stdint.h所提供的指定位长的整数类型。GCC对C99标准支持性较好,而VC的步骤很慢,貌似直到VC2010才支持stdint.h。而很多时候我们为了兼容旧代码,不得不使用VC等老版本的VC编译器。

      其次是Intrinsic函数的头文件问题,不同编译器所使用的头文件不同——
      对于早期版本VC,需要根据具体的指令集需求,手动引入mmintrin.h、xmmintrin.h等头文件。对于VC2005或更高版本,引入intrin.h就行了,它会自动引入当前编译器所支持的所有Intrinsic头文件。
      对于早期版本GCC,也是手动引入mmintrin.h、xmmintrin.h等头文件。而对于高版本的GCC,引入x86intrin.h就行了,它会自动引入当前编译环境所允许的Intrinsic头文件。

      再次是当前编译环境下的Intrinsic函数集支持性问题——
      对于VC来说,VC6支持MMX、3DNow!、SSE、SSE2,然后更高版本的VC支持更多的指令集。但是,VC没有提供检测Intrinsic函数集支持性的办法。例如你在VC2010上编写了一段使用了AVX Intrinsic函数的代码,但拿到VC2005上就不能通过编译了。其次,VC不支持64位下的MMX,这让一些老程序迁徙到64位版时遭来了一些麻烦。
      而对于GCC来说,它使用-mmmx、-msse等编译器开关来启用各种指令集,同时定义了对应的 __MMX__、__SSE__等宏,然后x86intrin.h会根据这些宏来声明相应的Intrinsic函数集。__MMX__、__SSE__等宏可以帮助我们判断Intrinsic函数集是否支持,但这只是GCC的专用功能。
      此外还有一些细节问题,例如某些Intrinsic函数仅在64下才能使用、有些老版本编译器的头文件缺少某个Intrinsic函数。所以我们希望有一种统一的方式来判断Intrinsic函数集的支持性。

      除了编译期间的问题外,还有运行期间的问题——
      在运行时,怎么检测当前处理器支持哪些指令集?
      虽然X86体系提供了用来检测处理器的CPUID指令,但它没有规范的Intrinsic函数,在不同的编译器上的用法不同。
      而且X86体系有很多种指令集,每种指令集具体的检测方法是略有区别的。尤其是SSE、AVX这样的SIMD指令集是需要操作系统配合才能正常使用的,所以在CPUID检查通过后,还需要进一步验证。


    二、范例讲解

    2.1 事先准备

      为了解决上面提到的问题,我编写了三个模块——
    stdint:智能支持C99的stdint.h,解决整数类型问题。最新版的地址是 http://www.cnblogs.com/zyl910/archive/2012/08/08/c99int.html
    zintrin:在编译时检测Intrinsic函数集支持性,并自动引入相关头文件、修正细节问题。最新版的地址是 http://www.cnblogs.com/zyl910/archive/2012/10/01/zintrin_v101.html
    ccpuid:在编译时检测指令集的支持性。最新版的地址是 http://www.cnblogs.com/zyl910/archive/2012/10/13/ccpuid_v103.html

      这三个模块的纯C版就是一个头文件,用起来很方便,将它们项目中,直接#include就行了。例如——

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


      因为stdint.h会被zintrin.h或ccpuid.h引用,所以不需要手动引入它。
      因为它们用到了C99整数范围常量,所以应该在程序的最前面定义__STDC_LIMIT_MACROS宏(或者可以在项目配置、编译器命令行等位置进行配置)。根据C99规范,纯C程序可以不用, 而C++程序必须定义该宏。本文为了演示,定义了该宏。


    2.2 C语言版

      我们先用C语言编写一个基本的单精度浮点数组求和函数——

    1. // 单精度浮点数组求和_基本版. 
    2. // 
    3. // result: 返回数组求和结果. 
    4. // pbuf: 数组的首地址. 
    5. // cntbuf: 数组长度. 
    6. float sumfloat_base(const float* pbuf, size_t cntbuf) 
    7.     float s = 0;    // 求和变量. 
    8.     size_t i; 
    9.     for(i=0; i<cntbuf; ++i) 
    10.     { 
    11.         s += pbuf[i]; 
    12.     } 
    13.     return s; 


      该函数很容易理解——先将返回值赋初值0,然后循环加上数组中每一项的值。


    2.3 SSE版

    2.3.1 SSE普通版

      SSE寄存器是128位的,对应__m128类型,它能一次能处理4个单精度浮点数。
      很多SSE指令要求内存地址按16字节对齐。本文为了简化,假定浮点数组的首地址是总是16字节对齐的,仅需要考虑数组长度不是4的整数倍问题。
      因使用了SSE Intrinsic函数,我们可以根据zintrin.h所提供的INTRIN_SSE宏进行条件编译。
      代码如下——

    1. #ifdef INTRIN_SSE 
    2. // 单精度浮点数组求和_SSE版. 
    3. float sumfloat_sse(const float* pbuf, size_t cntbuf) 
    4.     float s = 0;    // 求和变量. 
    5.     size_t i; 
    6.     size_t nBlockWidth = 4; // 块宽. SSE寄存器能一次处理4个float. 
    7.     size_t cntBlock = cntbuf / nBlockWidth; // 块数. 
    8.     size_t cntRem = cntbuf % nBlockWidth;   // 剩余数量. 
    9.     __m128 xfsSum = _mm_setzero_ps();   // 求和变量。[SSE] 赋初值0 
    10.     __m128 xfsLoad; // 加载. 
    11.     const float* p = pbuf;  // SSE批量处理时所用的指针. 
    12.     const float* q; // 将SSE变量上的多个数值合并时所用指针. 
    13.  
    14.     // SSE批量处理. 
    15.     for(i=0; i<cntBlock; ++i) 
    16.     { 
    17.         xfsLoad = _mm_load_ps(p);   // [SSE] 加载 
    18.         xfsSum = _mm_add_ps(xfsSum, xfsLoad);   // [SSE] 单精浮点紧缩加法 
    19.         p += nBlockWidth; 
    20.     } 
    21.     // 合并. 
    22.     q = (const float*)&xfsSum; 
    23.     s = q[0] + q[1] + q[2] + q[3]; 
    24.  
    25.     // 处理剩下的. 
    26.     for(i=0; i<cntRem; ++i) 
    27.     { 
    28.         s += p[i]; 
    29.     } 
    30.  
    31.     return s; 
    32.  
    33. #endif  // #ifdef INTRIN_SSE 


      上述代码大致可分为四个部分——
    1. 变量定义与初始化。
    2. SSE批量处理。即对前面能凑成4个一组的数据,利用SSE的128位宽度同时对4个数累加。
    3. 合并。将__m128上的多个数值合并到求和变量。因考虑某些编译器不能直接使用“.”来访问__m128变量中的数据,于是利用指针q来访问xfsSum中的数据。
    4. 处理剩下的。即对尾部不能凑成4个一组的数据,采用基本的逐项相加算法。

      上述代码总共用到了3个SSE Intrinsic函数——
    _mm_setzero_ps:对应XORPS指令。将__m128上的每一个单精度浮点数均赋0值,伪代码:for(i=0;i<4;++i) C[i]=0.0f。
    _mm_load_ps:对应MOVPS指令。从内存中对齐加载4个单精度浮点数到__m128变量,伪代码:for(i=0;i<4;++i) C[i]=_A[i]。
    _mm_add_ps:对应ADDPS指令。相加,即对2个__m128变量的4个单精度浮点数进行垂直相加,伪代码:for(i=0;i<4;++i) C[i]=A[i]+B[i]。


    2.3.2 SSE四路循环展开版

      循环展开可以降低循环开销,提高指令级并行性能。
      一般来说,四路循环展开就差不多够了。我们可以很方便的将上一节的代码改造为四路循环展开版——

    1. // 单精度浮点数组求和_SSE四路循环展开版. 
    2. float sumfloat_sse_4loop(const float* pbuf, size_t cntbuf) 
    3.     float s = 0;    // 返回值. 
    4.     size_t i; 
    5.     size_t nBlockWidth = 4*4;   // 块宽. SSE寄存器能一次处理4个float,然后循环展开4次. 
    6.     size_t cntBlock = cntbuf / nBlockWidth; // 块数. 
    7.     size_t cntRem = cntbuf % nBlockWidth;   // 剩余数量. 
    8.     __m128 xfsSum = _mm_setzero_ps();   // 求和变量。[SSE] 赋初值0 
    9.     __m128 xfsSum1 = _mm_setzero_ps(); 
    10.     __m128 xfsSum2 = _mm_setzero_ps(); 
    11.     __m128 xfsSum3 = _mm_setzero_ps(); 
    12.     __m128 xfsLoad; // 加载. 
    13.     __m128 xfsLoad1; 
    14.     __m128 xfsLoad2; 
    15.     __m128 xfsLoad3; 
    16.     const float* p = pbuf;  // SSE批量处理时所用的指针. 
    17.     const float* q; // 将SSE变量上的多个数值合并时所用指针. 
    18.  
    19.     // SSE批量处理. 
    20.     for(i=0; i<cntBlock; ++i) 
    21.     { 
    22.         xfsLoad = _mm_load_ps(p);   // [SSE] 加载. 
    23.         xfsLoad1 = _mm_load_ps(p+4); 
    24.         xfsLoad2 = _mm_load_ps(p+8); 
    25.         xfsLoad3 = _mm_load_ps(p+12); 
    26.         xfsSum = _mm_add_ps(xfsSum, xfsLoad);   // [SSE] 单精浮点紧缩加法 
    27.         xfsSum1 = _mm_add_ps(xfsSum1, xfsLoad1); 
    28.         xfsSum2 = _mm_add_ps(xfsSum2, xfsLoad2); 
    29.         xfsSum3 = _mm_add_ps(xfsSum3, xfsLoad3); 
    30.         p += nBlockWidth; 
    31.     } 
    32.     // 合并. 
    33.     xfsSum = _mm_add_ps(xfsSum, xfsSum1);   // 两两合并(0~1). 
    34.     xfsSum2 = _mm_add_ps(xfsSum2, xfsSum3); // 两两合并(2~3). 
    35.     xfsSum = _mm_add_ps(xfsSum, xfsSum2);   // 两两合并(0~3). 
    36.     q = (const float*)&xfsSum; 
    37.     s = q[0] + q[1] + q[2] + q[3]; 
    38.  
    39.     // 处理剩下的. 
    40.     for(i=0; i<cntRem; ++i) 
    41.     { 
    42.         s += p[i]; 
    43.     } 
    44.  
    45.     return s; 



    2.4 AVX版

    2.4.1 AVX普通版

      AVX寄存器是256位的,对应__m256类型,它能一次能处理8个单精度浮点数。
      很多AVX指令要求内存地址按32字节对齐。本文为了简化,假定浮点数组的首地址是总是32字节对齐的,仅需要考虑数组长度不是8的整数倍问题。
      因使用了AVX Intrinsic函数,我们可以根据zintrin.h所提供的INTRIN_AVX宏进行条件编译。

      代码如下——

    1. #ifdef INTRIN_AVX 
    2. // 单精度浮点数组求和_AVX版. 
    3. float sumfloat_avx(const float* pbuf, size_t cntbuf) 
    4.     float s = 0;    // 求和变量. 
    5.     size_t i; 
    6.     size_t nBlockWidth = 8; // 块宽. AVX寄存器能一次处理8个float. 
    7.     size_t cntBlock = cntbuf / nBlockWidth; // 块数. 
    8.     size_t cntRem = cntbuf % nBlockWidth;   // 剩余数量. 
    9.     __m256 yfsSum = _mm256_setzero_ps();    // 求和变量。[AVX] 赋初值0 
    10.     __m256 yfsLoad; // 加载. 
    11.     const float* p = pbuf;  // AVX批量处理时所用的指针. 
    12.     const float* q; // 将AVX变量上的多个数值合并时所用指针. 
    13.  
    14.     // AVX批量处理. 
    15.     for(i=0; i<cntBlock; ++i) 
    16.     { 
    17.         yfsLoad = _mm256_load_ps(p);    // [AVX] 加载 
    18.         yfsSum = _mm256_add_ps(yfsSum, yfsLoad);    // [AVX] 单精浮点紧缩加法 
    19.         p += nBlockWidth; 
    20.     } 
    21.     // 合并. 
    22.     q = (const float*)&yfsSum; 
    23.     s = q[0] + q[1] + q[2] + q[3] + q[4] + q[5] + q[6] + q[7]; 
    24.  
    25.     // 处理剩下的. 
    26.     for(i=0; i<cntRem; ++i) 
    27.     { 
    28.         s += p[i]; 
    29.     } 
    30.  
    31.     return s; 
    32.  
    33. #endif  // #ifdef INTRIN_AVX 


      由上可见,将SSE Intrinsic代码(sumfloat_sse)升级为 AVX Intrinsic代码(sumfloat_avx)是很容易的——
    1. 升级数据类型,将__m128升级成了__m256。
    2. 升级Intrinsic函数,在函数名中加入255。例如_mm_setzero_ps、_mm_load_ps、_mm_add_ps,对应的AVX版函数是 _mm256_setzero_ps、_mm256_load_ps、_mm256_add_ps。
    3. 因位宽翻倍,地址计算与数据合并的代码需稍加改动。

      当使用VC2010编译含有AVX的代码时,VC会提醒你——
    warning C4752: 发现 Intel(R) 高级矢量扩展;请考虑使用 /arch:AVX

      目前“/arch:AVX”尚未整合到项目属性的“C++\代码生成\启用增强指令集”中,需要手动在项目属性的“C++\命令行”的附加选项中加上“/arch:AVX”——

    详见MSDN——
    http://msdn.microsoft.com/zh-cn/library/7t5yh4fd(v=vs.100).aspx
    在 Visual Studio 中设置 /arch:AVX 编译器选项
    1.打开项目的“属性页”对话框。 有关更多信息,请参见 如何:打开项目属性页。
    2.单击“C/C++”文件夹。
    3.单击“命令行”属性页。
    4.在“附加选项”框中添加 /arch:AVX。


    2.4.2 AVX四路循环展开版

      同样的,我们可以编写AVX四路循环展开版——

    1. // 单精度浮点数组求和_AVX四路循环展开版. 
    2. float sumfloat_avx_4loop(const float* pbuf, size_t cntbuf) 
    3.     float s = 0;    // 求和变量. 
    4.     size_t i; 
    5.     size_t nBlockWidth = 8*4;   // 块宽. AVX寄存器能一次处理8个float,然后循环展开4次. 
    6.     size_t cntBlock = cntbuf / nBlockWidth; // 块数. 
    7.     size_t cntRem = cntbuf % nBlockWidth;   // 剩余数量. 
    8.     __m256 yfsSum = _mm256_setzero_ps();    // 求和变量。[AVX] 赋初值0 
    9.     __m256 yfsSum1 = _mm256_setzero_ps(); 
    10.     __m256 yfsSum2 = _mm256_setzero_ps(); 
    11.     __m256 yfsSum3 = _mm256_setzero_ps(); 
    12.     __m256 yfsLoad; // 加载. 
    13.     __m256 yfsLoad1; 
    14.     __m256 yfsLoad2; 
    15.     __m256 yfsLoad3; 
    16.     const float* p = pbuf;  // AVX批量处理时所用的指针. 
    17.     const float* q; // 将AVX变量上的多个数值合并时所用指针. 
    18.  
    19.     // AVX批量处理. 
    20.     for(i=0; i<cntBlock; ++i) 
    21.     { 
    22.         yfsLoad = _mm256_load_ps(p);    // [AVX] 加载. 
    23.         yfsLoad1 = _mm256_load_ps(p+8); 
    24.         yfsLoad2 = _mm256_load_ps(p+16); 
    25.         yfsLoad3 = _mm256_load_ps(p+24); 
    26.         yfsSum = _mm256_add_ps(yfsSum, yfsLoad);    // [AVX] 单精浮点紧缩加法 
    27.         yfsSum1 = _mm256_add_ps(yfsSum1, yfsLoad1); 
    28.         yfsSum2 = _mm256_add_ps(yfsSum2, yfsLoad2); 
    29.         yfsSum3 = _mm256_add_ps(yfsSum3, yfsLoad3); 
    30.         p += nBlockWidth; 
    31.     } 
    32.     // 合并. 
    33.     yfsSum = _mm256_add_ps(yfsSum, yfsSum1);    // 两两合并(0~1). 
    34.     yfsSum2 = _mm256_add_ps(yfsSum2, yfsSum3);  // 两两合并(2~3). 
    35.     yfsSum = _mm256_add_ps(yfsSum, yfsSum2);    // 两两合并(0~3). 
    36.     q = (const float*)&yfsSum; 
    37.     s = q[0] + q[1] + q[2] + q[3] + q[4] + q[5] + q[6] + q[7]; 
    38.  
    39.     // 处理剩下的. 
    40.     for(i=0; i<cntRem; ++i) 
    41.     { 
    42.         s += p[i]; 
    43.     } 
    44.  
    45.     return s; 


    2.5 测试框架

    2.5.1 测试所用的数组

      首先考虑一下测试所用的数组的长度应该是多少比较好。
      为了避免内存带宽问题,这个数组最好能放在L1 Data Cache中。现在的处理器的L1 Data Cache一般是32KB,为了保险最好再除以2,那么数组的长度应该是 32KB/(2*sizeof(float))=4096。
      其次考虑内存对齐问题,avx要求32字节对齐。我们可以定义一个ATTR_ALIGN宏来统一处理变量的内存对齐问题。
      该数组定义如下——

    1. // 变量对齐. 
    2. #ifndef ATTR_ALIGN 
    3. #  if defined(__GNUC__) // GCC 
    4. #    define ATTR_ALIGN(n)   __attribute__((aligned(n))) 
    5. #  else // 否则使用VC格式. 
    6. #    define ATTR_ALIGN(n)   __declspec(align(n)) 
    7. #  endif 
    8. #endif  // #ifndef ATTR_ALIGN 
    9.  
    10.  
    11. #define BUFSIZE 4096    // = 32KB{L1 Cache} / (2 * sizeof(float)) 
    12. ATTR_ALIGN(32) float buf[BUFSIZE]; 



    2.5.2 测试函数

      如果为每一个函数都编写一套测试代码,那不仅代码量大,而且不易维护。
      可以考虑利用函数指针来实现一套测试框架。
      因sumfloat_base等函数的签名是一致的,于是可以定义这样的一种函数指针——
    // 测试时的函数类型
    typedef float (*TESTPROC)(const float* pbuf, size_t cntbuf);

      然后再编写一个对TESTPROC函数指针进行测试的函数——

    1. // 进行测试 
    2. void runTest(const char* szname, TESTPROC proc) 
    3.     const int testloop = 4000;  // 重复运算几次延长时间,避免计时精度问题. 
    4.     const clock_t TIMEOUT = CLOCKS_PER_SEC/2;   // 最短测试时间. 
    5.     int i,j,k; 
    6.     clock_t tm0, dt;    // 存储时间. 
    7.     double mps; // M/s. 
    8.     double mps_good = 0;    // 最佳M/s. 因线程切换会导致的数值波动, 于是选取最佳值. 
    9.     volatile float n=0; // 避免内循环被优化. 
    10.     for(i=1; i<=3; ++i)  // 多次测试. 
    11.     { 
    12.         tm0 = clock(); 
    13.         // main 
    14.         k=0; 
    15.         do 
    16.         { 
    17.             for(j=1; j<=testloop; ++j)   // 重复运算几次延长时间,避免计时开销带来的影响. 
    18.             { 
    19.                 n = proc(buf, BUFSIZE); // 避免内循环被编译优化消掉. 
    20.             } 
    21.             ++k; 
    22.             dt = clock() - tm0; 
    23.         }while(dt<TIMEOUT); 
    24.         // show 
    25.         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 . 
    26.         if (mps_good<mps)    mps_good=mps;   // 选取最佳值. 
    27.         //printf("%s:\t%.0f M/s\t//%f\n", szname, mps, n); 
    28.     } 
    29.     printf("%s:\t%.0f M/s\t//%f\n", szname, mps_good, n); 


      j是最内层的循环,负责多次调用TESTPROC函数指针。如果每调用一次TESTPROC函数指针后又调用clock函数,那会带来较大的计时开销,影响评测成绩。
      k循环负责检测超时。当发现超过预定时限,便计算mps,即每秒钟处理了多少百万个单精度浮点数。然后存储最佳的mps。
      i是最外层循环的循环变量,循环3次然后报告最佳值。


    2.5.3 进行测试

      在进行测试之前,需要对buf数组进行初始化,将数组元素赋随机值——

    1. // init buf 
    2. srand( (unsigned)time( NULL ) ); 
    3. for (i = 0; i < BUFSIZE; i++) buf[i] = (float)(rand() & 0x3f);   // 使用&0x3f是为了让求和后的数值不会超过float类型的有效位数,便于观察结果是否正确. 


      然后可以开始测试了——

    1.     // test 
    2.     runTest("sumfloat_base", sumfloat_base);    // 单精度浮点数组求和_基本版. 
    3. #ifdef INTRIN_SSE 
    4.     if (simd_sse_level(NULL) >= SIMD_SSE_1) 
    5.     { 
    6.         runTest("sumfloat_sse", sumfloat_sse);  // 单精度浮点数组求和_SSE版. 
    7.         runTest("sumfloat_sse_4loop", sumfloat_sse_4loop);  // 单精度浮点数组求和_SSE四路循环展开版. 
    8.     } 
    9. #endif  // #ifdef INTRIN_SSE 
    10. #ifdef INTRIN_AVX 
    11.     if (simd_avx_level(NULL) >= SIMD_AVX_1) 
    12.     { 
    13.         runTest("sumfloat_avx", sumfloat_avx);  // 单精度浮点数组求和_SSE版. 
    14.         runTest("sumfloat_avx_4loop", sumfloat_avx_4loop);  // 单精度浮点数组求和_SSE四路循环展开版. 
    15.     } 
    16. #endif  // #ifdef INTRIN_AVX 



    2.6 杂项

      为了方便对比测试,可以在程序启动时显示程序版本、编译器名称、CPU型号信息。即在main函数中加上——

    1. char szBuf[64]; 
    2. int i; 
    3.  
    4. printf("simdsumfloat v1.00 (%dbit)\n", INTRIN_WORDSIZE); 
    5. printf("Compiler: %s\n", COMPILER_NAME); 
    6. cpu_getbrand(szBuf); 
    7. printf("CPU:\t%s\n", szBuf); 
    8. printf("\n"); 


      INTRIN_WORDSIZE 宏是 zintrin.h 提供的,为当前机器的字长。
      cpu_getbrand是 ccpuid.h 提供的,用于获得CPU型号字符串。
      COMPILER_NAME 是一个用来获得编译器名称的宏,它的详细定义是——

    1. // Compiler name 
    2. #define MACTOSTR(x) #x 
    3. #define MACROVALUESTR(x)    MACTOSTR(x) 
    4. #if defined(__ICL)  // Intel C++ 
    5. #  if defined(__VERSION__) 
    6. #    define COMPILER_NAME   "Intel C++ " __VERSION__ 
    7. #  elif defined(__INTEL_COMPILER_BUILD_DATE) 
    8. #    define COMPILER_NAME   "Intel C++ (" MACROVALUESTR(__INTEL_COMPILER_BUILD_DATE) ")" 
    9. #  else 
    10. #    define COMPILER_NAME   "Intel C++" 
    11. #  endif    // #  if defined(__VERSION__) 
    12. #elif defined(_MSC_VER) // Microsoft VC++ 
    13. #  if defined(_MSC_FULL_VER) 
    14. #    define COMPILER_NAME   "Microsoft VC++ (" MACROVALUESTR(_MSC_FULL_VER) ")" 
    15. #  elif defined(_MSC_VER) 
    16. #    define COMPILER_NAME   "Microsoft VC++ (" MACROVALUESTR(_MSC_VER) ")" 
    17. #  else 
    18. #    define COMPILER_NAME   "Microsoft VC++" 
    19. #  endif    // #  if defined(_MSC_FULL_VER) 
    20. #elif defined(__GNUC__) // GCC 
    21. #  if defined(__CYGWIN__) 
    22. #    define COMPILER_NAME   "GCC(Cygmin) " __VERSION__ 
    23. #  elif defined(__MINGW32__) 
    24. #    define COMPILER_NAME   "GCC(MinGW) " __VERSION__ 
    25. #  else 
    26. #    define COMPILER_NAME   "GCC " __VERSION__ 
    27. #  endif    // #  if defined(_MSC_FULL_VER) 
    28. #else 
    29. #  define COMPILER_NAME "Unknown Compiler" 
    30. #endif  // #if defined(__ICL)   // Intel C++ 


    三、全部代码

    3.1 simdsumfloat.c

      全部代码——

    1. #define __STDC_LIMIT_MACROS 1   // C99整数范围常量. [纯C程序可以不用, 而C++程序必须定义该宏.] 
    2.  
    3. #include <stdlib.h> 
    4. #include <stdio.h> 
    5. #include <time.h> 
    6.  
    7. #include "zintrin.h" 
    8. #include "ccpuid.h" 
    9.  
    10.  
    11. // Compiler name 
    12. #define MACTOSTR(x) #x 
    13. #define MACROVALUESTR(x)    MACTOSTR(x) 
    14. #if defined(__ICL)  // Intel C++ 
    15. #  if defined(__VERSION__) 
    16. #    define COMPILER_NAME   "Intel C++ " __VERSION__ 
    17. #  elif defined(__INTEL_COMPILER_BUILD_DATE) 
    18. #    define COMPILER_NAME   "Intel C++ (" MACROVALUESTR(__INTEL_COMPILER_BUILD_DATE) ")" 
    19. #  else 
    20. #    define COMPILER_NAME   "Intel C++" 
    21. #  endif    // #  if defined(__VERSION__) 
    22. #elif defined(_MSC_VER) // Microsoft VC++ 
    23. #  if defined(_MSC_FULL_VER) 
    24. #    define COMPILER_NAME   "Microsoft VC++ (" MACROVALUESTR(_MSC_FULL_VER) ")" 
    25. #  elif defined(_MSC_VER) 
    26. #    define COMPILER_NAME   "Microsoft VC++ (" MACROVALUESTR(_MSC_VER) ")" 
    27. #  else 
    28. #    define COMPILER_NAME   "Microsoft VC++" 
    29. #  endif    // #  if defined(_MSC_FULL_VER) 
    30. #elif defined(__GNUC__) // GCC 
    31. #  if defined(__CYGWIN__) 
    32. #    define COMPILER_NAME   "GCC(Cygmin) " __VERSION__ 
    33. #  elif defined(__MINGW32__) 
    34. #    define COMPILER_NAME   "GCC(MinGW) " __VERSION__ 
    35. #  else 
    36. #    define COMPILER_NAME   "GCC " __VERSION__ 
    37. #  endif    // #  if defined(_MSC_FULL_VER) 
    38. #else 
    39. #  define COMPILER_NAME "Unknown Compiler" 
    40. #endif  // #if defined(__ICL)   // Intel C++ 
    41.  
    42.  
    43. ////////////////////////////////////////////////// 
    44. // sumfloat: 单精度浮点数组求和的函数 
    45. ////////////////////////////////////////////////// 
    46.  
    47. // 单精度浮点数组求和_基本版. 
    48. // 
    49. // result: 返回数组求和结果. 
    50. // pbuf: 数组的首地址. 
    51. // cntbuf: 数组长度. 
    52. float sumfloat_base(const float* pbuf, size_t cntbuf) 
    53.     float s = 0;    // 求和变量. 
    54.     size_t i; 
    55.     for(i=0; i<cntbuf; ++i) 
    56.     { 
    57.         s += pbuf[i]; 
    58.     } 
    59.     return s; 
    60.  
    61. #ifdef INTRIN_SSE 
    62. // 单精度浮点数组求和_SSE版. 
    63. float sumfloat_sse(const float* pbuf, size_t cntbuf) 
    64.     float s = 0;    // 求和变量. 
    65.     size_t i; 
    66.     size_t nBlockWidth = 4; // 块宽. SSE寄存器能一次处理4个float. 
    67.     size_t cntBlock = cntbuf / nBlockWidth; // 块数. 
    68.     size_t cntRem = cntbuf % nBlockWidth;   // 剩余数量. 
    69.     __m128 xfsSum = _mm_setzero_ps();   // 求和变量。[SSE] 赋初值0 
    70.     __m128 xfsLoad; // 加载. 
    71.     const float* p = pbuf;  // SSE批量处理时所用的指针. 
    72.     const float* q; // 将SSE变量上的多个数值合并时所用指针. 
    73.  
    74.     // SSE批量处理. 
    75.     for(i=0; i<cntBlock; ++i) 
    76.     { 
    77.         xfsLoad = _mm_load_ps(p);   // [SSE] 加载 
    78.         xfsSum = _mm_add_ps(xfsSum, xfsLoad);   // [SSE] 单精浮点紧缩加法 
    79.         p += nBlockWidth; 
    80.     } 
    81.     // 合并. 
    82.     q = (const float*)&xfsSum; 
    83.     s = q[0] + q[1] + q[2] + q[3]; 
    84.  
    85.     // 处理剩下的. 
    86.     for(i=0; i<cntRem; ++i) 
    87.     { 
    88.         s += p[i]; 
    89.     } 
    90.  
    91.     return s; 
    92.  
    93. // 单精度浮点数组求和_SSE四路循环展开版. 
    94. float sumfloat_sse_4loop(const float* pbuf, size_t cntbuf) 
    95.     float s = 0;    // 返回值. 
    96.     size_t i; 
    97.     size_t nBlockWidth = 4*4;   // 块宽. SSE寄存器能一次处理4个float,然后循环展开4次. 
    98.     size_t cntBlock = cntbuf / nBlockWidth; // 块数. 
    99.     size_t cntRem = cntbuf % nBlockWidth;   // 剩余数量. 
    100.     __m128 xfsSum = _mm_setzero_ps();   // 求和变量。[SSE] 赋初值0 
    101.     __m128 xfsSum1 = _mm_setzero_ps(); 
    102.     __m128 xfsSum2 = _mm_setzero_ps(); 
    103.     __m128 xfsSum3 = _mm_setzero_ps(); 
    104.     __m128 xfsLoad; // 加载. 
    105.     __m128 xfsLoad1; 
    106.     __m128 xfsLoad2; 
    107.     __m128 xfsLoad3; 
    108.     const float* p = pbuf;  // SSE批量处理时所用的指针. 
    109.     const float* q; // 将SSE变量上的多个数值合并时所用指针. 
    110.  
    111.     // SSE批量处理. 
    112.     for(i=0; i<cntBlock; ++i) 
    113.     { 
    114.         xfsLoad = _mm_load_ps(p);   // [SSE] 加载. 
    115.         xfsLoad1 = _mm_load_ps(p+4); 
    116.         xfsLoad2 = _mm_load_ps(p+8); 
    117.         xfsLoad3 = _mm_load_ps(p+12); 
    118.         xfsSum = _mm_add_ps(xfsSum, xfsLoad);   // [SSE] 单精浮点紧缩加法 
    119.         xfsSum1 = _mm_add_ps(xfsSum1, xfsLoad1); 
    120.         xfsSum2 = _mm_add_ps(xfsSum2, xfsLoad2); 
    121.         xfsSum3 = _mm_add_ps(xfsSum3, xfsLoad3); 
    122.         p += nBlockWidth; 
    123.     } 
    124.     // 合并. 
    125.     xfsSum = _mm_add_ps(xfsSum, xfsSum1);   // 两两合并(0~1). 
    126.     xfsSum2 = _mm_add_ps(xfsSum2, xfsSum3); // 两两合并(2~3). 
    127.     xfsSum = _mm_add_ps(xfsSum, xfsSum2);   // 两两合并(0~3). 
    128.     q = (const float*)&xfsSum; 
    129.     s = q[0] + q[1] + q[2] + q[3]; 
    130.  
    131.     // 处理剩下的. 
    132.     for(i=0; i<cntRem; ++i) 
    133.     { 
    134.         s += p[i]; 
    135.     } 
    136.  
    137.     return s; 
    138. #endif  // #ifdef INTRIN_SSE 
    139.  
    140.  
    141. #ifdef INTRIN_AVX 
    142. // 单精度浮点数组求和_AVX版. 
    143. float sumfloat_avx(const float* pbuf, size_t cntbuf) 
    144.     float s = 0;    // 求和变量. 
    145.     size_t i; 
    146.     size_t nBlockWidth = 8; // 块宽. AVX寄存器能一次处理8个float. 
    147.     size_t cntBlock = cntbuf / nBlockWidth; // 块数. 
    148.     size_t cntRem = cntbuf % nBlockWidth;   // 剩余数量. 
    149.     __m256 yfsSum = _mm256_setzero_ps();    // 求和变量。[AVX] 赋初值0 
    150.     __m256 yfsLoad; // 加载. 
    151.     const float* p = pbuf;  // AVX批量处理时所用的指针. 
    152.     const float* q; // 将AVX变量上的多个数值合并时所用指针. 
    153.  
    154.     // AVX批量处理. 
    155.     for(i=0; i<cntBlock; ++i) 
    156.     { 
    157.         yfsLoad = _mm256_load_ps(p);    // [AVX] 加载 
    158.         yfsSum = _mm256_add_ps(yfsSum, yfsLoad);    // [AVX] 单精浮点紧缩加法 
    159.         p += nBlockWidth; 
    160.     } 
    161.     // 合并. 
    162.     q = (const float*)&yfsSum; 
    163.     s = q[0] + q[1] + q[2] + q[3] + q[4] + q[5] + q[6] + q[7]; 
    164.  
    165.     // 处理剩下的. 
    166.     for(i=0; i<cntRem; ++i) 
    167.     { 
    168.         s += p[i]; 
    169.     } 
    170.  
    171.     return s; 
    172.  
    173. // 单精度浮点数组求和_AVX四路循环展开版. 
    174. float sumfloat_avx_4loop(const float* pbuf, size_t cntbuf) 
    175.     float s = 0;    // 求和变量. 
    176.     size_t i; 
    177.     size_t nBlockWidth = 8*4;   // 块宽. AVX寄存器能一次处理8个float,然后循环展开4次. 
    178.     size_t cntBlock = cntbuf / nBlockWidth; // 块数. 
    179.     size_t cntRem = cntbuf % nBlockWidth;   // 剩余数量. 
    180.     __m256 yfsSum = _mm256_setzero_ps();    // 求和变量。[AVX] 赋初值0 
    181.     __m256 yfsSum1 = _mm256_setzero_ps(); 
    182.     __m256 yfsSum2 = _mm256_setzero_ps(); 
    183.     __m256 yfsSum3 = _mm256_setzero_ps(); 
    184.     __m256 yfsLoad; // 加载. 
    185.     __m256 yfsLoad1; 
    186.     __m256 yfsLoad2; 
    187.     __m256 yfsLoad3; 
    188.     const float* p = pbuf;  // AVX批量处理时所用的指针. 
    189.     const float* q; // 将AVX变量上的多个数值合并时所用指针. 
    190.  
    191.     // AVX批量处理. 
    192.     for(i=0; i<cntBlock; ++i) 
    193.     { 
    194.         yfsLoad = _mm256_load_ps(p);    // [AVX] 加载. 
    195.         yfsLoad1 = _mm256_load_ps(p+8); 
    196.         yfsLoad2 = _mm256_load_ps(p+16); 
    197.         yfsLoad3 = _mm256_load_ps(p+24); 
    198.         yfsSum = _mm256_add_ps(yfsSum, yfsLoad);    // [AVX] 单精浮点紧缩加法 
    199.         yfsSum1 = _mm256_add_ps(yfsSum1, yfsLoad1); 
    200.         yfsSum2 = _mm256_add_ps(yfsSum2, yfsLoad2); 
    201.         yfsSum3 = _mm256_add_ps(yfsSum3, yfsLoad3); 
    202.         p += nBlockWidth; 
    203.     } 
    204.     // 合并. 
    205.     yfsSum = _mm256_add_ps(yfsSum, yfsSum1);    // 两两合并(0~1). 
    206.     yfsSum2 = _mm256_add_ps(yfsSum2, yfsSum3);  // 两两合并(2~3). 
    207.     yfsSum = _mm256_add_ps(yfsSum, yfsSum2);    // 两两合并(0~3). 
    208.     q = (const float*)&yfsSum; 
    209.     s = q[0] + q[1] + q[2] + q[3] + q[4] + q[5] + q[6] + q[7]; 
    210.  
    211.     // 处理剩下的. 
    212.     for(i=0; i<cntRem; ++i) 
    213.     { 
    214.         s += p[i]; 
    215.     } 
    216.  
    217.     return s; 
    218.  
    219. #endif  // #ifdef INTRIN_AVX 
    220.  
    221.  
    222.  
    223. ////////////////////////////////////////////////// 
    224. // main 
    225. ////////////////////////////////////////////////// 
    226.  
    227. // 变量对齐. 
    228. #ifndef ATTR_ALIGN 
    229. #  if defined(__GNUC__) // GCC 
    230. #    define ATTR_ALIGN(n)   __attribute__((aligned(n))) 
    231. #  else // 否则使用VC格式. 
    232. #    define ATTR_ALIGN(n)   __declspec(align(n)) 
    233. #  endif 
    234. #endif  // #ifndef ATTR_ALIGN 
    235.  
    236.  
    237. #define BUFSIZE 4096    // = 32KB{L1 Cache} / (2 * sizeof(float)) 
    238. ATTR_ALIGN(32) float buf[BUFSIZE]; 
    239.  
    240. // 测试时的函数类型 
    241. typedef float (*TESTPROC)(const float* pbuf, size_t cntbuf); 
    242.  
    243. // 进行测试 
    244. void runTest(const char* szname, TESTPROC proc) 
    245.     const int testloop = 4000;  // 重复运算几次延长时间,避免计时精度问题. 
    246.     const clock_t TIMEOUT = CLOCKS_PER_SEC/2;   // 最短测试时间. 
    247.     int i,j,k; 
    248.     clock_t tm0, dt;    // 存储时间. 
    249.     double mps; // M/s. 
    250.     double mps_good = 0;    // 最佳M/s. 因线程切换会导致的数值波动, 于是选取最佳值. 
    251.     volatile float n=0; // 避免内循环被优化. 
    252.     for(i=1; i<=3; ++i)  // 多次测试. 
    253.     { 
    254.         tm0 = clock(); 
    255.         // main 
    256.         k=0; 
    257.         do 
    258.         { 
    259.             for(j=1; j<=testloop; ++j)   // 重复运算几次延长时间,避免计时开销带来的影响. 
    260.             { 
    261.                 n = proc(buf, BUFSIZE); // 避免内循环被编译优化消掉. 
    262.             } 
    263.             ++k; 
    264.             dt = clock() - tm0; 
    265.         }while(dt<TIMEOUT); 
    266.         // show 
    267.         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 . 
    268.         if (mps_good<mps)    mps_good=mps;   // 选取最佳值. 
    269.         //printf("%s:\t%.0f M/s\t//%f\n", szname, mps, n); 
    270.     } 
    271.     printf("%s:\t%.0f M/s\t//%f\n", szname, mps_good, n); 
    272.  
    273. int main(int argc, char* argv[]) 
    274.     char szBuf[64]; 
    275.     int i; 
    276.  
    277.     printf("simdsumfloat v1.00 (%dbit)\n", INTRIN_WORDSIZE); 
    278.     printf("Compiler: %s\n", COMPILER_NAME); 
    279.     cpu_getbrand(szBuf); 
    280.     printf("CPU:\t%s\n", szBuf); 
    281.     printf("\n"); 
    282.  
    283.     // init buf 
    284.     srand( (unsigned)time( NULL ) ); 
    285.     for (i = 0; i < BUFSIZE; i++) buf[i] = (float)(rand() & 0x3f);   // 使用&0x3f是为了让求和后的数值不会超过float类型的有效位数,便于观察结果是否正确. 
    286.  
    287.     // test 
    288.     runTest("sumfloat_base", sumfloat_base);    // 单精度浮点数组求和_基本版. 
    289. #ifdef INTRIN_SSE 
    290.     if (simd_sse_level(NULL) >= SIMD_SSE_1) 
    291.     { 
    292.         runTest("sumfloat_sse", sumfloat_sse);  // 单精度浮点数组求和_SSE版. 
    293.         runTest("sumfloat_sse_4loop", sumfloat_sse_4loop);  // 单精度浮点数组求和_SSE四路循环展开版. 
    294.     } 
    295. #endif  // #ifdef INTRIN_SSE 
    296. #ifdef INTRIN_AVX 
    297.     if (simd_avx_level(NULL) >= SIMD_AVX_1) 
    298.     { 
    299.         runTest("sumfloat_avx", sumfloat_avx);  // 单精度浮点数组求和_AVX版. 
    300.         runTest("sumfloat_avx_4loop", sumfloat_avx_4loop);  // 单精度浮点数组求和_AVX四路循环展开版. 
    301.     } 
    302. #endif  // #ifdef INTRIN_AVX 
    303.  
    304.     return 0; 



    3.2 makefile

      全部代码——

    1. # flags 
    2. CC = g++ 
    3. CFS = -Wall -msse 
    4.  
    5. # args 
    6. RELEASE =0 
    7. BITS = 
    8. CFLAGS = 
    9.  
    10. # [args] 生成模式. 0代表debug模式, 1代表release模式. make RELEASE=1. 
    11. ifeq ($(RELEASE),0) 
    12.     # debug 
    13.     CFS += -g 
    14. else 
    15.     # release 
    16.     CFS += -O3 -DNDEBUG 
    17.     //CFS += -O3 -g -DNDEBUG 
    18. endif 
    19.  
    20. # [args] 程序位数. 32代表32位程序, 64代表64位程序, 其他默认. make BITS=32. 
    21. ifeq ($(BITS),32) 
    22.     CFS += -m32 
    23. else 
    24.     ifeq ($(BITS),64) 
    25.         CFS += -m64 
    26.     else 
    27.     endif 
    28. endif 
    29.  
    30. # [args] 使用 CFLAGS 添加新的参数. make CFLAGS="-mavx". 
    31. CFS += $(CFLAGS) 
    32.  
    33.  
    34. .PHONY : all clean 
    35.  
    36. # files 
    37. TARGETS = simdsumfloat 
    38. OBJS = simdsumfloat.o 
    39.  
    40. all : $(TARGETS) 
    41.  
    42. simdsumfloat : $(OBJS) 
    43.     $(CC) $(CFS) -o $@ $^ 
    44.  
    45.  
    46. simdsumfloat.o : simdsumfloat.c zintrin.h ccpuid.h 
    47.     $(CC) $(CFS) -c $< 
    48.  
    49.  
    50. clean : 
    51.     rm -f $(OBJS) $(TARGETS) $(addsuffix .exe,$(TARGETS)) 



    四、编译测试

    4.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版。



    4.2 测试

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

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

      然后分别运行VC与GCC编译的Release版可执行文件,即以下4个程序——
    exe\simdsumfloat_vc32.exe:VC2010 SP1 编译的32位程序,/O2 /arch:SSE2。
    exe\simdsumfloat_vc64.exe:VC2010 SP1 编译的64位程序,/O2 /arch:AVX。
    exe\simdsumfloat_gcc32.exe:GCC 4.7.1(TDM-GCC(MinGW-w64)) 编译的32位程序,-O3 -mavx。
    exe\simdsumfloat_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


    源码下载——
    https://files.cnblogs.com/zyl910/simdsumfloat.rar

    原帖地址:http://blog.csdn.net/zyl910/article/details/8100744

  • 相关阅读:
    3.2.2.5 BRE运算符优先级
    随机场(Random field)
    D-Separation(D分离)-PRML-8.22-Graphical Model 五 18 by 小军
    CVPR 2013 关于图像/场景分类(classification)的文章paper list
    Introduction to One-class Support Vector Machines
    SVM学习资料
    MIT牛人解说数学体系
    牛顿法与拟牛顿法学习笔记(五)L-BFGS 算法
    牛顿法与拟牛顿法学习笔记(四)BFGS 算法
    牛顿法与拟牛顿法学习笔记(三)DFP 算法
  • 原文地址:https://www.cnblogs.com/catkins/p/5270731.html
Copyright © 2020-2023  润新知