快來操縱你的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），如下圖所示。
可以看到GPU包括更多的運算核心，其特別適合資料並行的計算密集型任務，如大型矩陣運算，而CPU的運算核心較少，但是其可以實現複雜的邏輯運算，因此其適合控制密集型任務。另外，CPU上的執行緒是重量級的，上下文切換開銷大，但是GPU由於存在很多核心，其執行緒是輕量級的。因此，基於CPU+GPU的異構計算平臺可以優勢互補，CPU負責處理邏輯複雜的序列程式，而GPU重點處理資料密集型的平行計算程式，從而發揮最大功效。
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
1CUDA程式設計模型基礎
在給出CUDA的程式設計例項之前，這裡先對CUDA程式設計模型中的一些概念及基礎知識做個簡單介紹。CUDA程式設計模型是一個異構模型，需要CPU和GPU協同工作。在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__符號宣告，在呼叫時需要用<<<grid, block>>>來指定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）<<<grid, block>>>來指定kernel所使用的執行緒數及結構。
dim3 grid(3, 2);
dim3 block(4, 3);
kernel_fun<<< grid, block >>>(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 < N && j < N)
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<<<numBlocks, threadsPerBlock>>>(A, B, C);
...
}此外這裡簡單介紹一下CUDA的記憶體模型，如下圖所示。可以看到，每個執行緒有自己的私有本地記憶體（Local Memory），而每個執行緒塊有包含共享記憶體（Shared Memory）,可以被執行緒塊中所有執行緒共享，其生命週期與執行緒塊一致。此外，所有的執行緒都可以訪問全域性記憶體（Global Memory）。還可以訪問一些只讀記憶體塊：常量記憶體（Constant Memory）和紋理記憶體（Texture Memory）。記憶體結構涉及到程式優化，這裡不深入探討它們。
還有重要一點，你需要對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程式設計前，可以先檢查一下自己的GPU的硬體配置，這樣才可以有的放矢，可以通過下面的程式獲得GPU的配置屬性：
int dev = 0;
cudaDeviceProp devProp;
CHECK(cudaGetDeviceProperties(&devProp, dev));
std::cout << "使用GPU device " << dev << ": " << devProp.name << std::endl;
std::cout << "SM的數量：" << devProp.multiProcessorCount << std::endl;
std::cout << "每個執行緒塊的共享記憶體大小：" << devProp.sharedMemPerBlock / 1024.0 << " KB" << std::endl;
std::cout << "每個執行緒塊的最大執行緒數：" << devProp.maxThreadsPerBlock << std::endl;
std::cout << "每個EM的最大執行緒數：" << devProp.maxThreadsPerMultiProcessor << std::endl;
std::cout << "每個EM的最大執行緒束數：" << devProp.maxThreadsPerMultiProcessor / 32 << std::endl;
// 輸出如下
使用GPU device 0: GeForce GT 730
SM的數量：2
每個執行緒塊的共享記憶體大小：48 KB
每個執行緒塊的最大執行緒數：1024
每個EM的最大執行緒數：2048
每個EM的最大執行緒束數：64好吧，GT 730顯示卡確實有點渣，只有2個SM，嗚嗚......2
向量加法例項
知道了CUDA程式設計基礎，我們就來個簡單的實戰，利用CUDA程式設計實現兩個向量的加法，在實現之前，先簡單介紹一下CUDA程式設計中記憶體管理API。首先是在device上分配記憶體的cudaMalloc函式：
cudaError_t cudaMalloc(void** devPtr, size_t size);這個函式和C語言中的malloc類似，但是在device上申請一定位元組大小的視訊記憶體，其中devPtr是指向所分配記憶體的指標。同時要釋放分配的記憶體使用cudaFree函式，這和C語言中的free函式對應。另外一個重要的函式是負責host和device之間資料通訊的cudaMemcpy函式：
cudaError_t cudaMalloc(void** devPtr, size_t size);其中src指向資料來源，而dst是目標區域，count是複製的位元組數，其中kind控制複製的方向：cudaMemcpyHostToHost, cudaMemcpyHostToDevice,
cudaMemcpyDeviceToHost及cudaMemcpyDeviceToDevice，如cudaMemcpyHostToDevice將host上資料拷貝到device上。
現在我們來實現一個向量加法的例項，這裡grid和block都設計為1-dim，首先定義kernel如下：
// 兩個向量加法kernel，grid和block均為一維
__global__ void add(float* x, float * y, float* z, int n)
{
// 獲取全域性索引
int index = threadIdx.x + blockIdx.x * blockDim.x;
// 步長
int stride = blockDim.x * gridDim.x;
for (int i = index; i < n; i += stride)
{
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裡面的迴圈是不執行的。下面我們具體實現向量加法：
int main()
{
int N = 1 << 20;
int nBytes = N * sizeof(float);
// 申請host記憶體
float *x, *y, *z;
x = (float*)malloc(nBytes);
y = (float*)malloc(nBytes);
z = (float*)malloc(nBytes);
// 初始化資料
for (int i = 0; i < N; ++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的執行配置
dim3 blockSize(256);
dim3 gridSize((N + blockSize.x - 1) / blockSize.x);
// 執行kernel
add << < gridSize, blockSize >> >(d_x, d_y, d_z, N);
// 將device得到的結果拷貝到host
cudaMemcpy((void*)z, (void*)d_z, nBytes, cudaMemcpyHostToDevice);
// 檢查執行結果
float maxError = 0.0;
for (int i = 0; i < N; i++)
maxError = fmax(maxError, fabs(z[i] - 30.0));
std::cout << "最大誤差: " << maxError << std::endl;
// 釋放device記憶體
cudaFree(d_x);
cudaFree(d_y);
cudaFree(d_z);
// 釋放host記憶體
free(x);
free(y);
free(z);
return 0;
}這裡我們的向量大小為1<<20，而block大小為256，那麼grid大小是4096，
kernel的執行緒層級結構如下圖所示：
使用nvprof工具可以分析kernel執行情況，結果如下所示，可以看到kernel函式費時約1.5ms。
nvprof cuda9.exe
==7244== NVPROF is profiling process 7244, 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.2256ms 2 1.6128ms 1.6017ms 1.6239ms [CUDA memcpy HtoD]
32.43% 1.5478ms 1 1.5478ms 1.5478ms 1.5478ms add(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_t cudaMallocManaged(void **devPtr, size_t size, unsigned int flag=0);利用統一記憶體，可以將上面的程式簡化如下：
int main()
{
int N = 1 << 20;
int nBytes = N * sizeof(float);
// 申請託管記憶體
float *x, *y, *z;
cudaMallocManaged((void**)&x, nBytes);
cudaMallocManaged((void**)&y, nBytes);
cudaMallocManaged((void**)&z, nBytes);
// 初始化資料
for (int i = 0; i < N; ++i)
{
x[i] = 10.0;
y[i] = 20.0;
}
// 定義kernel的執行配置
dim3 blockSize(256);
dim3 gridSize((N + blockSize.x - 1) / blockSize.x);
// 執行kernel
add << < gridSize, blockSize >> >(x, y, z, N);
// 同步device 保證結果能正確訪問
cudaDeviceSynchronize();
// 檢查執行結果
float maxError = 0.0;
for (int i = 0; i < N; i++)
maxError = fmax(maxError, fabs(z[i] - 30.0));
std::cout << "最大誤差: " << maxError << std::endl;
// 釋放記憶體
cudaFree(x);
cudaFree(y);
cudaFree(z);
return 0;
}相比之前的程式碼，使用統一記憶體更簡潔了，值得注意的是kernel執行是與host非同步的，由於託管記憶體自動進行資料傳輸，這裡要用cudaDeviceSynchronize()函式保證device和host同步，這樣後面才可以正確訪問kernel計算的結果。
3
矩陣乘法例項
最後我們再實現一個稍微複雜一些的例子，就是兩個矩陣的乘法，設輸入矩陣為A和B,要得到C=A*B。實現思路是每個執行緒計算C的一個元素值，對於矩陣運算，應該選用grid和block為2-D的。首先定義矩陣的結構體：
// 矩陣型別，行優先，M(row, col) = *(M.elements + row * M.width + col)
struct Matrix
{
int width;
int height;
float *elements;
Matrix(int w, int h, float* e = NULL)
{
width = w;
height = h;
elements = e;
}
};然後實現矩陣乘法的核函式，這裡我們定義了兩個輔助的__device__函式分別用於獲取矩陣的元素值和為矩陣元素賦值，具體程式碼如下：
// 獲取矩陣A的(row, col)元素
__device__ float getElement(const Matrix A, int row, int col)
{
return A.elements[row * A.width + col];
}
// 為矩陣A的(row, col)元素賦值
__device__ void setElement(Matrix A, int row, int col, float value)
{
A.elements[row * A.width + col] = value;
}
// 矩陣相乘kernel，2-D，每個執行緒計算一個元素
__global__ void matMulKernel(const Matrix A, const Matrix B, Matrix C)
{
float Cvalue = 0.0;
int row = threadIdx.y + blockIdx.y * blockDim.y;
int col = threadIdx.x + blockIdx.x * blockDim.x;
for (int i = 0; i < A.width; ++i)
{
Cvalue += getElement(A, row, i) * getElement(B, i, col);
}
setElement(C, row, col, Cvalue);
}最後我們採用統一記憶體編寫矩陣相乘的測試例項：
int main()
{
int width = 1 << 10;
int height = 1 << 10;
Matrix A(width, height, NULL);
Matrix B(width, height, NULL);
Matrix C(width, height, NULL);
int nBytes = width * height * sizeof(float);
// 申請託管記憶體
cudaMallocManaged((void**)&A.elements, nBytes);
cudaMallocManaged((void**)&B.elements, nBytes);
cudaMallocManaged((void**)&C.elements, nBytes);
// 初始化資料
for (int i = 0; i < width * height; ++i)
{
A.elements[i] = 1.0;
B.elements[i] = 2.0;
}
// 定義kernel的執行配置
dim3 blockSize(32, 32);
dim3 gridSize((width + blockSize.x - 1) / blockSize.x,
(height + blockSize.y - 1) / blockSize.y);
// 執行kernel
matMulKernel << < gridSize, blockSize >> >(A, B, C);
// 同步device 保證結果能正確訪問
cudaDeviceSynchronize();
// 檢查執行結果
float maxError = 0.0;
for (int i = 0; i < width * height; i++)
maxError = fmax(maxError, fabs(C.elements[i] -