遇見C++ AMP：GPU的執行緒模式和記憶體模型

遇見C++ AMP：GPU的執行緒模式和記憶體模型

 

Written by Allen Lee

 

I don't care where the enemies are / Can't be stopped / All I know / Go hard
– Linkin Park, Lost In The Echo

 

C++ AMP、CUDA和OpenCL，選擇哪個？

      在《遇見C++ AMP：在GPU上做並行計算》發布之後，我曾被多次問及為何選擇C++ AMP，以及它與CUDA、OpenCL等相比有何優勢，看來有必要在進入正題之前就這個問題發表一下看法了。

      在眾多可以影響決策的因素之中，平台種類的支援和GPU種類的支援是兩個非常重要的因素，它們聯合起來足以直接否決某些選擇。如果我們把這兩個因素看作兩個維度，可以把平面分成四個象限，C++ AMP、CUDA和OpenCL分別位於第二象限、第四象限和第一象限，1所示。如果你想通吃所有平台和所有GPU，OpenCL是目前唯一的選擇，當然，你也需要為此承擔相當的複雜性。CUDA是一個有趣的選擇，緊貼最新的硬體技術、數量可觀的行業應用和類庫支援使之成為一個無法忽視的選擇，但是，它只能用於NVIDIA的GPU極大地限制了它在商業應用上的採用，我想你不會為了運行我的應用程式特意把顯卡換成NVIDIA的。C++ AMP的情況剛好相反，它適用於各種支援DirectX 11的GPU，但只能在Windows上運行。

圖 1

      這些技術都有自己的特點和位置，你應該根據項目的具體情況選擇合適的解決方案。如果你正在從事的工作需要進行大量計算，你想儘可能利用硬體特性對演算法進行最佳化，而你的機器剛好有一塊NVIDIA的顯卡，並且你不需要在其他機器上重複執行這些計算，那麼CUDA將是你的不二之選。儘管NVIDIA已經開源CUDA編譯器，並且歡迎其他廠商通過CUDA編譯器SDK添加新的語言/處理器，但AMD不太可能會為它提供在AMD的GPU上啟動並執行擴充，畢竟它也有自己的基於OpenCL的AMD APP技術。如果你正在從事Windows應用程式的開發工作，熟悉C++和Visual Studio，並且希望藉助GPU進一步提升應用程式的效能，那麼C++ AMP將是你的不二之選。儘管微軟已經開放C++ AMP規範，Intel的Dillon Sharlet也通過Shevlin Park項目驗證了在Clang/LLVM上使用OpenCL實現C++ AMP是可行的，但這不是一個產品層級的商用編譯器，Intel也沒有宣布任何發布計劃。如果你確實需要同時相容Windows、Mac OS X和Linux等多個作業系統，並且需要同時支援NVIDIA和AMD的GPU，那麼OpenCL將是你的不二之選。

 

GPU線程的執行

      在《遇見C++ AMP：在GPU上做並行計算》裡，我們通過extent對象告訴parallel_for_each函數建立多少個GPU線程，那麼，這些GPU線程又是如何組織、分配和執行的呢？

      首先，我們建立的GPU線程會被分組，分組的規格並不固定，但必須滿足兩個條件：對應的維度必須能被整除，分組的大小不能超過1024。假設我們的GPU線程是一維的，共8個，2所示，則可以選擇每2個GPU線程為1組或者每4個GPU線程為1組，但不能選擇每3個GPU線程為1組，因為剩下的2個GPU線程不足1組。

圖 2

      假設我們建立的GPU線程是二維的，3 x 4，共12個，3所示，則可以選擇3 x 1或者3 x 2作為分組的規格，但不能選擇2 x 2作為分組的規格，因為剩下的4個GPU線程雖然滿足分組的大小，但不滿足分組的形狀。每個分組必須完全相同，包括大小和形狀。

圖 3

      為了便於解釋，我們的GPU線程只有寥寥數個，但真實案例的GPU線程往往是幾十萬甚至幾百萬個，這個時候，分組的規格會有大量選擇，我們必須仔細判斷它們是否滿足條件。假設我們的GPU線程是640 x 480，那麼16 x 48、32 x 16和32 x 32都可以選擇，它們分別產生40 x 10、20 x 30和20 x 15個分組，但32 x 48不能選擇，因為它的大小已經超過1024了。

      接著，這些分組會被分配到GPU的流多處理器（streaming multiprocessor），每個流多處理器根據資源的使用方式可能分得一組或多組GPU線程。在執行的過程中，同一組的GPU線程可以同步，不同組的GPU線程無法同步。你可能會覺得這種有限同步的做法會極大地限制GPU的作為，但正因為組與組之間是相互獨立的，GPU才能隨意決定這些分組的執行順序。這有什麼好處呢？假設低端的GPU每次只能同時執行2個分組，那麼執行8個分組需要4個執行循環，假設高端的GPU每次可以同時執行4個分組，執行8個分組只需2個執行循環，4所示，這意味著我們寫出來的程式具備延展性，能夠自動適應GPU的計算資源。

