GCC中SIMD指令的应用方法 at GCCFeli.cn

GCC中SIMD指令的应用方法

评论186次阅读2009.04.14 22:44; 作者：Felicia　

最近做图形学方面的编程，对SIMD指令比较感兴趣，于是转载了这篇文章。文章格式我稍微修正了一下。

X86的SIMD指令 – SIMD instrucitons in X86

IA-32 Intel体系结构的指令主要分为以下几类 [1]：

通用
x87 FPU
MMX技术
SSE/SSE2/SSE3扩展

MMX/SSE类扩展引入了SIMD（单指令多数据）的执行模式，可用于加速多媒体应用。下面简要介绍一下这些指令的执行环境和特征。

8个32位通用寄存器可为各个SIMD扩展所使用；
MMX：8个64位MMX寄存器（mm0 – mm7），也可为各SSE扩展所使用；
- 数据为整数，最多支持两个32位
- 运算中没有寄存器能够进行溢出指示
SSE：8个128位xmm寄存器，MXSCR寄存器，EFLAGS寄存器
- 支持单精度浮点
- MXSCR含有rounding, overflow标志
- 支持64位SIMD整数
SSE2：执行环境同sse
- 双精度浮点
- 128位整数
- 双—单精度转换
SSE3：与Inte Prescott处理器一同发布不久，共13条指令
- 主要增强了视频解码、3D图形优化和超线程性能

MMX技术出现最早，目前几乎所有的X86处理器都提供支持，包括嵌入式X86，所以下面的讨论主要基于MMX，但方法完全适用于SSEn，包括像AMD的3D Now等其它SIMD扩展。

MMX指令又分为以下几种：

数据传送：movd, movq
数据转换：packsswb, packssdw, packuswb, punpckhbw, punpckhwd, punpckhdq, punpcklbw, punpcklwd, punpckldq
并行算术：paddb, paddw, paddd, paddsb, paddsw, paddusb, paddusw, psubb, psubw, psubd, psubsb, psubsw, psubusb, psubusb, psubusw, pmulhw, pmullw, pmaddwd
并行比较：pcmpeqb, pcmpeqw, pcmpeqd, pcmpgtb, pcmpgtw, pcmpgtd
并行逻辑：pand, pandn, por, pxor
移位与旋转：psllw, pslld, psllq, psrlw, psrld, psrlq, psraw, psrad
状态管理：emms

这些指令除了需要注意功能外，还需要注意处理的数据类型。以上内容为背景介绍，细节请参考手册。

性能优化 – Performance Optimization

当使用C/C++完成了一个嵌入式应用的所有功能，性能问题常摆在面前，这时可以使用profile工具(如gprof)找出产生瓶颈的函数，将这些函数使用汇编彻底重写，例如MPEG-4编解码器xvid项目 [4]就使用了这种方法，而且针对不同处理器／指令集分别给出了不同的优化，正是如此该项目无论功能、还是性能均为一流，显然这是深度优化的目标所在。

在使用流水线、VLIW以及SIMD的体系结构（比如某些DSP）上，整个函数的手工优化可以带来几倍到几十倍的性能提升。不过，性能允许，对于函数内关键部分使用一些特定的实现，既突出重点提高性能，又可以尽多地利用C/C++的高级特征，相对缩短开发周期。下面给出使用GCC时，应用MMX指令的几种混合编程方法：

Intel C/C++ 编译器intrinsics
GCC builtin操作
嵌入汇编asm construct

Intel C/C++ 编译器intrinsics – Intel C/C++ Compiler Intrinsics

查看IA-32 Intel指令集手册[2]时，部分指令的解释中会有一项“Intel C/C++ Compiler Intrinsic Equivalent”，会指出该指令对等的intrinsic。 intrinsic在C/C++程序中的语法是以函数形式出现，编译时可以直接翻译为一条MMX指令（复合情况会生成最直接的几条），换言之，如果不使用intrinsic，可能需要多条C/C++语句完成，而编译器却并不能保证将这几条语句能够生成这条最高效的MMX指令。并不是每条MMX指令都有对等的intrinsic，手册的附录中列出了所有的，它们分为简单型（simple）和复合型（composite）两种，每个简单型的就是对应一条指令，而复合型则对应多条指令。

GCC支持Intel C/C++ Compiler Intrinsics。用法如下示例：

#include <stdio.h>
#include <xmmintrin.h> /*一定需要包括此头文件*/
/*gcc -Wall -march=pentium4 -mmmx -o ins mmx_ins.c*/
int main(int argc,char *argv[]) {
   /*使用MMX做以下向量的点积*/
   short in1[] = {1, 2, 3, 4};
   short in2[] = {2, 3, 4, 5};
   int out1;
   int out2;
   __m64 m1; /* MMX支持64位整数的mm寄存器 */
   __m64 m2; /* MMX操作需要使用mm寄存器 */
   __m128 m128; /* for SSEn only*/
   /*每次往mm寄存器装入两个short型的数，注意是两个*/
   m1 = _mm_cvtsi32_si64(((int*)in1)[0]);
   m2 = _mm_cvtsi32_si64(((int*)in2)[0]);
   /*一条指令进行4个16位整数的乘加*/
   /*生成两个32位整数*/
   m2 = _mm_madd_pi16(m1, m2);
   /*将低32位整数放入通用寄存器*/
   out1 = _mm_cvtsi64_si32(m2);
   /*将高32位整数右移后，放入通用寄存器*/
   m2 = _mm_slli_pi32(m2, 32);
   out2 = _mm_cvtsi64_si32(m2);
   /*清除MMX状态*/
   _mm_empty();
   /*将两个32位数相加，结果为8*/
   out1 += out2;
   printf("a: %d\n", out1);
   return 0;
}

