SIMD 编程的优势 TickerTape Part 2

简介

Ticker Tape 是一种技术演示，旨在鼓励开发人员在粒子系统中执行更为复杂的操作。参与该演示的开发人员会运用大量技术，来提高包括多线程和针对英特尔® SIMD 流指令扩展（SSE）的优化等在内的性能。请访问：http://software.intel.com/zh-cn/articles/tickertape/ ，查看本文概述，下载本演示。本文将重点谈论通过在 Ticker Tape 演示中引入 SSE 指令所获得的性能提升。但在此之前，我们将首先介绍 SSE 编程，指导您规划您的数据，使之能为 SSE 带来最大优势。最后，我们将示范如何采用 SSE 计算点积.

背景

SSE 是一套专门为 SIMD（单指令多数据）架构设计的指令集。通过它，用户可以同时在多个数据片段上执行运算，实现数据并行（有时又称矢量处理）。例如，我们可以利用这套指令集使两个数组各自相乘：

float a[nElements], b[nElements], c[nElements];

for (unsigned int i = 0; i < nElements; i++) {
    c[i] = a[i] * b[i];
}

列表 1：两个数组相乘的标量法，其中每次迭代处理一个元素

一般而言，如列表 1 所示，存在一个让所有元素进行迭代的循环，在这个循环中每个元素会相互相乘，然后保存乘积。现在，我们除了可以在每次循环迭代时执行一个乘法运算外，还可以执行多个乘法运算。下面是可以执行多个乘法运算的函数 MultiplyFourElements()。

float a[nElements], b[nElements], c[nElements];

void MultiplyFourElements(float *a, float *b, float *c) {
    c[0] = a[0] * b[0];
    c[1] = a[1] * b[1];
    c[2] = a[2] * b[2];
    c[3] = a[3] * b[3];
}

for (unsigned int i = 0; i < nElements; i += 4) {
    MultiplyFourElements(&a[i], &b[i], &c[i]);
}

列表 2：两个数组相乘的标量法，其中每次迭代处理四个元素

如列表 2 所示，我们创建了一个可以同时处理四个元素的函数。利用该函数，我们执行循环迭代的次数减少了四倍。尽管如此，由于每次迭代所需执行的数学运算相同，我们在效率上并未获得较大提升。然而有了 SSE 指令，我们不再需要通过一个函数连续执行四次乘法运算，只需借助一条指令即可同时执行四个乘法运算。SSE 会使用处理器上的 128 位宽专用寄存器。这些寄存器可以保存任何 128 位数据，如两个双精度数字、四个单精度数字或 16 字节数字等。采用 SSE 进行编程的方式有两种：一种是直接编写汇编指令代码, 另一种是使用 intrinsics 函数编程。Ticker Tape 演示以及本文都将只侧重于使用 intrinsic 函数进行编程。使用 intrinsics 而非汇编指令编程是一种更加直接的方法，与标准 C/C++ 编程类似。此外，使用 intrinsics 编程还有助于编译器更好地优化代码，对此本文将稍后阐释。

// Assembly SSE Instruction
/// xmm0 and xmm1 are actual registers, not variables
mulps xmm0,xmm1

// Intrinsic SSE Instruction
__m128 a, b, c;
c = _mm_mul_ps(a, b);

列表 3：汇编指令与 Intrinsic SSE 指令 
type __m128 是针对一个可映射至其中一个 SIMD 寄存器的 16 字节对齐变量的定义。该程序首先需要将已完成 _mm_load_ps() intrinsic（未在列表 3 中显示）运算的数据明确加载到 SIMD 寄存器中。_mm_mul_ps() 是实际执行运算的指令。一旦我们获得运算结果，我们可以利用 _mm_store_ps()intrinsic 将其作为输出数组保存下来。

#include <xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h>

float a[nElements], b[nElements], c[nElements];
__m128 A, B, C;

//... load data in a, b

for (unsigned int i = 0; i < nElements; i += 4) {
    A = _mm_load_ps(&a[i]);
    B = _mm_load_ps(&b[i]);
    C = _mm_mul_ps(A, B);
    _mm_store_ps(&c[i], C);
}
<xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h><xmmintrin.h></xmmintrin.h>

列表 4：两个数组相乘的 SIMD 方法 
显然，运算差异已在图 1 中显示：

 
图 1：对比标量循环运算法与 SIMD 运算法

数据布局

许多应用的瓶颈并非在于算法的运算部分，而在于数据读写上。在上一个示例中，75% 的指令被用于加载和保存数据。一旦您想访问的数据保存在了不同的区域中，则需花费大量时间才能将该数据加载至高速缓存中并继而加载至 SIMD 寄存器中。例如，如果上一示例中的数组实际上是我们所拥有的数组类的成员变量，则将该数据加载至寄存器中会是一件十分困难且耗时的工作。各个元素将不得不被逐个加载至 SIMD 寄存器中。

class foo {
    float a;
    float b;
    ... other data ...
};