圖 4

      說了這麼多，是時候看看代碼了。parallel_for_each函數有兩種模式，一種是簡單模式，我們通過extent對象告訴它建立多少GPU線程，C++ AMP負責對GPU線程進行分組，另一種是分組模式，我們通過tiled_extent對象告訴它建立多少GPU線程以及如何進行分組。建立tiled_extent對象非常簡單，只需在現有的extent對象上調用tile方法，並告知分組的規格就行了，如代碼1所示。值得提醒的是，分組的規格是通過模板參數告訴tile方法的，這意味著分組的規格必須在編譯時間確定下來，C++ AMP目前無法做到運行時動態分組。

代碼 1

      既然C++ AMP不支援運行時動態分組，肯定會為簡單模式預先定義一些分組的規格，那麼C++ AMP又是如何確保它們能被整除？假設我們建立的GPU線程是一維的，共10000個，C++ AMP會選擇每256個GPU線程為1組，把前面9984個GPU線程分成39個分組，然後補充240個GPU線程和剩下的16個GPU線程湊夠1組，執行的時候會通過邊界測試確保只有前10000個GPU線程執行我們的代碼。對於二維和三維的情況，C++ AMP也會採取這種補充GPU線程的策略，只是分組的規格不同，必要時還會重新排列GPU線程，以便分組能夠順利完成。需要說明的是，簡單模式背後採取的策略屬於實現細節，在這裡提及是為了滿足部分讀者的好奇心，你的演算法不該對它有所依賴。

 

共用記憶體的訪問

      既然簡單模式可以自動分組，為何還要大費周章使用分組模式？為了回答這個問題，我們先要瞭解一下GPU的記憶體模型。在Kernel裡，我們可以訪問全域記憶體、共用記憶體和寄存器，5所示。當我們通過array_view對象把資料從主機記憶體複製到顯卡記憶體時，這些資料會被儲存在全域記憶體，直到應用程式退出，所有GPU線程都能訪問全域記憶體，不過訪問速度很慢，大概需要1000個GPU刻度，大量的GPU線程反覆執行這種高延遲的操作將會導致GPU計算資源的閑置，從而降低整體的計算效能。

圖 5

      為了避免反覆從全域記憶體訪問相同的資料，我們可以把這些資料緩衝到寄存器或者共用記憶體，因為它們整合在GPU晶片裡，所以訪問速度很快。當我們在Kernel裡聲明一個基本類型的變數時，它的資料會被儲存在寄存器，直到GPU線程執行完畢，每個GPU線程只能訪問自己的寄存器，寄存器的容量非常小，不過訪問速度非常快，只需1個GPU刻度。當我們在Kernel裡通過tile_static關鍵字聲明一個變數時，它的資料會被儲存在共用記憶體（也叫tile_static記憶體），直到分組裡的所有GPU線程都執行完畢，同一組的GPU線程都能訪問相同的共用記憶體，共用記憶體的容量很小，不過訪問速度很快，大概需要10個GPU刻度。tile_static關鍵字只能在分組模式裡使用，因此，如果我們想使用共用記憶體，就必須使用分組模式。

      如果資料只在單個GPU線程裡反覆使用，可以考慮把資料緩衝到寄存器。如果資料會在多個GPU線程裡反覆使用，可以考慮把資料緩衝到共用記憶體。共用記憶體的緩衝策略是對全域記憶體的資料進行分組，然後把這些分組從全域記憶體複製到共用記憶體。假設我們需要緩衝4 x 4的資料，可以選擇2 x 2作為分組的規格把資料分成4組，6所示。以右上方的分組為例，我們需要4個GPU線程分別把這4個資料從全域記憶體複製到共用記憶體。複製的過程涉及兩種不同的索引，一種是相對於所有資料的全域索引，用於從全域記憶體訪問資料，另一種是相對於單個分組的本地索引，用於從共用記憶體訪問資料，比如說，全域索引(1, 2)對應本地索引(1, 0)。

