原文鏈接：https://zhuanlan.zhihu.com/p/77307505
Python是當(dāng)前最流行的編程語言俏险，被廣泛應(yīng)用在深度學(xué)習(xí)、金融建模纽甘、科學(xué)和工程計算上充石。作為一門解釋型語言廷臼，它運(yùn)行速度慢也常常被用戶詬病。著名Python發(fā)行商Anaconda公司開發(fā)的Numba庫為程序員提供了Python版CPU和GPU編程工具，速度比原生Python快數(shù)十倍甚至更多彼棍。使用Numba進(jìn)行GPU編程灭忠，你可以享受：

Python簡單易用的語法；

極快的開發(fā)速度座硕；

成倍的硬件加速弛作。

為了既保證Python語言的易用性和開發(fā)速度，又達(dá)到并行加速的目的华匾，本系列主要從Python的角度給大家分享GPU編程方法映琳。關(guān)于Numba的入門可以參考我的另一篇文章。更加令人興奮的是蜘拉，Numba提供了一個GPU模擬器萨西，即使你手頭暫時沒有GPU機(jī)器，也可以先使用這個模擬器來學(xué)習(xí)GPU編程旭旭！

本系列為英偉達(dá)GPU入門介紹的第二篇谎脯，主要介紹CUDA編程的基本流程和核心概念，并使用Python Numba編寫GPU并行程序持寄。為了更好地理解GPU的硬件架構(gòu)穿肄，建議讀者先閱讀我的第一篇文章。

GPU硬件知識和基礎(chǔ)概念：包括CPU與GPU的區(qū)別际看、GPU架構(gòu)咸产、CUDA軟件棧簡介。

GPU編程入門：主要介紹CUDA核函數(shù)仲闽，Thread脑溢、Block和Grid概念，并使用Python Numba進(jìn)行簡單的并行計算赖欣。

GPU編程進(jìn)階：主要介紹一些優(yōu)化方法屑彻。

GPU編程實(shí)踐：使用Python Numba解決復(fù)雜問題。

初識GPU編程

兵馬未動顶吮，糧草先行社牲。在開始GPU編程前，需要明確一些概念悴了，并準(zhǔn)備好相關(guān)工具搏恤。

CUDA是英偉達(dá)提供給開發(fā)者的一個GPU編程框架，程序員可以使用這個框架輕松地編寫并行程序湃交。本系列第一篇文章提到熟空，CPU和主存被稱為主機(jī)（Host），GPU和顯存（顯卡內(nèi)存）被稱為設(shè)備（Device）搞莺，CPU無法直接讀取顯存數(shù)據(jù)息罗，GPU無法直接讀取主存數(shù)據(jù)，主機(jī)與設(shè)備必須通過總線（Bus）相互通信才沧。

GPU和CPU架構(gòu)

在進(jìn)行GPU編程前迈喉，需要先確認(rèn)是否安裝了CUDA工具箱绍刮，可以使用echo $CUDA_HOME檢查CUDA環(huán)境變量，返回值不為空說明已經(jīng)安裝好CUDA挨摸。也可以直接用Anaconda里的conda命令安裝CUDA：

$ conda install cudatoolkit

然后可以使用nvidia-smi命令查看顯卡情況孩革，比如這臺機(jī)器上幾張顯卡，CUDA版本油坝，顯卡上運(yùn)行的進(jìn)程等嫉戚。我這里是一臺32GB顯存版的Telsa V100機(jī)器刨裆。

nvidia-smi命令返回結(jié)果

安裝Numba庫：

$ conda install numba

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

fromnumbaimportcudaprint(cuda.gpus)

如果上述步驟沒有問題澈圈，可以得到結(jié)果：<Managed Device 0>...。如果機(jī)器上沒有GPU或沒安裝好上述包帆啃，會有報錯瞬女。CUDA程序執(zhí)行時會獨(dú)霸一張卡，如果你的機(jī)器上有多張GPU卡努潘，CUDA默認(rèn)會選用0號卡诽偷。如果你與其他人共用這臺機(jī)器，最好協(xié)商好誰在用哪張卡疯坤。一般使用CUDA_VISIBLE_DEVICES這個環(huán)境變量來選擇某張卡报慕。如選擇5號GPU卡運(yùn)行你的程序。

CUDA_VISIBLE_DEVICES='5'python example.py

