2010年11月20日 星期六

[轉錄] CUDA速成篇(中)

延續上一篇...




※ [本文轉錄自PTT VideoCard 看板]

作者: a5000ml (咖啡裏的海洋藍) 看板: VideoCard
標題: [分享] CUDA 程式設計(10) -- 速成篇(中)


============================================================================
第五招  記憶體 (主機、裝置、共享、暫存器)
============================================================================
基本的記憶體部份,最重要的是區分主機、裝置、共享記憶體以及暫存器,硬體位置
和存取速度列於下表:

    -------------------------------------------------------------------------------------------------
    中文名稱       英文名稱               硬體位置                    存取速度
    -------------------------------------------------------------------------------------------------
    主機記憶體   Host Memory         PC 上的 DRAM        透過 PCI-E, 很慢
    裝置記憶體   Device Memory    顯示卡上的 DRAM    400~600 cycles
    共享記憶體   Shared Memory    顯示晶片                    4 cycles issue
    暫存器           Register                 顯示晶片                     立即
   -------------------------------------------------------------------------------------------------

這些記憶體的大小、功能和使用方式如下

    --------------------------------------------------------------------------------------------------
    名稱               大小                           功能                                 使用方式
    --------------------------------------------------------------------------------------------------
    主機記憶體   0.5~10GB                 存放傳統的主機資料     透過 API
    裝置記憶體   0.5~10GB                 存放給 GPU 使用資料   kernel 直接存取
    共享記憶體   16KB/BLOCK          執行緒之間交換資料      kernel 直接存取
    暫存器            32~64KB/BLOCK  執行緒的區域變數          kernel 直接存取
    ---------------------------------------------------------------------------------------------------

相關細節如下:

(1) 【主機記憶體】我們都很熟了,就是傳統 C/C++ 使用的變數,可以透過 malloc()、
    free()、new、delete 等來配置或釋放。

(2) 【裝置記憶體】的地位和主機記憶體很像,隻是主機記憶體是對付 CPU,而
    裝置記憶體是對付 GPGPU,兩者間的資料傳送要透過 CUDA 的 API 達成,
    可以透過 cudaMalloc()、cudaFree() 等函式來配置或釋放 (詳見第二招),
    在 CUDA 中,這類的記憶體稱為【全域記憶體 (global memory)】。

(3) 【共享記憶體】是比較特殊的,在傳統的序列化程式設計裏沒有直接的對應,
    初學者可能要花多一點的時間在這上面,它用在『區塊內執行緒間交換資料』,
    隻能在 kernel 中使用,宣告時使用 __shared__ 這個標籤,使用時必需同步,
    以確保資料讀寫時序上的正確性,現階段在每個區塊上最大的容量為 16KB
    超過便無法執行或編譯,因為 on chip 的關係,它屬於快速記憶體。

(4) 【暫存器】kernel 中大部份的區域變數都是以暫存器的型式存放,不需做額外的
    宣告手序,這些暫存器是區塊中執行緒共享的,也就是如果每個執行緒使用到 8 個
    暫存器,呼叫這個 kernel 時區塊大小指定為 10,則整個區塊使用到 8*10=80 個
    暫存器,若呼叫時指定區塊大小為 50,則整個區塊使用 8*50=400 個暫存器。

(5) 當使用過多的暫存器時 (>32~64KB/BLOCK,看是那個世代的 GPU),系統會自動把
    一些資料置換到全域記憶體中,導致執行緒變多,但效率反而變慢 (類似作業系統
    虛擬記憶體的 swap);另一個會引發這種 swap 的情況是在使用動態索引存取陣列,
    因為此時需要陣列的順序性,而暫存器本身是沒有所謂的順序的,所以系統會自動
    把陣列置於全域記憶體中,再按索引存取,這種情況建議使用共享記憶體手動避免。

(6) 【暫存器】和【共享記憶體】的使用量會限制執行緒的數目,在開發複雜 kernel 時
    宜註意,可使用 nvcc --ptxas-options=-v 這個選項在設計時期監控,或使用
    nvcc --maxrregcount 選項限制每個執行緒的暫存器使用量。


//----------------------------------------------------------------------------
//範例(5): 平滑處理 (使用相鄰的三點做加權平均,使資料變平滑)[081119-smooth.cu]
//        執行緒同步 __syncthreads() 和 cudaThreadSynchronize(),詳見第六招
//---------------------------------------------------------------------------
#include<stdio.h>
#include<time.h>
#include<cuda.h>

//設定區塊大小 (shared 版本會用到, 所以先宣告).
#define BLOCK 512