圖 6

      在分組模式裡，我們可以通過tiled_index對象訪問索引資訊，它的global屬性返回全域索引，local屬性返回本地索引，tile屬性返回分組索引，它是分組作為一個整體相對於其他分組的索引，tile_origin屬性返回分組原點的全域索引，它是分組裡的(0, 0)位置上的元素的全域索引。還是以右上方的分組為例，(1, 2)位置的global屬性的值是(1, 2)，local屬性的值是(1, 0)，tile屬性的值是(0, 1)，tile_origin屬性的值是(0, 2)。tiled_index對象將會通過Lambda的參數傳給我們，我們將會在Kernel裡通過它的屬性訪問全域記憶體和共用記憶體。

      說了這麼多，是時候看看代碼了。正如extent對象搭配index對象用於簡單模式，tiled_extent對象搭配tiled_index對象用於分組模式，使用的時候，兩者的模板參數必須完全符合，如代碼2所示。parallel_for_each函數將會建立16個GPU線程，每4個GPU線程為1組，同一組的GPU線程共用一個2 x 2的陣列變數，每個元素由一個GPU線程負責複製，每個GPU線程通過tiled_index對象的global屬性獲知從全域記憶體的哪個位置讀取資料，通過local屬性獲知向共用記憶體的哪個位置寫入資料。

代碼 2

      因為緩衝的資料會在多個GPU線程裡使用，所以每個GPU線程必須等待其他GPU線程緩衝完畢才能繼續執行後面的代碼，否則，一些GPU線程還沒開始快取資料，另一些GPU線程就開始使用資料了，這樣計算出來的結果肯定是錯的。為了避免這種情況的發生，我們需要在代碼2後面加上一句idx.barrier.wait();，加上之後的效果就像設了一道閘門，7所示，它把整個代碼分成兩個階段，第一階段快取資料，第二階段計算結果，緩衝完畢的GPU線程會在閘門前面等待，當所有GPU線程都緩衝完畢時，就會開啟閘門讓它們進入第二階段。

圖 7

      總的來說，使用分組模式是為了藉助共用記憶體減少全域記憶體的訪問，緩衝的過程已經包含了一次全域記憶體的訪問，因此，如果我們的演算法只需訪問全域記憶體一次，比如《遇見C++ AMP：在GPU上做並行計算》的"並行計算矩陣之和"，那麼快取資料不會帶來任何改善，反而增加了代碼的複雜性。

 

並行計算矩陣之積

      矩陣的乘法需要反覆訪問相同的元素，非常適合用來示範分組模式。接下來，我們將會分別使用簡單模式和分組模式實現矩陣的乘法，然後通過對比瞭解這兩種實現的區別。

      設矩陣

 

求AB。設C = AB，根據定義，，其中，。你可以把這個公式想象成矩陣A的第i行和矩陣B的第j列兩個數組對應位置的元素相乘，然後相加。

      如何把這些數學描述翻譯成代碼呢？第一步，定義A、B和C三個矩陣，如代碼3所示，iota函數可以在指定的起止位置之間填充連續的數字，正好滿足這裡的需求。

代碼 3

      第二步，計算矩陣C的元素，如代碼4所示，整個Kernel就是計算的求和公式， 因為每個元素的計算都是獨立的，所以非常適合并行執行。

代碼 4

      在執行代碼4的時候，parallel_for_each函數將會建立36個GPU線程，每個GPU線程計算矩陣C的一個元素，因為這36個GPU線程會同時執行，所以計算矩陣C的時間就是計算一個元素的時間。這聽起來已經很好，還能更好嗎？仔細想想，計算需要訪問矩陣A的第i行一次，那麼，計算矩陣C的第i行將會訪問矩陣A的第i行M次，M是矩陣C的列數，在這裡是6；同理，計算矩陣C的第j列將會訪問矩陣B的第j列M次，M是矩陣C的行數，在這裡也是6。因為A、B和C三個矩陣的資料是儲存在全域記憶體的，所以最佳化的關鍵就是減少全域記憶體的訪問。

      根據上一節的討論，我們將會使用分組模式，並把需要反覆訪問的資料從全域記憶體緩衝到共用記憶體，那麼，使用分組模式會對效能帶來多少改善，又對演算法造成多少影響呢，這正是我們接下來需要探討的。

      第一步，選擇2 x 2作為分塊的規格對A、B兩個矩陣進行分塊處理

 

分塊矩陣的乘法和普通矩陣的乘法是一樣的，設C = AB，根據定義，分塊矩陣，其中，。

      第二步，把parallel_for_each函數改成分組模式，如代碼5所示。T是子塊的邊長，W是分塊矩陣A的列數，也是分塊矩陣B的行數。

