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