GPU加速02:超詳細Python Cuda零基礎入門教程，沒有顯示卡也能學！

  • 2019 年 12 月 26 日
  • 筆記

Python是當前最流行的程式語言，被廣泛應用在深度學習、金融建模、科學和工程計算上。作為一門解釋型語言，它運行速度慢也常常被用戶詬病。著名Python發行商Anaconda公司開發的Numba庫為程式設計師提供了Python版CPU和GPU編程工具，速度比原生Python快數十倍甚至更多。使用Numba進行GPU編程，你可以享受：

  1. Python簡單易用的語法；
  2. 極快的開發速度；
  3. 成倍的硬體加速。

為了既保證Python語言的易用性和開發速度，又達到並行加速的目的，本系列主要從Python的角度給大家分享GPU編程方法。關於Numba的入門可以參考我的Numba入門文章。更加令人興奮的是，Numba提供了一個GPU模擬器，即使你手頭暫時沒有GPU機器，也可以先使用這個模擬器來學習GPU編程！

本系列為NVIDIA GPU入門介紹的第二篇，主要介紹CUDA編程的基本流程和核心概念，並使用Python Numba編寫GPU並行程式。為了更好地理解GPU的硬體架構，建議讀者先閱讀我的第一篇文章。

  1. GPU硬體知識和基礎概念：包括CPU與GPU的區別、GPU架構、CUDA軟體棧簡介。
  2. GPU編程入門：主要介紹CUDA核函數，Thread、Block和Grid概念，並使用Python Numba進行簡單的並行計算。
  3. GPU編程進階：主要介紹一些優化方法。
  4. GPU編程實踐：使用Python Numba解決複雜問題。

初識GPU編程

兵馬未動，糧草先行。在開始GPU編程前，需要明確一些概念，並準備好相關工具。

CUDA是NVIDIA 提供給開發者的一個GPU編程框架，程式設計師可以使用這個框架輕鬆地編寫並行程式。本系列第一篇文章提到，CPU和主存被稱為主機（Host），GPU和顯示記憶體（顯示卡記憶體）被稱為設備（Device），CPU無法直接讀取顯示記憶體數據，GPU無法直接讀取主存數據，主機與設備必須通過匯流排（Bus）相互通訊。

GPU和CPU架構

在進行GPU編程前，需要先確認是否安裝了CUDA工具箱，可以使用echo $CUDA_HOME檢查CUDA環境變數，返回值不為空說明已經安裝好CUDA。也可以直接用Anaconda里的conda命令安裝CUDA：

$ conda install cudatoolkit

然後可以使用nvidia-smi命令查看顯示卡情況，比如這台機器上幾張顯示卡，CUDA版本，顯示卡上運行的進程等。我這裡是一台32GB顯示記憶體版的Telsa V100機器。

nvidia-smi命令返回結果

安裝Numba庫： 

$ conda install numba

然後檢查一下CUDA和Numba是否安裝成功：

from numba import cuda  print(cuda.gpus)

如果上述步驟沒有問題，可以得到結果：<Managed Device 0>...。如果機器上沒有GPU或沒安裝好上述包，會有報錯。CUDA程式執行時會獨霸一張卡，如果你的機器上有多張GPU卡，CUDA默認會選用0號卡。如果你與其他人共用這台機器，最好協商好誰在用哪張卡。一般使用CUDA_VISIBLE_DEVICES這個環境變數來選擇某張卡。如選擇5號GPU卡運行你的程式。

CUDA_VISIBLE_DEVICES='5' python example.py

如果手頭暫時沒有GPU設備，Numba提供了一個模擬器，供用戶學習和調試，只需要在命令行里添加一個環境變數。

Mac/Linux:

export NUMBA_ENABLE_CUDASIM=1

Windows:

SET NUMBA_ENABLE_CUDASIM=1

需要注意的是，模擬器只是一個調試的工具，在模擬器中使用Numba並不能加速程式，有可能速度更慢，而且在模擬器能夠運行的程式，並不能保證一定能在真正的GPU上運行，最終還是要以GPU為準。

有了以上的準備工作，我們就可以開始我們的GPU編程之旅了！

GPU程式與CPU程式的區別

一個傳統的CPU程式的執行順序如下圖所示：

CPU程式執行流程

CPU程式是順序執行的，一般需要：

  1. 初始化。
  2. CPU計算。
  3. 得到計算結果。

在CUDA編程中，CPU和主存被稱為主機（Host），GPU被稱為設備（Device）。

GPU程式執行流程