代碼 5

      第三步，分別緩衝和，如代碼6所示。因為它們都是2 x 2的矩陣，所以緩衝它們的工作需要4個GPU線程協同完成。正確緩衝的關鍵在於弄清每個GPU線程負責全域記憶體和共用記憶體的哪些位置，共用記憶體的位置可以通過tiled_index對象的local屬性獲知，而全域記憶體的位置則需要換算一下，因為i和j是針對矩陣C而不是矩陣A和矩陣B的。每個GPU線程只是分別從矩陣A和矩陣B緩衝一個元素，根據定義，從矩陣A緩衝的元素必定位於第i行，而從矩陣B緩衝的元素必定位於第j列。當我們緩衝時，子塊位於分塊矩陣A的左上方，tiled_index對象的local屬性和global屬性指向相同的列，因此，目標元素位於矩陣A的第h列，當我們緩衝時，我們已經從左至右跨過了w個子塊，因此，目標元素位於矩陣A的第h + w * T列。同理，當我們緩衝時，子塊位於分塊矩陣B的左上方，目標元素位於矩陣B的第g行，當我們緩衝時，我們已經從上到下跨過了w個子塊，因此，目標元素位於矩陣B的第g + w * T行。4個GPU線程都緩衝完畢就會進入第二階段。

代碼 6

      第四步，計算，如代碼7所示。這是兩個2 x 2的普通矩陣相乘，需要4個GPU線程協同完成，每個GPU線程計算結果矩陣的一個元素，然後加到變數sum上。4個GPU線程都計算完畢就會重複第三、四步，緩衝和，計算。如果還有其他子塊，那麼這個過程會一直重複下去。最終，每個GPU線程匯總矩陣的一個元素。

代碼 7

      最後一步，把匯總的結果儲存到矩陣C，如代碼8所示。

代碼 8

      至此，我相信你已經深刻地體會到分組模式的複雜性。CPU擁有更強的控制組件和更大的快取區域，可以預測和決定應該緩衝哪些資料，而GPU則把原本屬於它們的空間留給更多的運算組件，把緩衝的控制權交給程式員，這意味著緩衝的邏輯將會滲透到業務的邏輯，從而增加了代碼的複雜性。

      那麼，這樣做是否值得？我們可以算一下，在代碼4裡，av、bv和cv都是位於全域記憶體，每次訪問都要1000個GPU刻度， 讀取N次av和bv，寫入一次cv，總共耗時2000 * N + 1000個GPU刻度，當N = 4時，總共耗時9000個GPU刻度。在代碼6、7、8裡，at和bt都是位於共用記憶體，每次訪問只要10個GPU刻度，讀取W次av和bv，寫入W次at和bt，讀取W * T次at和bt，寫入一次cv，總共耗時 (2000 + 20 + 20 * T) * W + 1000，當N = 4，T = 2時，W = 2，總共耗時5120個GPU刻度，約為簡單模式的56.89%，效能的改善非常明顯。如果我們增加矩陣和子塊的大小，這個差距就會更加明顯，令N = 1024，T = 16，則簡單模式總共耗時2049000個GPU刻度，而分組模式總共耗時150760個GPU刻度，後者是前者的7.36%。

      當然，這並不是簡單模式和分組模式在效能上的精確差距，因為我們還沒考慮訪問寄存器和算術運算的耗時，但這些操作的耗時和訪問全域記憶體的相比簡直就是小巫見大巫，即使把它們考慮進去也不會對結果造成太大影響。

 

你可能會問的問題

1. 什麼時候不能使用分組模式？

      分組模式最高只能處理三維的資料結構，四維或者更高的資料結構必須使用簡單模式，事實上，簡單模式會把四維或者更高的資料結構換算成三維的。如果演算法沒有反覆從全域記憶體訪問相同的資料，那也不必使用分組模式。

 

2. 分組的限制有哪些？

      分組的大小不能超過1024，對應的維度必須能被整除，第一、二維度大小不能超過1024，第三維度大小不能超過64，分組的總數不能超過65535。

 

3. 分組的大小越大越好嗎？

      不是的，我們使用分組模式主要是為了使用共用記憶體，一般情況下，分組的大小和它使用的共用記憶體成正比，每個流多處理器的共用記憶體是有限的，比如說，NVIDIA的GK110 GPU的最大規格是48K，每個流多處理器最多可以同時容納16個分組，這意味著每個分組最多隻能使用3K，如果每個分組使用4K，那麼每個流多處理器最多隻能同時容納12個分組，這意味著流多處理器的計算能力沒被最大限度使用。

 

4. 在設定分組的大小時需要考慮Warp Size嗎？

      在計算域允許的情況下儘可能考慮，NVIDIA的Warp Size是32，AMD的Wavefront Size是64，因此，分組的大小最好是64的倍數。

 

*聲明：本文已經首發於InfoQ中文站，著作權，《遇見C++ AMP：GPU的執行緒模式和記憶體模型》，如需轉載，請務必附帶本聲明，謝謝。