CUDA程式設計之快速入門

阿新 • • 發佈：2018-12-10

CUDA（Compute Unified Device Architecture）的中文全稱為計算統一裝置架構。做影象視覺領域的同學多多少少都會接觸到CUDA，畢竟要做效能速度優化，CUDA是個很重要的工具，CUDA是做視覺的同學難以繞過的一個坑，必須踩一踩才踏實。CUDA程式設計真的是入門容易精通難，具有計算機體系結構和C語言程式設計知識儲備的同學上手CUDA程式設計應該難度不會很大。本文章將通過以下五個方面幫助大家比較全面地瞭解CUDA程式設計最重要的知識點，做到快速入門：

GPU架構特點
CUDA執行緒模型
CUDA記憶體模型
CUDA程式設計模型
CUDA應用小例子

1. GPU架構特點

首先我們先談一談序列計算和平行計算。我們知道，高效能運算的關鍵利用多核處理器進行平行計算。

當我們求解一個計算機程式任務時，我們很自然的想法就是將該任務分解成一系列小任務，把這些小任務一一完成。在序列計算時，我們的想法就是讓我們的處理器每次處理一個計算任務，處理完一個計算任務後再計算下一個任務，直到所有小任務都完成了，那麼這個大的程式任務也就完成了。如下圖所示，就是我們怎麼用序列程式設計思想求解問題的步驟。

但是序列計算的缺點非常明顯，如果我們擁有多核處理器，我們可以利用多核處理器同時處理多個任務時，而且這些小任務並沒有關聯關係（不需要相互依賴，比如我的計算任務不需要用到你的計算結果），那我們為什麼還要使用序列程式設計呢？為了進一步加快大任務的計算速度，我們可以把一些獨立的模組分配到不同的處理器上進行同時計算（這就是並行），最後再將這些結果進行整合，完成一次任務計算。下圖就是將一個大的計算任務分解為小任務，然後將獨立的小任務分配到不同處理器進行平行計算，最後再通過序列程式把結果彙總完成這次的總的計算任務。

所以，一個程式可不可以進行平行計算，關鍵就在於我們要分析出該程式可以拆分出哪幾個執行模組，這些執行模組哪些是獨立的，哪些又是強依賴強耦合的，獨立的模組我們可以試著設計平行計算，充分利用多核處理器的優勢進一步加速我們的計算任務，強耦合模組我們就使用序列程式設計，利用序列+並行的程式設計思路完成一次高效能運算。

接下來我們談談CPU和GPU有什麼區別，他們倆各自有什麼特點，我們在談並行、序列計算時多次談到“多核”的概念，現在我們先從“核”的角度開始這個話題。首先CPU是專為順序序列處理而優化的幾個核心組成。而GPU則由數以千計的更小、更高效的核心組成，這些核心專門為同時處理多工而設計，可高效地處理並行任務。也就是，CPU雖然每個核心自身能力極強，處理任務上非常強悍，無奈他核心少，在平行計算上表現不佳；反觀GPU，雖然他的每個核心的計算能力不算強，但他勝在核心非常多，可以同時處理多個計算任務，在平行計算的支援上做得很好。

GPU和CPU的不同硬體特點決定了他們的應用場景，CPU是計算機的運算和控制的核心，GPU主要用作圖形影象處理。影象在計算機呈現的形式就是矩陣，我們對影象的處理其實就是操作各種矩陣進行計算，而很多矩陣的運算其實可以做並行化，這使得影象處理可以做得很快，因此GPU在圖形影象領域也有了大展拳腳的機會。下圖表示的就是一個多GPU計算機硬體系統，可以看出，一個GPU記憶體就有很多個SP和各類記憶體，這些硬體都是GPU進行高效平行計算的基礎。

現在再從資料處理的角度來對比CPU和GPU的特點。CPU需要很強的通用性來處理各種不同的資料型別，比如整型、浮點數等，同時它又必須擅長處理邏輯判斷所導致的大量分支跳轉和中斷處理，所以CPU其實就是一個能力很強的夥計，他能把很多事處理得妥妥當當，當然啦我們需要給他很多資源供他使用（各種硬體），這也導致了CPU不可能有太多核心（核心總數不超過16）。而GPU面對的則是型別高度統一的、相互無依賴的大規模資料和不需要被打斷的純淨的計算環境，GPU有非常多核心（費米架構就有512核），雖然其核心的能力遠沒有CPU的核心強，但是勝在多，在處理簡單計算任務時呈現出“人多力量大”的優勢，這就是平行計算的魅力。