几点说明：

即使你不是P4平台，编译时也请使用以下选项，
gcc -Wall -march=pentium4 -mmmx -o ins mmx_ins.c
否则，会出现如下类似信息：
...xmmintrin.h:34:3: #error "SSE instruction set not enabled"
最终结果实际并没有求得四对乘积的和，只是前两对的， instrinsic _mm_cvtsi32_si64只向mm寄存器放入了低32位，高32位为零，但mmx有指令movq可以做到64位的数据传送，intrinsic没有对应，这也说明并不是所有的指令有等价的intrinsic。
当计算的向量为两对0×8000, 0×8000时，即(-2^15)*(-2^15) + (-2^15)*(-2^15) ，结果应该为 2^31，但计算出来的值是
-2^31，因为发生了溢出，可程序无从知道。这是使用MMX时，应特别注意的，计算溢出没有任何标志位指示，一个极大的值变为极小，SSE对此做了改善。
程序不再使用MMX之时，注意使用emms指令清除MMX状态。

使用built-in操作 – GCC built-in Operation

什么是built-in操作？就是对待MMX操作数，就如int, float等基本数据类型一般，有相应定义的操作，如加(+)、减(-)，或者数据类型之间的转换。详细内容参考GNU GCC Manual
[5] Extensions to the C Language Family4#4Built-in Functions4#4 X86 Built-in Functions一节。

一些MMX指令有其相应的built-in操作，下面一段代码为例：

#include <stdio.h>
/*无需特别的头文件，built-in嘛*/
/* gcc -Wall -o bins builtinmmx.c*/
/*定义了一个vector数据类型，hi表示16位，4表示4个*/
/*typedef int v4hi __attribute__ ((mode(V4HI)));*/
/*新版的gcc认为这么定义更好，vector_size(8)表示8byte长度的vector，short表示按照short方式存储*/
typedef short v4hi __attribute__ ((vector_size(8)));
/*定义了2个32位的vector类型，si表示32位*/
typedef int v2si __attribute__ ((mode(V2SI)));

int main(int argc,char *argv[]) {
   short pa[4] = {0x8000, 0x8000, 1, -1};
   short pb[4] = {0x8000, 0x7FFF, -1, -2};

   v4hi va, vb;
   v4hi vsum;

   va = ((v4hi*)pa)[0];
   vb = ((v4hi*)pb)[0];

   /* 4个16位进行饱和加 */
   //vsum = __builtin_ia32_paddsw(va, vb);
   /* 4个16位还可以直接进行加法，但不同于两个long long相加 */
   vsum = va + vb;

   /*vector的输出还需要强制转换为long long*/
   printf("...with MMX instructions...to compute vec_add: %llx \n", (long long)vsum);

   //结果1：0xfffd0000ffff8000
   //结果2：0xfffd0000ffff0000

   return 0;
}

几点说明：

是的，这里built-in vector及其操作，随着GCC的发展正在加强。如果需要使用以上范例，应使用GCC 3.4以上版本；
使用builtin函数时，与intrinsic相似；但本质却是不同，这里两个向量使用‘+’操作就说明了vector也如其它数据类型一样，编译器直接支持，只不过这里的加法就是指四个单元数分别相加，低位单元的进位不会影响相邻高位单元的数据；
vector还可以强制转换为通用数据。

嵌入汇编 – Inline asm

GCC一开始就允许C代码中嵌入asm指令，并不只是针对MMX指令，不过对于MMX技术，显然也是一个很好的利用方法，详细的语法请参考GNU GCC手册[5]，或者GCC: The Complete Reference[6]”Inline Assembly”一节。
如下是一个点积的例子：

#include <stdio.h>
/** GCC -o ins inlinemmx.c **/
int main(int argc,char *argv[]) {
   int i;
   int result;
   short a[] = {1, 2, 3, 4, 5, 6, 7, 8};
   short b[] = {1, 1, 1, 1, 1, 1, 1, 1};
   printf("...with MMX instructions...\n");

   /*首先，将点积合累积寄存器清零，实际缺省就为0？*/
   asm("pandn %%mm5,%%mm5;"::);
   /*读入a, b，每四对数相乘后分两组相加，形成两组和*/
   /*这里的循环控制是C在做*/
   for (i = 0; i < sizeof(a)/sizeof(short); i += 4) {
   asm("movq %0,%%mm0;\
   movq %1,%%mm1;\
   pmaddwd %%mm1,%%mm0;\
   paddd %%mm0,%%mm5; #相乘后相加 "
   :
   : "m" (a[i]), "m" (b[i]));
   }
   /*将两组和分离，并相加*/
   asm("movq %%mm5, %%mm0;\
   psrlq $32,%%mm5;\
   paddd %%mm0, %%mm5;\
   movd %%mm5,%0;\
   emms"
   :"=r" (result)
   :);
   printf("result: 0x%x\n", result);
   //这里结果为0x24
   return 0;
}

