作者：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。而很多時候我們為了兼容舊代碼，不得不使用VC6等老版本的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就行了。例如——
#define __STDC_LIMIT_MACROS 1 // C99整數范圍常量. [純C程序可以不用, 而C++程序必須定義該宏.] #include "zintrin.h" #include "ccpuid.h"
因為stdint.h會被zintrin.h或ccpuid.h引用，所以不需要手動引入它。
因為它們用到了C99整數范圍常量，所以應該在程序的最前面定義__STDC_LIMIT_MACROS宏（或者可以在項目配置、編譯器命令行等位置進行配置）。根據C99規范，純C程序可以不用, 而C++程序必須定義該宏。本文為了演示，定義了該宏。
2.2 C語言版
我們先用C語言編寫一個基本的單精度浮點數組求和函數——
// 單精度浮點數組求和_基本版. // // result: 返回數組求和結果. // pbuf: 數組的首地址. // cntbuf: 數組長度. float sumfloat_base(const float* pbuf, size_t cntbuf) { float s = 0; // 求和變量. size_t i; for(i=0; i<cntbuf; ++i) { s += pbuf[i]; } return s; }
該函數很容易理解——先將返回值賦初值0，然后循環加上數組中每一項的值。
2.3 SSE版
2.3.1 SSE普通版
SSE寄存器是128位的，對應__m128類型，它能一次能處理4個單精度浮點數。
很多SSE指令要求內存地址按16字節對齊。本文為了簡化，假定浮點數組的首地址是總是16字節對齊的，僅需要考慮數組長度不是4的整數倍問題。
因使用了SSE Intrinsic函數，我們可以根據zintrin.h所提供的INTRIN_SSE宏進行條件編譯。
代碼如下——
#ifdef INTRIN_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; } #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四路循環展開版
循環展開可以降低循環開銷，提高指令級並行性能。
一般來說，四路循環展開就差不多夠了。我們可以很方便的將上一節的代碼改造為四路循環展開版——
// 單精度浮點數組求和_SSE四路循環展開版. float sumfloat_sse_4loop(const float* pbuf, size_t cntbuf) { float s = 0; // 返回值. size_t i; size_t nBlockWidth = 4*4; // 塊寬. SSE寄存器能一次處理4個float，然后循環展開4次. size_t cntBlock = cntbuf / nBlockWidth; // 塊數. size_t cntRem = cntbuf % nBlockWidth; // 剩余數量. __m128 xfsSum = _mm_setzero_ps(); // 求和變量。[SSE] 賦初值0 __m128 xfsSum1 = _mm_setzero_ps(); __m128 xfsSum2 = _mm_setzero_ps(); __m128 xfsSum3 = _mm_setzero_ps(); __m128 xfsLoad; // 加載. __m128 xfsLoad1; __m128 xfsLoad2; __m128 xfsLoad3; const float* p = pbuf; // SSE批量處理時所用的指針. const float* q; // 將SSE變量上的多個數值合並時所用指針. // SSE批量處理. for(i=0; i<cntBlock; ++i) { xfsLoad = _mm_load_ps(p); // [SSE] 加載. xfsLoad1 = _mm_load_ps(p+4); xfsLoad2 = _mm_load_ps(p+8); xfsLoad3 = _mm_load_ps(p+12); xfsSum = _mm_add_ps(xfsSum, xfsLoad); // [SSE] 單精浮點緊縮加法 xfsSum1 = _mm_add_ps(xfsSum1, xfsLoad1); xfsSum2 = _mm_add_ps(xfsSum2, xfsLoad2); xfsSum3 = _mm_add_ps(xfsSum3, xfsLoad3); p += nBlockWidth; } // 合並. xfsSum = _mm_add_ps(xfsSum, xfsSum1); // 兩兩合並(0~1). xfsSum2 = _mm_add_ps(xfsSum2, xfsSum3); // 兩兩合並(2~3). xfsSum = _mm_add_ps(xfsSum, xfsSum2); // 兩兩合並(0~3). q = (const float*)&xfsSum; s = q[0] + q[1] + q[2] + q[3]; // 處理剩下的. for(i=0; i<cntRem; ++i) { s += p[i]; } return s; }
2.4 AVX版
2.4.1 AVX普通版
AVX寄存器是256位的，對應__m256類型，它能一次能處理8個單精度浮點數。
很多AVX指令要求內存地址按32字節對齊。本文為了簡化，假定浮點數組的首地址是總是32字節對齊的，僅需要考慮數組長度不是8的整數倍問題。
因使用了AVX Intrinsic函數，我們可以根據zintrin.h所提供的INTRIN_AVX宏進行條件編譯。
代碼如下——
#ifdef INTRIN_AVX // 單精度浮點數組求和_AVX版. float sumfloat_avx(const float* pbuf, size_t cntbuf) { float s = 0; // 求和變量. size_t i; size_t nBlockWidth = 8; // 塊寬. AVX寄存器能一次處理8個float. size_t cntBlock = cntbuf / nBlockWidth; // 塊數. size_t cntRem = cntbuf % nBlockWidth; // 剩余數量. __m256 yfsSum = _mm256_setzero_ps(); // 求和變量。[AVX] 賦初值0 __m256 yfsLoad; // 加載. const float* p = pbuf; // AVX批量處理時所用的指針. const float* q; // 將AVX變量上的多個數值合並時所用指針. // AVX批量處理. for(i=0; i<cntBlock; ++i) { yfsLoad = _mm256_load_ps(p); // [AVX] 加載 yfsSum = _mm256_add_ps(yfsSum, yfsLoad); // [AVX] 單精浮點緊縮加法 p += nBlockWidth; } // 合並. q = (const float*)&yfsSum; s = q[0] + q[1] + q[2] + q[3] + q[4] + q[5] + q[6] + q[7]; // 處理剩下的. for(i=0; i<cntRem; ++i) { s += p[i]; } return s; } #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四路循環展開版——
// 單精度浮點數組求和_AVX四路循環展開版. float sumfloat_avx_4loop(const float* pbuf, size_t cntbuf) { float s = 0; // 求和變量. size_t i; size_t nBlockWidth = 8*4; // 塊寬. AVX寄存器能一次處理8個float，然后循環展開4次. size_t cntBlock = cntbuf / nBlockWidth; // 塊數. size_t cntRem = cntbuf % nBlockWidth; // 剩余數量. __m256 yfsSum = _mm256_setzero_ps(); // 求和變量。[AVX] 賦初值0 __m256 yfsSum1 = _mm256_setzero_ps(); __m256 yfsSum2 = _mm256_setzero_ps(); __m256 yfsSum3 = _mm256_setzero_ps(); __m256 yfsLoad; // 加載. __m256 yfsLoad1; __m256 yfsLoad2; __m256 yfsLoad3; const float* p = pbuf; // AVX批量處理時所用的指針. const float* q; // 將AVX變量上的多個數值合並時所用指針. // AVX批量處理. for(i=0; i<cntBlock; ++i) { yfsLoad = _mm256_load_ps(p); // [AVX] 加載. yfsLoad1 = _mm256_load_ps(p+8); yfsLoad2 = _mm256_load_ps(p+16); yfsLoad3 = _mm256_load_ps(p+24); yfsSum = _mm256_add_ps(yfsSum, yfsLoad); // [AVX] 單精浮點緊縮加法 yfsSum1 = _mm256_add_ps(yfsSum1, yfsLoad1); yfsSum2 = _mm256_add_ps(yfsSum2, yfsLoad2); yfsSum3 = _mm256_add_ps(yfsSum3, yfsLoad3); p += nBlockWidth; } // 合並. yfsSum = _mm256_add_ps(yfsSum, yfsSum1); // 兩兩合並(0~1). yfsSum2 = _mm256_add_ps(yfsSum2, yfsSum3); // 兩兩合並(2~3). yfsSum = _mm256_add_ps(yfsSum, yfsSum2); // 兩兩合並(0~3). q = (const float*)&yfsSum; s = q[0] + q[1] + q[2] + q[3] + q[4] + q[5] + q[6] + q[7]; // 處理剩下的. for(i=0; i<cntRem; ++i) { s += p[i]; } return s; }
2.5 測試框架
2.5.1 測試所用的數組
首先考慮一下測試所用的數組的長度應該是多少比較好。
為了避免內存帶寬問題，這個數組最好能放在L1 Data Cache中。現在的處理器的L1 Data Cache一般是32KB，為了保險最好再除以2，那么數組的長度應該是 32KB/(2*sizeof(float))=4096。
其次考慮內存對齊問題，avx要求32字節對齊。我們可以定義一個ATTR_ALIGN宏來統一處理變量的內存對齊問題。
該數組定義如下——
// 變量對齊. #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 4096 // = 32KB{L1 Cache} / (2 * sizeof(float)) ATTR_ALIGN(32) float buf[BUFSIZE];
2.5.2 測試函數
如果為每一個函數都編寫一套測試代碼，那不僅代碼量大，而且不易維護。
可以考慮利用函數指針來實現一套測試框架。
因sumfloat_base等函數的簽名是一致的，於是可以定義這樣的一種函數指針——
// 測試時的函數類型
typedef float (*TESTPROC)(const float* pbuf, size_t cntbuf);
然后再編寫一個對TESTPROC函數指針進行測試的函數——
// 進行測試 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 float 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); }
j是最內層的循環，負責多次調用TESTPROC函數指針。如果每調用一次TESTPROC函數指針后又調用clock函數，那會帶來較大的計時開銷，影響評測成績。
k循環負責檢測超時。當發現超過預定時限，便計算mps，即每秒鍾處理了多少百萬個單精度浮點數。然后存儲最佳的mps。
i是最外層循環的循環變量，循環3次然后報告最佳值。
2.5.3 進行測試
在進行測試之前，需要對buf數組進行初始化，將數組元素賦隨機值——
// init buf srand( (unsigned)time( NULL ) ); for (i = 0; i < BUFSIZE; i++) buf[i] = (float)(rand() & 0x3f); // 使用&0x3f是為了讓求和后的數值不會超過float類型的有效位數，便於觀察結果是否正確.
然后可以開始測試了——
// test runTest("sumfloat_base", sumfloat_base); // 單精度浮點數組求和_基本版. #ifdef INTRIN_SSE if (simd_sse_level(NULL) >= SIMD_SSE_1) { runTest("sumfloat_sse", sumfloat_sse); // 單精度浮點數組求和_SSE版. runTest("sumfloat_sse_4loop", sumfloat_sse_4loop); // 單精度浮點數組求和_SSE四路循環展開版. } #endif // #ifdef INTRIN_SSE #ifdef INTRIN_AVX if (simd_avx_level(NULL) >= SIMD_AVX_1) { runTest("sumfloat_avx", sumfloat_avx); // 單精度浮點數組求和_AVX版. runTest("sumfloat_avx_4loop", sumfloat_avx_4loop); // 單精度浮點數組求和_AVX四路循環展開版. } #endif // #ifdef INTRIN_AVX
INTRIN_SSE、INTRIN_AVX 宏是 zintrin.h 提供的，用於在編譯時檢測編譯器是否支持SSE、AVX指令集。
simd_sse_level、simd_avx_level函數是 ccpuid.h 提供的，用於在運行時檢測當前系統環境是否支持SSE、AVX指令集。
2.6 雜項
為了方便對比測試，可以在程序啟動時顯示程序版本、編譯器名稱、CPU型號信息。即在main函數中加上——
char szBuf[64]; int i; printf("simdsumfloat v1.00 (%dbit)\n", INTRIN_WORDSIZE); printf("Compiler: %s\n", COMPILER_NAME); cpu_getbrand(szBuf); printf("CPU:\t%s\n", szBuf); printf("\n");
INTRIN_WORDSIZE 宏是 zintrin.h 提供的，為當前機器的字長。
cpu_getbrand是 ccpuid.h 提供的，用於獲得CPU型號字符串。
COMPILER_NAME 是一個用來獲得編譯器名稱的宏，它的詳細定義是——
// 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++
三、全部代碼
3.1 simdsumfloat.c
全部代碼——
simdsumfloat.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++ ////////////////////////////////////////////////// // sumfloat: 單精度浮點數組求和的函數 ////////////////////////////////////////////////// // 單精度浮點數組求和_基本版. // // result: 返回數組求和結果. // pbuf: 數組的首地址. // cntbuf: 數組長度. float sumfloat_base(const float* pbuf, size_t cntbuf) { float s = 0; // 求和變量. size_t i; for(i=0; i<cntbuf; ++i) { s += pbuf[i]; } return s; } #ifdef INTRIN_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四路循環展開版. float sumfloat_sse_4loop(const float* pbuf, size_t cntbuf) { float s = 0; // 返回值. size_t i; size_t nBlockWidth = 4*4; // 塊寬. SSE寄存器能一次處理4個float，然后循環展開4次. size_t cntBlock = cntbuf / nBlockWidth; // 塊數. size_t cntRem = cntbuf % nBlockWidth; // 剩余數量. __m128 xfsSum = _mm_setzero_ps(); // 求和變量。[SSE] 賦初值0 __m128 xfsSum1 = _mm_setzero_ps(); __m128 xfsSum2 = _mm_setzero_ps(); __m128 xfsSum3 = _mm_setzero_ps(); __m128 xfsLoad; // 加載. __m128 xfsLoad1; __m128 xfsLoad2; __m128 xfsLoad3; const float* p = pbuf; // SSE批量處理時所用的指針. const float* q; // 將SSE變量上的多個數值合並時所用指針. // SSE批量處理. for(i=0; i<cntBlock; ++i) { xfsLoad = _mm_load_ps(p); // [SSE] 加載. xfsLoad1 = _mm_load_ps(p+4); xfsLoad2 = _mm_load_ps(p+8); xfsLoad3 = _mm_load_ps(p+12); xfsSum = _mm_add_ps(xfsSum, xfsLoad); // [SSE] 單精浮點緊縮加法 xfsSum1 = _mm_add_ps(xfsSum1, xfsLoad1); xfsSum2 = _mm_add_ps(xfsSum2, xfsLoad2); xfsSum3 = _mm_add_ps(xfsSum3, xfsLoad3); p += nBlockWidth; } // 合並. xfsSum = _mm_add_ps(xfsSum, xfsSum1); // 兩兩合並(0~1). xfsSum2 = _mm_add_ps(xfsSum2, xfsSum3); // 兩兩合並(2~3). xfsSum = _mm_add_ps(xfsSum, xfsSum2); // 兩兩合並(0~3). q = (const float*)&xfsSum; s = q[0] + q[1] + q[2] + q[3]; // 處理剩下的. for(i=0; i<cntRem; ++i) { s += p[i]; } return s; } #endif // #ifdef INTRIN_SSE #ifdef INTRIN_AVX // 單精度浮點數組求和_AVX版. float sumfloat_avx(const float* pbuf, size_t cntbuf) { float s = 0; // 求和變量. size_t i; size_t nBlockWidth = 8; // 塊寬. AVX寄存器能一次處理8個float. size_t cntBlock = cntbuf / nBlockWidth; // 塊數. size_t cntRem = cntbuf % nBlockWidth; // 剩余數量. __m256 yfsSum = _mm256_setzero_ps(); // 求和變量。[AVX] 賦初值0 __m256 yfsLoad; // 加載. const float* p = pbuf; // AVX批量處理時所用的指針. const float* q; // 將AVX變量上的多個數值合並時所用指針. // AVX批量處理. for(i=0; i<cntBlock; ++i) { yfsLoad = _mm256_load_ps(p); // [AVX] 加載 yfsSum = _mm256_add_ps(yfsSum, yfsLoad); // [AVX] 單精浮點緊縮加法 p += nBlockWidth; } // 合並. q = (const float*)&yfsSum; s = q[0] + q[1] + q[2] + q[3] + q[4] + q[5] + q[6] + q[7]; // 處理剩下的. for(i=0; i<cntRem; ++i) { s += p[i]; } return s; } // 單精度浮點數組求和_AVX四路循環展開版. float sumfloat_avx_4loop(const float* pbuf, size_t cntbuf) { float s = 0; // 求和變量. size_t i; size_t nBlockWidth = 8*4; // 塊寬. AVX寄存器能一次處理8個float，然后循環展開4次. size_t cntBlock = cntbuf / nBlockWidth; // 塊數. size_t cntRem = cntbuf % nBlockWidth; // 剩余數量. __m256 yfsSum = _mm256_setzero_ps(); // 求和變量。[AVX] 賦初值0 __m256 yfsSum1 = _mm256_setzero_ps(); __m256 yfsSum2 = _mm256_setzero_ps(); __m256 yfsSum3 = _mm256_setzero_ps(); __m256 yfsLoad; // 加載. __m256 yfsLoad1; __m256 yfsLoad2; __m256 yfsLoad3; const float* p = pbuf; // AVX批量處理時所用的指針. const float* q; // 將AVX變量上的多個數值合並時所用指針. // AVX批量處理. for(i=0; i<cntBlock; ++i) { yfsLoad = _mm256_load_ps(p); // [AVX] 加載. yfsLoad1 = _mm256_load_ps(p+8); yfsLoad2 = _mm256_load_ps(p+16); yfsLoad3 = _mm256_load_ps(p+24); yfsSum = _mm256_add_ps(yfsSum, yfsLoad); // [AVX] 單精浮點緊縮加法 yfsSum1 = _mm256_add_ps(yfsSum1, yfsLoad1); yfsSum2 = _mm256_add_ps(yfsSum2, yfsLoad2); yfsSum3 = _mm256_add_ps(yfsSum3, yfsLoad3); p += nBlockWidth; } // 合並. yfsSum = _mm256_add_ps(yfsSum, yfsSum1); // 兩兩合並(0~1). yfsSum2 = _mm256_add_ps(yfsSum2, yfsSum3); // 兩兩合並(2~3). yfsSum = _mm256_add_ps(yfsSum, yfsSum2); // 兩兩合並(0~3). q = (const float*)&yfsSum; s = q[0] + q[1] + q[2] + q[3] + q[4] + q[5] + q[6] + q[7]; // 處理剩下的. 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 4096 // = 32KB{L1 Cache} / (2 * sizeof(float)) ATTR_ALIGN(32) float buf[BUFSIZE]; // 測試時的函數類型 typedef float (*TESTPROC)(const float* 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 float 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("simdsumfloat 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] = (float)(rand() & 0x3f); // 使用&0x3f是為了讓求和后的數值不會超過float類型的有效位數，便於觀察結果是否正確. // test runTest("sumfloat_base", sumfloat_base); // 單精度浮點數組求和_基本版. #ifdef INTRIN_SSE if (simd_sse_level(NULL) >= SIMD_SSE_1) { runTest("sumfloat_sse", sumfloat_sse); // 單精度浮點數組求和_SSE版. runTest("sumfloat_sse_4loop", sumfloat_sse_4loop); // 單精度浮點數組求和_SSE四路循環展開版. } #endif // #ifdef INTRIN_SSE #ifdef INTRIN_AVX if (simd_avx_level(NULL) >= SIMD_AVX_1) { runTest("sumfloat_avx", sumfloat_avx); // 單精度浮點數組求和_AVX版. runTest("sumfloat_avx_4loop", sumfloat_avx_4loop); // 單精度浮點數組求和_AVX四路循環展開版. } #endif // #ifdef INTRIN_AVX return 0; }
3.2 makefile
全部代碼——
makefile
# flags CC = g++ CFS = -Wall -msse # 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 = simdsumfloat OBJS = simdsumfloat.o all : $(TARGETS) simdsumfloat : $(OBJS) $(CC) $(CFS) -o $@ $^ simdsumfloat.o : simdsumfloat.c zintrin.h ccpuid.h $(CC) $(CFS) -c $< clean : 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