用SIMD指令优化程序之抛砖引玉(ZT)

小谈 CPU 缓存体系

　　现在的 CPU 依旧采用冯诺伊曼体系，喜欢像傻子一样从头执行到尾，中途没有任何的跳转停顿等待。可是现实情况是，大部分程序里面还是少不了 IF ELSE 之类的判断，循环就更加得多了。如何优化循环大家可以自己琢磨，其实不难，可以参考一下《高质量 C\C++ 编程指南》

　　现在 CPU 上都有 Level 1 指令缓存（又叫做 L1 Trace ）与 Level 1 数据缓存（ L1 Data Cache ）。 PMMX ， P2 ， P3 为二者都准备了 16kb ，我的 P4 Northwood （以下简称 P4NW ）有 8kbL1 数据缓存和 12kb 指令缓存。 CPU 读取 L1 Data Cache 中的数据只需要 1 个时钟周期，速度非常快，应该是仅次于寄存器了。数据缓存是由 256 或者 512 行 32bytes 组成的，也就是 32bytes 对齐的，而 P4NW 是 64bytes 字节对齐的，并行 4 路，总共 128 行。当你处理的数据没有载入缓存的时候， CPU 将从内存读取缓存行大小的数据，所以缓存行总是对齐到能被 32 整除的物理地址。 CPU 对 L1 数据缓存中的数据进行操作是最快速的。所以推荐内存地址最起码是 32byte 对齐的。目前编译器在这个地方的优化已经非常好了，一般都是 4byte 对齐，当然也都是 32 对齐的。在后面你将会看到， SSE2 要求数据是 16 字节对齐的。

  　　缓存类似一个 C++ set 容器，但是不能赋值到一个任意的内存地址。每行本身都有 1 个 7bit 大小的关联值（ set value ）要和目标内存地址的 5 到 11 位对应（ 0-4 位已经忽略了），也可以理解为，关联值是内存段地址的一部分。 PPro 中，有 128 个关联值对应到 2 行，所以最多可以为任意的内存单元准备 2 个缓存行。 PMMX P2 P3 P4NW 有 4 个。由于内存是分段的，所以说 CPU 只能为， 5-11 位地址相同的内存准备 2 或者 4 个不同的缓存行。如何为两个内存地址赋予相同的关联值呢？把 2 个地址的低 5bit 去掉，这样就能被 32 整除了。如果这 2 个截断了的地址都是 4096 （ 1000H ）的倍数，那么这两个地址就有了相同的关联值。

  　　让我们用汇编加深一下印象，假设 ESI 中是 32 对齐的地址。

                                          AGAIN:  MOV  EAX,  [ESI]

MOV  EBX,  [ESI+13*4096+4]

MOV  ECX,  [ESI+20*4096+28]

DEC   EDX

JNZ   AGAIN

　　 Oh Year ，这里 3 个地址都有相同的关联值，而且地址跨度都超过了数据缓存的大小，可这个循环在 PPro 上效率会相当低。当你想读取 ECX 的值的时候，将没有空闲的缓存行了 —— 因为共享一个关联值，而且 2 行已经被使用了。此时 CPU 将腾出最近使用的 2 个缓存行，一个已经被 EAX 使用。然后 CPU 把这个缓存行用 [ESI+20*4096] 到 [ESI+20*4096+31] 的内存数据填充，然后从缓存中读取 ECX 。听起来好象相当的烦琐。更加糟糕的是，当又需要读取 EAX 的时候，还需要重复上述的过程，需要对内存缓存来回操作，效率相当的低，甚至不如不用缓存。可是，如果我们把第三行改成：

MOV  ECX,  [ESI+20*4096+32]