整理一下兩者特點就是：

CPU：擅長流程控制和邏輯處理，不規則資料結構，不可預測儲存結構，單執行緒程式，分支密集型演算法
GPU：擅長資料平行計算，規則資料結構，可預測儲存模式

現在的計算機體系架構中，要完成CUDA平行計算，單靠GPU一人之力是不能完成計算任務的，必須藉助CPU來協同配合完成一次高效能的平行計算任務。

一般而言，並行部分在GPU上執行，序列部分在CPU執行，這就是異構計算。具體一點，異構計算的意思就是不同體系結構的處理器相互協作完成計算任務。CPU負責總體的程式流程，而GPU負責具體的計算任務，當GPU各個執行緒完成計算任務後，我們就將GPU那邊計算得到的結果拷貝到CPU端，完成一次計算任務。

所以應用程式利用GPU實現加速的總體分工就是：密集計算程式碼（約佔5%的程式碼量）由GPU負責完成，剩餘序列程式碼由CPU負責執行。

2. CUDA執行緒模型

下面我們介紹CUDA的執行緒組織結構。首先我們都知道，執行緒是程式執行的最基本單元，CUDA的平行計算就是通過成千上萬個執行緒的並行執行來實現的。下面的機構圖說明了GPU的不同層次的結構。

CUDA的執行緒模型從小往大來總結就是：

Thread：執行緒，並行的基本單位
Thread Block：執行緒塊，互相合作的執行緒組，執行緒塊有如下幾個特點：

允許彼此同步
可以通過共享記憶體快速交換資料
以1維、2維或3維組織

Grid：一組執行緒塊

以1維、2維組織
共享全域性記憶體

Kernel：在GPU上執行的核心程式，這個kernel函式是執行在某個Grid上的。

One kernel <-> One Grid

每一個block和每個thread都有自己的ID，我們通過相應的索引找到相應的執行緒和執行緒塊。

threadIdx，blockIdx
Block ID: 1D or 2D
Thread ID: 1D, 2D or 3D

理解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結構，kernel呼叫時也必須通過執行配置<<<grid, block>>>來指定kernel所使用的網格維度和執行緒塊維度。舉個例子，我們以上圖為例，分析怎麼通過<<<grid,block>>>>這種標記方式索引到我們想要的那個執行緒。CUDA的這種<<<grid,block>>>其實就是一個多級索引的方法，第一級索引是(grid.xIdx, grid.yIdy)，對應上圖例子就是(1, 1)，通過它我們就能找到了這個執行緒塊的位置，然後我們啟動二級索引(block.xIdx, block.yIdx, block.zIdx)來定位到指定的執行緒。這就是我們CUDA的執行緒組織結構。

這裡想談談SP和SM（流處理器），很多人會被這兩個專業名詞搞得暈頭轉向。

SP：最基本的處理單元，streaming processor，也稱為CUDA core。最後具體的指令和任務都是在SP上處理的。GPU進行平行計算，也就是很多個SP同時做處理。
SM：多個SP加上其他的一些資源組成一個streaming multiprocessor。也叫GPU大核，其他資源如：warp scheduler，register，shared memory等。SM可以看做GPU的心臟（對比CPU核心），register和shared memory是SM的稀缺資源。CUDA將這些資源分配給所有駐留在SM中的threads。因此，這些有限的資源就使每個SM中active warps有非常嚴格的限制，也就限制了並行能力。

需要指出，每個SM包含的SP數量依據GPU架構而不同，Fermi架構GF100是32個，GF10X是48個，Kepler架構都是192個，Maxwell都是128個。

簡而言之，SP是執行緒執行的硬體單位，SM中包含多個SP，一個GPU可以有多個SM（比如16個），最終一個GPU可能包含有上千個SP。這麼多核心“同時執行”，速度可想而知，這個引號只是想表明實際上，軟體邏輯上是所有SP是並行的，但是物理上並不是所有SP都能同時執行計算（比如我們只有8個SM卻有1024個執行緒塊需要排程處理），因為有些會處於掛起，就緒等其他狀態，這有關GPU的執行緒排程。

