SSE指令集入门

原文地址：http://blog.csdn.net/wensishuai/article/details/8439418

Intel公司的单指令多数据流式扩展（SSE，Streaming SIMD Extensions）技术能够有效增强CPU浮点运算的能力。Visual Studio .NET 2003提供了对SSE指令集的编程支持，从而允许用户在C++代码中不用编写汇编代码就可直接使用SSE指令的功能。MSDN中有关SSE技术的主题[1]有可能会使不熟悉使用SSE汇编指令编程的初学者感到困惑，但是在阅读MSDN有关文档的同时，参考一下Intel软件说明书（Intel Software manuals）[2]会使你更清楚地理解使用SSE指令编程的要点。

SIMD（single-instruction, multiple-data）是一种使用单道指令处理多道数据流的CPU执行模式，即在一个CPU指令执行周期内用一道指令完成处理多个数据的操作。考虑一下下面这个任务：计算一个很长的浮点型数组中每一个元素的平方根。实现这个任务的算法可以这样写：

for each  f in array        //对数组中的每一个元素
    f = sqrt(f)             //计算它的平方根

为了了解实现的细节，我们把上面的代码这样写：

for each  f in array
{
    把f从内存加载到浮点寄存器
    计算平方根
    再把计算结果从寄存器中取出放入内存
}

具有Intel SSE指令集支持的处理器有8个128位的寄存器，每一个寄存器可以存放4个（32位）单精度的浮点数。SSE同时提供了一个指令集，其中的指令可以允许把浮点数加载到这些128位的寄存器之中，这些数就可以在这些寄存器中进行算术逻辑运算，然后把结果放回内存。采用SSE技术后，算法可以写成下面的样子：

for each  4 members in array  //对数组中的每4个元素
{
    把数组中的这4个数加载到一个128位的SSE寄存器中
    在一个CPU指令执行周期中完成计算这4个数的平方根的操作
    把所得的4个结果取出写入内存
}

下面是一个演示的例子

使用纯C＋＋

[cpp]  view plain copy

1. void CSSETestDlg::ComputeArrayCPlusPlus(  
2.           float* pArray1,                   // [in] first source array  
3.           float* pArray2,                   // [in] second source array  
4.           float* pResult,                   // [out] result array  
5.           int nSize)                        // [in] size of all arrays  
6. {  
7.   
8.     int i;  
9.   
10.     float* pSource1 = pArray1;  
11.     float* pSource2 = pArray2;  
12.     float* pDest = pResult;  
13.   
14.     for ( i = 0; i < nSize; i++ )  
15.     {  
16.         *pDest = (float)sqrt((*pSource1) * (*pSource1) + (*pSource2)  
17.                  * (*pSource2)) + 0.5f;  
18.   
19.         pSource1++;  
20.         pSource2++;  
21.         pDest++;  
22.     }  
23. }  

使用SSE内嵌原语

[cpp]  view plain copy

1. void CSSETestDlg::ComputeArrayCPlusPlusSSE(  
2.           float* pArray1,                   // [in] first source array  
3.           float* pArray2,                   // [in] second source array  
4.           float* pResult,                   // [out] result array  
5.           int nSize)                        // [in] size of all arrays  
6. {  
7.     int nLoop = nSize/ 4;  
8.   
9.     __m128 m1, m2, m3, m4;  
10.   
11.     __m128* pSrc1 = (__m128*) pArray1;  
12.     __m128* pSrc2 = (__m128*) pArray2;  
13.     __m128* pDest = (__m128*) pResult;  
14.   
15.   
16.     __m128 m0_5 = _mm_set_ps1(0.5f);        // m0_5[0, 1, 2, 3] = 0.5  
17.   
18.     for ( int i = 0; i < nLoop; i++ )  
19.     {  
20.         m1 = _mm_mul_ps(*pSrc1, *pSrc1);        // m1 = *pSrc1 * *pSrc1  
21.         m2 = _mm_mul_ps(*pSrc2, *pSrc2);        // m2 = *pSrc2 * *pSrc2  
22.         m3 = _mm_add_ps(m1, m2);                // m3 = m1 + m2  
23.         m4 = _mm_sqrt_ps(m3);                   // m4 = sqrt(m3)  
24.         *pDest = _mm_add_ps(m4, m0_5);          // *pDest = m4 + 0.5  
25.           
26.         pSrc1++;  
27.         pSrc2++;  
28.         pDest++;  
29.     }  
30. }  