如果手頭暫時沒有GPU設(shè)備压怠，Numba提供了一個模擬器眠冈，供用戶學(xué)習(xí)和調(diào)試，只需要在命令行里添加一個環(huán)境變量菌瘫。

Mac/Linux:

exportNUMBA_ENABLE_CUDASIM=1

Windows:

SETNUMBA_ENABLE_CUDASIM=1

需要注意的是蜗顽，模擬器只是一個調(diào)試的工具，在模擬器中使用Numba并不能加速程序雨让，有可能速度更慢雇盖，而且在模擬器能夠運(yùn)行的程序，并不能保證一定能在真正的GPU上運(yùn)行栖忠，最終還是要以GPU為準(zhǔn)崔挖。

有了以上的準(zhǔn)備工作，我們就可以開始我們的GPU編程之旅了庵寞！

GPU程序與CPU程序的區(qū)別

一個傳統(tǒng)的CPU程序的執(zhí)行順序如下圖所示：

CPU程序執(zhí)行流程

CPU程序是順序執(zhí)行的虚汛，一般需要：

初始化。

CPU計算皇帮。

得到計算結(jié)果卷哩。

在CUDA編程中，CPU和主存被稱為主機(jī)（Host）属拾，GPU被稱為設(shè)備（Device）将谊。

GPU程序執(zhí)行流程

當(dāng)引入GPU后冷溶，計算流程變?yōu)椋?/p>

初始化，并將必要的數(shù)據(jù)拷貝到GPU設(shè)備的顯存上尊浓。

CPU調(diào)用GPU函數(shù)逞频，啟動GPU多個核心同時進(jìn)行計算。

CPU與GPU異步計算栋齿。

將GPU計算結(jié)果拷貝回主機(jī)端苗胀，得到計算結(jié)果。

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

fromnumbaimportcudadefcpu_print():print("print by cpu.")@cuda.jitdefgpu_print():# GPU核函數(shù)print("print by gpu.")defmain():gpu_print[1,2]()cuda.synchronize()cpu_print()if__name__=="__main__":main()

使用CUDA_VISIBLE_DEVICES='0' python gpu_print.py執(zhí)行這段代碼瓦堵，得到的結(jié)果為：

print by gpu.

print by gpu.

print by cpu.

與傳統(tǒng)的Python CPU代碼不同的是：

使用from numba import cuda引入cuda庫

在GPU函數(shù)上添加@cuda.jit裝飾符基协，表示該函數(shù)是一個在GPU設(shè)備上運(yùn)行的函數(shù)，GPU函數(shù)又被稱為核函數(shù)菇用。

主函數(shù)調(diào)用GPU核函數(shù)時澜驮，需要添加如[1, 2]這樣的執(zhí)行配置，這個配置是在告知GPU以多大的并行粒度同時進(jìn)行計算惋鸥。gpu_print[1, 2]()表示同時開啟2個線程并行地執(zhí)行g(shù)pu_print函數(shù)杂穷，函數(shù)將被并行地執(zhí)行2次。下文會深入探討如何設(shè)置執(zhí)行配置卦绣。

GPU核函數(shù)的啟動方式是異步的：啟動GPU函數(shù)后耐量，CPU不會等待GPU函數(shù)執(zhí)行完畢才執(zhí)行下一行代碼。必要時滤港，需要調(diào)用cuda.synchronize()廊蜒，告知CPU等待GPU執(zhí)行完核函數(shù)后，再進(jìn)行CPU端后續(xù)計算蜗搔。這個過程被稱為同步劲藐，也就是GPU執(zhí)行流程圖中的紅線部分。如果不調(diào)用cuda.synchronize()函數(shù)樟凄，執(zhí)行結(jié)果也將改變聘芜，"print by cpu.將先被打印。雖然GPU函數(shù)在前缝龄，但是程序并沒有等待GPU函數(shù)執(zhí)行完汰现，而是繼續(xù)執(zhí)行后面的cpu_print函數(shù)，由于CPU調(diào)用GPU有一定的延遲叔壤，反而后面的cpu_print先被執(zhí)行瞎饲，因此cpu_print的結(jié)果先被打印了出來。

Thread層次結(jié)構(gòu)

