出處：http://www.cnblogs.com/allenlooplee/archive/2013/01/17/2863866.html
遇見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的倍數。
