[平行運算] CUDA入門

最近因為興趣使然, 開始碰了CUDA, 因此參考惹Nvidia的資料, 來寫個筆記

首先是CPU和GPU的比較, 這裡先簡略帶過, 之後有機會再來寫個一篇(是要挖多少坑....), 簡單來說GPU和CPU比較起來, GPU有非常多的核心(core)可以一次做很多事情, 因此如果需要大量計算, 且每次計算結果間彼此獨立的話, 可以使用GPU來同時進行運算, 這個概念被稱為是平行運算(parallel programming)

這樣講還是有點抽象, 因此用向量加法來作例子, 舉例有兩個向量, A=(1,2,3), B=(3,4,5), 要把兩個矩陣相加, 如果一一相加的話, 要加三次, 但如果分給GPU做, 雖然還是要做三次加法, 但三次加法可以同時做, 因此完成的時間變短, 這就是平行運算的精神

CUDA就是Nvidia拿來整合大家熟悉的程式語言(C,C++,java, fortan....), 和GPU間的平台

由於小弟寫比較多的是C/C++語言, 而且官方文件也是用C/C++語言, 因此以下用C/C++語言語法來介紹,

CUDA裡面有所謂的 kernel 概念, kernel內的程式碼是由GPU執行, 用 __global__ 做宣告(注意底標的橫槓有兩個), 這裡參考一下官方文件, 寫法如下:

//目標是在GPU裡面把兩個陣列相加

__global__ void add( int *a, int *b, int *c ) 
{
    int i = threadIdx.x;//threadIdx是CUDA內建用來指向特定thread的變數
    c[i]=a[i]+b[i];  
}

P.S. 在這裡可能有些人會想說什麼是thread, thread 中文譯為執行緒, 在作業系統的課本會花一大堆時間講,不過這裡簡單來說, 可以把thread想成電腦執行一件事情,就稱為thread

看完了GPU的code, 接著來看主程式怎麼和GPU互動

這裡其實有點複雜, 因為牽扯到兩個部分, 第一個是CUDA如何在主程式裡面表示thread, 第二個部分是CPU和GPU彼此間記憶體如何操作, 先看第一部分:

前面提到的threadIdx其實是一個三維的向量, 每組數字代表一個thread; 一群thread可以組成一個block(目前技術一個block最多可以有但1024個thread), 一群block可以組成一個grid; CUDA有一個特殊的語法, add<<<numBlocks, threadsPerBlock>>>(a, b, c); 如果我們不看<<< >>>裡面的內容的話,會發現就和一般的C語言函數一樣, 而加上的部份就是代表CUDA指定要用這麼多的執行緒去執行程式, 寫出來如下:

#include <iostream>
#define N 100
#define numBlocks 1
#define threadsPerBlock 100

__global__ void add( int *a, int *b, int *c ) 
{
    int i = threadIdx.x;
    c[i]=a[i]+b[i];  
}

int main(void )
{
    //處理記憶體的部分,待補
    add<<<numBlocks, threadsPerBlock>>>(a, b, c);
    //處理記憶體的部分,待補
}

寫到這裡, 其實程式還是不會動, 要讓程式會動還要處理記憶體的部份,沒錯大家害怕的指標(pointer)要出現惹, 當然還有超級煩人的動態記憶體配置, malloc!!!

在CUDA裡面, 分為host和device, host一般指得是用CPU處理, device一般指得是用GPU處理, 彼此間的記憶體並不共用, 因此在運算的時候, 要把host記憶體的東西複製給device, 算完後再複製回來, 在這裡為了方便辨識, host的變數都加上h_, device的變數都加上d_, 程式碼如下:

#include <iostream>
#include <cstdlib> 
#include <ctime>   
#define N 100
#define numBlocks 1
#define threadsPerBlock 100

using namespace std;

__global__ void add( int *a, int *b, int *c ) 
{
    int i = threadIdx.x;
    c[i]=a[i]+b[i];  
}

int main( void ) {
    
    int *h_a, *h_b, *h_c;   
    int *d_a, *d_b, *d_c;   
    
    // 在host分配記憶體空間
    h_a = (int*)malloc( N * sizeof(int) ); 
    h_b = (int*)malloc( N * sizeof(int) );
    h_c = (int*)malloc( N * sizeof(int) );

    // 亂數指定值給陣列
    for (int i=0; i<N; i++) {
        h_a[i] = rand();
        h_b[i] = rand();
    }

    // 在device上分配記憶體
     cudaMalloc( (void**)&d_a, N * sizeof(int) );
     cudaMalloc( (void**)&d_b, N * sizeof(int) );
     cudaMalloc( (void**)&d_c, N * sizeof(int) );

    // 把host的東西複製到device上
     cudaMemcpy( d_a, h_a, N * sizeof(int), cudaMemcpyHostToDevice );
     cudaMemcpy( d_b, h_b, N * sizeof(int), cudaMemcpyHostToDevice );

    
    add<<<numBlocks, threadsPerBlock>>>( d_a, d_b, d_c );

    // 把device算出來的結果複製回CPU
    cudaMemcpy( h_c, d_c, N * sizeof(int),cudaMemcpyDeviceToHost );

    //驗算
    int ok=1;
    for(int i=0;i<N;i++)
    {
        if(h_c[i]!=h_a[i]+h_b[i])
        {
            cout<<"fail";
            ok=0;
            break;
        }
    }
    if(ok==1)
        printf("success");
   
    // 釋放host的記憶體
    free( h_a );
    free( h_b );
    free( h_c );

    // 釋放device的記憶體
     cudaFree( d_a );
     cudaFree( d_b );
     cudaFree( d_c );

    return 0;
}

reference:
1. Nvidia developer zone
2. http://selkie.macalester.edu/csinparallel/modules/ConceptDataDecomposition/build/html/Decomposition/CUDA_VecAdd.html