void bar() {
    foo fooArray[nElements];
    float c[nElements];
    __m128 A, B, C;

    // Non_SIMD Method
    for (unsigned int i = 0; i < nElements; i++) {
        c[i] = fooArray[i].a * fooArray[i].b;
    }

    // SIMD Method
    for (unsigned int i = 0; i < nElements; i += 4) {
        A = _mm_load_ps(&fooArray[i].a); // What will this do?
        B = _mm_load_ps(&fooArray[i].b); // Load incorrect data
        C = _mm_mul_ps(A, B);
        _mm_store_ps(&c[i], C);
    }
}

列表 5：支持 SIMD 不适用数据布局（结构数组）的类 
在列表 5 中，数据加载方式出现了错误。该程序打算将连续内存加载至 SIMD 寄存器中，但这将会形成错误数据。要将数据加载至寄存器中不是不可能，但这样做涉及到重新安排数据的格式，而这需要一定的步骤才能完成。无论重新安排数据实际生成的处理成本如何，需要转移的数据很多，加载至内存的高速缓存行实际上也很多。然而，尽管如此，使用 SSE 指令仍然存在诸多优势，尤其是当相同的数据正在执行多项运算时。重新安排数据布局是指按照顺序安排类似数据片断，便于高效加载和保存这些数据。除此以外，安排数据布局的另一个优势是在 SIMD 寄存器尺寸增加时，改写代码会更轻松。这样一来，您便可以一次加载八个而非四个浮点数据。在上一示例的基础上采用更高效的数据布局进行改编后的程序如下：

class foo {
    float *a;
    float *b;
    ... other data ...
};

foo::foo(unsigned int nElements) {
    a = (float *)_mm_malloc(nElements * sizeof(float), 16);
    b = (float *)_mm_malloc(nElements * sizeof(float), 16);
}

void bar() {
    foo fooVariable(nElements);
    float c[nElements];
    __m128 A, B, C;

    // Non_SIMD Method
    for (unsigned int i = 0; i < nElements; i++) {
        c[i] = fooVariable.a[i] * fooVariable.b[i];
    }

    // SIMD Method
    for (unsigned int i = 0; i < nElements; i += 4) {
        A = _mm_load_ps(&fooVariable.a[i]);
        B = _mm_load_ps(&fooVariable.b[i]);
        C = _mm_mul_ps(A, B);
        _mm_store_ps(&c[i], C);
    }
}

列表 6：支持 SIMD 适用数据布局（结构数组）的类

 
图 2：SoA 与 AoS 内存布局

在 Ticker Tape 中实施 SIMD 优化

当我们通过优化 Ticker Tape 以发挥由两个重要部分组成的 SIMD 优势、重新安排数据使之更适用于 SIMD，以及重新编写部分代码以使用 SSE 指令时，我们对 Ticker Tape 演示进行了改造，最后发现我们的大部分时间花费在了 Newton() 函数上。这个函数主要计算各种力如何影响粒子。在原始架构中，每个粒子都会调用一次该函数，并会在每个角执行四次内部运算。在较高层面上进行的概述如下：

class RigidBody
{
    void Newton();
    // ... Bunch more functions ...

    D3DXVECTOR3 Position;
    D3DXVECTOR3 Rotation;

    // ... Lots of other data ...
};

void RigidBody::Newton()
{
    for (unsigned int = 0; i < 4; i++)
    {
        // ... some math goes on ...

        // Velocity_Ground
        D3DXVECTOR3 Vel_Ang_CM_global;
        D3DXVec3TransformCoord(&Vel_Ang_CM_global,
            &Velocity_Angular_CM, &this->LocalGlobal);
        D3DXVec3Cross(&tmp, &Radial_Vec, &(Vel_Ang_CM_global));
        Velocity_Ground = (this->Velocity_Linear_CM) + tmp;

        // ... more math ...
    }

    return;
}

列表 7：原始 Ticker Tape 布局 
每个粒子都代表着不同的对象并包含着各自的操作方法。此外，程序中还存在一个通过每个调用其成员方法的粒子进行迭代的循环。在向 SSE 进行迁移时首先应创建另外一个名为 NewtonArray() 的 Newton() 函数。这与原始布局并无二致，但却是在整个粒子数组而非一个粒子上创建的新函数。同时，还应创建包含整个数组的 RigidBody 类，以便程序从拥有 RigidBody 对象的数组迁移至一个拥有数组的 RigidBody。根据针对不同数据布局的测试，数据进一步分解成了各自的组成部分（如列表 8 所示）。

// Original
D3DXVECTOR3 rotation;

// Half way
D3DXVECTOR3 *rotation;

// Final
float *rotationX;
float *rotationY;
float *rotationZ;

列表 8：改变数据布局

待数据重组后，NewtonArray() 已经被修改，可以使用 SSE 指令。实施这一修改的基本前提是存在嵌套循环结构。每个采用内部循环索引的变量一般都会被加载到 SIMD 寄存器中，但每个采用外部循环索引的变量却都会被重复加载到寄存器中。事实上，我们删除了内部循环（如列表 9 所示）。_mm_set1_ps() 指令与 _mm_load_ps() 指令类似，唯一的不同之处在于前者可加载一个浮点值并能将其复制到寄存器的四个分区中。