下面這個圖將從硬體角度和軟體角度解釋CUDA的執行緒模型。

每個執行緒由每個執行緒處理器（SP）執行
執行緒塊由多核處理器（SM）執行
一個kernel其實由一個grid來執行，一個kernel一次只能在一個GPU上執行

block是軟體概念，一個block只會由一個sm排程，程式設計師在開發時，通過設定block的屬性，告訴GPU硬體，我有多少個執行緒，執行緒怎麼組織。而具體怎麼排程由sm的warps scheduler負責，block一旦被分配好SM，該block就會一直駐留在該SM中，直到執行結束。一個SM可以同時擁有多個blocks，但需要序列執行。下圖顯示了GPU內部的硬體架構：

3. CUDA記憶體模型

CUDA中的記憶體模型分為以下幾個層次：

每個執行緒都用自己的registers（暫存器）
每個執行緒都有自己的local memory（區域性記憶體）
每個執行緒塊內都有自己的shared memory（共享記憶體），所有執行緒塊內的所有執行緒共享這段記憶體資源
每個grid都有自己的global memory（全域性記憶體），不同執行緒塊的執行緒都可使用
每個grid都有自己的constant memory（常量記憶體）和texture memory（紋理記憶體），），不同執行緒塊的執行緒都可使用

執行緒訪問這幾類儲存器的速度是register > local memory >shared memory > global memory

下面這幅圖表示就是這些記憶體在計算機架構中的所在層次。

4. CUDA程式設計模型

上面講了這麼多硬體相關的知識點，現在終於可以開始說說CUDA是怎麼寫程式的了。

我們先捋一捋常見的CUDA術語：

第一個要掌握的程式設計要點：我們怎麼寫一個能在GPU跑的程式或函式呢？

通過關鍵字就可以表示某個程式在CPU上跑還是在GPU上跑！如下表所示，比如我們用__global__定義一個kernel函式，就是CPU上呼叫，GPU上執行，注意__global__函式的返回值必須設定為void。

第二個程式設計要點：CPU和GPU間的資料傳輸怎麼寫？

首先介紹在GPU記憶體分配回收記憶體的函式介面：

cudaMalloc(): 在裝置端分配global memory
cudaFree(): 釋放儲存空間

CPU的資料和GPU端資料做資料傳輸的函式介面是一樣的，他們通過傳遞的函式實參（列舉型別）來表示傳輸方向：

cudaMemcpy(void dst, void src, size_t nbytes, enum cudaMemcpyKind direction)

enum cudaMemcpyKind:

cudaMemcpyHostToDevice（CPU到GPU）
cudaMemcpyDeviceToHost（GPU到CPU）
cudaMemcpyDeviceToDevice（GPU到GPU）

第三個程式設計要點是：怎麼用程式碼表示執行緒組織模型？我們可以用dim3類來表示網格和執行緒塊的組織方式，網格grid可以表示為一維和二維格式，執行緒塊block可以表示為一維、二維和三維的資料格式。

dim3 DimGrid(100, 50);  //5000個執行緒塊，維度是100*50
dim3 DimBlock(4, 8, 8);  //每個線層塊內包含256個執行緒，執行緒塊內的維度是4*8*8

接下來介紹一個非常重要又很難懂的一個知識點，我們怎麼計算執行緒號呢？

1.使用N個執行緒塊，每一個執行緒塊只有一個執行緒，即

dim3 dimGrid(N);
dim3 dimBlock(1);

此時的執行緒號的計算方式就是

threadId = blockIdx.x;

其中threadId的取值範圍為0到N-1。對於這種情況，我們可以將其看作是一個列向量，列向量中的每一行對應一個執行緒塊。列向量中每一行只有1個元素，對應一個執行緒。

2.使用M×N個執行緒塊，每個執行緒塊1個執行緒

由於執行緒塊是2維的，故可以看做是一個M*N的2維矩陣，其執行緒號有兩個維度，即：

dim3 dimGrid(M,N);
dim3 dimBlock(1);

其中

blockIdx.x 取值0到M-1
blcokIdx.y 取值0到N-1

