作者： 葉 虎
編輯：李雪冬
2006年，NVIDIA公司發(fā)布了CUDA（http://docs.nvidia.com/cuda/），CUDA是建立在NVIDIA的CPUs上的一個通用并行計算平臺和編程模型，基于CUDA編程可以利用GPUs的并行計算引擎來更加高效地解決比較復雜的計算難題。近年來，GPU最成功的一個應用就是深度學習領域，基于GPU的并行計算已經(jīng)成為訓練深度學習模型的標配。目前，最新的CUDA版本為CUDA 9。
GPU并不是一個獨立運行的計算平臺，而需要與CPU協(xié)同工作，可以看成是CPU的協(xié)處理器，因此當我們在說GPU并行計算時，其實是指的基于CPU+GPU的異構(gòu)計算架構(gòu)。在異構(gòu)計算架構(gòu)中，GPU與CPU通過PCIe總線連接在一起來協(xié)同工作，CPU所在位置稱為為主機端（host），而GPU所在位置稱為設備端（device），如下圖所示。
基于CPU+GPU的異構(gòu)計算. 來源：Preofessional CUDA? C Programming
可以看到GPU包括更多的運算核心，其特別適合數(shù)據(jù)并行的計算密集型任務，如大型矩陣運算，而CPU的運算核心較少，但是其可以實現(xiàn)復雜的邏輯運算，因此其適合控制密集型任務。另外，CPU上的線程是重量級的，上下文切換開銷大，但是GPU由于存在很多核心，其線程是輕量級的。因此，基于CPU+GPU的異構(gòu)計算平臺可以優(yōu)勢互補，CPU負責處理邏輯復雜的串行程序，而GPU重點處理數(shù)據(jù)密集型的并行計算程序，從而發(fā)揮最大功效。
基于CPU+GPU的異構(gòu)計算應用執(zhí)行邏輯. 來源：Preofessional CUDA? C Programming
CUDA是NVIDIA公司所開發(fā)的GPU編程模型，它提供了GPU編程的簡易接口，基于CUDA編程可以構(gòu)建基于GPU計算的應用程序。CUDA提供了對其它編程語言的支持，如C/C++，Python，F(xiàn)ortran等語言，這里我們選擇CUDA C/C++接口對CUDA編程進行講解。開發(fā)平臺為Windows 10 + VS 2013，Windows系統(tǒng)下的CUDA安裝教程可以參考這里http://docs.nvidia.com/cuda/cuda-installation-guide-microsoft-windows/index.html
在給出CUDA的編程實例之前，這里先對CUDA編程模型中的一些概念及基礎知識做個簡單介紹。CUDA編程模型是一個異構(gòu)模型，需要CPU和GPU協(xié)同工作。在CUDA中，host和device是兩個重要的概念，我們用host指代CPU及其內(nèi)存，而用device指代GPU及其內(nèi)存。CUDA程序中既包含host程序，又包含device程序，它們分別在CPU和GPU上運行。同時，host與device之間可以進行通信，這樣它們之間可以進行數(shù)據(jù)拷貝。典型的CUDA程序的執(zhí)行流程如下：
分配host內(nèi)存，并進行數(shù)據(jù)初始化；
分配device內(nèi)存，并從host將數(shù)據(jù)拷貝到device上；
調(diào)用CUDA的核函數(shù)在device上完成指定的運算；
將device上的運算結(jié)果拷貝到host上；
釋放device和host上分配的內(nèi)存。
上面流程中最重要的一個過程是調(diào)用CUDA的核函數(shù)來執(zhí)行并行計算，kernel（http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#kernels）是CUDA中一個重要的概念，kernel是在device上線程中并行執(zhí)行的函數(shù)，核函數(shù)用__global__符號聲明，在調(diào)用時需要用<<<grid, block>>>來指定kernel要執(zhí)行的線程數(shù)量，在CUDA中，每一個線程都要執(zhí)行核函數(shù)，并且每個線程會分配一個唯一的線程號thread ID，這個ID值可以通過核函數(shù)的內(nèi)置變量threadIdx來獲得。
由于GPU實際上是異構(gòu)模型，所以需要區(qū)分host和device上的代碼，在CUDA中是通過函數(shù)類型限定詞開區(qū)別host和device上的函數(shù)，主要的三個函數(shù)類型限定詞如下：
__global__：在device上執(zhí)行，從host中調(diào)用（一些特定的GPU也可以從device上調(diào)用），返回類型必須是void，不支持可變參數(shù)參數(shù)，不能成為類成員函數(shù)。注意用__global__定義的kernel是異步的，這意味著host不會等待kernel執(zhí)行完就執(zhí)行下一步。
__device__：在device上執(zhí)行，僅可以從device中調(diào)用，不可以和__global__同時用。
__host__：在host上執(zhí)行，僅可以從host上調(diào)用，一般省略不寫，不可以和__global__同時用，但可和__device__，此時函數(shù)會在device和host都編譯。
要深刻理解kernel，必須要對kernel的線程層次結(jié)構(gòu)有一個清晰的認識。首先GPU上很多并行化的輕量級線程。kernel在device上執(zhí)行時實際上是啟動很多線程，一個kernel所啟動的所有線程稱為一個網(wǎng)格（grid），同一個網(wǎng)格上的線程共享相同的全局內(nèi)存空間，grid是線程結(jié)構(gòu)的第一層次，而網(wǎng)格又可以分為很多線程塊（block），一個線程塊里面包含很多線程，這是第二個層次。線程兩層組織結(jié)構(gòu)如下圖所示，這是一個gird和block均為2-dim的線程組織。grid和block都是定義為dim3類型的變量，dim3可以看成是包含三個無符號整數(shù)（x，y，z）成員的結(jié)構(gòu)體變量，在定義時，缺省值初始化為1。因此grid和block可以靈活地定義為1-dim，2-dim以及3-dim結(jié)構(gòu)，對于圖中結(jié)構(gòu)（主要水平方向為x軸），定義的grid和block如下所示，kernel在調(diào)用時也必須通過執(zhí)行配置（http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#execution-configuration）<<<grid, block>>>來指定kernel所使用的線程數(shù)及結(jié)構(gòu)。
dim3 grid(3, 2); dim3 block(4, 3); kernel_fun<<< grid, block >>>(prams...);
所以，一個線程需要兩個內(nèi)置的坐標變量（blockIdx，threadIdx）來唯一標識，它們都是dim3類型變量，其中blockIdx指明線程所在grid中的位置，而threaIdx指明線程所在block中的位置，如圖中的Thread (1,1)滿足：
threadIdx.x = 1threadIdx.y = 1blockIdx.x = 1blockIdx.y = 1
一個線程塊上的線程是放在同一個流式多處理器（SM)上的，但是單個SM的資源有限，這導致線程塊中的線程數(shù)是有限制的，現(xiàn)代GPUs的線程塊可支持的線程數(shù)可達1024個。有時候，我們要知道一個線程在blcok中的全局ID，此時就必須還要知道block的組織結(jié)構(gòu)，這是通過線程的內(nèi)置變量blockDim來獲得。它獲取線程塊各個維度的大小。對于一個2-dim的block，線程（x,y)的ID值為,如果是3-dim的block,線程(x,y,z)的ID值為。另外線程還有內(nèi)置變量gridDim，用于獲得網(wǎng)格塊各個維度的大小。
kernel的這種線程組織結(jié)構(gòu)天然適合vector,matrix等運算，如我們將利用上圖2-dim結(jié)構(gòu)實現(xiàn)兩個矩陣的加法，每個線程負責處理每個位置的兩個元素相加，代碼如下所示。線程塊大小為(16, 16)，然后將N*N大小的矩陣均分為不同的線程塊來執(zhí)行加法運算。
// 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調(diào)用
MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
...
}
此外這里簡單介紹一下CUDA的內(nèi)存模型，如下圖所示?？梢钥吹?，每個線程有自己的私有本地內(nèi)存（Local Memory），而每個線程塊有包含共享內(nèi)存（Shared Memory）,可以被線程塊中所有線程共享，其生命周期與線程塊一致。此外，所有的線程都可以訪問全局內(nèi)存（Global Memory）。還可以訪問一些只讀內(nèi)存塊：常量內(nèi)存（Constant Memory）和紋理內(nèi)存（Texture Memory）。內(nèi)存結(jié)構(gòu)涉及到程序優(yōu)化，這里不深入探討它們。
CUDA內(nèi)存模型
還有重要一點，你需要對GPU的硬件實現(xiàn)有一個基本的認識。上面說到了kernel的線程組織層次，那么一個kernel實際上會啟動很多線程，這些線程是邏輯上并行的，但是在物理層卻并不一定。這其實和CPU的多線程有類似之處，多線程如果沒有多核支持，在物理層也是無法實現(xiàn)并行的。但是好在GPU存在很多CUDA核心，充分利用CUDA核心可以充分發(fā)揮GPU的并行計算能力。GPU硬件的一個核心組件是SM，前面已經(jīng)說過，SM是英文名是 Streaming Multiprocessor，翻譯過來就是流式多處理器。SM的核心組件包括CUDA核心，共享內(nèi)存，寄存器等，SM可以并發(fā)地執(zhí)行數(shù)百個線程，并發(fā)能力就取決于SM所擁有的資源數(shù)。當一個kernel被執(zhí)行時，它的gird中的線程塊被分配到SM上，一個線程塊只能在一個SM上被調(diào)度。SM一般可以調(diào)度多個線程塊，這要看SM本身的能力。那么有可能一個kernel的各個線程塊被分配多個SM，所以grid只是邏輯層，而SM才是執(zhí)行的物理層。SM采用的是SIMT（鏈接：http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#simt-architecture）(Single-Instruction, Multiple-Thread，單指令多線程)架構(gòu)，基本的執(zhí)行單元是線程束（wraps)，線程束包含32個線程，這些線程同時執(zhí)行相同的指令，但是每個線程都包含自己的指令地址計數(shù)器和寄存器狀態(tài)，也有自己獨立的執(zhí)行路徑。所以盡管線程束中的線程同時從同一程序地址執(zhí)行，但是可能具有不同的行為，比如遇到了分支結(jié)構(gòu)，一些線程可能進入這個分支，但是另外一些有可能不執(zhí)行，它們只能死等，因為GPU規(guī)定線程束中所有線程在同一周期執(zhí)行相同的指令，線程束分化會導致性能下降。當線程塊被劃分到某個SM上時，它將進一步劃分為多個線程束，因為這才是SM的基本執(zhí)行單元，但是一個SM同時并發(fā)的線程束數(shù)是有限的。這是因為資源限制，SM要為每個線程塊分配共享內(nèi)存，而也要為每個線程束中的線程分配獨立的寄存器。所以SM的配置會影響其所支持的線程塊和線程束并發(fā)數(shù)量?？傊褪蔷W(wǎng)格和線程塊只是邏輯劃分，一個kernel的所有線程其實在物理層是不一定同時并發(fā)的。所以kernel的grid和block的配置不同，性能會出現(xiàn)差異，這點是要特別注意的。還有，由于SM的基本執(zhí)行單元是包含32個線程的線程束，所以block大小一般要設置為32的倍數(shù)。
CUDA編程的邏輯層和物理層
在進行CUDA編程前，可以先檢查一下自己的GPU的硬件配置，這樣才可以有的放矢，可以通過下面的程序獲得GPU的配置屬性：
int dev = 0;
cudaDeviceProp devProp;
CHECK(cudaGetDeviceProperties(&devProp, dev));
std::cout << "使用GPU device " << dev << ": " << devProp.name << std::endl;
std::cout << "SM的數(shù)量：" << devProp.multiProcessorCount << std::endl;
std::cout << "每個線程塊的共享內(nèi)存大?。? << devProp.sharedMemPerBlock / 1024.0 << " KB" << std::endl;
std::cout << "每個線程塊的最大線程數(shù)：" << devProp.maxThreadsPerBlock << std::endl;
std::cout << "每個EM的最大線程數(shù)：" << devProp.maxThreadsPerMultiProcessor << std::endl;
std::cout << "每個EM的最大線程束數(shù)：" << devProp.maxThreadsPerMultiProcessor / 32 << std::endl;
// 輸出如下
使用GPU device 0: GeForce GT 730
SM的數(shù)量：2
每個線程塊的共享內(nèi)存大?。?8 KB
每個線程塊的最大線程數(shù)：1024
每個EM的最大線程數(shù)：2048
每個EM的最大線程束數(shù)：64
好吧，GT 730顯卡確實有點渣，只有2個SM，嗚嗚......
知道了CUDA編程基礎，我們就來個簡單的實戰(zhàn)，利用CUDA編程實現(xiàn)兩個向量的加法，在實現(xiàn)之前，先簡單介紹一下CUDA編程中內(nèi)存管理API。首先是在device上分配內(nèi)存的cudaMalloc函數(shù)：
cudaError_t cudaMalloc(void** devPtr, size_t size);
這個函數(shù)和C語言中的malloc類似，但是在device上申請一定字節(jié)大小的顯存，其中devPtr是指向所分配內(nèi)存的指針。同時要釋放分配的內(nèi)存使用cudaFree函數(shù)，這和C語言中的free函數(shù)對應。另外一個重要的函數(shù)是負責host和device之間數(shù)據(jù)通信的cudaMemcpy函數(shù)：
cudaError_t cudaMalloc(void** devPtr, size_t size);
其中src指向數(shù)據(jù)源，而dst是目標區(qū)域，count是復制的字節(jié)數(shù)，其中kind控制復制的方向：cudaMemcpyHostToHost, cudaMemcpyHostToDevice,
cudaMemcpyDeviceToHost及cudaMemcpyDeviceToDevice，如cudaMemcpyHostToDevice將host上數(shù)據(jù)拷貝到device上。
現(xiàn)在我們來實現(xiàn)一個向量加法的實例，這里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的線程數(shù)，有時候向量的元素數(shù)很多，這時候可以將在每個線程實現(xiàn)多個元素（元素總數(shù)/線程總數(shù)）的加法，相當于使用了多個grid來處理，這是一種grid-stride loop（鏈接：https://devblogs.nvidia.com/cuda-pro-tip-write-flexible-kernels-grid-stride-loops/）方式，不過下面的例子一個線程只處理一個元素，所以kernel里面的循環(huán)是不執(zhí)行的。下面我們具體實現(xiàn)向量加法：
int main()
{
int N = 1 << 20;
int nBytes = N * sizeof(float);
// 申請host內(nèi)存
float *x, *y, *z;
x = (float*)malloc(nBytes);
y = (float*)malloc(nBytes);
z = (float*)malloc(nBytes);
// 初始化數(shù)據(jù)
for (int i = 0; i < N; ++i)
{
x[i] = 10.0;
y[i] = 20.0;
}
// 申請device內(nèi)存
float *d_x, *d_y, *d_z;
cudaMalloc((void**)&d_x, nBytes);
cudaMalloc((void**)&d_y, nBytes);
cudaMalloc((void**)&d_z, nBytes);
// 將host數(shù)據(jù)拷貝到device
cudaMemcpy((void*)d_x, (void*)x, nBytes, cudaMemcpyHostToDevice);
cudaMemcpy((void*)d_y, (void*)y, nBytes, cudaMemcpyHostToDevice);
// 定義kernel的執(zhí)行配置
dim3 blockSize(256);
dim3 gridSize((N + blockSize.x - 1) / blockSize.x);
// 執(zhí)行kernel
add << < gridSize, blockSize >> >(d_x, d_y, d_z, N);
// 將device得到的結(jié)果拷貝到host
cudaMemcpy((void*)z, (void*)d_z, nBytes, cudaMemcpyHostToDevice);
// 檢查執(zhí)行結(jié)果
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內(nèi)存
cudaFree(d_x);
cudaFree(d_y);
cudaFree(d_z);
// 釋放host內(nèi)存
free(x);
free(y);
free(z);
return0;
}
這里我們的向量大小為1<<20，而block大小為256，那么grid大小是4096，kernel的線程層級結(jié)構(gòu)如下圖所示：
kernel的線程層次結(jié)構(gòu). 來源：https://devblogs.nvidia.com/even-easier-introduction-cuda/
使用nvprof工具可以分析kernel運行情況，結(jié)果如下所示，可以看到kernel函數(shù)費時約1.5ms。
你調(diào)整block的大小，對比不同配置下的kernel運行情況，我這里測試的是當block為128時，kernel費時約1.6ms，而block為512時kernel費時約1.7ms，當block為64時，kernel費時約2.3ms?？磥聿皇莃lock越大越好，而要適當選擇。
在上面的實現(xiàn)中，我們需要單獨在host和device上進行內(nèi)存分配，并且要進行數(shù)據(jù)拷貝，這是很容易出錯的。好在CUDA 6.0引入統(tǒng)一內(nèi)存（Unified Memory）（鏈接：http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#um-unified-memory-programming-hd）來避免這種麻煩，簡單來說就是統(tǒng)一內(nèi)存使用一個托管內(nèi)存來共同管理host和device中的內(nèi)存，并且自動在host和device中進行數(shù)據(jù)傳輸。CUDA中使用cudaMallocManaged函數(shù)分配托管內(nèi)存：
cudaError_t cudaMallocManaged(void **devPtr, size_t size, unsigned int flag=0);
利用統(tǒng)一內(nèi)存，可以將上面的程序簡化如下：
int main()
{
int N = 1 << 20;
int nBytes = N * sizeof(float);
// 申請托管內(nèi)存
float *x, *y, *z;
cudaMallocManaged((void**)&x, nBytes);
cudaMallocManaged((void**)&y, nBytes);
cudaMallocManaged((void**)&z, nBytes);
// 初始化數(shù)據(jù)
for (int i = 0; i < N; ++i)
{
x[i] = 10.0;
y[i] = 20.0;
}
// 定義kernel的執(zhí)行配置
dim3 blockSize(256);
dim3 gridSize((N + blockSize.x - 1) / blockSize.x);
// 執(zhí)行kernel
add << < gridSize, blockSize >> >(x, y, z, N);
// 同步device 保證結(jié)果能正確訪問
cudaDeviceSynchronize();
// 檢查執(zhí)行結(jié)果
float maxError = 0.0;
for (int i = 0; i < N; i++)
maxError = fmax(maxError, fabs(z[i] - 30.0));
std::cout << "最大誤差: " << maxError << std::endl;
// 釋放內(nèi)存
cudaFree(x);
cudaFree(y);
cudaFree(z);
return0;
}
相比之前的代碼，使用統(tǒng)一內(nèi)存更簡潔了，值得注意的是kernel執(zhí)行是與host異步的，由于托管內(nèi)存自動進行數(shù)據(jù)傳輸，這里要用cudaDeviceSynchronize()函數(shù)保證device和host同步，這樣后面才可以正確訪問kernel計算的結(jié)果。
最后我們再實現(xiàn)一個稍微復雜一些的例子，就是兩個矩陣的乘法，設輸入矩陣為A和B,要得到C=A*B。實現(xiàn)思路是每個線程計算C的一個元素值：
，對于矩陣運算，應該選用grid和block為2-D的。首先定義矩陣的結(jié)構(gòu)體：
矩陣乘法實現(xiàn)模式
然后實現(xiàn)矩陣乘法的核函數(shù)，這里我們定義了兩個輔助的__device__函數(shù)分別用于獲取矩陣的元素值和為矩陣元素賦值，具體代碼如下：
最后我們采用統(tǒng)一內(nèi)存編寫矩陣相乘的測試實例：
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);
// 申請托管內(nèi)存
cudaMallocManaged((void**)&A.elements, nBytes);
cudaMallocManaged((void**)&B.elements, nBytes);
cudaMallocManaged((void**)&C.elements, nBytes);
// 初始化數(shù)據(jù)
for (int i = 0; i < width * height; ++i)
{
A.elements[i] = 1.0;
B.elements[i] = 2.0;
}
// 定義kernel的執(zhí)行配置
dim3 blockSize(32, 32);
dim3 gridSize((width + blockSize.x - 1) / blockSize.x,
(height + blockSize.y - 1) / blockSize.y);
// 執(zhí)行kernel
matMulKernel << < gridSize, blockSize >> >(A, B, C);
// 同步device 保證結(jié)果能正確訪問
cudaDeviceSynchronize();
// 檢查執(zhí)行結(jié)果
float maxError = 0.0;
for (int i = 0; i < width * height; i++)
maxError = fmax(maxError, fabs(C.elements[i] - 2 * width));
std::cout << C.elements[0] << std::endl;
std::cout << "最大誤差: " << maxError << std::endl;
return0;
}
這里矩陣大小為1024*1024，設計的線程的block大小為(32, 32)，那么grid大小為(32, 32)，最終測試結(jié)果如下：
當然，這不是最高效的實現(xiàn)，后面可以繼續(xù)優(yōu)化...
最后只有一句話：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/）
