CUDA編程入門極簡教程（轉）

前言

2006年，NVIDIA公司發布了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安裝教程可以參考這里。

CUDA編程模型支持的編程語言

CUDA編程模型基礎

在給出CUDA的編程實例之前，這里先對CUDA編程模型中的一些概念及基礎知識做個簡單介紹。CUDA編程模型是一個異構模型，需要CPU和GPU協同工作。在CUDA中，host和device是兩個重要的概念，我們用host指代CPU及其內存，而用device指代GPU及其內存。CUDA程序中既包含host程序，又包含device程序，它們分別在CPU和GPU上運行。同時，host與device之間可以進行通信，這樣它們之間可以進行數據拷貝。典型的CUDA程序的執行流程如下：

1. 分配host內存，並進行數據初始化；
2. 分配device內存，並從host將數據拷貝到device上；
3. 調用CUDA的核函數在device上完成指定的運算；
4. 將device上的運算結果拷貝到host上；
5. 釋放device和host上分配的內存。

上面流程中最重要的一個過程是調用CUDA的核函數來執行並行計算，kernel是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在調用時也必須通過執行配置<<<grid, block>>>來指定kernel所使用的線程數及結構。

dim3 grid(3, 2);
dim3 block(5, 3);
kernel_fun<<< grid, block >>>(prams...);

Kernel上的兩層線程組織結構（2-dim）

所以，一個線程需要兩個內置的坐標變量（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  ，線程  的ID值為  ，如果是3-dim的block  ，線程  的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）。內存結構涉及到程序優化，這里不深入探討它們。

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 (Single-Instruction, Multiple-Thread，單指令多線程)架構，基本的執行單元是線程束（wraps)，線程束包含32個線程，這些線程同時執行相同的指令，但是每個線程都包含自己的指令地址計數器和寄存器狀態，也有自己獨立的執行路徑。所以盡管線程束中的線程同時從同一程序地址執行，但是可能具有不同的行為，比如遇到了分支結構，一些線程可能進入這個分支，但是另外一些有可能不執行，它們只能死等，因為GPU規定線程束中所有線程在同一周期執行相同的指令，線程束分化會導致性能下降。當線程塊被划分到某個SM上時，它將進一步划分為多個線程束，因為這才是SM的基本執行單元，但是一個SM同時並發的線程束數是有限的。這是因為資源限制，SM要為每個線程塊分配共享內存，而也要為每個線程束中的線程分配獨立的寄存器。所以SM的配置會影響其所支持的線程塊和線程束並發數量。總之，就是網格和線程塊只是邏輯划分，一個kernel的所有線程其實在物理層是不一定同時並發的。所以kernel的grid和block的配置不同，性能會出現差異，這點是要特別注意的。還有，由於SM的基本執行單元是包含32個線程的線程束，所以block大小一般要設置為32的倍數。

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的數量：" << 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，嗚嗚......

向量加法實例

知道了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 cudaMemcpy(void* dst, const void* src, size_t count, cudaMemcpyKind kind)

其中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方式，不過下面的例子一個線程只處理一個元素，所以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的線程層級結構如下圖所示：

kernel的線程層次結構. 來源：https://devblogs.nvidia.com/even-easier-introduction-cuda/

使用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）來避免這種麻煩，簡單來說就是統一內存使用一個托管內存來共同管理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計算的結果。

矩陣乘法實例

最后我們再實現一個稍微復雜一些的例子，就是兩個矩陣的乘法，設輸入矩陣為  和  ，要得到  。實現思路是每個線程計算  的一個元素值  ，對於矩陣運算，應該選用grid和block為2-D的。首先定義矩陣的結構體：

// 矩陣類型，行優先，M(row, col) = *(M.elements + row * M.width + col)
struct Matrix
{
    int width;
    int height;
    float *elements;
};

矩陣乘法實現模式

然后實現矩陣乘法的核函數，這里我們定義了兩個輔助的__device__函數分別用於獲取矩陣的元素值和為矩陣元素賦值，具體代碼如下：

// 獲取矩陣A的(row, col)元素
__device__ float getElement(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(Matrix *A, 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, *B, *C;
    // 申請托管內存
    cudaMallocManaged((void**)&A, sizeof(Matrix));
    cudaMallocManaged((void**)&B, sizeof(Matrix));
    cudaMallocManaged((void**)&C, sizeof(Matrix));
    int nBytes = width * height * sizeof(float);
    cudaMallocManaged((void**)&A->elements, nBytes);
    cudaMallocManaged((void**)&B->elements, nBytes);
    cudaMallocManaged((void**)&C->elements, nBytes);

    // 初始化數據
    A->height = height;
    A->width = width;
    B->height = height;
    B->width = width;
    C->height = height;
    C->width = width;
    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] - 2 * width));
    std::cout << "最大誤差: " << maxError << std::endl;

    return 0;
}

這里矩陣大小為，設計的線程的block大小為(32, 32)，那么grid大小為(32, 32)，最終測試結果如下：

nvprof cuda9.exe
==16304== NVPROF is profiling process 16304, command: cuda9.exe
最大誤差: 0
==16304== Profiling application: cuda9.exe
==16304== Profiling result:
            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
 GPU activities:  100.00%  1.32752s         1  1.32752s  1.32752s  1.32752s  matMulKernel(Matrix*, Matrix*, Matrix*)
      API calls:   83.11%  1.32762s         1  1.32762s  1.32762s  1.32762s  cudaDeviceSynchronize
                   13.99%  223.40ms         6  37.233ms  37.341us  217.66ms  cudaMallocManaged
                    2.81%  44.810ms         1  44.810ms  44.810ms  44.810ms  cudaLaunch
                    0.08%  1.3300ms        94  14.149us       0ns  884.64us  cuDeviceGetAttribute
                    0.01%  199.03us         1  199.03us  199.03us  199.03us  cuDeviceGetName
                    0.00%  10.009us         1  10.009us  10.009us  10.009us  cuDeviceTotalMem
                    0.00%  6.5440us         1  6.5440us  6.5440us  6.5440us  cudaConfigureCall
                    0.00%  3.0800us         3  1.0260us     385ns  1.5400us  cudaSetupArgument
                    0.00%  2.6940us         3     898ns     385ns  1.5390us  cuDeviceGetCount
                    0.00%  1.9250us         2     962ns     385ns  1.5400us  cuDeviceGet

==16304== Unified Memory profiling result:
Device "GeForce GT 730 (0)"
   Count  Avg Size  Min Size  Max Size  Total Size  Total Time  Name
    2051  4.0000KB  4.0000KB  4.0000KB  8.011719MB  21.20721ms  Host To Device
     270  45.570KB  4.0000KB  1.0000MB  12.01563MB  7.032508ms  Device To Host

當然，這不是最高效的實現，后面可以繼續優化...

小結

最后只有一句話：CUDA入門容易，但是深入難！希望不是從入門到放棄...

 

參考資料

1. John Cheng, Max Grossman, Ty McKercher. Professional CUDA C Programming, 2014.
2. CUDA docs.
3. An Even Easier Introduction to CUDA.
4. Unified Memory in CUDA 6.
5. Maximizing Unified Memory Performance in CUDA.

轉自：https://zhuanlan.zhihu.com/p/34587739