1. 程式人生 > >《GPU高效能程式設計CUDA實戰》學習筆記(九)

《GPU高效能程式設計CUDA實戰》學習筆記(九)

第9章 原子性

在某些情況下,對於單執行緒應用程式來說非常簡單的任務,或許使用大規模的並行架構實現卻會變成一個複雜的問題。這裡我們將在這些情況中使用特殊的原語從而確保安全地完成傳統單執行緒應用程式中的簡單任務。

9.1 本章目標

  • 瞭解不同NVIDIA GPU的計算功能集。
  • 瞭解原子操作以及為什麼需要使用它們。
  • 瞭解如何在CUDA C核函式中執行帶有原子操作的運算。

9.2 計算功能集

1-8章,介紹的是所以支援CUDA 的GPU的通用功能。如,啟動核函式、訪問全域性記憶體、讀取常量記憶體和紋理記憶體。 但不同架構有不同功能,NVIDIA 將GPU支援的各種功能統稱為計算功能集(Compute Capability)。

9.2.1 NVIDIA GPU的計算功能集

目前支援1.0、1.1、1.2、1.3以及2.0,以後可能支援更高版本,需要從官網來了解。高版本計算功能集是低版本計算功能集的超集。如支援1.2版本的GPU,同樣會支援1.0和1.1版本的所有功能。NVIDIA CUDA程式設計指南中包含所以最新列表及計算功能集。

9.2.2 基於最小計算功能集的編譯

假設在編寫的程式碼中要求計算功能集的版本最低不能低於某個版本。例如,假設你閱讀完本章,並開始編寫一個需要使用全域性記憶體原子操作的應用程式。你知道要支援全域性記憶體原子操作,計算功能集的最低版本為1.1。當編譯程式碼時,你需要告訴編譯器,如果硬體支援的計算功能集低於1.1,那麼將無法執行這個核函式。而且,當告訴編譯器這個要求時,還可以指定一些只在1.1或者更高版本的計算功能集中才支援的編譯優化。要將這個資訊告訴編譯器,只需在呼叫 nvcc 時增加一個命令選項: nvcc -arch = sm_11 同樣的,在編譯需要使用共享記憶體原子操作的核函式時,你要告訴編譯器程式碼需要1.2版本或者更高的計算功能集。 nvcc -arch = sm_12

9.2.3 原子操作簡介

在編寫傳統的單執行緒應用程式時,程式設計師通常不需要使用原子操作。但有時也是需要的。 C++ 遞增運算子:
x++;
步驟:1)讀取x值;2)增加1;3)遞增後結果寫回x。讀取-修改-寫入(Read-Modify-Write) 多執行緒時會出現混亂,結果不可預測。所以我們需要將這3步變成1個不可分割為更小的操作,滿足這種條件限制的操作稱為原子操作。 CUDA C支援多種原子操作,當有數千個執行緒在記憶體訪問上發生競爭時,這些操作能夠確保在記憶體上實現安全的操作。 現在我們已經看到一個只有使用原子操作才能計算出正確結果的示例。

9.4 計算直方圖

直方圖(Histogram)又稱質量分佈圖。是一種統計報告圖,由一系列高度不等的縱向條紋或線段表示資料分佈的情況。 一般用橫軸表示資料型別,縱軸表示分佈情況。
給定一個包含一組元素的資料集,直方圖表示每個元素的出現頻率。例如,“Programming with CUDA C”中字元頻率的直方圖,結果如下 2 2 1 2 1 2 2 1 1 1 2 1 1 1 a c d g h i m n o p r t u w 直方圖定義簡單,但應用廣泛,包括影象處理、資料壓縮、計算機視覺、機器學習、音訊編碼等等。

9.4.1 在 CPU 上計算直方圖

#include "../common/book.h"

#define SIZE    (100*1024*1024)