前面的程序中炼绘，核函數(shù)被GPU并行地執(zhí)行了2次嗅战。在進(jìn)行GPU并行編程時需要定義執(zhí)行配置來告知以怎樣的方式去并行計算，比如上面打印的例子中，是并行地執(zhí)行2次驮捍，還是8次疟呐，還是并行地執(zhí)行20萬次，或者2000萬次东且。2000萬的數(shù)字太大启具，遠(yuǎn)遠(yuǎn)多于GPU的核心數(shù)，如何將2000萬次計算合理分配到所有GPU核心上珊泳。解決這些問題就需要弄明白CUDA的Thread層次結(jié)構(gòu)鲁冯。

并行執(zhí)行8次的執(zhí)行配置

CUDA將核函數(shù)所定義的運(yùn)算稱為線程（Thread），多個線程組成一個塊（Block）色查，多個塊組成網(wǎng)格（Grid）薯演。這樣一個grid可以定義成千上萬個線程，也就解決了并行執(zhí)行上萬次操作的問題综慎。例如涣仿，把前面的程序改為并行執(zhí)行8次：可以用2個block勤庐，每個block中有4個thread示惊。原來的代碼可以改為gpu_print[2, 4]()，其中方括號中第一個數(shù)字表示整個grid有多少個block愉镰，方括號中第二個數(shù)字表示一個block有多少個thread米罚。

實(shí)際上，線程（thread）是一個編程上的軟件概念丈探。從硬件來看录择，thread運(yùn)行在一個CUDA核心上，多個thread組成的block運(yùn)行在Streaming Multiprocessor（SM的概念詳見本系列第一篇文章）碗降，多個block組成的grid運(yùn)行在一個GPU顯卡上隘竭。

軟硬件對應(yīng)關(guān)系

CUDA提供了一系列內(nèi)置變量，以記錄thread和block的大小及索引下標(biāo)讼渊。以[2, 4]這樣的配置為例：blockDim.x變量表示block的大小是4动看，即每個block有4個thread，threadIdx.x變量是一個從0到blockDim.x - 1（4-1=3）的索引下標(biāo)爪幻，記錄這是第幾個thread菱皆；gridDim.x變量表示grid的大小是2，即每個grid有2個block挨稿，blockIdx.x變量是一個從0到gridDim.x - 1（2-1=1）的索引下標(biāo)仇轻，記錄這是第幾個block。

CUDA內(nèi)置變量示意圖

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

使用內(nèi)置變量計算某個thread編號

利用上述變量篷店，我們可以把之前的代碼豐富一下：

fromnumbaimportcudadefcpu_print(N):foriinrange(0,N):print(i)@cuda.jitdefgpu_print(N):idx=cuda.threadIdx.x+cuda.blockIdx.x*cuda.blockDim.xif(idx<N):print(idx)defmain():print("gpu print:")gpu_print[2,4](8)cuda.synchronize()print("cpu print:")cpu_print(8)if__name__=="__main__":main()

執(zhí)行結(jié)果為：

gpu print:

0

1

2

3

4

5

6

7

cpu print:

0

1

2

3

4

5

6

7

這里的GPU函數(shù)在每個CUDA thread中打印了當(dāng)前thread的編號，起到了CPU函數(shù)for循環(huán)同樣的作用臭家。因?yàn)閒or循環(huán)中的計算內(nèi)容互相不依賴疲陕，也就是說吭产，某次循環(huán)只是專心做自己的事情，循環(huán)第i次不影響循環(huán)第j次的計算鸭轮，所以這樣互相不依賴的for循環(huán)非常適合放到CUDA thread里做并行計算臣淤。在實(shí)際使用中，我們一般將CPU代碼中互相不依賴的的for循環(huán)適當(dāng)替換成CUDA代碼窃爷。

這份代碼打印了8個數(shù)字邑蒋，核函數(shù)有一個參數(shù)N，N = 8按厘，假如我們只想打印5個數(shù)字呢医吊？當(dāng)前的執(zhí)行配置共2 * 4 = 8個線程，線程數(shù)8與要執(zhí)行的次數(shù)5不匹配逮京，不過我們已經(jīng)在代碼里寫好了if (idx < N)的判斷語句卿堂，判斷會幫我們過濾不需要的計算。我們只需要把N = 5傳遞給gpu_print函數(shù)中就好懒棉，CUDA仍然會啟動8個thread草描，但是大于等于N的thread不進(jìn)行計算。注意策严，當(dāng)線程數(shù)與計算次數(shù)不一致時穗慕，一定要使用這樣的判斷語句，以保證某個線程的計算不會影響其他線程的數(shù)據(jù)妻导。