當引入GPU後，計算流程變為：

  1. 初始化，並將必要的數據拷貝到GPU設備的顯示記憶體上。
  2. CPU調用GPU函數，啟動GPU多個核心同時進行計算。
  3. CPU與GPU非同步計算。
  4. 將GPU計算結果拷貝回主機端，得到計算結果。

一個名為gpu_print.py的GPU程式如下所示：

from numba import cuda    def cpu_print():      print("print by cpu.")    @cuda.jit  def gpu_print():      # GPU核函數      print("print by gpu.")    def main():      gpu_print[1, 2]()      cuda.synchronize()      cpu_print()    if __name__ == "__main__":      main()

使用CUDA_VISIBLE_DEVICES='0' python gpu_print.py執行這段程式碼，得到的結果為：

print by gpu.  print by gpu.  print by cpu.

與傳統的Python CPU程式碼不同的是：

  • 使用from numba import cuda引入cuda庫
  • 在GPU函數上添加@cuda.jit裝飾符，表示該函數是一個在GPU設備上運行的函數，GPU函數又被稱為核函數。
  • 主函數調用GPU核函數時，需要添加如[1, 2]這樣的執行配置，這個配置是在告知GPU以多大的並行粒度同時進行計算。gpu_print[1, 2]()表示同時開啟2個執行緒並行地執行gpu_print函數，函數將被並行地執行2次。下文會深入探討如何設置執行配置。
  • GPU核函數的啟動方式是非同步的：啟動GPU函數後，CPU不會等待GPU函數執行完畢才執行下一行程式碼。必要時，需要調用cuda.synchronize()，告知CPU等待GPU執行完核函數後，再進行CPU端後續計算。這個過程被稱為同步，也就是GPU執行流程圖中的紅線部分。如果不調用cuda.synchronize()函數，執行結果也將改變，"print by cpu."將先被列印。雖然GPU函數在前，但是程式並沒有等待GPU函數執行完，而是繼續執行後面的cpu_print函數，由於CPU調用GPU有一定的延遲，反而後面的cpu_print先被執行，因此cpu_print的結果先被列印了出來。

Thread層次結構

前面的程式中，核函數被GPU並行地執行了2次。在進行GPU並行編程時需要定義執行配置來告知以怎樣的方式去並行計算，比如上面列印的例子中，是並行地執行2次，還是8次，還是並行地執行20萬次，或者2000萬次。2000萬的數字太大，遠遠多於GPU的核心數，如何將2000萬次計算合理分配到所有GPU核心上。解決這些問題就需要弄明白CUDA的Thread層次結構。

並行執行8次的執行配置

CUDA將核函數所定義的運算稱為執行緒（Thread），多個執行緒組成一個塊（Block），多個塊組成網格（Grid）。這樣一個grid可以定義成千上萬個執行緒，也就解決了並行執行上萬次操作的問題。例如，把前面的程式改為並行執行8次：可以用2個block，每個block中有4個thread。原來的程式碼可以改為gpu_print[2, 4]()，其中方括弧中第一個數字表示整個grid有多少個block，方括弧中第二個數字表示一個block有多少個thread。

實際上，執行緒（thread）是一個編程上的軟體概念。從硬體來看，thread運行在一個CUDA核心上，多個thread組成的block運行在Streaming Multiprocessor（SM的概念詳見本系列第一篇文章），多個block組成的grid運行在一個GPU顯示卡上。

軟硬體對應關係

CUDA提供了一系列內置變數，以記錄thread和block的大小及索引下標。以[2, 4]這樣的配置為例：blockDim.x變數表示block的大小是4，即每個block有4個thread，threadIdx.x變數是一個從0到blockDim.x - 1（4-1=3）的索引下標，記錄這是第幾個thread；gridDim.x變數表示grid的大小是2，即每個grid有2個block，blockIdx.x變數是一個從0到gridDim.x - 1（2-1=1）的索引下標，記錄這是第幾個block。

CUDA內置變數示意圖

某個thread在整個grid中的位置編號為：threadIdx.x + blockIdx.x * blockDim.x。

使用內置變數計算某個thread編號

利用上述變數，我們可以把之前的程式碼豐富一下：

from numba import cuda    def cpu_print(N):      for i in range(0, N):          print(i)    @cuda.jit  def gpu_print(N):      idx = cuda.threadIdx.x + cuda.blockIdx.x * cuda.blockDim.x      if (idx < N):          print(idx)    def main():      print("gpu print:")      gpu_print[2, 4](8)      cuda.synchronize()      print("cpu print:")      cpu_print(8)    if __name__ == "__main__":      main()