　　哦，不好，看起来，我们的地址超过了 32 ，不能被整除了。可是这样有了不同的关联值，也就意味着有了 1 个新行，不再共享可怜的 2 个行。这样一来，对三个寄存器的操作就不需要反复的用 2 个缓存行进行调度了，各有一个了。嘿嘿，这次只需要 3 个时钟周期了，而上一个要 60 个周期。这是在 PPro 上的，在后来的 CPU 中都是 4 路的，也就不存在上面的问题了。搞笑的是， Intel 的文档却错误的说 P2 的缓存是 2 路的。虽然说很少人在用那么古老的 CPU ，可是其中的道理大家应该明白。

　　可是判断要访问的部分数据是否有相同的关联值，也就是关于缓存是否能够命中的问题，是相当困难的，汇编还好，用高等级语言编译过的程序鬼知道是否对缓存做过优化呢。所以么，推荐，在程序的核心部分，对性能要求最高的部分，先对齐数据，然后确保使用的单个数据块不要超过缓存大小， 2 个数据块，单个不要超过缓存大小的一半（仔细想想为什么，因为关联值的问题，可以缓存分为两部分处理两块）。可是大部分情况下，我们都是使用远比数据缓存大的多的结构，以及编译器自己返回的指针，然后为了优化你可能希望把所有频繁使用的变量放到一个连续的数据块中以充分利用缓存。我们可以这样做，把静态变量数值拷贝到栈中的局部变量中，等子函数或者循环结束后再拷贝回来。这样一来就相当于把静态变量放入了连续的地址空间中去。

当读取的数据不在 L1 Cache 内时， CPU 将要从 L2 Cache 读取 L1 缓存行大小的数据到 L1 里去，大概需要 200ns 的时间（也就是 100Mhz 系统的 20 个时钟周期），但是直到你能够使用这些数据前，又需要有 50-100ns 的延迟。最糟糕的是，如果数据也不在 L2 Cache 中，那么就只能从最慢速的内存里读取了，内存的龟速哪能和全速的缓存相比。

好了，关于缓存的知识可以就此打住了，下面开始讲如何优化缓存。无非就是 3 种方法，硬件预取（ Prefetch ）、软件预取、使用缓存指令。关于预取的注意事项主要有这些：

1、  合理安排内存的数据，使用块结构，提高缓存命中率。

2、  使用编译器提供的预取指令。比如ICC中的_mm_prefetch _mm_stream，甚至_mm_load等比较“传统”的指令。

3、  尽可能少的使用全局的变量或者指针。

4、  程序尽可能少的进行判断跳转循环。

5、  使用const标记，不要在代码中混合register声明。

不过要提醒一句，真正提高程序效率的方法不是那种，从头到尾由于外科手术般的解剖，一个一个地方的优化，请抓住程序最核心的部分进行优化，记住 80-20 规则。

 

使用 SIMD