int main( void ) {
    unsigned char *buffer =
                     (unsigned char*)big_random_block( SIZE );

    // capture the start time
    clock_t         start, stop;
    start = clock();

    unsigned int    histo[256];
    for (int i=0; i<256; i++)
        histo[i] = 0;

    for (int i=0; i<SIZE; i++)
        histo[buffer[i]]++;

    stop = clock();
    float   elapsedTime = (float)(stop - start) /
                          (float)CLOCKS_PER_SEC * 1000.0f;
    printf( "Time to generate:  %3.1f ms\n", elapsedTime );

    long histoCount = 0;
    for (int i=0; i<256; i++) {
        histoCount += histo[i];
    }
    printf( "Histogram Sum:  %ld\n", histoCount );

    free( buffer );
    return 0;
}

big_random_block() 生成隨機的位元組流。
1位元組 = 8位元 思想: 每當buffer中出現某個z時,就遞增直方圖陣列中索引為z的元素。這樣就計算出z的出現次數。 後邊的程式碼是驗證個數是否正確。

9.4.2 在GPU上計算直方圖

#include "../common/book.h"

#define SIZE    (100*1024*1024)


__global__ void histo_kernel( unsigned char *buffer,
                              long size,
                              unsigned int *histo ) {
    // calculate the starting index and the offset to the next
    // block that each thread will be processing
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    int stride = blockDim.x * gridDim.x;
    while (i < size) {
        atomicAdd( &histo[buffer[i]], 1 );
        i += stride;
    }
}

int main( void ) {
    unsigned char *buffer =
                     (unsigned char*)big_random_block( SIZE );

    // capture the start time
    // starting the timer here so that we include the cost of
    // all of the operations on the GPU.
    cudaEvent_t     start, stop;
    HANDLE_ERROR( cudaEventCreate( &start ) );
    HANDLE_ERROR( cudaEventCreate( &stop ) );
    HANDLE_ERROR( cudaEventRecord( start, 0 ) );

    // allocate memory on the GPU for the file's data
    unsigned char *dev_buffer;
    unsigned int *dev_histo;
    HANDLE_ERROR( cudaMalloc( (void**)&dev_buffer, SIZE ) );
    HANDLE_ERROR( cudaMemcpy( dev_buffer, buffer, SIZE,
                              cudaMemcpyHostToDevice ) );

    HANDLE_ERROR( cudaMalloc( (void**)&dev_histo,
                              256 * sizeof( int ) ) );
    HANDLE_ERROR( cudaMemset( dev_histo, 0,
                              256 * sizeof( int ) ) );

    // kernel launch - 2x the number of mps gave best timing
    cudaDeviceProp  prop;
    HANDLE_ERROR( cudaGetDeviceProperties( &prop, 0 ) );
    int blocks = prop.multiProcessorCount;
    histo_kernel<<<blocks*2,256>>>( dev_buffer, SIZE, dev_histo );
    
    unsigned int    histo[256];
    HANDLE_ERROR( cudaMemcpy( histo, dev_histo,
                              256 * sizeof( int ),
                              cudaMemcpyDeviceToHost ) );

    // get stop time, and display the timing results
    HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
    HANDLE_ERROR( cudaEventSynchronize( stop ) );
    float   elapsedTime;
    HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime,
                                        start, stop ) );
    printf( "Time to generate:  %3.1f ms\n", elapsedTime );

    long histoCount = 0;
    for (int i=0; i<256; i++) {
        histoCount += histo[i];
    }
    printf( "Histogram Sum:  %ld\n", histoCount );

    // verify that we have the same counts via CPU
    for (int i=0; i<SIZE; i++)
        histo[buffer[i]]--;
    for (int i=0; i<256; i++) {
        if (histo[i] != 0)
            printf( "Failure at %d!  Off by %d\n", i, histo[i] );
    }

    HANDLE_ERROR( cudaEventDestroy( start ) );
    HANDLE_ERROR( cudaEventDestroy( stop ) );
    cudaFree( dev_histo );
    cudaFree( dev_buffer );
    free( buffer );
    return 0;
}

經驗告訴我們: 當執行緒塊數量為GPU中處理器數量的2倍時,將達到最優效能。 1. 使用全域性記憶體原子操作的直方圖核函式 atomicAdd(addr, y)將生成一個原子的操作系列。 包括,讀取地址addr處的值,將y增加到這個值,以及將結果儲存回地址addr。底層硬體確保addr安全。
不過此時的效率比cpu更低,因為對相同記憶體位置的操作都將被硬體序列化,這導致儲存未完成操作的佇列非常長,抵消了並行效果。 2. 使用共享記憶體原子操作和全域性記憶體原子操作的直方圖核函式 解決上述問題,全部程式碼如下:
#include "../common/book.h"