這種情況一般用於處理2維資料結構，比如2維影象。每一個畫素用一個執行緒來處理，此時需要執行緒號來對映影象畫素的對應位置，如

pos = blockIdx.y * blcokDim.x + blockIdx.x; //其中gridDim.x等於M

3.使用一個執行緒塊，該執行緒具有N個執行緒，即

dim3 dimGrid(1);
dim3 dimBlock(N);

此時執行緒號的計算方式為

threadId = threadIdx.x;

其中threadId的範圍是0到N-1，對於這種情況，可以看做是一個行向量，行向量中的每一個元素的每一個元素對應著一個執行緒。

4.使用M個執行緒塊，每個執行緒塊內含有N個執行緒，即

dim3 dimGrid(M);
dim3 dimBlock(N);

這種情況，可以把它想象成二維矩陣，矩陣的行與執行緒塊對應，矩陣的列與執行緒編號對應，那執行緒號的計算方式為

threadId = threadIdx.x + blcokIdx*blockDim.x;

上面其實就是把二維的索引空間轉換為一維索引空間的過程。

5.使用M×N的二維執行緒塊，每一個執行緒塊具有P×Q個執行緒，即

dim3 dimGrid(M, N);
dim3 dimBlock(P, Q);

這種情況其實是我們遇到的最多情況，特別適用於處理具有二維資料結構的演算法，比如影象處理領域。

其索引有兩個維度

threadId.x = blockIdx.x*blockDim.x+threadIdx.x;
threadId.y = blockIdx.y*blockDim.y+threadIdx.y;

上述公式就是把執行緒和執行緒塊的索引對映為影象畫素座標的計算方法。

CUDA應用例子

我們已經掌握了CUDA程式設計的基本語法，現在我們開始以一些小例子來真正上手CUDA。

首先我們編寫一個程式，檢視我們GPU的一些硬體配置情況。

#include "device_launch_parameters.h"
#include <iostream>

int main()
{
    int deviceCount;
    cudaGetDeviceCount(&deviceCount);
    for(int i=0;i<deviceCount;i++)
    {
        cudaDeviceProp devProp;
        cudaGetDeviceProperties(&devProp, i);
        std::cout << "使用GPU device " << i << ": " << devProp.name << std::endl;
        std::cout << "裝置全域性記憶體總量： " << devProp.totalGlobalMem / 1024 / 1024 << "MB" << 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 << "裝置上一個執行緒塊（Block）種可用的32位暫存器數量： " << devProp.regsPerBlock << std::endl;
        std::cout << "每個EM的最大執行緒數：" << devProp.maxThreadsPerMultiProcessor << std::endl;
        std::cout << "每個EM的最大執行緒束數：" << devProp.maxThreadsPerMultiProcessor / 32 << std::endl;
        std::cout << "裝置上多處理器的數量： " << devProp.multiProcessorCount << std::endl;
        std::cout << "======================================================" << std::endl;     
        
    }
    return 0;
}

我們利用nvcc來編譯程式。

nvcc test1.cu -o test1

輸出結果：因為我的伺服器是8個TITAN GPU，為了省略重複資訊，下面只顯示兩個GPU結果

使用GPU device 0: TITAN X (Pascal)
裝置全域性記憶體總量： 12189MB
SM的數量：28
每個執行緒塊的共享記憶體大小：48 KB
每個執行緒塊的最大執行緒數：1024
裝置上一個執行緒塊（Block）種可用的32位暫存器數量： 65536
每個EM的最大執行緒數：2048
每個EM的最大執行緒束數：64
裝置上多處理器的數量： 28
======================================================
使用GPU device 1: TITAN X (Pascal)
裝置全域性記憶體總量： 12189MB
SM的數量：28
每個執行緒塊的共享記憶體大小：48 KB
每個執行緒塊的最大執行緒數：1024
裝置上一個執行緒塊（Block）種可用的32位暫存器數量： 65536
每個EM的最大執行緒數：2048
每個EM的最大執行緒束數：64
裝置上多處理器的數量： 28
======================================================

.......

第一個計算任務：將兩個元素數目為1024×1024的float陣列相加。

首先我們思考一下如果只用CPU我們怎麼序列完成這個任務。

