快來操縱你的 GPU:CUDA 編程入門極簡教程
前言
2006年,NVIDIA公司發布了CUDA(http://docs.nvidia.com/cuda/),CUDA是建立在NVIDIA的CPUs上的一個通用並行計算平台和編程模型,基於CUDA編程可以利用GPUs的並行計算引擎來更加高效地解決比較複雜的計算難題。近年來,GPU最成功的一個應用就是深度學習領域,基於GPU的並行計算已經成為訓練深度學習模型的標配。目前,最新的CUDA版本為CUDA 9。
GPU並不是一個獨立運行的計算平台,而需要與CPU協同工作,可以看成是CPU的協處理器,因此當我們在說GPU並行計算時,其實是指的基於CPU+GPU的異構計算架構。在異構計算架構中,GPU與CPU通過PCIe匯流排連接在一起來協同工作,CPU所在位置稱為為主機端(host),而GPU所在位置稱為設備端(device),如下圖所示。
基於CPU+GPU的異構計算. 來源:Preofessional CUDA C Programming
可以看到GPU包括更多的運算核心,其特別適合數據並行的計算密集型任務,如大型矩陣運算,而CPU的運算核心較少,但是其可以實現複雜的邏輯運算,因此其適合控制密集型任務。另外,CPU上的線程是重量級的,上下文切換開銷大,但是GPU由於存在很多核心,其線程是輕量級的。因此,基於CPU+GPU的異構計算平台可以優勢互補,CPU負責處理邏輯複雜的串列程序,而GPU重點處理數據密集型的並行計算程序,從而發揮最大功效。
基於CPU+GPU的異構計算應用執行邏輯. 來源:Preofessional CUDA C Programming
CUDA是NVIDIA公司所開發的GPU編程模型,它提供了GPU編程的簡易介面,基於CUDA編程可以構建基於GPU計算的應用程序。CUDA提供了對其它編程語言的支持,如C/C++,Python,Fortran等語言,這裡我們選擇CUDA C/C++介面對CUDA編程進行講解。開發平台為Windows 10 + VS 2013,Windows系統下的CUDA安裝教程可以參考這裡 http://docs.nvidia.com/cuda/cuda-installation-guide-microsoft-windows/index.html。
CUDA編程模型基礎
在給出CUDA的編程實例之前,這裡先對CUDA編程模型中的一些概念及基礎知識做個簡單介紹。CUDA編程模型是一個異構模型,需要CPU和GPU協同工作。在CUDA中,host和device是兩個重要的概念,我們用host指代CPU及其內存,而用device指代GPU及其內存。CUDA程序中既包含host程序,又包含device程序,它們分別在CPU和GPU上運行。同時,host與device之間可以進行通信,這樣它們之間可以進行數據拷貝。典型的CUDA程序的執行流程如下:
分配host內存,並進行數據初始化;
分配device內存,並從host將數據拷貝到device上;
調用CUDA的核函數在device上完成指定的運算;
將device上的運算結果拷貝到host上;
釋放device和host上分配的內存。
上面流程中最重要的一個過程是調用CUDA的核函數來執行並行計算,kernel(http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#kernels)是CUDA中一個重要的概念,kernel是在device上線程中並行執行的函數,核函數用__global__符號聲明,在調用時需要用>>來指定kernel要執行的線程數量,在CUDA中,每一個線程都要執行核函數,並且每個線程會分配一個唯一的線程號thread ID,這個ID值可以通過核函數的內置變數threadIdx來獲得。
由於GPU實際上是異構模型,所以需要區分host和device上的代碼,在CUDA中是通過函數類型限定詞開區別host和device上的函數,主要的三個函數類型限定詞如下:
__global__:在device上執行,從host中調用(一些特定的GPU也可以從device上調用),返回類型必須是void,不支持可變參數參數,不能成為類成員函數。注意用__global__定義的kernel是非同步的,這意味著host不會等待kernel執行完就執行下一步。
__device__:在device上執行,僅可以從device中調用,不可以和__global__同時用。
__host__:在host上執行,僅可以從host上調用,一般省略不寫,不可以和__global__同時用,但可和__device__,此時函數會在device和host都編譯。
要深刻理解kernel,必須要對kernel的線程層次結構有一個清晰的認識。首先GPU上很多並行化的輕量級線程。kernel在device上執行時實際上是啟動很多線程,一個kernel所啟動的所有線程稱為一個網格(grid),同一個網格上的線程共享相同的全局內存空間,grid是線程結構的第一層次,而網格又可以分為很多線程塊(block),一個線程塊裡面包含很多線程,這是第二個層次。線程兩層組織結構如下圖所示,這是一個gird和block均為2-dim的線程組織。grid和block都是定義為dim3類型的變數,dim3可以看成是包含三個無符號整數(x,y,z)成員的結構體變數,在定義時,預設值初始化為1。因此grid和block可以靈活地定義為1-dim,2-dim以及3-dim結構,對於圖中結構(主要水平方向為x軸),定義的grid和block如下所示,kernel在調用時也必須通過執行配置(http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#execution-configuration)>>來指定kernel所使用的線程數及結構。
dim3 grid(3, 2);
dim3 block(4, 3);
kernel_fun>>(prams...);
所以,一個線程需要兩個內置的坐標變數(blockIdx,threadIdx)來唯一標識,它們都是dim3類型變數,其中blockIdx指明線程所在grid中的位置,而threaIdx指明線程所在block中的位置,如圖中的Thread (1,1)滿足:
threadIdx.x = 1
threadIdx.y = 1
blockIdx.x = 1
blockIdx.y = 1
一個線程塊上的線程是放在同一個流式多處理器(SM)上的,但是單個SM的資源有限,這導致線程塊中的線程數是有限制的,現代GPUs的線程塊可支持的線程數可達1024個。有時候,我們要知道一個線程在blcok中的全局ID,此時就必須還要知道block的組織結構,這是通過線程的內置變數blockDim來獲得。它獲取線程塊各個維度的大小。對於一個2-dim的block,線程(x,y)的ID值為,如果是3-dim的block,線程(x,y,z)的ID值為。另外線程還有內置變數gridDim,用於獲得網格塊各個維度的大小。
kernel的這種線程組織結構天然適合vector,matrix等運算,如我們將利用上圖2-dim結構實現兩個矩陣的加法,每個線程負責處理每個位置的兩個元素相加,代碼如下所示。線程塊大小為(16, 16),然後將N*N大小的矩陣均分為不同的線程塊來執行加法運算。
// Kernel定義
__global__void MatAdd(float A[N][N], float B[N][N], float C[N][N])
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
if (i
C[i][j] = A[i][j] + B[i][j];
}
int main()
{
...
// Kernel 線程配置
dim3 threadsPerBlock(16, 16);
dim3 numBlocks(N / threadsPerBlock.x, N / threadsPerBlock.y);
// kernel調用
MatAdd>>(A, B, C);
...
}
此外這裡簡單介紹一下CUDA的內存模型,如下圖所示。可以看到,每個線程有自己的私有本地內存(Local Memory),而每個線程塊有包含共享內存(Shared Memory),可以被線程塊中所有線程共享,其生命周期與線程塊一致。此外,所有的線程都可以訪問全局內存(Global Memory)。還可以訪問一些只讀內存塊:常量內存(Constant Memory)和紋理內存(Texture Memory)。內存結構涉及到程序優化,這裡不深入探討它們。
CUDA內存模型
還有重要一點,你需要對GPU的硬體實現有一個基本的認識。上面說到了kernel的線程組織層次,那麼一個kernel實際上會啟動很多線程,這些線程是邏輯上並行的,但是在物理層卻並不一定。這其實和CPU的多線程有類似之處,多線程如果沒有多核支持,在物理層也是無法實現並行的。但是好在GPU存在很多CUDA核心,充分利用CUDA核心可以充分發揮GPU的並行計算能力。GPU硬體的一個核心組件是SM,前面已經說過,SM是英文名是 Streaming Multiprocessor,翻譯過來就是流式多處理器。SM的核心組件包括CUDA核心,共享內存,寄存器等,SM可以並發地執行數百個線程,並發能力就取決於SM所擁有的資源數。當一個kernel被執行時,它的gird中的線程塊被分配到SM上,一個線程塊只能在一個SM上被調度。SM一般可以調度多個線程塊,這要看SM本身的能力。那麼有可能一個kernel的各個線程塊被分配多個SM,所以grid只是邏輯層,而SM才是執行的物理層。SM採用的是SIMT(鏈接:http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#simt-architecture)(Single-Instruction, Multiple-Thread,單指令多線程)架構,基本的執行單元是線程束(wraps),線程束包含32個線程,這些線程同時執行相同的指令,但是每個線程都包含自己的指令地址計數器和寄存器狀態,也有自己獨立的執行路徑。所以儘管線程束中的線程同時從同一程序地址執行,但是可能具有不同的行為,比如遇到了分支結構,一些線程可能進入這個分支,但是另外一些有可能不執行,它們只能死等,因為GPU規定線程束中所有線程在同一周期執行相同的指令,線程束分化會導致性能下降。當線程塊被劃分到某個SM上時,它將進一步劃分為多個線程束,因為這才是SM的基本執行單元,但是一個SM同時並發的線程束數是有限的。這是因為資源限制,SM要為每個線程塊分配共享內存,而也要為每個線程束中的線程分配獨立的寄存器。所以SM的配置會影響其所支持的線程塊和線程束並發數量。總之,就是網格和線程塊只是邏輯劃分,一個kernel的所有線程其實在物理層是不一定同時並發的。所以kernel的grid和block的配置不同,性能會出現差異,這點是要特別注意的。還有,由於SM的基本執行單元是包含32個線程的線程束,所以block大小一般要設置為32的倍數。
CUDA編程的邏輯層和物理層
在進行CUDA編程前,可以先檢查一下自己的GPU的硬體配置,這樣才可以有的放矢,可以通過下面的程序獲得GPU的配置屬性:
intdev =;
cudaDeviceProp devProp;
CHECK(cudaGetDeviceProperties(&devProp, dev));
std::cout
std::cout
std::cout
std::cout
std::cout
std::cout
// 輸出如下
使用GPU device: GeForce GT730
SM的數量:2
每個線程塊的共享內存大小:48KB
每個線程塊的最大線程數:1024
每個EM的最大線程數:2048
每個EM的最大線程束數:64
好吧,GT 730顯卡確實有點渣,只有2個SM,嗚嗚......
向量加法實例
知道了CUDA編程基礎,我們就來個簡單的實戰,利用CUDA編程實現兩個向量的加法,在實現之前,先簡單介紹一下CUDA編程中內存管理API。首先是在device上分配內存的cudaMalloc函數:
cudaError_tcudaMalloc(void** devPtr,size_tsize);
這個函數和C語言中的malloc類似,但是在device上申請一定位元組大小的顯存,其中devPtr是指向所分配內存的指針。同時要釋放分配的內存使用cudaFree函數,這和C語言中的free函數對應。另外一個重要的函數是負責host和device之間數據通信的cudaMemcpy函數:
cudaError_tcudaMalloc(void** devPtr,size_tsize);
其中src指向數據源,而dst是目標區域,count是複製的位元組數,其中kind控制複製的方向:cudaMemcpyHostToHost, cudaMemcpyHostToDevice,
cudaMemcpyDeviceToHost及cudaMemcpyDeviceToDevice,如cudaMemcpyHostToDevice將host上數據拷貝到device上。
現在我們來實現一個向量加法的實例,這裡grid和block都設計為1-dim,首先定義kernel如下:
// 兩個向量加法kernel,grid和block均為一維
__global__voidadd(float* x,float* y,float* z,intn)
{
// 獲取全局索引
intindex = threadIdx.x + blockIdx.x * blockDim.x;
// 步長
intstride = blockDim.x * gridDim.x;
for(inti = index; i
{
z[i] = x[i] + y[i];
}
}
其中stride是整個grid的線程數,有時候向量的元素數很多,這時候可以將在每個線程實現多個元素(元素總數/線程總數)的加法,相當於使用了多個grid來處理,這是一種grid-stride loop(鏈接:https://devblogs.nvidia.com/cuda-pro-tip-write-flexible-kernels-grid-stride-loops/)方式,不過下面的例子一個線程只處理一個元素,所以kernel裡面的循環是不執行的。下面我們具體實現向量加法:
intmain()
{
intN =1
intnBytes = N * sizeof(float);
// 申請host內存
float*x, *y, *z;
x = (float*)malloc(nBytes);
y = (float*)malloc(nBytes);
z = (float*)malloc(nBytes);
// 初始化數據
for(inti =; i
{
x[i] =10.0;
y[i] =20.0;
}
// 申請device內存
float*d_x, *d_y, *d_z;
cudaMalloc((void**)&d_x, nBytes);
cudaMalloc((void**)&d_y, nBytes);
cudaMalloc((void**)&d_z, nBytes);
// 將host數據拷貝到device
cudaMemcpy((void*)d_x, (void*)x, nBytes, cudaMemcpyHostToDevice);
cudaMemcpy((void*)d_y, (void*)y, nBytes, cudaMemcpyHostToDevice);
// 定義kernel的執行配置
dim3blockSize(256);
dim3gridSize((N + blockSize.x -1)/ blockSize.x);
// 執行kernel
add > >(d_x, d_y, d_z, N);
// 將device得到的結果拷貝到host
cudaMemcpy((void*)z, (void*)d_z, nBytes, cudaMemcpyHostToDevice);
// 檢查執行結果
floatmaxError =0.0;
for(inti =; i
maxError = fmax(maxError, fabs(z[i] -30.0));
std::cout
// 釋放device內存
cudaFree(d_x);
cudaFree(d_y);
cudaFree(d_z);
// 釋放host內存
free(x);
free(y);
free(z);
return;
}
這裡我們的向量大小為1
kernel的線程層級結構如下圖所示:
kernel的線程層次結構. 來源:https://devblogs.nvidia.com/even-easier-introduction-cuda/
使用nvprof工具可以分析kernel運行情況,結果如下所示,可以看到kernel函數費時約1.5ms。
nvprof cuda9.exe
==7244== NVPROFisprofiling process7244,command: cuda9.exe
最大誤差:4.31602e+008
==7244== Profiling application: cuda9.exe
==7244== Profiling result:
Type Time(%) Time Calls Avg Min Max Name
GPU activities:67.57%3.2256ms21.6128ms1.6017ms1.6239ms [CUDA memcpy HtoD]
32.43%1.5478ms11.5478ms1.5478ms1.5478msadd(float*, float*, float*,int)
你調整block的大小,對比不同配置下的kernel運行情況,我這裡測試的是當block為128時,kernel費時約1.6ms,而block為512時kernel費時約1.7ms,當block為64時,kernel費時約2.3ms。看來不是block越大越好,而要適當選擇。
在上面的實現中,我們需要單獨在host和device上進行內存分配,並且要進行數據拷貝,這是很容易出錯的。好在CUDA 6.0引入統一內存(Unified Memory)(鏈接:http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#um-unified-memory-programming-hd)來避免這種麻煩,簡單來說就是統一內存使用一個託管內存來共同管理host和device中的內存,並且自動在host和device中進行數據傳輸。CUDA中使用cudaMallocManaged函數分配託管內存:
cudaError_tcudaMallocManaged(void**devPtr,size_tsize,unsignedintflag=);
利用統一內存,可以將上面的程序簡化如下:
intmain()
{
intN =1
intnBytes = N * sizeof(float);
// 申請託管內存
float*x, *y, *z;
cudaMallocManaged((void**)&x, nBytes);
cudaMallocManaged((void**)&y, nBytes);
cudaMallocManaged((void**)&z, nBytes);
// 初始化數據
for(inti =; i
{
x[i] =10.0;
y[i] =20.0;
}
// 定義kernel的執行配置
dim3blockSize(256);
dim3gridSize((N + blockSize.x -1)/ blockSize.x);
// 執行kernel
add > >(x, y, z, N);
// 同步device 保證結果能正確訪問
cudaDeviceSynchronize();
// 檢查執行結果
floatmaxError =0.0;
for(inti =; i
maxError = fmax(maxError, fabs(z[i] -30.0));
std::cout
// 釋放內存
cudaFree(x);
cudaFree(y);
cudaFree(z);
return;
}
相比之前的代碼,使用統一內存更簡潔了,值得注意的是kernel執行是與host非同步的,由於託管內存自動進行數據傳輸,這裡要用cudaDeviceSynchronize()函數保證device和host同步,這樣後面才可以正確訪問kernel計算的結果。
矩陣乘法實例
最後我們再實現一個稍微複雜一些的例子,就是兩個矩陣的乘法,設輸入矩陣為A和B,要得到C=A*B。實現思路是每個線程計算C的一個元素值,對於矩陣運算,應該選用grid和block為2-D的。首先定義矩陣的結構體:
// 矩陣類型,行優先,M(row, col) = *(M.elements + row * M.width + col)
structMatrix
{
intwidth;
intheight;
float*elements;
Matrix(intw,inth,float* e =NULL)
{
width = w;
height = h;
elements = e;
}
};
矩陣乘法實現模式
然後實現矩陣乘法的核函數,這裡我們定義了兩個輔助的__device__函數分別用於獲取矩陣的元素值和為矩陣元素賦值,具體代碼如下:
// 獲取矩陣A的(row, col)元素
__device__floatgetElement(constMatrix A,introw,intcol)
{
returnA.elements[row * A.width + col];
}
// 為矩陣A的(row, col)元素賦值
__device__voidsetElement(Matrix A,introw,intcol,floatvalue)
{
A.elements[row * A.width + col] =value;
}
// 矩陣相乘kernel,2-D,每個線程計算一個元素
__global__voidmatMulKernel(constMatrix A,constMatrix B, Matrix C)
{
floatCvalue =0.0;
introw = threadIdx.y + blockIdx.y * blockDim.y;
intcol = threadIdx.x + blockIdx.x * blockDim.x;
for(inti =; i
{
Cvalue += getElement(A, row, i) * getElement(B, i, col);
}
setElement(C, row, col, Cvalue);
}
最後我們採用統一內存編寫矩陣相乘的測試實例:
intmain()
{
intwidth =1
intheight =1
MatrixA(width, height, NULL);
MatrixB(width, height, NULL);
MatrixC(width, height, NULL);
intnBytes = width * height * sizeof(float);
// 申請託管內存
cudaMallocManaged((void**)&A.elements, nBytes);
cudaMallocManaged((void**)&B.elements, nBytes);
cudaMallocManaged((void**)&C.elements, nBytes);
// 初始化數據
for(inti =; i
{
A.elements[i] =1.0;
B.elements[i] =2.0;
}
// 定義kernel的執行配置
dim3blockSize(32,32);
dim3gridSize((width + blockSize.x -1)/ blockSize.x,
(height + blockSize.y -1)/ blockSize.y);
// 執行kernel
matMulKernel > >(A, B, C);
// 同步device 保證結果能正確訪問
cudaDeviceSynchronize();
// 檢查執行結果
floatmaxError =0.0;
for(inti =; i
maxError = fmax(maxError, fabs(C.elements[i] -2* width));
std::cout
std::cout
return;
}
這裡矩陣大小為1024*1024,設計的線程的block大小為(32, 32),那麼grid大小為(32, 32),最終測試結果如下:
nvprof cuda9.exe
==2456== NVPROFisprofiling process2456,command: cuda9.exe
最大誤差:
==2456== Profiling application: cuda9.exe
==2456== Profiling result:
Type Time(%) Time Calls Avg Min Max Name
GPU activities:100.00%2.67533s12.67533s2.67533s2.67533s matMulKernel(Matrix, Matrix, Matrix)
APIcalls:92.22%2.67547s12.67547s2.67547s2.67547s cudaDeviceSynchronize
6.06%175.92ms358.640ms2.3933ms170.97ms cudaMallocManaged
1.65%47.845ms147.845ms47.845ms47.845ms cudaLaunch
0.05%1.4405ms9415.324usns938.54us cuDeviceGetAttribute
0.01%371.49us1371.49us371.49us371.49us cuDeviceGetName
0.00%13.474us113.474us13.474us13.474us cuDeviceTotalMem
0.00%6.9300us16.9300us6.9300us6.9300us cudaConfigureCall
0.00%3.8500us31.2830us385ns1.9250us cuDeviceGetCount
0.00%3.4650us31.1550usns2.3100us cudaSetupArgument
0.00%2.3100us21.1550us385ns1.9250us cuDeviceGet
==2456== Unified Memory profiling result:
Device"GeForce GT 730 (0)"
Count Avg Size Min Size Max Size Total Size Total Time Name
20484.0000KB4.0000KB4.0000KB8.000000MB22.70431ms Host To Device
26646.195KB32.000KB1.0000MB12.00000MB7.213048ms Device To Host
當然,這不是最高效的實現,後面可以繼續優化...
小結
最後只有一句話:CUDA入門容易,但是深入難!希望不是從入門到放棄.
參考資料:
John Cheng, Max Grossman, Ty McKercher. Professional CUDA C Programming, 2014.
(http://www.wrox.com/WileyCDA/WroxTitle/Professional-CUDA-C-Programming.productCd-1118739329.html)
CUDA docs.
(http://docs.nvidia.com/cuda/)
An Even Easier Introduction to CUDA
(https://devblogs.nvidia.com/even-easier-introduction-cuda/)
4 月 AI 求職季
8 大明星企業
10 場分享盛宴
20 小時獨門秘籍
4.10-4.19,我們準時相約!
新人福利
關注 AI 研習社(okweiwu),回復1領取
【超過 1000G 神經網路 / AI / 大數據資料】
用 GPU 加速深度學習: Windows 安裝 CUDA+TensorFlow 教程
![](https://pic.pimg.tw/zzuyanan/1488615166-1259157397.png)
![](https://pic.pimg.tw/zzuyanan/1482887990-2595557020.jpg)
※野外動物監測圖像挑戰賽:預測捕捉到的野外圖像是否包含動物
※用 Hinton 的膠囊神經網路來識別空間關係 Part1:CNNs及其缺點
TAG:AI研習社 |