線程數(shù)與計算次數(shù)不匹配

Block大小設(shè)置

不同的執(zhí)行配置會影響GPU程序的速度逛绵，一般需要多次調(diào)試才能找到較好的執(zhí)行配置，在實(shí)際編程中倔韭，執(zhí)行配置[gridDim, blockDim]應(yīng)參考下面的方法：

block運(yùn)行在SM上术浪，不同硬件架構(gòu)（Turing、Volta寿酌、Pascal...）的CUDA核心數(shù)不同胰苏，一般需要根據(jù)當(dāng)前硬件來設(shè)置block的大小blockDim（執(zhí)行配置中第二個參數(shù)）。一個block中的thread數(shù)最好是32份名、128碟联、256的倍數(shù)。==注意僵腺，限于當(dāng)前硬件的設(shè)計鲤孵，block大小不能超過1024。==

grid的大小gridDim（執(zhí)行配置中第一個參數(shù)）辰如，即一個grid中block的個數(shù)可以由總次數(shù)N除以blockDim普监，并向上取整。

例如，我們想并行啟動1000個thread凯正，可以將blockDim設(shè)置為128毙玻，1000 ÷ 128 = 7.8，向上取整為8廊散。使用時桑滩，執(zhí)行配置可以寫成gpuWork[8, 128]()，CUDA共啟動8 * 128 = 1024個thread允睹，實(shí)際計算時只使用前1000個thread运准，多余的24個thread不進(jìn)行計算。

注意缭受，這幾個變量比較容易混淆胁澳，再次明確一下：blockDim是block中thread的個數(shù)，一個block中的threadIdx最大不超過blockDim米者；gridDim是grid中block的個數(shù)韭畸，一個grid中的blockIdx最大不超過gridDim。

以上討論中蔓搞，block和grid大小均是一維胰丁，實(shí)際編程使用的執(zhí)行配置常常更復(fù)雜，block和grid的大小可以設(shè)置為二維甚至三維败明，如下圖所示隘马。這部分內(nèi)容將在下篇文章中討論太防。

Thread Block Grid

內(nèi)存分配

前文提到妻顶，GPU計算時直接從顯存中讀取數(shù)據(jù)，因此每當(dāng)計算時要將數(shù)據(jù)從主存拷貝到顯存上蜒车，用CUDA的術(shù)語來說就是要把數(shù)據(jù)從主機(jī)端拷貝到設(shè)備端讳嘱。CUDA強(qiáng)大之處在于它能自動將數(shù)據(jù)從主機(jī)和設(shè)備間相互拷貝，不需要程序員在代碼中寫明酿愧。這種方法對編程者來說非常方便沥潭，不必對原有的CPU代碼做大量改動。

我們以一個向量加法為例嬉挡，編寫一個向量加法的核函數(shù)如下：

@cuda.jitdefgpu_add(a,b,result,n):# a, b為輸入向量钝鸽，result為輸出向量# 所有向量都是n維# 得到當(dāng)前thread的索引idx=cuda.threadIdx.x+cuda.blockDim.x*cuda.blockIdx.xifidx<n:result[idx]=a[idx]+b[idx]

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

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

把上述代碼整合起來庞钢，與CPU代碼做對比拔恰，并驗(yàn)證結(jié)果正確性：

fromnumbaimportcudaimportnumpyasnpimportmathfromtimeimporttime@cuda.jitdefgpu_add(a,b,result,n):idx=cuda.threadIdx.x+cuda.blockDim.x*cuda.blockIdx.xifidx<n:result[idx]=a[idx]+b[idx]defmain():n=20000000x=np.arange(n).astype(np.int32)y=2*xgpu_result=np.zeros(n)cpu_result=np.zeros(n)threads_per_block=1024blocks_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()

運(yùn)行結(jié)果，GPU代碼竟然比CPU代碼慢10+倍基括！

gpu vector add time 13.589356184005737

cpu vector add time 1.2823548316955566

result correct

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

向量加法的這個計算比較簡單，CPU的numpy已經(jīng)優(yōu)化到了極致，無法突出GPU的優(yōu)勢河爹，我們要解決實(shí)際問題往往比這個復(fù)雜得多匠璧，當(dāng)解決復(fù)雜問題時，優(yōu)化后的GPU代碼將遠(yuǎn)快于CPU代碼咸这。