#include <iostream>
#include <stdlib.h>
#include <sys/time.h>
#include <math.h>

using namespace std;

int main()
{
    struct timeval start, end;
    gettimeofday( &start, NULL );
    float*A, *B, *C;
    int n = 1024 * 1024;
    int size = n * sizeof(float);
    A = (float*)malloc(size);
    B = (float*)malloc(size);
    C = (float*)malloc(size);

    for(int i=0;i<n;i++)
    {
        A[i] = 90.0;
        B[i] = 10.0;
    }
    
    for(int i=0;i<n;i++)
    {
        C[i] = A[i] + B[i];
    }

    float max_error = 0.0;
    for(int i=0;i<n;i++)
    {
        max_error += fabs(100.0-C[i]);
    }
    cout << "max_error is " << max_error << endl;
    gettimeofday( &end, NULL );
    int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
    cout << "total time is " << timeuse/1000 << "ms" <<endl;
    return 0;
}

CPU方式輸出結果

max_error is 0
total time is 22ms

如果我們使用GPU來做平行計算，速度將會如何呢？

程式設計要點：

每個Block中的Thread數最大不超過512；
為了充分利用SM，Block數儘可能多，>100。

#include "cuda_runtime.h"
#include <stdlib.h>
#include <iostream>
#include <sys/time.h>

using namespace std;

__global__ void Plus(float A[], float B[], float C[], int n)
{
    int i = blockDim.x * blockIdx.x + threadIdx.x;
    C[i] = A[i] + B[i];
}

int main()
{
    struct timeval start, end;
    gettimeofday( &start, NULL );
    float*A, *Ad, *B, *Bd, *C, *Cd;
    int n = 1024 * 1024;
    int size = n * sizeof(float);

    // CPU端分配記憶體
    A = (float*)malloc(size);
    B = (float*)malloc(size);
    C = (float*)malloc(size);

    // 初始化陣列
    for(int i=0;i<n;i++)
    {
        A[i] = 90.0;
        B[i] = 10.0;
    }

    // GPU端分配記憶體
    cudaMalloc((void**)&Ad, size);
    cudaMalloc((void**)&Bd, size);
    cudaMalloc((void**)&Cd, size);

    // CPU的資料拷貝到GPU端
    cudaMemcpy(Ad, A, size, cudaMemcpyHostToDevice);
    cudaMemcpy(Bd, B, size, cudaMemcpyHostToDevice);
    cudaMemcpy(Bd, B, size, cudaMemcpyHostToDevice);

    // 定義kernel執行配置，4個block，每個block裡面有256個執行緒
    dim3 dimBlock(512);
    dim3 dimGrid(n/512);

    // 執行kernel
    Plus<<<dimGrid, dimBlock>>>(Ad, Bd, Cd, n);

    // 將在GPU端計算好的結果拷貝回CPU端
    cudaMemcpy(C, Cd, size, cudaMemcpyDeviceToHost);

    // 校驗誤差
    float max_error = 0.0;
    for(int i=0;i<n;i++)
    {
        max_error += fabs(100.0 - C[i]);
    }

    cout << "max error is " << max_error << endl;

    // 釋放CPU端、GPU端的記憶體
    free(A);
    free(B);
    free(C);
    cudaFree(Ad);
    cudaFree(Bd);
    cudaFree(Cd);
    gettimeofday( &end, NULL );
    int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
    cout << "total time is " << timeuse/1000 << "ms" <<endl;
    return 0;
}

GPU方式輸出結果

max error is 0
total time is 1278ms

由上面的例子看出，使用CUDA程式設計時我們看不到for迴圈了，因為CPU程式設計的迴圈已經被分散到各個thread上做了，所以我們也就看到不到for一類的語句。從結果上看，CPU的迴圈計算的速度比GPU計算快多了，原因就在於CUDA中有大量的記憶體拷貝操作（資料傳輸花費了大量時間，而計算時間卻非常少），如果計算量比較小的話，CPU計算會更合適一些。

下面計算一個稍微複雜的例子，矩陣加法，即對兩個矩陣對應座標的元素相加後的結果儲存在第三個的對應位置的元素上。

值得注意的是，這個計算任務我採用了二維陣列的計算方式，注意一下二維陣列在CUDA程式設計中的寫法。