執行結果為： 

gpu print:  0  1  2  3  4  5  6  7  cpu print:  0  1  2  3  4  5  6  7

這裡的GPU函數在每個CUDA thread中列印了當前thread的編號，起到了CPU函數for循環同樣的作用。因為for循環中的計算內容互相不依賴，也就是說，某次循環只是專心做自己的事情，循環第i次不影響循環第j次的計算，所以這樣互相不依賴的for循環非常適合放到CUDA thread里做並行計算。在實際使用中，我們一般將CPU程式碼中互相不依賴的的for循環適當替換成CUDA程式碼。

這份程式碼列印了8個數字，核函數有一個參數N，N = 8，假如我們只想列印5個數字呢？當前的執行配置共2 * 4 = 8個執行緒，執行緒數8與要執行的次數5不匹配，不過我們已經在程式碼里寫好了if (idx < N)的判斷語句，判斷會幫我們過濾不需要的計算。我們只需要把N = 5傳遞給gpu_print函數中就好，CUDA仍然會啟動8個thread，但是大於等於N的thread不進行計算。注意，當執行緒數與計算次數不一致時，一定要使用這樣的判斷語句，以保證某個執行緒的計算不會影響其他執行緒的數據。

執行緒數與計算次數不匹配

Block大小設置

不同的執行配置會影響GPU程式的速度，一般需要多次調試才能找到較好的執行配置，在實際編程中，執行配置[gridDim, blockDim]應參考下面的方法：

  • block運行在SM上，不同硬體架構（Turing、Volta、Pascal…）的CUDA核心數不同，一般需要根據當前硬體來設置block的大小blockDim（執行配置中第二個參數）。一個block中的thread數最好是32、128、256的倍數。注意，限於當前硬體的設計，block大小不能超過1024。
  • grid的大小gridDim（執行配置中第一個參數），即一個grid中block的個數可以由總次數N除以blockDim，並向上取整。

例如，我們想並行啟動1000個thread，可以將blockDim設置為128，1000 ÷ 128 = 7.8，向上取整為8。使用時，執行配置可以寫成gpuWork[8, 128]()，CUDA共啟動8 * 128 = 1024個thread，實際計算時只使用前1000個thread，多餘的24個thread不進行計算。

注意，這幾個變數比較容易混淆，再次明確一下：blockDim是block中thread的個數，一個block中的threadIdx最大不超過blockDim；gridDim是grid中block的個數，一個grid中的blockIdx最大不超過gridDim。

以上討論中，block和grid大小均是一維，實際編程使用的執行配置常常更複雜，block和grid的大小可以設置為二維甚至三維，如下圖所示。這部分內容將在下篇文章中討論。

記憶體分配

前文提到，GPU計算時直接從顯示記憶體中讀取數據，因此每當計算時要將數據從主存拷貝到顯示記憶體上，用CUDA的術語來說就是要把數據從主機端拷貝到設備端。CUDA強大之處在於它能自動將數據從主機和設備間相互拷貝，不需要程式設計師在程式碼中寫明。這種方法對編程者來說非常方便，不必對原有的CPU程式碼做大量改動。

我們以一個向量加法為例，編寫一個向量加法的核函數如下：

@cuda.jit  def gpu_add(a, b, result, n):      # a, b為輸入向量，result為輸出向量      # 所有向量都是n維      # 得到當前thread的索引      idx = cuda.threadIdx.x + cuda.blockDim.x * cuda.blockIdx.x      if idx < n:          result[idx] = a[idx] + b[idx]

初始化兩個2千萬維的向量，作為參數傳遞給核函數：

n = 20000000  x = np.arange(n).astype(np.int32)  y = 2 * x  gpu_result = np.zeros(n)    # CUDA執行配置  threads_per_block = 1024  blocks_per_grid = math.ceil(n / threads_per_block)    gpu_add[blocks_per_grid, threads_per_block](x, y, gpu_result, n)

把上述程式碼整合起來，與CPU程式碼做對比，並驗證結果正確性：

from numba import cuda  import numpy as np  import math  from time import time    @cuda.jit  def gpu_add(a, b, result, n):      idx = cuda.threadIdx.x + cuda.blockDim.x * cuda.blockIdx.x      if idx < n:          result[idx] = a[idx] + b[idx]    def main():      n = 20000000      x = np.arange(n).astype(np.int32)      y = 2 * x        gpu_result = np.zeros(n)      cpu_result = np.zeros(n)        threads_per_block = 1024      blocks_per_grid = math.ceil(n / threads_per_block)      start = time()      gpu_add[blocks_per_grid, threads_per_block](x, y, gpu_result, n)      cuda.synchronize()      print("gpu vector add time " + str(time() - start))      start = time()      cpu_result = np.add(x, y)      print("cpu vector add time " + str(time() - start))        if (np.array_equal(cpu_result, gpu_result)):          print("result correct")    if __name__ == "__main__":      main()