// Original code example
for (unsigned int i = 0; i < nParticles; i++) {
    for (unsigned int j = 0; j < 4; j++) {
        c = a[i] * b[j]; // Notice the different indexers
    }
}

// SSE version, inner loop removed
for (unsigned int i = 0; i < nParticles; i++) {
    A = _mm_set1_ps(&a[i]);  // copy the value a[i] into all 4 slots
    B = _mm_load_ps(&b[i * 4]);
    C = _mm_mul_ps(A, B);
    _mm_store_ps(&c[i], C);
}

列表 9：删除内部循环

即使只修改该函数的一小部分代码，实施上述优化的优势仍然十分明显。在完成所有优化操作后，我们计算了 NewtonArray() 函数的执行时间。通过使用支持 Visual Studio 2008 的微软编译器（MSVCC），函数执行速度提高了 1.8 倍。由于具备自动矢量化性能，采用英特尔® C++ 编译器（ICC）编写原始代码能够带来巨大优势。采用 ICC 编写手写 SSE 可使函数执行速度提高达 4.5 倍。

Compiler:	MSVCC	MSVCC + SSE	ICC	ICC + SSE
Time (ms):	16.3	8.9	9.9	3.6
Improvement:	1.0x	1.8x	1.6x	4.5x

表 1：函数 Newton 的执行时间（以毫秒为单位）

点积示例

本章节将演示如何使用 SSE 指令计算点积（实际上是四个点积），并简要介绍具体算法。所涉及的变量全部被命名为 xmm#，因为八个 SSE 寄存器的名称从 xmm0 到 xmm7 不等。尽管如此，您没有必要采用这种方式命名变量，也不用只保留八个变量。

// Dot Product
// Computes dot products on two arrays of vectors, 1 at a time
for (unsigned int i = 0; i < nElements; i++) {
        result[i] = v1[i].x * v2[i].x +
                    v1[i].y * v2[i].y +
                    v1[i].z * v2[i].z;
}

// Dot Product
// Computes dot products on two arrays of vectors, 4 at a time
for (unsigned int i = 0; i < nElements; i += 4) {
    xmm0 = _mm_load_ps(X1 + i);     // Load data into SIMD registers
    xmm1 = _mm_load_ps(Y1 + i);
    xmm2 = _mm_load_ps(Z1 + i);

    xmm3 = _mm_load_ps(X2 + i);
    xmm4 = _mm_load_ps(Y2 + i);
    xmm5 = _mm_load_ps(Z2 + i);

    xmm6 = _mm_mul_ps(xmm0, xmm3);  // Multiply x's together
    xmm7 = _mm_mul_ps(xmm1, xmm4);  // Multiply y's together
    xmm8 = _mm_mul_ps(xmm2, xmm5);  // Multiply z's together

    xmm0 = _mm_add_ps(xmm6, xmm7);  // Add all the values together
    xmm7 = _mm_add_ps(xmm0, xmm8);

    _mm_store_ps(result + i, xmm7); // Save the results
}

列表 10：SSE 版点积示例

采用 _mm_load_ps() 指令运行上述算法首先要将所有数据加载到 SIMD 寄存器中。这样做会获得作为一项参数的数据加载地址。所有拥有 ps 后缀的指令都是单精度版本的指令。许多指令都拥有多个版本，如 _mm_load_pd() 指令还可以用于加载两个双精度数字。同样至关重要的是，数据必须为 16 字节对齐数据。如果不是 16 字节对齐数据，则必须采用一条不同的指令——_mm_loadu_ps() 来运行函数，但这样做不会带来相同的性能提升优势。 
数据加载完毕后，应采用 _mm_mul_ps() 指令完成数字相乘运算。所相乘的元素来自两个寄存器的相匹配元素。

图 3：_mm_mul_ps() 图示

在两个数组的相对应元素互相相乘后，应将各个乘积相加。执行求和运算应采用 _mm_add_ps() 函数。它与求积函数的工作原理一样，唯一的不同是它解决数字相加问题。最后，应采用 _mm_store_ps() 将运算结果写入内存。需要再次提醒的是，写入地址必须可容纳 16 字节数据。请访问：http://software.intel.com/zh-cn/articles/tickertape/，获取包含该函数的代码样本以及一个用于比较不同方法性能的简单测试框架。

未来展望

在 Ticker Tape 演示中显然存在 SSE 需要改进的方面。其中一个应该是它需要进一步支持对代码进行矢量化，而不应仅针对一个函数实施矢量化。此外，SSE 还应能与即将推出的英特尔® 高级矢量扩展指令集（英特尔® AVX）进行协作。另外一个值得探索的有趣领域是使用英特尔® C++ 编译器自动矢量化代码。目前存在一些可帮助编译器执行矢量化操作的代码模式和结构。经修改后的 Ticker Tape 代码可采用这些代码模式和结构，便于 SSE 仍然保持幕后运行状态。

转自：http://blog.csdn.net/qiaoliang90/archive/2010/06/10/5660406.aspx