CPU版本

#include <stdlib.h>
#include <iostream>
#include <sys/time.h>
#include <math.h>

#define ROWS 1024
#define COLS 1024

using namespace std;

int main()
{
    struct timeval start, end;
    gettimeofday( &start, NULL );
    int *A, **A_ptr, *B, **B_ptr, *C, **C_ptr;
    int total_size = ROWS*COLS*sizeof(int);
    A = (int*)malloc(total_size);
    B = (int*)malloc(total_size);
    C = (int*)malloc(total_size);
    A_ptr = (int**)malloc(ROWS*sizeof(int*));
    B_ptr = (int**)malloc(ROWS*sizeof(int*));
    C_ptr = (int**)malloc(ROWS*sizeof(int*));
    
    //CPU一維陣列初始化
    for(int i=0;i<ROWS*COLS;i++)
    {
        A[i] = 80;
        B[i] = 20;
    }
    
    for(int i=0;i<ROWS;i++)
    {
        A_ptr[i] = A + COLS*i;
        B_ptr[i] = B + COLS*i;
        C_ptr[i] = C + COLS*i;
    }
    
    for(int i=0;i<ROWS;i++)
        for(int j=0;j<COLS;j++)
        {
            C_ptr[i][j] = A_ptr[i][j] + B_ptr[i][j];
        }
        
    //檢查結果
    int max_error = 0;
    for(int i=0;i<ROWS*COLS;i++)
    {
        //cout << C[i] << endl;
        max_error += abs(100-C[i]);
    }
    
    cout << "max_error is " << max_error <<endl;     
    gettimeofday( &end, NULL );
    int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
    cout << "total time is " << timeuse/1000 << "ms" <<endl;
    
    return 0;
}

CPU方式輸出

max_error is 0
total time is 29ms

GPU版本

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <sys/time.h> 
#include <stdio.h>
#include <math.h>
#define Row  1024
#define Col 1024
 
 
__global__ void addKernel(int **C,  int **A, int ** B)
{
    int idx = threadIdx.x + blockDim.x * blockIdx.x;
    int idy = threadIdx.y + blockDim.y * blockIdx.y;
    if (idx < Col && idy < Row) {
        C[idy][idx] = A[idy][idx] + B[idy][idx];
    }
}
 
int main()
{
    struct timeval start, end;
    gettimeofday( &start, NULL );

    int **A = (int **)malloc(sizeof(int*) * Row);
    int **B = (int **)malloc(sizeof(int*) * Row);
    int **C = (int **)malloc(sizeof(int*) * Row);
    int *dataA = (int *)malloc(sizeof(int) * Row * Col);
    int *dataB = (int *)malloc(sizeof(int) * Row * Col);
    int *dataC = (int *)malloc(sizeof(int) * Row * Col);
    int **d_A;
    int **d_B;
    int **d_C;
    int *d_dataA;
    int *d_dataB;
    int *d_dataC;
    //malloc device memory
    cudaMalloc((void**)&d_A, sizeof(int **) * Row);
    cudaMalloc((void**)&d_B, sizeof(int **) * Row);
    cudaMalloc((void**)&d_C, sizeof(int **) * Row);
    cudaMalloc((void**)&d_dataA, sizeof(int) *Row*Col);
    cudaMalloc((void**)&d_dataB, sizeof(int) *Row*Col);
    cudaMalloc((void**)&d_dataC, sizeof(int) *Row*Col);
    //set value
    for (int i = 0; i < Row*Col; i++) {
        dataA[i] = 90;
        dataB[i] = 10;
    }
    //將主機指標A指向裝置資料位置，目的是讓裝置二級指標能夠指向裝置資料一級指標
    //A 和  dataA 都傳到了裝置上，但是二者還沒有建立對應關係
    for (int i = 0; i < Row; i++) {
        A[i] = d_dataA + Col * i;
        B[i] = d_dataB + Col * i;
        C[i] = d_dataC + Col * i;
    }
                                                                
    cudaMemcpy(d_A, A, sizeof(int*) * Row, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, B, sizeof(int*) * Row, cudaMemcpyHostToDevice);
    cudaMemcpy(d_C, C, sizeof(int*) * Row, cudaMemcpyHostToDevice);
    cudaMemcpy(d_dataA, dataA, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);
    cudaMemcpy(d_dataB, dataB, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);
    dim3 threadPerBlock(16, 16);
    dim3 blockNumber( (Col + threadPerBlock.x - 1)/ threadPerBlock.x, (Row + threadPerBlock.y - 1) / threadPerBlock.y );
    printf("Block(%d,%d)   Grid(%d,%d).\n", threadPerBlock.x, threadPerBlock.y, blockNumber.x, blockNumber.y);
    addKernel << <blockNumber, threadPerBlock >> > (d_C, d_A, d_B);
    //拷貝計算資料-一級資料指標
    cudaMemcpy(dataC, d_dataC, sizeof(int) * Row * Col, cudaMemcpyDeviceToHost);
                                                                                             
    int max_error = 0;
    for(int i=0;i<Row*Col;i++)
    {
        //printf("%d\n", dataC[i]);
        max_error += abs(100-dataC[i]);
    }

    //釋放記憶體
    free(A);
    free(B);
    free(C);
    free(dataA);
    free(dataB);
    free(dataC);
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    cudaFree(d_dataA);
    cudaFree(d_dataB);
    cudaFree(d_dataC);

    printf("max_error is %d\n", max_error);
    gettimeofday( &end, NULL );
    int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
    printf("total time is %d ms\n", timeuse/1000);

    return 0;
}

