之前我寫了一篇用GPU繪制曼德勃羅特（Mandelbrot）集圖像的文章，里面使用的技術是與DirectX 11繼承在一起的DirectCompute。DirectCompute執行在GPU上的kernel代碼，必須用一種特殊的HLSL語言來編寫。雖然這種語言有些類似於C，但一些特殊的細節使得沒接觸過DirectX的開發人員很不適應。相比於kernel代碼，驅動HLSL所要進行的准備工作那簡直麻煩得要命，所以我在那篇博客里索性略去了。如果要想要體會一下DirectCompute那麻煩又不直觀的API，可以參考我翻譯的另一篇博客。后來我轉向了OpenCL，這是一種異構並行計算的行業標准。它的的kernel語言和API相比DirectCompute那可是簡單了好幾倍，又有.NET封裝類庫，無疑是個方便的選擇。但是經過我的實踐，至少在我的AMD 5850顯卡上，OpenCL的性能沒有DirectCompute那么好。能否魚和熊掌兼得呢？現在這個選擇出現了：C++ AMP
C++ AMP全名C++ Accelerated Massive Parallelism（加速大規模並行計算）。是微軟提出的基於C++的異構化並行計算平台。它將隨Visual Studio 11一起發布，目前為預覽版本。所謂異構並行計算，主要的需求就來自於GPU通用計算的崛起。GPU非常適合大規模數據並行算法，即同一程序應多多組不同的數據進行並行運算。然而GPU的架構與主流CPU不同，而且常常更新換代，沒法采用傳統編程語言來編程。現有的GPU多數編程方案，如DirectCompute和OpenCL，都要使用不同的語言或編譯器來編寫運行於GPU上的kernel部分和運行在CPU上的host部分。C++ AMP統一了這兩部分，可以用同一個編譯器，同一種語法來編寫kernel代碼；無需任何編譯器選項或設置。更甚至C++ AMP的API簡單到了極致，比OpenCL的方便程度更上了一個層次。下面我們就用C++ AMP來實現經典的曼德勃羅特集圖形繪制。
先來回顧一下曼德勃羅特集，這個是很多並行計算書都用來做例子的經典問題，因為它是一個典型的易並行算法。選任意復數c，讓z從0開始采用以下公式
進行迭代。不同的c會產生很不一樣的效果。有的c可以一直迭代下去（z的模都在半徑為2的圓以內），而有的c會導致迭代若干次之后逃逸（其模會變得大於2，進而趨近於無窮）。不同的c迭代次數相差很大，如果把平面上的每個點坐標都當作c來進行迭代，然后把迭代次數繪制成不同的顏色，就會形成絢麗的圖形。而不同c可以完全獨立地進行獨立，不會影響其他c值的計算結果。因此算法可以是完全並行的。
首先我們要開辟出一塊空間來表示最后生成的圖像。圖像上的每一個像素將會映射到一個復數c上，這樣a x b的圖像，就表示a x b個不同的復數，分別進行上述迭代。C++ AMP提供了兩種表示數據的方式：array模板類和array_view模板類。array模板類就像C++ 11的STL容器array一樣，是一個固定尺寸、固定維度的數組。和STL array唯一的不同是C++ AMP的array直接分配在GPU的顯存中。聲明array時要傳入兩個模板參數——array元素的類型和維度。比如我們要構建一個unsigned int類型的二維數組，就可以這樣寫：
int a = 100, b = 100; array<unsigned int, 2> buffer(a, b); |
注意這個array模板是定義在Concurrency命名空間里的，需要提前using。
有時我們需要把CPU上的數據傳送到GPU上，或者希望GPU計算的結果直接寫入指定的CPU內存中，這時我們就可以用第二個模板類array_view。array_view本身不會分配任何內存，它只是一個CPU數據buffer的包裝，其中數據buffer可以是一個指針、數組，也可以是諸如std::vector等各種常用容器。array_view僅僅當數據訪問發生時才進行必要的數據拷貝。聲明array_view和聲明array一樣，但必須傳入它的基礎容器：
int a = 100, b = 100; std::vector<unsigned int> myarray(a * b); array_view<unsigned int, 2> view(a, b, myarray); |
接下來最重要的是讓代碼在GPU上執行，完成這個任務的入口方法是parallel_for_each方法。它非常像.NET 4里面的Parallel.ForEach方法，是一個通過lambda函數來分派任務的方法。我們先來完成一個簡單的任務，現在有兩個array_view，要生成一個array，其內容是兩個array_view中相應元素的和。
std::vector<unsigned int> arr1(10); std::vector<unsigned int> arr2(10); array_view<unsigned int, 1> view1(10, arr1); array_view<unsigned int, 1> view2(10, arr2); array<unsigned int, 1> r(10); parallel_for_each(r.grid, [view1, view2, &r](index<1> i) restrict(direct3d) { r[i] = view1[i] + view2[i]; }); |
觀察這里的parallel_for_each，有許多有趣的信息。首先它的第一個參數是一個grid類型的對象。grid是C++ AMP用來描述並行計算拓撲結構的對象。因為經過parallel_for_each方法分派的計算是在顯卡計算的，我們需要事先指定有多少GPU線程來完成計算，以及這些線程的排列方式。這個任務當中，我們想要把兩個array_view表示的數據想加，為了達到最大並行性，我們給數組的每一個元素都分配一個線程。代碼里的r.grid就表示線程的排列方式由數組r自身的拓撲結構決定——在這個例子中就是一個1 x 10的一維結構。注意，GPU上通常可以開成百上千的線程，幾乎沒有任何代價，所以線程的拓撲結構應當盡量達到最大的並行度。
接下來那個方括號語法是C++的lambda表達式語法。C++ 11開始支持lambda表達式，和其他語言的lambda表達式一樣，這是一種就地聲明匿名函數的語法。lambda表達式在C++ 11中具有極其重要的地位。一開始的方括號是lambda表達式的捕獲列表。與C#等語言的lambda表達式不同，C++得lambda表達式需要顯式進行變量捕獲，只有列在捕獲列表中的變量才能在lambda函數體內使用。注意，捕獲的時候直接寫名字的變量（如上面例子中的view1和view2）是按值捕獲的變量，在lambda函數體內修改它們的值，不會改變原始變量。而名字前加一個&符號的捕獲變量是按引用捕獲的變量，lambda函數對其修改會反應到原始變量中。在C++ AMP的parallel_for_each調用時有一個規定：所有變量都必須按值捕獲，除了array對象，array對象必須按引用捕獲。這就是上面例子為何要這樣寫的原因。
捕獲列表之后是lambda函數的參數，定義方法和普通函數的參數方法是一樣的。parallel_for_each要求傳入的函數接受一個index類型的特殊參數。這個參數是實際運行時的線程編號。因為線程的拓撲結構已經由grid參數決定了，所以這個index運行時就會是grid表示范圍中的一個值，維度和grid一樣。比如上面的例子，index就會表示0-9號線程中的某一個，這些線程都是並行執行的。
參數表之后，出現了一個新的修飾符restrict(direct3d)。這是專為C++ AMP設計的新修飾符，用來表示函數的語法約束。眾所周知GPU和CPU是有很多不同的，在GPU上執行的方法無法進行某些操作。比如GPU代碼無法調用Windows API，也無法遞歸。restrict(direct3d)告訴編譯器檢查函數體中的代碼，並且在編譯時就能檢測出所有不滿足約束的語法。restrict修飾符還是函數簽名的一部分，可以針對不同的restrict修飾符進行函數重載，這和const等C++原有的修飾符是一樣的。除了restrict(direct3d)，還有一種restrict(cpu)修飾符。所有現有C++函數都默認帶有restrict(cpu)修飾符。一個函數還可以同時約束為CPU和GPU代碼，只要這么寫即可：restrict(cpu, direct3d)。在restrict(direct3d)函數中調用其它的函數，會自動選擇restrict(direct3d)的重載。稍后我們會看到，C++ AMP定義了一組restrict(direct3d)版本的數學函數，如log和sin等。在restrict(direct3d)代碼中調用這些數學函數就會自動使用GPU版本，而普通的CPU代碼則會去調用不帶修飾的CPU版本。
最后是lambda函數的函數體，里面的代碼非常直白，就是直接對當前線程index所對應的數組元素進行相加。當你運行上述代碼，該代碼就會在GPU上以並行計算的方式完成。沒有任何額外的初始化代碼，隱藏代碼或者編譯器選項設置——只要這么寫就可以運行了。是不是很方便？
下面我們就來完成GPU上的曼德勃羅特集生成算法。先看頭文件mandelbrot.h，里面include了amp.h——C++ AMP的主要頭文件。還聲明了一個fp_t類型，根據宏來控制它表示float或double。目前Windows 7上的C++ AMP還不支持雙精度浮點運算。
#include "amp.h" //#define FP64 #if defined FP64 #define _F(x) x typedef double fp_t; #else #define _F(x) x##f typedef float fp_t; #endif void generate_mandelbrot( Concurrency::array_view<unsigned int, 2> result, unsigned int max_iter, fp_t real_min, fp_t imag_min, fp_t real_max, fp_t imag_max ); |
接下來是實現的代碼，還順便生成了一個HSB（色調、飽和度、亮度）到RGB轉換的GPU函數：
#include "stdafx.h" #include "mandelbrot.h" using namespace Concurrency; unsigned int set_hsb (float hue, float saturate, float bright) restrict (direct3d); void generate_mandelbrot( array_view<unsigned int, 2> result, unsigned int max_iter, fp_t real_min, fp_t imag_min, fp_t real_max, fp_t imag_max ) { int width = result.extent.get_x(); int height = result.extent.get_y(); fp_t scale_real = (real_max - real_min) / width; fp_t scale_imag = (imag_max - imag_min) / height; parallel_for_each(result.grid, [=](index<2> i) restrict(direct3d) { int gx = i.get_x(); int gy = i.get_y(); fp_t cx = real_min + gx * scale_real; fp_t cy = imag_min + (height - gy) * scale_imag; fp_t zx = _F(0.0); fp_t zy = _F(0.0); fp_t temp; fp_t length_sqr; unsigned int count = 0; do { count++; temp = zx * zx - zy * zy + cx; zy = 2 * zx * zy + cy; zx = temp; length_sqr = zx * zx + zy * zy; } while((length_sqr < _F(4.0)) && (count < max_iter)); //faster using multiplication than division float n = count * 0.0078125f; // n = count / 128.0f; float h = 1.0f - 2.0f * fabs(0.5f - n + floor(n)); //turn points at maximum iteration to black float bfactor = direct3d::clamp((float)(max_iter - count), 0.0f, 1.0f); result[i] = set_hsb(h, 0.7f, (1.0f - h * h * 0.83f) * bfactor); }); } unsigned int set_hsb (float hue, float saturate, float bright) restrict (direct3d) { float red, green, blue; float h = (hue * 256) / 60; float p = bright * (1 - saturate); float q = bright * (1 - saturate * (h - (int)h)); float t = bright * (1 - saturate * (1 - (h - (int)h))); switch ((int)h) { case 0: red = bright, green = t, blue = p; break; case 1: red = q, green = bright, blue = p; break; case 2: red = p, green = bright, blue = t; break; case 3: red = p, green = q, blue = bright; break; case 4: red = t, green = p, blue = bright; break; case 5: case 6: red = bright, green = p, blue = q; break; } unsigned int ired, igreen, iblue; ired = (unsigned int)(red * 255.0f); igreen = (unsigned int)(green * 255.0f); iblue = (unsigned int)(blue * 255.0f); return 0xff000000 | (ired << 16) | (igreen << 8) | iblue; } |
注意，在這段代碼中我們使用的array_view是二維的，所以index也相應是一個二維的結構。我們可以從index的get_x()，get_y()等方法得到線程的編號。代碼中還使用了若干小技巧，比如判斷復數模的平方大於4而不是模大於2。因為開平方運算即使在GPU上也是個比較慢的操作，需要盡量避免。方法的最后部分使用一個特殊的算法將迭代數轉化成了顏色。這個公式是我反復試驗得到的，只是為了視覺上好看。大家自己實踐可以根據自己喜好隨便改。
C++ AMP的parallel_for_each調用是一種“准同步”的調用，當parallel_for_each完成時，GPU計算並不一定結束了，但代碼會繼續執行。只有當你訪問計算的結果——比如包含結果的array_view時，它就會阻塞線程等待GPU計算的完成。我們可以調用array_view的synchornize方法來顯式地等待。C++ AMP也提供真正異步接口，但用起來比較麻煩，在這個例子中就不采用了。下面的代碼就是調用以及同步的代碼。這里展示代碼將結果保存在一個位圖中：
int _tmain(int argc, _TCHAR* argv[]) { GdiplusStartupInput gdiplusStartupInput; ULONG_PTR gdiplusToken; GdiplusStartup(&gdiplusToken, &gdiplusStartupInput, nullptr); const int edge = 1024; std::vector<unsigned int> a(edge * edge); array_view<unsigned int, 2> av(edge, edge, a); generate_mandelbrot(av, 512, -2, -2, 2, 2); av.synchronize(); //construct a image { Bitmap output(edge, edge, edge * 4, PixelFormat32bppARGB, reinterpret_cast<BYTE*>(a.data())); CLSID pngClsid; GetEncoderClsid(L"image/png", &pngClsid); output.Save(L"test.png", &pngClsid, nullptr); } GdiplusShutdown(gdiplusToken); return 0; } |
以上代碼繪制了實部從-2到2，虛部從-2i到2i范圍的曼德勃羅特集圖形，繪制到一個邊長1024的位圖上，最大迭代次數512次。注意其中synchronize方法的調用。所得結果（縮小圖）為：
曼德勃羅特集圖形每一點放大都有無窮無盡美麗的花紋，為了能夠任意探索圖形的細節，我還編寫了一個Windows界面的程序，使用Direct2D將所繪圖像展示出來，而且加入了拖拽和鼠標滾輪的交互。實際運行效果非常酷，因為GPU並行計算的強大性能，所有交互都是平滑實時進行的。就好像是一張可以無限放大的圖片一樣。美中不足的是當前C++ AMP的Windows 7實現只支持單精度浮點數，所以放大到一定倍數就會模糊了。等到Windows 8的穩定版本出現，大家才能體驗到雙精度浮點的美妙結果。
我已經將所有源代碼上傳到了Github上，地址為：https://github.com/Ninputer/AMP-Demo
主要，要編譯和運行此例子，需要以下軟硬件環境：
1. Windows 7 （不支持XP和Vista）
2. 硬件支持DirectX 11所有特性的顯卡。如Nvidia Geforce 400,500系列顯卡，AMD Radeon HD5000,6000,7000系列顯卡。無法運行於不支持DirectX 11的顯卡上。
3. 需要Visual Studio 11 Developer Preview。下載地址
4. 需要Windows SDK。下載地址
我還翻譯了一個Build大會關於C++ AMP的視頻，是快速了解C++ AMP的最佳途徑。觀看地址：http://www.tudou.com/playlist/p/l13690258i105735621.html
歡迎關注我的微博http://weibo.com/ninputer，一起討論C++ AMP和其他有趣的技術~