這份代碼使用CUDA默認(rèn)的統(tǒng)一內(nèi)存管理機(jī)制夷恍，沒有對數(shù)據(jù)的拷貝做優(yōu)化。CUDA的統(tǒng)一內(nèi)存系統(tǒng)是當(dāng)GPU運(yùn)行到某塊數(shù)據(jù)發(fā)現(xiàn)不在設(shè)備端時媳维，再去主機(jī)端中將數(shù)據(jù)拷貝過來裁厅，當(dāng)執(zhí)行完核函數(shù)后，又將所有的內(nèi)存拷貝回主存侨艾。在上面的代碼中执虹，輸入的兩個向量是只讀的，沒必要再拷貝回主存唠梨。

這份代碼沒有做流水線優(yōu)化袋励。CUDA并非同時計算2千萬個數(shù)據(jù)，一般分批流水線工作：一邊對2000萬中的某批數(shù)據(jù)進(jìn)行計算当叭，一邊將下一批數(shù)據(jù)從主存拷貝過來茬故。計算占用的是CUDA核心，數(shù)據(jù)拷貝占用的是總線蚁鳖，所需資源不同磺芭，互相不存在競爭關(guān)系。這種機(jī)制被稱為流水線醉箕。這部分內(nèi)容將在下篇文章中討論钾腺。

原因2中本該程序員動腦思考的問題交給了CUDA解決，增加了時間開銷讥裤，所以CUDA非常方便的統(tǒng)一內(nèi)存模型缺點(diǎn)是計算速度慢放棒。針對原因2，我們可以繼續(xù)優(yōu)化這個程序己英，告知GPU哪些數(shù)據(jù)需要拷貝到設(shè)備间螟，哪些需要拷貝回主機(jī)。

fromnumbaimportcudaimportnumpyasnpimportmathfromtimeimporttime@cuda.jitdefgpu_add(a,b,result,n):idx=cuda.threadIdx.x+cuda.blockDim.x*cuda.blockIdx.xifidx<n:result[idx]=a[idx]+b[idx]defmain():n=20000000x=np.arange(n).astype(np.int32)y=2*x# 拷貝數(shù)據(jù)到設(shè)備端x_device=cuda.to_device(x)y_device=cuda.to_device(y)# 在顯卡設(shè)備上初始化一塊用于存放GPU計算結(jié)果的空間gpu_result=cuda.device_array(n)cpu_result=np.empty(n)threads_per_block=1024blocks_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()

這段代碼的運(yùn)行結(jié)果為：

gpu vector add time 0.19940638542175293

cpu vector add time 1.132070541381836

result correct!

至此损肛，可以看到GPU速度終于比CPU快了很多厢破。

Numba對Numpy的比較友好，編程中一定要使用Numpy的數(shù)據(jù)類型治拿。用到的比較多的內(nèi)存分配函數(shù)有：

cuda.device_array()： 在設(shè)備上分配一個空向量摩泪，類似于numpy.empty()

cuda.to_device()：將主機(jī)的數(shù)據(jù)拷貝到設(shè)備

ary=np.arange(10)device_ary=cuda.to_device(ary)

cuda.copy_to_host()：將設(shè)備的數(shù)據(jù)拷貝回主機(jī)

host_ary=device_ary.copy_to_host()

總結(jié)

Python Numba庫可以調(diào)用CUDA進(jìn)行GPU編程，CPU端被稱為主機(jī)忍啤，GPU端被稱為設(shè)備加勤，運(yùn)行在GPU上的函數(shù)被稱為核函數(shù)仙辟，調(diào)用核函數(shù)時需要有執(zhí)行配置，以告知CUDA以多大的并行粒度來計算鳄梅。使用GPU編程時要合理地將數(shù)據(jù)在主機(jī)和設(shè)備間互相拷貝叠国。

GPU程序執(zhí)行流程

CUDA編程的基本流程為：

初始化，并將必要的數(shù)據(jù)拷貝到GPU設(shè)備的顯存上戴尸。

使用某個執(zhí)行配置粟焊，以一定的并行粒度調(diào)用CUDA核函數(shù)。

CPU和GPU異步計算孙蒙。

將GPU計算結(jié)果拷貝回主機(jī)项棠。