[CUDA]CUDA編程實戰二——向量加法

CUDA C是一種在GPU上支持多線程並行化的語言，有了CUDA，很多需要多線程運行的程序變得簡單起來，今天我們將從CUDA的的向量加法說起。

問題定義

向量加法是十分常見的操作，對於一個長度為n的向量，其運算規則如下：

\[{c[i] = a[i] + b[i] for i < n} \]

即將對應位置上的元素依次進行相加。

C++實現

有了上述的算法，我們可以很快地寫出一個C++版本的實現，其實就是一個循環的事情。

#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;
}

很明顯，遍歷相加，幾乎沒啥代碼量，這里為了對比，我們加上了時間測量。

測試結果

最終的運行結果為16ms。

CUDA版本

在CUDA中，我們稱CPU為host，GPU為device，稱在device上運行的函數為核(kernel)函數,需要使用__global__來修飾。
還有兩個其他的修飾符號__device__和__host__，三者區別在於globa可以被cpu函數調用，device只可以被cuda代碼調用，host和device可以同時使用，以便在某個函數中可以同時兼容使用GPU或CPU。
在運行核函數時，可以將其放入多個blocks和多個threads中運行，所以在每次運行核函數需要定義每個block中的threads和需要的blocks數。
其函數形式如下所示：Kernel_fun<<<Blocks, ThreadsPreBlock>>>(...);
在邏輯上grid>block>thread。

於是我們可以實現一個GPU版本的向量加法：

#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執行配置，（1024*1024/512）個block，每個block里面有512個線程
    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;
}

cuda的運行分為以下幾個步驟：

1. cudaMalloc為在顯存上開辟一段內存空間，具體用法和malloc類似。
2. cudaMemcpy為內存拷貝函數，需要將Host上的數據拷貝到device上，不然無法運行。
3. 經過kernel函數運行，計算對應的結果，結果保存在顯存中。
4. 最后將算好的結果拷貝回Host，不要忘了free掉內存和顯存。

這里配置了（1024*1024/512）個block，每個block里面有512個線程，在計算時，這些線程理論上將同時運行。

測試結果

結果為179ms，確實出乎意外，這是由於kernel函數計算過於簡單，而GPU的調度同樣需要時間，使得GPU的時間實際上要高於cpu。