GPU加速03:多流和共享內存—讓你的CUDA程序如虎添翼的優化技術！

  • 2019 年 12 月 26 日
  • 筆記

本文為英偉達GPU計算加速系列的第三篇，前兩篇文章為：

  1. AI時代人人都應該了解的GPU知識：主要介紹了CPU與GPU的區別、GPU架構、CUDA軟件棧簡介。
  2. 超詳細Python Cuda零基礎入門教程：主要介紹了CUDA核函數，Thread、Block和Grid概念，內存分配，並使用Python Numba進行簡單的並行計算。

閱讀完前兩篇文章後，相信讀者應該能夠將一些簡單的CPU代碼修改成GPU並行代碼，但是對計算密集型任務，僅僅使用前文的方法還是遠遠不夠的，GPU的並行計算能力未能充分利用。本文將主要介紹一些常用性能優化的進階技術，這部分對編程技能和硬件知識都有更高的要求，建議讀者先閱讀本系列的前兩篇文章，甚至閱讀英偉達官方的編程手冊，熟悉CUDA編程的底層知識。當然，將這些優化技巧應用之後，程序將獲得更大的加速比，這對於需要跑數小時甚至數天的程序來說，收益非常之大。

本文仍然使用Python版的Numba庫調用CUDA，有更複雜需求的朋友可以直接使用C/C++調用CUDA，並閱讀英偉達的官方文檔。C/C++對數據的控制更細緻，是英偉達官方推薦的編程語言，所能提供的編程接口更全面。

CUDA C Programming Guide ：https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html

下一篇文章將提供實戰案例，包括金融領域期權定價的GPU實現。

CUDA優化方向

我之前的文章中提到，CPU + GPU 是一種異構計算的組合，各有獨立的內存，GPU的優勢是更多的計算核心。該架構在並行計算上有很大優勢，但是數據需要從主機和設備間相互拷貝，會造成一定的延遲。因此，要從下面兩個方面來優化GPU程序：

  • 充分利用GPU的多核心，最大化並行執行度
  • 優化內存使用，最大化數據吞吐量，減少不必要的數據拷貝

哪個方向有更大收益，最終還是要看具體的計算場景。英偉達提供了非常強大的性能分析器nvprof和可視化版nvvp，使用性能分析器能監控到當前程序的瓶頸。據我了解，分析器只支持C/C++編譯後的可執行文件，Python Numba目前應該不支持。

Profiler User's Guide : https://docs.nvidia.com/cuda/profiler-users-guide/index.html

並行計算優化

網格跨度

在上一篇文章中，我曾提到，CUDA的執行配置：[gridDim, blockDim]中的blockDim最大只能是1024，但是並沒提到gridDim的最大限制。英偉達給出的官方回復是gridDim最大為一個32位整數的最大值，也就是2,147,483,648，大約二十億。這個數字已經非常大了，足以應付絕大多數的計算，但是如果對並行計算的維度有更高需求呢？網格跨度有更好的並行計算效率。

並行計算數大於線程數

這裡仍然以[2, 4]的執行配置為例，該執行配置中整個grid只能並行啟動8個線程，假如我們要並行計算的數據是32，會發現後面8號至31號數據共計24個數據無法被計算。

網格跨度

我們可以在0號線程中，處理第0、8、16、24號數據，就能解決數據遠大於執行配置中的線程總數的問題，用程序表示，就是在核函數里再寫個for循環。以打印為例，代碼如下：

from numba import cuda

注意，跨步大小為網格中線程總數，用gridDim.x * blockDim.x來計算。for循環的step是網格中線程總數，這也是為什麼將這種方式稱為網格跨步。如果網格總線程數為1024，那麼0號線程將計算第0、1024、2048…號的數據。這裡我們也不用再明確使用if (idx < N)來判斷是否越界，因為for循環也有這個判斷。