//--------------------------------------------------------
//(1) 對照組 (host 版).
//--------------------------------------------------------
void smooth_host(float* b, float* a, int n){
        for(int k=1; k<n-1; k++){
                b[k]=(a[k-1]+2*a[k]+a[k+1])*0.25;
        }

        //邊界為0
        b[0]=(2*a[0]+a[1])*0.25;
        b[n-1]=(a[n-2]+2*a[n-1])*0.25;
}

//--------------------------------------------------------
//(2) 裝置核心(global 版).
//--------------------------------------------------------
__global__ void smooth_global(float* b, float* a, int n){
        int k = blockIdx.x*blockDim.x+threadIdx.x;

        if(k==0){
                b[k]=(2*a[0]+a[1])*0.25;
        }
        else if(k==n-1){
                b[k]=(a[n-2]+2*a[n-1])*0.25;
        }
        else if(k<n){
                b[k]=(a[k-1]+2*a[k]+a[k+1])*0.25;
        }

}


//--------------------------------------------------------
//(3) 裝置核心(shared 版).
//--------------------------------------------------------
__global__ void smooth_shared(float* b, float* a, int n){
    //----------------------------------------
    //計算區塊的基底
    //----------------------------------------
        int base = blockIdx.x*blockDim.x;
        int t = threadIdx.x;

    //----------------------------------------
    //宣告共享記憶體.
    //----------------------------------------
        __shared__ float s[BLOCK+2];

    //----------------------------------------
    //載入主要資料 s[1]~s[BLOCK]
    //----------------------------------------
    // s[0] <-- a[base-1]  (左邊界)
    // s[1] <-- a[base]
    // s[2] <-- a[base+1]
    // s[3] <-- a[base+2]
    //      ...
    // s[BLOCK]   <-- a[base+BLOCK-1]
    // s[BLOCK+1] <-- a[base+BLOCK]  (右邊界)
    //----------------------------------------
        if(base+t<n){
                s[t+1]=a[base+t];
        }

    //----------------------------------------
    //載入邊界資料 s[0] & s[BLOCK+1] (隻用兩個執行緒處理)
    //----------------------------------------
        if(t==0){
                //左邊界.
                if(base==0){
                        s[0]=0;
                }
                else{
                        s[0]=a[base-1];
                }
        }

        //*** 使用獨立的 warp 讓 branch 更快 ***
        if(t==32){
                //右邊界.
                if(base+BLOCK>=n){
                        s[n-base+1]=0;
                }
                else{
                        s[BLOCK+1] = a[base+BLOCK];
                }
        }

    //----------------------------------------
    //同步化 (確保共享記憶體已寫入)
    //----------------------------------------
        __syncthreads();

    //----------------------------------------
    //輸出三點加權平均值
    //----------------------------------------
        if(base+t<n){
                b[base+t]=(s[t]+2*s[t+1]+s[t+2])*0.25;
        }

};


