SIMD 編程的優勢 –TickerTape Part 2

引自：http://software.intel.com/zh-cn/articles/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>

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>

列表 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 倍。

 

 

表 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 仍然保持幕後運行狀態。