使用網格跨步的優勢主要有：

  1. 擴展性：可以解決數據量比線程數大的問題
  2. 線程復用：CUDA線程啟動和銷毀都有開銷，主要是線程內存空間初始化的開銷；不使用網格跨步，CUDA需要啟動大於計算數的線程，每個線程內只做一件事情，做完就要被銷毀；使用網格跨步，線程內有for循環，每個線程可以干更多事情，所有線程的啟動銷毀開銷更少。
  3. 方便調試：我們可以把核函數的執行配置寫為[1, 1]，如下所示，那麼核函數的跨步大小就成為了1，核函數里的for循環與CPU函數中順序執行的for循環的邏輯一樣，非常方便驗證CUDA並行計算與原來的CPU函數計算邏輯是否一致。
kernel_function[1,1](...)

多流

之前我們討論的並行，都是線程級別的，即CUDA開啟多個線程，並行執行核函數內的代碼。GPU最多就上千個核心，同一時間只能並行執行上千個任務。當我們處理千萬級別的數據，整個大任務無法被GPU一次執行，所有的計算任務需要放在一個隊列中，排隊順序執行。CUDA將放入隊列順序執行的一系列操作稱為流（Stream）。

由於異構計算的硬件特性，CUDA中以下操作是相互獨立的，通過編程，是可以操作他們並發地執行的：

  • 主機端上的計算
  • 設備端的計算（核函數）
  • 數據從主機和設備間相互拷貝
  • 數據從設備內拷貝或轉移
  • 數據從多個GPU設備間拷貝或轉移

數據拷貝和計算的重疊

針對這種互相獨立的硬件架構，CUDA使用多流作為一種高並發的方案：把一個大任務中的上述幾部分拆分開，放到多個流中，每次只對一部分數據進行拷貝、計算和回寫，並把這個流程做成流水線。因為數據拷貝不佔用計算資源，計算不佔用數據拷貝的總線（Bus）資源，因此計算和數據拷貝完全可以並發執行。如圖所示，將數據拷貝和函數計算重疊起來的，形成流水線，能獲得非常大的性能提升。實際上，流水線作業的思想被廣泛應用於CPU和GPU等計算機芯片設計上，以加速程序。

默認流與多流

以向量加法為例，上圖中第一行的Stream 0部分是我們之前的邏輯，沒有使用多流技術，程序的三大步驟是順序執行的：先從主機拷貝初始化數據到設備（Host To Device）；在設備上執行核函數（Kernel）；將計算結果從設備拷貝回主機（Device To Host）。當數據量很大時，每個步驟的耗時很長，後面的步驟必須等前面執行完畢才能繼續，整體的耗時相當長。以2000萬維的向量加法為例，向量大約有幾十M大小，將整個向量在主機和設備間拷貝將佔用佔用上百毫秒的時間，有可能遠比核函數計算的時間多得多。將程序改為多流後，每次只計算一小部分，流水線並發執行，會得到非常大的性能提升。

默認情況下，CUDA使用0號流，又稱默認流。不使用多流時，所有任務都在默認流中順序執行，效率較低。在使用多流之前，必須先了解多流的一些規則：

  • 給定流內的所有操作會按序執行。
  • 非默認流之間的不同操作，無法保證其執行順序。
  • 所有非默認流執行完後，才能執行默認流；默認流執行完後，才能執行其他非默認流。

多流

參照上圖，可將這三個規則解釋為：

  • 非默認流1中，根據進流的先後順序，核函數1和2是順序執行的。
  • 無法保證核函數2與核函數4的執行先後順序，因為他們在不同的流中。他們執行的開始時間依賴於該流中前一個操作結束時間，例如核函數2的開始依賴於核函數1的結束，與核函數3、4完全不相關。
  • 默認流有阻塞的作用。如圖中紅線所示，如果調用默認流，那麼默認流會等非默認流都執行完才能執行；同樣，默認流執行完，才能再次執行其他非默認流。

可見，某個流內的操作是順序的，非默認流之間是異步的，默認流有阻塞作用。

如果想使用多流時，必須先定義流：