運行結果，GPU程式碼竟然比CPU程式碼慢10+倍！ 

gpu vector add time 13.589356184005737  cpu vector add time 1.2823548316955566  result correct

說好的GPU比CPU快幾十倍上百倍呢？這裡GPU比CPU慢很多原因主要在於：

  1. 向量加法的這個計算比較簡單，CPU的numpy已經優化到了極致，無法突出GPU的優勢，我們要解決實際問題往往比這個複雜得多，當解決複雜問題時，優化後的GPU程式碼將遠快於CPU程式碼。
  2. 這份程式碼使用CUDA默認的統一記憶體管理機制，沒有對數據的拷貝做優化。CUDA的統一記憶體系統是當GPU運行到某塊數據發現不在設備端時，再去主機端中將數據拷貝過來，當執行完核函數後，又將所有的記憶體拷貝回主存。在上面的程式碼中，輸入的兩個向量是只讀的，沒必要再拷貝回主存。
  3. 這份程式碼沒有做流水線優化。CUDA並非同時計算2千萬個數據，一般分批流水線工作：一邊對2000萬中的某批數據進行計算，一邊將下一批數據從主存拷貝過來。計算佔用的是CUDA核心，數據拷貝佔用的是匯流排，所需資源不同，互相不存在競爭關係。這種機制被稱為流水線。這部分內容將在下篇文章中討論。

原因2中本該程式設計師動腦思考的問題交給了CUDA解決，增加了時間開銷，所以CUDA非常方便的統一記憶體模型缺點是計算速度慢。針對原因2，我們可以繼續優化這個程式，告知GPU哪些數據需要拷貝到設備，哪些需要拷貝回主機。

from numba import cuda  import numpy as np  import math  from time import time    @cuda.jit  def gpu_add(a, b, result, n):      idx = cuda.threadIdx.x + cuda.blockDim.x * cuda.blockIdx.x      if idx < n :          result[idx] = a[idx] + b[idx]    def main():      n = 20000000      x = np.arange(n).astype(np.int32)      y = 2 * x        # 拷貝數據到設備端      x_device = cuda.to_device(x)      y_device = cuda.to_device(y)      # 在顯示卡設備上初始化一塊用於存放GPU計算結果的空間      gpu_result = cuda.device_array(n)      cpu_result = np.empty(n)        threads_per_block = 1024      blocks_per_grid = math.ceil(n / threads_per_block)      start = time()      gpu_add[blocks_per_grid, threads_per_block](x_device, y_device, gpu_result, n)      cuda.synchronize()      print("gpu vector add time " + str(time() - start))      start = time()      cpu_result = np.add(x, y)      print("cpu vector add time " + str(time() - start))        if (np.array_equal(cpu_result, gpu_result.copy_to_host())):          print("result correct!")    if __name__ == "__main__":      main()

這段程式碼的運行結果為： 

gpu vector add time 0.19940638542175293  cpu vector add time 1.132070541381836  result correct!

至此，可以看到GPU速度終於比CPU快了很多。

Numba對Numpy的比較友好，編程中一定要使用Numpy的數據類型。用到的比較多的記憶體分配函數有：

  • cuda.device_array()：在設備上分配一個空向量，類似於numpy.empty()
  • cuda.to_device()：將主機的數據拷貝到設備
ary = np.arange(10)  device_ary = cuda.to_device(ary)
  • cuda.copy_to_host()：將設備的數據拷貝回主機
host_ary = device_ary.copy_to_host()

總結

Python Numba庫可以調用CUDA進行GPU編程，CPU端被稱為主機，GPU端被稱為設備，運行在GPU上的函數被稱為核函數，調用核函數時需要有執行配置，以告知CUDA以多大的並行粒度來計算。使用GPU編程時要合理地將數據在主機和設備間互相拷貝。

GPU程式執行流程

CUDA編程的基本流程為：

  1. 初始化，並將必要的數據拷貝到GPU設備的顯示記憶體上。
  2. 使用某個執行配置，以一定的並行粒度調用CUDA核函數。
  3. CPU和GPU非同步計算。
  4. 將GPU計算結果拷貝回主機。

VirMach 便宜 VPS

QNews

QNews