先复习一下对齐指令， __declspec(aliagn(#)) ， # 替换为字节数。比如想声明一个 16 字结对齐的浮点数组， __declspec(aliagn(16)) float Array[128] 。需要注意的是，最好充分了解你 CPU 的类型，支持哪些指令集。 SIMD 主要使用在需要同时操作大量数据的工作领域，比如 3D 图形处理（游戏），物理建模（ CAD ），加密，以及科学计算领域。据我所知，目前 GPGPU 也是使用 SIMD 的代表之一。

MMX

主要特性： 57 条指令， 64bit 的 FP 寄存器 MM0-MM7 ，对齐到 8 个 80bit 的 FP 寄存器 ST0-ST7 。需要数据 8 字节对齐，也就是使用 Packed 数字。

PS ：这里冒出了一个问题，为什么 Intel 要把 MMX 的寄存器和 FPU 的寄存器混合起来使用呢？因为这里牵涉到一个 FPU 状态切换问题，后面会提到，当你在一段代码中又要用到 MMX 指令又要用到传统的 FPU 指令，那么需要保存 FPU 状态，或者退出 MMX 。可是这种操作对于 FPU 来说非常昂贵，而且对于多任务操作系统来说，近乎于不可能完成的任务 —— 同时有许多程序，有些需要 MMX ，有些不需要，而正确地进行调度会变得非常困难。所以 Intel 将保存状态的工作完全交给了 CPU 自己，软件人员无须作太多这方面的工作，这样一来，就向前向后兼容了多任务操作系统，比如 Windows 和 Linux 。后来随着操作系统和 CPU 的不断升级，操作系统开发人员发布了一个补丁包，就可以让操作系统使用新的寄存器。这时人们都发现 Intel 的这种做法是相当短视的，这可以当作一个重大的失误。后来 Intel 通过引入了新的浮点指令集，这时才加入 XMM 寄存器。可造成这段故事的原因却根本不是技术问题，保证兼容性也是一个方面，总之真的说不清楚。你只要记得无法同时使用 MMX 与 FPU 就可以了， CPU 要进行模式切换。

SSE1

主要特性： 128bit 的 FP 寄存器 XMM0-XMM7 。增加了数据预取指令。额外的 64bit 整数支持。支持同时处理 4 个单精度浮点数，也就是 C\C++ 里的 float 。

适用范围：多媒体信号处理

SSE2

主要特性： 128bit 的 FP 寄存器支持处理同时处理 2 个双精度 double 浮点数，以及 16byte 8word 4dword 2quadword 整数。

适用范围： 3D 处理 语音识别 视频编码解码

SSE3

主要特性：增加支持非对称 asymmetric 和水平 horizontal 计算的 SIMD 指令。为 SIMD 提供了一条特殊的寄存器 load 指令。线程同步指令。

适用范围：科学计算 多线程程序

手头工具

1 、选择一个合适的编译器，推荐用 Intel C++ Compiler （以下简称 ICC ），以及 Visual Studio .NET 2003 及以上 IDE 附带的 C++ 编译器。同时， Microsoft C++ Compiler 也支持 AMD 的 3DNow 。 GCC C++ Compiler 没有测试。

2 、 Intel 以及 AMD 的汇编指令集手册。这个是必需的，强烈建议每个C++ Coder人手准备一份。

  所有的都用 C++ 混合变成的方式实现

使用范例：

向量乘法在 3D 处理中非常非常多，多半用于计算单位矢量的夹角。

我们先定义一个顶点结构。

__declspec(align( 16 ))  struct  Vertex{
     float  x,y,z,w;
};

    16字节对齐的结构，其实本身也是16字节的东西。如果没有对齐，运行时会报错。

w是其次坐标系的参数，处理向量的时候不需要用到。我的函数是这样的：

float Dot(Vertex* v1,Vertex* v2)
{
    Vertex tmp;
    __asm{
        MOV EAX,[v1];
        MOVAPS XMM0,[EAX];
        MOV EAX,[v2];
        MOVAPS XMM1,[EAX];
        MULPS XMM0,XMM1;
        MOVAPS tmp,XMM0;
    };
    return tmp.x + tmp.y + tmp.z;
};

    VC中反汇编之：

 1 float Dot(Vertex* v1,Vertex* v2)
 2 {
 3 0041C690  push        ebx  
 4 0041C691  mov         ebx,esp 
 5 0041C693  sub         esp,8 
 6 0041C696  and         esp,0FFFFFFF0h 
 7 0041C699  add         esp,4 
 8 0041C69C  push        ebp  
 9 0041C69D  mov         ebp,dword ptr [ebx+4] 
10 0041C6A0  mov         dword ptr [esp+4],ebp 
11 0041C6A4  mov         ebp,esp 
12 0041C6A6  sub         esp,0E8h 
13 0041C6AC  push        esi  
14 0041C6AD  push        edi  
15 0041C6AE  lea         edi,[ebp-0E8h] 
16 0041C6B4  mov         ecx,3Ah 
17 0041C6B9  mov         eax,0CCCCCCCCh 
18 0041C6BE  rep stos    dword ptr [edi] 
19     Vertex tmp;
20     __asm{
21         MOV EAX,[v1];
22 0041C6C0  mov         eax,dword ptr [v1] 
23         MOVAPS XMM0,[EAX];
24 0041C6C3  movaps      xmm0,xmmword ptr [eax] 
25         MOV EAX,[v2];
26 0041C6C6  mov         eax,dword ptr [v2] 
27         MOVAPS XMM1,[EAX];
28 0041C6C9  movaps      xmm1,xmmword ptr [eax] 
29         MULPS XMM0,XMM1;
30 0041C6CC  mulps       xmm0,xmm1 
31         MOVAPS tmp,XMM0;
32 0041C6CF  movaps      xmmword ptr [tmp],xmm0 
33     };
34     return tmp.x + tmp.y + tmp.z;
35 0041C6D3  fld         dword ptr [tmp] 
36 0041C6D6  fadd        dword ptr [ebp-1Ch] 
37 0041C6D9  fadd        dword ptr [ebp-18h] 
38 };

    前面都是保护现场入Stack的代码，没有必要管。我之所以这样，在Stack中声明了一个零时变量返回之，是为了减少代码的行数。有兴趣地可以参考本文后面引用资料中的Intel范例，代码多的多，功能却一样。这样就可以利用SIMD计算点乘了。图示：
    这种顶点格式称为AoS（Array of structure），这种结构的好处是，能够和现有的程序结构，比如D3D中的FVF顶点格式，和GL中的顶点格式。但是，由于许多情况下，并没有使用第四各浮点数，这就让SIMD指令浪费了25%的性能。于是有了SoA格式，让我们重新来过。
    我借用了一下上面一个结构的指令，还是没有用_mm_128格式，让大家看得清楚一些：

__declspec(align(16)) struct Vertex_soa{
     float x[4],y[4],z[4],w[4];
};

    依旧16字节对齐。计算函数如下：

 1 void Dot(Vertex_soa* v1,Vertex* v2,float* result)
 2 {
 3     Vertex tmp1,tmp2;
 4     __asm{
 5         MOV ECX,v1;
 6         MOV EDX,v2;
 7 
 8         MOVAPS XMM7,[ECX];
 9         MOVAPS XMM6,[ECX+16];
10         MOVAPS XMM5,[ECX+32];
11         MOVAPS XMM4,[ECX+48];
12         MOVAPS XMM0,XMM7;
13         UNPCKLPS XMM7,XMM6;
14         MOVLPS [EDX],XMM7;
15         MOVHPS [EDX+16],XMM7;
16         UNPCKHPS XMM0,XMM6;
17         MOVLPS [EDX+32],XMM0;
18         MOVHPS [EDX+48],XMM0;
19 
20         MOVAPS XMM0,XMM5;
21         UNPCKLPS XMM5,XMM4;
22         UNPCKHPS XMM0,XMM4;
23         MOVLPS [EDX+8],XMM5;
24         MOVHPS [EDX+24],XMM5;
25         MOVLPS [EDX+40],XMM0;
26         MOVHPS [EDX+56],XMM0;
27 
28         MOVAPS XMM3,[EDX];
29         MOVAPS XMM2,[EDX+16];
30         MOVAPS XMM1,[EDX+32];
31         MOVAPS XMM0,[EDX+48];
32 
33         MULPS XMM3,XMM2;
34         MULPS XMM1,XMM0;
35         MOVAPS tmp2,XMM1;
36         MOVAPS tmp1,XMM3;
37     };
38     result[0] = tmp1.x + tmp1.y + tmp1.z;
39     result[1] = tmp2.x + tmp2.y + tmp2.z;
40 };

    Oh Yeah，就是这样了，同时计算了1对乘法。我在代码中借用了一下前面的顶点结构，这样方便一些。至于SOA格式，请看前面的声明。很多代码都是转换Stack中的内存格式，转换成AOS格式，这样才能使用SIMD指令计算。

    通过上面的演示，想必大家已经对SIMD有了个直观地认识，其实在自己的代码中加入这些是非常方便与容易的。虽然说现在的CPU性能已经提高了许多，性能也强了许多，可是在诸多对性能要求高的地方，还是非常烤烟程序员的水平的。