#define SIZE    (100*1024*1024)


__global__ void histo_kernel( unsigned char *buffer,
                              long size,
                              unsigned int *histo ) {

    // clear out the accumulation buffer called temp
    // since we are launched with 256 threads, it is easy
    // to clear that memory with one write per thread
    __shared__  unsigned int temp[256];
    temp[threadIdx.x] = 0;
    __syncthreads();

    // calculate the starting index and the offset to the next
    // block that each thread will be processing
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    int stride = blockDim.x * gridDim.x;
    while (i < size) {
        atomicAdd( &temp[buffer[i]], 1 );
        i += stride;
    }
    // sync the data from the above writes to shared memory
    // then add the shared memory values to the values from
    // the other thread blocks using global memory
    // atomic adds
    // same as before, since we have 256 threads, updating the
    // global histogram is just one write per thread!
    __syncthreads();
    atomicAdd( &(histo[threadIdx.x]), temp[threadIdx.x] );
}

int main( void ) {
    unsigned char *buffer =
                     (unsigned char*)big_random_block( SIZE );

    // capture the start time
    // starting the timer here so that we include the cost of
    // all of the operations on the GPU.  if the data were
    // already on the GPU and we just timed the kernel
    // the timing would drop from 74 ms to 15 ms.  Very fast.
    cudaEvent_t     start, stop;
    HANDLE_ERROR( cudaEventCreate( &start ) );
    HANDLE_ERROR( cudaEventCreate( &stop ) );
    HANDLE_ERROR( cudaEventRecord( start, 0 ) );

    // allocate memory on the GPU for the file's data
    unsigned char *dev_buffer;
    unsigned int *dev_histo;
    HANDLE_ERROR( cudaMalloc( (void**)&dev_buffer, SIZE ) );
    HANDLE_ERROR( cudaMemcpy( dev_buffer, buffer, SIZE,
                              cudaMemcpyHostToDevice ) );

    HANDLE_ERROR( cudaMalloc( (void**)&dev_histo,
                              256 * sizeof( int ) ) );
    HANDLE_ERROR( cudaMemset( dev_histo, 0,
                              256 * sizeof( int ) ) );

    // kernel launch - 2x the number of mps gave best timing
    cudaDeviceProp  prop;
    HANDLE_ERROR( cudaGetDeviceProperties( &prop, 0 ) );
    int blocks = prop.multiProcessorCount;
    histo_kernel<<<blocks*2,256>>>( dev_buffer,
                                    SIZE, dev_histo );
    
    unsigned int    histo[256];
    HANDLE_ERROR( cudaMemcpy( histo, dev_histo,
                              256 * sizeof( int ),
                              cudaMemcpyDeviceToHost ) );

    // get stop time, and display the timing results
    HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
    HANDLE_ERROR( cudaEventSynchronize( stop ) );
    float   elapsedTime;
    HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime,
                                        start, stop ) );
    printf( "Time to generate:  %3.1f ms\n", elapsedTime );

    long histoCount = 0;
    for (int i=0; i<256; i++) {
        histoCount += histo[i];
    }
    printf( "Histogram Sum:  %ld\n", histoCount );

    // verify that we have the same counts via CPU
    for (int i=0; i<SIZE; i++)
        histo[buffer[i]]--;
    for (int i=0; i<256; i++) {
        if (histo[i] != 0)
            printf( "Failure at %d!\n", i );
    }

    HANDLE_ERROR( cudaEventDestroy( start ) );
    HANDLE_ERROR( cudaEventDestroy( stop ) );
    cudaFree( dev_histo );
    cudaFree( dev_buffer );
    free( buffer );
    return 0;
}
使用共享記憶體原子操作和全域性原子性,效能將數量級的提升。  這裡主要使用了重構演算法,分成兩個階段計算,降低了記憶體訪問上的競爭程度,帶來了不錯的效果。以後用得到哦,要記住這種策略。