使用SSE汇编

[cpp]  view plain copy

1. void CSSETestDlg::ComputeArrayAssemblySSE(  
2.           float* pArray1,                   // [输入] 源数组1  
3.           float* pArray2,                   // [输入] 源数组2  
4.           float* pResult,                   // [输出] 用来存放结果的数组  
5.           int nSize)                        // [输入] 数组的大小  
6. {  
7.     int nLoop = nSize/4;  
8.     float f = 0.5f;  
9.   
10.     _asm  
11.     {  
12.         movss   xmm2, f                         // xmm2[0] = 0.5  
13.         shufps  xmm2, xmm2, 0                   // xmm2[1, 2, 3] = xmm2[0]  
14.   
15.         mov         esi, pArray1                // 输入的源数组1的地址送往esi  
16.         mov         edx, pArray2                // 输入的源数组2的地址送往edx  
17.   
18.         mov         edi, pResult                // 输出结果数组的地址保存在edi  
19.         mov         ecx, nLoop                  //循环次数送往ecx  
20.   
21. start_loop:  
22.         movaps      xmm0, [esi]                 // xmm0 = [esi]  
23.         mulps       xmm0, xmm0                  // xmm0 = xmm0 * xmm0  
24.   
25.         movaps      xmm1, [edx]                 // xmm1 = [edx]  
26.         mulps       xmm1, xmm1                  // xmm1 = xmm1 * xmm1  
27.   
28.         addps       xmm0, xmm1                  // xmm0 = xmm0 + xmm1  
29.         sqrtps      xmm0, xmm0                  // xmm0 = sqrt(xmm0)  
30.   
31.         addps       xmm0, xmm2                  // xmm0 = xmm1 + xmm2  
32.   
33.         movaps      [edi], xmm0                 // [edi] = xmm0  
34.   
35.         add         esi, 16                     // esi += 16  
36.         add         edx, 16                     // edx += 16  
37.         add         edi, 16                     // edi += 16  
38.   
39.         dec         ecx                         // ecx--  
40.         jnz         start_loop                //如果不为0则转向start_loop  
41.     }  
42. }  

在信号处理中的实际应用（sse2）：

获得信号能量

[cpp]  view plain copy

1. /* 
2. * Compute Energy of a complex signal vector, removing the DC component!  
3. * input  : points to vector 
4. * length : length of vector in complex samples 
5. */  
6.   
7. #define shift 4  
8. #define shift_DC 0  
9.   
10. int signal_energy(int *input, unsigned int length)  
11. {  
12.     int i;  
13.     int temp, temp2;  
14.     register __m64 mm0, mm1, mm2, mm3;  
15.     __m64 *in;  
16.   
17.     in = (__m64 *)input;  
18.   
19.     mm0 = _m_pxor(mm0,mm0);  
20.     mm3 = _m_pxor(mm3,mm3);  
21.   
22.     for (i = 0; i < length >> 1; i++) {  
23.         mm1 = in[i];  
24.         mm2 = mm1;  
25.         mm1 = _m_pmaddwd(mm1, mm1);  
26.         mm1 = _m_psradi(mm1, shift);  
27.         mm0 = _m_paddd(mm0, mm1);  
28.         mm2 = _m_psrawi(mm2, shift_DC);  
29.         mm3 = _m_paddw(mm3, mm2);  
30.     }  
31.   
32.     mm1 = mm0;  
33.     mm0 = _m_psrlqi(mm0, 32);  
34.     mm0 = _m_paddd(mm0, mm1);  
35.     temp = _m_to_int(mm0);  
36.     temp /= length;  
37.     temp <<= shift;   
38.   
39.     /*now remove the DC component*/  
40.     mm2 = _m_psrlqi(mm3, 32);  
41.     mm2 = _m_paddw(mm2, mm3);  
42.     mm2 = _m_pmaddwd(mm2, mm2);  
43.     temp2 = _m_to_int(mm2);  
44.     temp2 /= (length * length);  
45.     temp2 <<= (2 * shift_DC);  
46.     temp -= temp2;  
47.     _mm_empty();  
48.     _m_empty();  
49.   
50.     return((temp > 0) ? temp : 1);  
51. }