//--------------------------------------------------------
//主函式.
//--------------------------------------------------------
int main(){
    //--------------------------------------------------
    //參數.
    //--------------------------------------------------
        int num=10*1000*1000;
        int loop=130;  //測試迴圈次數 (量時間用)

    //--------------------------------------------------
    //配置主機記憶體.
    //--------------------------------------------------
        float* a=new float[num];
        float* b=new float[num];
        float* bg=new float[num];
        float* bs=new float[num];

    //--------------------------------------------------
    //配置裝置記憶體.
    //--------------------------------------------------
        float *GA, *GB;
        cudaMalloc((void**) &GA, sizeof(float)*num);
        cudaMalloc((void**) &GB, sizeof(float)*num);

    //--------------------------------------------------
    //初始化(亂數) & 複制資料到顯示卡的 DRAM.
    //--------------------------------------------------
        for(int k=0; k<num; k++){
                a[k]=(float)rand()/RAND_MAX;
                b[k]=bg[k]=bs[k]=0;
        }
        cudaMemcpy(GA, a, sizeof(float)*num, cudaMemcpyHostToDevice);

    //--------------------------------------------------
    //Test(1): smooth_host
    //--------------------------------------------------
        double t_host=(double)clock()/CLOCKS_PER_SEC;
        for(int k=0; k<loop; k++){
                smooth_host(b,a,num);
        }
        t_host=((double)clock()/CLOCKS_PER_SEC-t_host)/loop;


    //--------------------------------------------------
    //Test(2): smooth_global (GRID*BLOCK 必需大於 num).
    //--------------------------------------------------
        double t_global=(double)clock()/CLOCKS_PER_SEC;
        cudaThreadSynchronize();
        for(int k=0; k<loop; k++){
                smooth_global<<<num/BLOCK+1,BLOCK>>>(GB,GA,num);
        }
        cudaThreadSynchronize();
        t_global=((double)clock()/CLOCKS_PER_SEC-t_global)/loop;

        //下載裝置記憶體內容到主機上.
        cudaMemcpy(bg, GB, sizeof(float)*num, cudaMemcpyDeviceToHost);

    //--------------------------------------------------
    //Test(3): smooth_shared (GRID*BLOCK 必需大於 num).
    //--------------------------------------------------
        double t_shared=(double)clock()/CLOCKS_PER_SEC;
        cudaThreadSynchronize();
        for(int k=0; k<loop; k++){
                smooth_shared<<<num/BLOCK+1,BLOCK>>>(GB,GA,num);
        }
        cudaThreadSynchronize();
        t_shared=((double)clock()/CLOCKS_PER_SEC-t_shared)/loop;

        //下載裝置記憶體內容到主機上.
        cudaMemcpy(bs, GB, sizeof(float)*num, cudaMemcpyDeviceToHost);

    //--------------------------------------------------
    //比較正確性
    //--------------------------------------------------
        double sum_dg2=0, sum_ds2=0, sum_b2=0;
        for(int k=0; k<num; k++){
                double dg=bg[k]-b[k];
                double ds=bs[k]-b[k];

                sum_b2+=b[k]*b[k];
                sum_dg2+=dg*dg;
                sum_ds2+=ds*ds;
        }

    //--------------------------------------------------
    //報告
    //--------------------------------------------------
        //組態.
        printf("vector size: %d \n", num);
        printf("\n");

        //時間.
        printf("Smooth_Host:   %g ms\n", t_host*1000);
        printf("Smooth_Global: %g ms\n", t_global*1000);
        printf("Smooth_Shared: %g ms\n", t_shared*1000);
        printf("\n");

        //相對誤差.
        printf("Diff(Smooth_Global): %g \n", sqrt(sum_dg2/sum_b2));
        printf("Diff(Smooth_Shared): %g \n", sqrt(sum_ds2/sum_b2));
        printf("\n");

    //--------------------------------------------------
    //釋放裝置記憶體.
    //--------------------------------------------------
        cudaFree(GA);
        cudaFree(GB);
        delete [] a;
        delete [] b;
        delete [] bg;
        delete [] bs;
        return 0;
}


//--------------------------------------------------------
//範例(5): 執行結果 (測試 10M 個 float)
//--------------------------------------------------------


          vector size: 10000000

          Smooth_Host:   36.9231 ms
          Smooth_Global: 14.4615 ms
          Smooth_Shared: 5.07692 ms

          Diff(Smooth_Global): 3.83862e-08
          Diff(Smooth_Shared): 3.83862e-08


(1) 這次測試的機器比較爛: P4-3.2 (prescott) vs. 9600GT
    不過我們仍可看到共享記憶體使載入的資料量變為 1/3 所得到的增速

(2) 在 smooth_shared() 裏我們用 2 個 warp 使得條件判斷可以獨立,
    如果第 2 個 if(t==32) 改成 if(t==16) 或其他小於 32 的值,
    也就是和第 1 個 if 使用同一個 warp, 則速度會變慢, 有興趣的朋友
    可以去試試看,warp 不打算在新手篇講,之後硬體時才詳細討論。

(3) 測試效能時使用 cudaThreadSynchronize() 同步主機和裝置核心,
    以免量到錯誤的時間

(4) 在 smooth_shared() 裏使用 __syncthreads() 同步化執行緒,
    以免在計算 output 時仍有共享記憶體還沒完成寫入動作,
    卻有執行緒已經需要使用它的資料。


============================================================================
第六招  執行緒同步 (網格、區塊)
============================================================================
同步執行緒有兩個函式,分別是 __syncthreads() 和 cudaThreadSynchronize()

-----------------------------------------------------------------
同步化函式                 使用地點     功能
-----------------------------------------------------------------
__syncthreads()            核心程序中   同步化【區塊內的執行緒】
cudaThreadSynchronize()    主機程序中   同步化【核心和主機程序】
-----------------------------------------------------------------

(1) 在 kernel 中,使用 __syncthreads() 來進行區塊內的執行緒的同步,
    避免資料時序上的問題 (來自不同 threads),時常和共享記憶體一起使用,
    在範例(5)中示範了使用 __syncthreads() 來隔開共享記憶體的【寫入週期】
    和【讀取週期】,避免 WAR 錯誤 (write after read)。

(2) 在主機程序中,使用 cudaThreadSynchronize() 來進行核心和主機程序的同步,
    範例(5)中示範了用它來避免量到不正確的主機時間 (kernel仍未完成就量時間),
    因為主機的程序和裝置程序預設是不同步的 (直到下載結果資料之前),這個 API
    可以強迫它們同步。
--