stream = numba.cuda.stream()

CUDA的數據拷貝以及核函數都有專門的stream參數來接收流，以告知該操作放入哪個流中執行：

  • numba.cuda.to_device(obj, stream=0, copy=True, to=None)
  • numba.cuda.copy_to_host(self, ary=None, stream=0)

核函數調用的地方除了要寫清執行配置，還要加一項stream參數：

  • kernel[blocks_per_grid, threads_per_block, stream=0]

根據這些函數定義也可以知道，不指定stream參數時，這些函數都使用默認的0號流。

對於程序員來說，需要將數據和計算做拆分，分別放入不同的流里，構成一個流水線操作。

將之前的向量加法的例子改為多流處理，完整的代碼為：

from numba import cuda

是否使用多流的計算時間差距非常大：

gpu vector add time 9.33862018585205

在上面的程序中，我將向量分拆成了5份，同時也創建了5個流，每個流執行1/5的「拷貝、計算、回寫」操作，多個流之間異步執行，最終得到非常大的性能提升。

多流不僅需要程序員掌握流水線思想，還需要用戶對數據和計算進行拆分，並編寫更多的代碼，但是收益非常明顯。對於計算密集型的程序，這種技術非常值得認真研究。

內存優化

我在本系列第一篇文章提到，CPU和GPU組成異構計算架構，如果想從內存上優化程序，我們必須盡量減少主機與設備間的數據拷貝，並將更多計算從主機端轉移到設備端。盡量在設備端初始化數據，並計算中間數據，並盡量不做無意義的數據回寫。

GPU內存硬件結構

GPU的內存結構如圖所示：GPU的計算核心都在Streaming Multiprocessor（SM）上，Multiprocessor里有計算核心可直接訪問的寄存器（Register）和共享內存（Shared Memory）；多個SM可以讀取顯卡上的顯存，包括全局內存（Global Memory）。每個Multiprocessor上的Shared Memory相當於該Multiprocessor上的一個緩存，一般都很小，當前最強的GPU Telsa V100的Shared Memory也只有96KB。注意，Shared Memory和Global Memory的字面上都有共享的意思，但是不要將兩者的概念混淆，Shared Memory離計算核心更近，延遲很低；Global Memory是整個顯卡上的全局內存，延遲高。

英偉達GPU存儲結構

從軟件角度來看，CUDA的線程可以訪問不同級別的存儲，每個Thread有獨立的私有內存；每個Block中多個Thread都可以在該Block的Shared Memory中讀寫數據；整個Grid中所有Thread都可以讀寫Global Memory。Shared Memory的讀寫訪問速度會遠高於Global Memory。內存優化一般主要利用Shared Memory技術。下文將以矩陣乘法為例，展示如何使用Shared Memory來優化程序。

二維和三維執行配置

在解釋內存優化前，先填一下之前埋下的多維執行配置的坑。我們之前使用的threadIdx 和blockIdx變量都是一維的，實際上，CUDA允許這兩個變量最多為三維，一維、二維和三維的大小配置可以適應向量、矩陣和張量等不同的場景。

二維block size

一個二維的執行配置如上圖所示，其中，每個block有(3 * 4)個Thread，每個grid有(2 * 3)個Block。二維塊大小為 (Dx, Dy)，某個線程號 (x, y) 的公式為 (x + y Dx)；三維塊大小為 (Dx, Dy, Dz)，某個線程號(x, y, z) 的公式為 (x + y Dx + z Dx Dy)。各個內置變量中.x .y和.z為不同維度下的值。

例如，一個二維配置，某個線程在矩陣中的位置可以表示為：

col = cuda.threadIdx.y + cuda.blockDim.y * cuda.blockIdx.y

如何將二維Block映射到自己的數據上並沒有固定的映射方法，一般情況將.x映射為矩陣的行，將.y映射為矩陣的列。Numba提供了一個更簡單的方法幫我們計算線程的編號：

row, col = cuda.grid(2)