GPU輸出

Block(16,16)   Grid(64,64).
max_error is 0
total time is 442 ms

從結果看出，CPU計算時間還是比GPU的計算時間短。這裡需要指出的是，這種二維陣列的程式寫法的效率並不高（雖然比較符合我們的思維方式），因為我們做了兩次訪存操作。所以一般而言，做高效能運算一般不會採取這種程式設計方式。

最後一個例子我們將計算一個更加複雜的任務，矩陣乘法

回顧一下矩陣乘法：兩矩陣相乘，左矩陣第一行乘以右矩陣第一列（分別相乘，第一個數乘第一個數），乘完之後相加，即為結果的第一行第一列的數，依次往下算，直到計算完所有矩陣元素。

CPU版本

#include <iostream>
#include <stdlib.h>
#include <sys/time.h>

#define ROWS 1024
#define COLS 1024

using namespace std;

void matrix_mul_cpu(float* M, float* N, float* P, int width)
{
    for(int i=0;i<width;i++)
        for(int j=0;j<width;j++)
        {
            float sum = 0.0;
            for(int k=0;k<width;k++)
            {
                float a = M[i*width+k];
                float b = N[k*width+j];
                sum += a*b;
            }
            P[i*width+j] = sum;
        }
}

int main()
{
    struct timeval start, end;
    gettimeofday( &start, NULL );
    float *A, *B, *C;
    int total_size = ROWS*COLS*sizeof(float);
    A = (float*)malloc(total_size);
    B = (float*)malloc(total_size);
    C = (float*)malloc(total_size);

    //CPU一維陣列初始化
    for(int i=0;i<ROWS*COLS;i++)
    {
        A[i] = 80.0;
        B[i] = 20.0;
    }

    matrix_mul_cpu(A, B, C, COLS);

    gettimeofday( &end, NULL );
    int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
    cout << "total time is " << timeuse/1000 << "ms" <<endl;

    return 0;
}

CPU輸出

total time is 7617ms

梳理一下CUDA求解矩陣乘法的思路：因為C=A×B，我們利用每個執行緒求解C矩陣每個(x, y)的元素，每個執行緒載入A的一行和B的一列，遍歷各自行列元素，對A、B對應的元素做一次乘法和一次加法。

GPU版本

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <sys/time.h> 
#include <stdio.h>
#include <math.h>
#define Row  1024
#define Col 1024

 
__global__ void matrix_mul_gpu(int *M, int* N, int* P, int width)
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;
    int j = threadIdx.y + blockDim.y * blockIdx.y;
                
    int sum = 0;
    for(int k=0;k<width;k++)
    {
        int a = M[j*width+k];
        int b = N[k*width+i];
        sum += a*b;
    }
    P[j*width+i] = sum;
}
 