几点说明：

这里是典型的在函数中C和汇编混合编程；
注意汇编指令中操作数的顺序；
这里可以直接使用movq等没有intrinsics/built-in对应的指令；
注意在asm指令序列中间不要加杂注释，可能导致生成的代码不正确。

MMX实用一例：合成滤波器 – Synthesis Filter in X86 SIMD INSTRUCTIONS

下面是合成滤波器(Synthesis Filter)的一个优化过程，合成滤波器在语音编解码中有广泛应用，运行时也占用了整个算法中较高比例的时间。

for (i = 0; i < lg; i++) {
   s = L_mult(x[i], a[0]); /*L_mult是相乘后左移*/
   for (j = 1; j <= M; j++) { /*M这里固定为10*/
   s = L_msu(s, a[j], yy[-j]); /*L_msu是乘减后左移操作*/
   }

   s = L_shl(s, 3); /*左移三位*/
   *yy++ = g729round(s);
}

上面的代码，因为内存循环为10，可以考虑展开，并统一操作为乘加指令。

/*为了使用乘加操作，需要调整10个系数的顺序*/
for (i = 0; i < M; i++)
   ta[i] = -a[M - i];

ta[11] = 0;
ta[10] = a[0];

for (i = 0; i < lg; i++) {
   *yy = x[i];
   yy[1] = 0;
   s = L_mac(s, ta[11], yy[1]);
   s = L_mac(s, ta[10], yy[0]);
   s = L_mac(s, ta[9], yy[-1]);
   s = L_mac(s, ta[8], yy[-2]);
   s = L_mac(s, ta[7], yy[-3]);
   s = L_mac(s, ta[6], yy[-4]);
   s = L_mac(s, ta[5], yy[-5]);
   s = L_mac(s, ta[4], yy[-6]);
   s = L_mac(s, ta[3], yy[-7]);
   s = L_mac(s, ta[2], yy[-8]);
   s = L_mac(s, ta[1], yy[-9]);
   s = L_mac(s, ta[0], yy[-10]);

   s = L_shl(s, 3);
   *yy++ = g729round(s);
}

以上循环内核正好可以将MMX的8个寄存器全部利用。

/*为了使用乘加操作，需要调整10个系数的顺序*/
for (i = 0; i < M; i++)
   ta[i] = -a[M - i];

ta[11] = 0;
ta[10] = a[0];

/*11个系数分别放入3个MMX寄存器，0作填充*/
asm("movq %0,%%mm0;\
   movq %1,%%mm1;\
   movq %2,%%mm2"
   :
   : "m" (ta[0]), "m" (ta[4]), "m"(ta[8]));

/*利用MMX技术进行滤波器核心操作*/
for (i = 0; i < lg; i++) {
   *yy = x[i];
   yy[1] = 0;
   asm("pandn %%mm6,%%mm6;\
   movq %1,%%mm3;\
   movq %2,%%mm4;\
   movq %3,%%mm5;\
   pmaddwd %%mm0,%%mm3;\
   pmaddwd %%mm1,%%mm4;\
   pmaddwd %%mm2,%%mm5;\
   paddd %%mm3, %%mm6;\
   paddd %%mm4, %%mm6;\
   paddd %%mm5, %%mm6;\
   movq %%mm6, %%mm7;\
   psrlq $32, %%mm6;\
   paddd %%mm7, %%mm6;\
   movd %%mm6,%0;\
   emms"
   :
   :"r"(s), "m" (yy[-10]), "m" (yy[-6]), "m"(yy[-2]));
   /*因为指令结果饱和属性的限制，s还没有左移，所以下面多做一位饱和左移*/
   s = L_shl(s, 4);
   *yy++ = g729round(s);
}

几点说明：

注意：以上嵌入的汇编代码输出结果s放在了输入处，属于实践中的个案；
MMX没有乘左移之类的DSP指令，甚至还没有加饱和之类的操作，SSE中有一定增强；
以上操作，理论上存在溢出可能，所以最后使用原有的饱和左移操作，减少了一定风险；
上面的部分代码操作显然允许并行，这在VLIW系统中十分有用；
这已经形成了该滤波器全面优化的核心。

总结 – Conclusion

如果愿意尽多地利用SIMD技术，可能需要更多地使用汇编级的编码，不过也有一些高级语言和汇编的混合编程技术能够帮助你，它们有的提高性能更大一些，有的形式上更优雅些，本质上效率也不错，都不失好的方法，建议尝试。

正是如此，一方面CPU上支持越来越多的SIMD指令集扩展，另一方面GCC也正在加紧支持这些扩展的易用，对，正在，碰到一些问题，先想办法绕过去，这里使用GCC 3.4.1，根据经验效果还是不错的。

关于文档