其中，參數2表示這是一個2維的執行配置。1維或3維的時候，可以將參數改為1或3。

對應的執行配置也要改為二維：

threads_per_block = (16, 16)

(16, 16)的二維Block是一個常用的配置，共256個線程。本系列第二篇文章也提到，每個Block的Thread個數最好是128、256或512，這與GPU的硬件架構高度相關。

矩陣運算

一個C = AB的矩陣乘法運算，需要我們把A的某一行與B的某一列的所有元素一一相乘，求和後，將結果存儲到結果矩陣C的(row, col)上。在這種實現中，每個線程都要讀取A的一整行和B的一整列，共計算M行*P列。以計算第row行為例，計算C[row, 0]、C[row, 1]…C[row, p-1]這些點時都需要從Global Memory中把整個第row行讀取一遍。可以算到，A矩陣中的每個點需要被讀 B.width 次，B矩陣中的每個點需要被讀 A.height 次。這樣比較浪費時間。因此，可以將多次訪問的數據放到Shared Memory中，減少重複讀取的次數，並充分利用Shared Memory的延遲低的優勢。

from numba import cuda

Shared Memory

接下來的程序利用了Shared Memory來做矩陣乘法。這個實現中，跟未做優化的版本相同的是，每個Thread計算結果矩陣中的一個元素，不同的是，每個CUDA Block會以一個 BLOCK_SIZE * BLOCK_SIZE 子矩陣為基本的計算單元。具體而言，需要聲明Shared Memory區域，數據第一次會從Global Memory拷貝到Shared Memory上，接下來可多次重複利用Shared Memory上的數據。

使用Shared Memory的矩陣乘法

from numba import cuda, float32

進行Shared Memory優化後，計算部分的耗時減少了近一半：

matmul time :1.4370720386505127

在上面的實現過程中，有些地方也比較容易讓人迷惑。

  1. 聲明Shared Memory。這裡使用了cuda.shared.array(shape,type)，shape為這塊數據的向量維度大小，type為Numba數據類型，例如是int32還是float32。這個函數只能在設備端使用。定義好後，這塊數據可被同一個Block的所有Thread共享。需要注意的是，這塊數據雖然在核函數中定義，但它不是單個Thread的私有數據， 它可被同Block中的所有Thread讀寫。
  2. 數據加載。每個Thread會將A中的一個元素加載到sA中，一個Block的 BLOCK_SIZE x BLOCK_SIZE 個Thread可以把sA填充滿。cuda.syncthreads()會等待Block中所有Thread執行完之後才執行下一步。所以，當執行完這個函數的時候，sA和sB的數據已經拷貝好了。
  3. 數據復用。A中的某個點，只會被讀取 B.width / BLOCK_SIZE 次；B中的某個點，只會被讀 A.height / BLOCK_SIZE 次。for n in range(BLOCK_SIZE)這個循環做子矩陣向量乘法時，可多次復用sA和sB的數據。
  4. 子矩陣的數據匯總。我們以一個 BLOCK_SIZE x BLOCK_SIZE 的子矩陣為單位分別對A從左到右，對B從上到下平移並計算，共循環 A.width / BLOCK_SIZE 次。在某一步平移，會得到子矩陣的點積。for m in range(math.ceil(A.shape[1] / BLOCK_SIZE))這個循環起到了計算A從左到右與B從上到下點積的過程。

總結

一般情況下，我們主要從「增大並行度」和「充分利用內存」兩個方向對CUDA來進行優化。本文針對這兩種方向，分別介紹了多流和共享內存技術。這兩種技術有一定的學習成本，但收益非常大，建議有計算密集型任務的朋友花一些時間了解一下這兩種技術和背景知識。本文展示的CUDA接口均為Python Numba版封裝，其他CUDA優化技巧可能還沒完全被Numba支持。CUDA C/C++的接口更豐富，可優化粒度更細，對於有更複雜需求的朋友，建議使用C/C++進行CUDA編程。

VirMach 便宜 VPS

QNews

QNews