int main()
{
    struct timeval start, end;
    gettimeofday( &start, NULL );

    int *A = (int *)malloc(sizeof(int) * Row * Col);
    int *B = (int *)malloc(sizeof(int) * Row * Col);
    int *C = (int *)malloc(sizeof(int) * Row * Col);
    //malloc device memory
    int *d_dataA, *d_dataB, *d_dataC;
    cudaMalloc((void**)&d_dataA, sizeof(int) *Row*Col);
    cudaMalloc((void**)&d_dataB, sizeof(int) *Row*Col);
    cudaMalloc((void**)&d_dataC, sizeof(int) *Row*Col);
    //set value
    for (int i = 0; i < Row*Col; i++) {
        A[i] = 90;
        B[i] = 10;
    }
                                                                
    cudaMemcpy(d_dataA, A, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);
    cudaMemcpy(d_dataB, B, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);
    dim3 threadPerBlock(16, 16);
    dim3 blockNumber((Col+threadPerBlock.x-1)/ threadPerBlock.x, (Row+threadPerBlock.y-1)/ threadPerBlock.y );
    printf("Block(%d,%d)   Grid(%d,%d).\n", threadPerBlock.x, threadPerBlock.y, blockNumber.x, blockNumber.y);
    matrix_mul_gpu << <blockNumber, threadPerBlock >> > (d_dataA, d_dataB, d_dataC, Col);
    //拷貝計算資料-一級資料指標
    cudaMemcpy(C, d_dataC, sizeof(int) * Row * Col, cudaMemcpyDeviceToHost);
                                                                                             
    //釋放記憶體
    free(A);
    free(B);
    free(C);
    cudaFree(d_dataA);
    cudaFree(d_dataB);
    cudaFree(d_dataC);

    gettimeofday( &end, NULL );
    int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
    printf("total time is %d ms\n", timeuse/1000);

    return 0;
}

GPU輸出

Block(16,16)   Grid(64,64).
total time is 506 ms

從這個矩陣乘法任務可以看出，我們通過GPU進行平行計算的方式僅花費了0.5秒，但是CPU序列計算方式卻花費了7.6秒，計算速度提升了十多倍，可見平行計算的威力！

CUDA程式設計之快速入門

1. GPU架構特點

2. CUDA執行緒模型

3. CUDA記憶體模型

4. CUDA程式設計模型

1.使用N個執行緒塊，每一個執行緒塊只有一個執行緒，即

2.使用M×N個執行緒塊，每個執行緒塊1個執行緒

3.使用一個執行緒塊，該執行緒具有N個執行緒，即

4.使用M個執行緒塊，每個執行緒塊內含有N個執行緒，即

5.使用M×N的二維執行緒塊，每一個執行緒塊具有P×Q個執行緒，即

CUDA應用例子

CUDA程式設計之快速入門

CUDA編程之快速入門

openvpn部署之快速入門實戰+一鍵部署openvpn腳本

php之快速入門學習-5(常量)

php之快速入門學習-6(常量)

php之快速入門學習-13(PHP 循環 - While 循環)

php之快速入門學習-16（PHP 魔術變量）

MyBatis 框架之快速入門程序

14.spark mllib之快速入門

10.spark sql之快速入門

正則表示式之快速入門篇

Python GUI程式設計之tkinter入門案例--簡單的自動檔案下載器

WebService一之快速入門

大資料Hadoop學習之快速入門方法

帶你走進Spark程式設計之Scala入門

Spark 2.4.0 程式設計指南--快速入門

Xamarin.IOS之快速入門

RockerMQ實戰之快速入門

React Native開發之快速入門React

react-native入門之快速入門---返回鍵實現

CUDA程式設計之快速入門

1. GPU架構特點

2. CUDA執行緒模型

3. CUDA記憶體模型

4. CUDA程式設計模型

1.使用N個執行緒塊，每一個執行緒塊只有一個執行緒，即

2.使用M×N個執行緒塊，每個執行緒塊1個執行緒

3.使用一個執行緒塊，該執行緒具有N個執行緒，即

4.使用M個執行緒塊，每個執行緒塊內含有N個執行緒，即

5.使用M×N的二維執行緒塊，每一個執行緒塊具有P×Q個執行緒，即

CUDA應用例子

相關推薦