CUDA學習筆記(5) 原子操作

原子操作是指，當一個線程（Thread）要對同一個顯存變量依次進行“讀-計算-寫”的操作時，這個“讀-計算-寫”的操作必須連貫地執行，中間不能插入任何其他操作。

舉個例子，假設我們想要用GPU統計“char data_0[32] = {1,0, … ,1}”這個數組中“0”和“1”的個數并寫入“int counter[2]”中。

如果我們不使用原子操作，直接在核函數中這樣寫：

extern "C" __global__ void kernel_func(int * counter, char * data_0)
{
    // 計算線程号
    unsigned int block_index = blockIdx.x + blockIdx.y * gridDim.x + blockIdx.z * gridDim.x * gridDim.y;
    unsigned int thread_index = block_index * blockDim.x * blockDim.y * blockDim.z + \
        threadIdx.x + threadIdx.y * blockDim.x + threadIdx.z * blockDim.x * blockDim.y;

    // 統計結果
    int value = data_0[thread_index];
    counter[value] ++;
}

我們會發現結果是“counter[2] = {1, 1}”，這顯然不是正确的結果。

CUDA學習筆記(5) 原子操作

在GPU中，線程在并行時，每個線程（Thread）都會有一組專供自己的用的寄存器。例如上面的代碼中，我們在核函數中定義了“block_index”、“thread_index”、“value”三個寄存器，此外編譯器還會自動生成一個“counter”寄存器。也就是說，當這32個線程（Thread）并行執行的時候，會占用 4×32 個32位寄存器。

*注1：自動生成一個“counter”寄存器的原因和CPU中彙編的的原理相同，由于“counter[value]”是存放在顯存中的變量，對他進行“++”操作時會先将它讀取到寄存器（Register）中，對這個寄存器變量進行操作後再将其寫入原來的顯存變量。如此便會在執行過程中占用一個寄存器。

CUDA學習筆記(5) 原子操作

原因是這樣的，在GPU中每32個線程（Thread）作為一個線程束（Warp）整體執行一系列操作。上圖中的執行過程是這樣的：

線程束0（即線程0到31），從全局顯存中讀取了數組“data_0”的數值放在每個線程對應的寄存器“value”中。
線程0到31幾乎同時分别占用一個SP，然後第 i 條線程根據各自的寄存器“value”中的數值，準備讀取顯存變量“counter[value]”的值。（此時“counter[2] = {0, 0}”）
線程束0讀取顯存變量“counter[value]”到各個線程的寄存器“counter”中。
線程0到31同時對自己的“value”寄存器中的數值執行“++”。
線程束0将寄存器“counter”值寫入顯存變量“counter[value]”中。（此時“counter[2] = {1, 1}”）

是以，得到的結果是“counter[2] = {1, 1}”。這樣，我們可以看到，錯誤實際上是出在了并行上——線程0還沒有将自己計算的“counter”寫回顯存變量“counter[value]”，其他線程就已經讀取了顯存變量“counter[value]”的值。

*注2：上面的例子作為簡單情況分析，例中隻有32條線程（Thread），其數量小于空閑的流處理器（SP）數量（我的電腦上384個）時是這樣的結果，每個線程都由一個流處理器（SP）來處理。線上程較多時可能多個線程都由一個流處理器（SP）處理。

正确的方法是使用原子操作，在核函數中這樣寫：

extern "C" __global__ void kernel_func(int * counter, char * data_0)
{
    // 計算線程号
    unsigned int block_index = blockIdx.x + blockIdx.y * gridDim.x + blockIdx.z * gridDim.x * gridDim.y;
    unsigned int thread_index = block_index * blockDim.x * blockDim.y * blockDim.z + \
        threadIdx.x + threadIdx.y * blockDim.x + threadIdx.z * blockDim.x * blockDim.y;

    // 統計結果
    int value = data_0[thread_index];
    atomicAdd(&counter[value], );
}

注意“atomicAdd”函數，它的作用在于當線程k要操作顯存變量“counter[0]”的值時，其他的線程若需要操作“counter[0]”就必須等待線程 k 操作結束，才能再對“counter[0]”執行操作。程式的執行步驟如下：

線程束0（即線程0到31），從全局顯存中讀取了數組“data_0”的數值放在每個線程對應的寄存器“value”中。
某個線程（記作線程k）檢查顯存變量“counter[value]”是否被其他線程的原子操作鎖定，如果已經被鎖定則等待，否則讀取顯存變量“counter[value]”到線程 k 的寄存器“counter”中并鎖定顯存變量“counter[value]”。僞代碼如下：

如果（顯存變量“counter[value]”被鎖定）

等待顯存變量“counter[value]”解鎖。

否則

讀取顯存變量“counter[value]”到寄存器“counter”中，并鎖定顯存變量“counter[value]”。
線程k對自己的“value”寄存器中的數值執行“++”。
線程 k <script type="math/tex" id="MathJax-Element-1600">k</script>将“value”寄存器中的數值寫回顯存變量“counter[value]”中，并解除對顯存變量“counter[value]”的鎖定。

這樣，我的就可以正确地統計得到結果了。然而，我們的程式并行程度也是以有所下降。假如“data_0”中的值全部為同一個值，那麼這個程式就變成了一個幾乎串行的程式，大大降低了程式的運作效率，當然這也是有辦法改進的。後面我會使用共享顯存（Shared memory）進行改進。

*注3：僞代碼步驟2中的判斷是原子操作在硬體層自動實作的，在編寫程式時，“atomicAdd(&counter[value], 1);”即可達到效果。

下面是各種原子操作的清單：

函數名	作用
atomicAdd(&value, add_num)	加法：value = value + add_num
atomicSub(&value, sub_num)	減法：value = value + sub_num
atomicExch(&value, num)	指派：value = num
atomicMax(&value, num)	求最大：value = max(value, num)
atomicMin(&value, num)	求最小：value = min(value, num)
atomicInc(&value, compare)	向上計數：如果(value <= compare)則 `value++` ,否則 `value = 0`
atomicDec(&value, compare)	向下計數：如果(value > compare或value == 0), 則 `value--` ,否則 `value = 0`
atomicCAS(&value, compare)	比較并交換：如果(value != compare),則 `value = compare`
atomicAnd(&value, add_num)	與運算：value = value & num
atomicOr(&value, add_num)	或運算 value = value
atomicXor(&value, add_num)	異或運算 value = value ^ num

以上這些原子操作函數的傳回值均為被修改前的“value”值。

我的“kernel.cu”代碼：

// This is kernel function file !

extern "C" __global__ void kernel_func(unsigned int * thread_index_array, unsigned char *source_data_array, \
    unsigned int * result_array, unsigned int *clock_counter_array)
{
    // 計算線程号
    unsigned int block_index = blockIdx.x + blockIdx.y * gridDim.x + blockIdx.z * gridDim.x * gridDim.y;
    unsigned int thread_index = block_index * blockDim.x * blockDim.y * blockDim.z + \
        threadIdx.x + threadIdx.y * blockDim.x + threadIdx.z * blockDim.x * blockDim.y;
    // 總線程數
    unsigned int thread_num = gridDim.x * gridDim.y * gridDim.z * blockDim.x * blockDim.y * blockDim.z;

    // 記錄線程号
    thread_index_array[thread_index] = thread_index;

    clock_counter_array[thread_index] = clock();

    unsigned int value = source_data_array[thread_index];

    // 正确地統計方法
    atomicAdd(&result_array[value], );
    // 錯誤的統計方法
    //result_array[value] ++;

    clock_counter_array[thread_index + thread_num] = clock();
}

我的“learn_CUDA_05.cpp”代碼：

// CUDA Atomic

// C/C++ IO
#include <stdio.h>
#include <iostream>
using namespace std;

// For the CUDA runtime routines (prefixed with "cuda_")
#include <cuda_runtime.h>
#include <cuda.h>

// helper functions and utilities to work with CUDA
#include <helper_cuda.h>
#include <helper_functions.h>
// NVRTC
#include <nvrtc_helper.h>

// timer
#include <timer.h>
#include <windows.h>


// GPU核心時鐘頻率
int GPU_clock_rate;
// Block最大線程數
int max_thread_per_block;
// 檢查顯示卡硬體屬性
void check_Cuda_information(int main_argc, char ** main_argv);
// 輸出PTX代碼
void print_PTX(char * PTX_ptr, int size_of_PTX);


int main(int argc, char **argv)
{
    // Host變量（記憶體變量）
    const int ARRAY_LENGTH = , STATIS_size = UCHAR_MAX + ;
    unsigned int *thread_index_array;
    unsigned char *source_data_array;
    unsigned int *clock_counter_array, *result_array;
    float *cost_time_array;

    // Device變量（顯存變量）
    CUdeviceptr dev_thread_index_array, dev_source_data_array, dev_clock_counter_array, dev_result_array;

    // 檢查顯示卡硬體屬性
    check_Cuda_information(argc, &argv[]);

#pragma region(CPU硬體計時器 開始計時)
    //用QueryPerformanceCounter()來計時  微秒
    LARGE_INTEGER  large_interger;
    double CPU_counter_frequency;
    __int64
        CPU_c1_start, CPU_c2_compile, CPU_c3_alloc, CPU_c4_caculate, CPU_c5_copyout, CPU_c6_end;

    QueryPerformanceFrequency(&large_interger);
    CPU_counter_frequency = large_interger.QuadPart;
    printf("CPU硬體計數器頻率：\t%.2lf Hz\r\n\r\n", CPU_counter_frequency);

    // 計數
    QueryPerformanceCounter(&large_interger);
    CPU_c1_start = large_interger.QuadPart;
#pragma endregion


#pragma region(讀取、編譯CUDA核函數)
    // 讀取、編譯、加載CUDA核函數
    char *ptx, *kernel_file;
    size_t ptxSize;

    // 打開核函數的 .cu 檔案 并編譯為PTX
    kernel_file = sdkFindFilePath("kernel.cu", argv[]);
    // 如果找不到 kernel_file 會報錯：error: unable to open ./xxxx.cu for reading!
    compileFileToPTX(kernel_file, NULL, NULL, &ptx, &ptxSize, );
    // 輸出PTX代碼
    print_PTX(ptx, (int)ptxSize);


    // 選擇GPU裝置，加載核函數到GPU裝置
    CUmodule module = loadPTX(ptx, argc, argv);
    CUfunction kernel_addr;
    checkCudaErrors(cuModuleGetFunction(&kernel_addr, module, "kernel_func"));

    // CPU 計數
    QueryPerformanceCounter(&large_interger);
    CPU_c2_compile = large_interger.QuadPart;
#pragma endregion


#pragma region(記憶體、顯存配置設定/初始化)
    thread_index_array = (unsigned int*)malloc(ARRAY_LENGTH * sizeof(unsigned int));
    source_data_array = (unsigned char*)malloc(ARRAY_LENGTH * sizeof(unsigned char));
    result_array = (unsigned int*)malloc(STATIS_size * sizeof(unsigned int));
    clock_counter_array = (unsigned int*)malloc( * ARRAY_LENGTH * sizeof(unsigned int));
    cost_time_array = (float*)malloc(ARRAY_LENGTH * sizeof(float));

    checkCudaErrors(cuMemAlloc(&dev_thread_index_array, ARRAY_LENGTH * sizeof(unsigned int)));
    checkCudaErrors(cuMemAlloc(&dev_source_data_array,  * ARRAY_LENGTH * sizeof(unsigned char)));
    checkCudaErrors(cuMemAlloc(&dev_result_array, STATIS_size * sizeof(unsigned int)));
    checkCudaErrors(cuMemAlloc(&dev_clock_counter_array,  * ARRAY_LENGTH * sizeof(unsigned int)));

    for (int i = ; i < ARRAY_LENGTH; i++)
    {
        source_data_array[i] = (unsigned char)(rand() % );//STATIS_size);
    }
    memset(result_array, , STATIS_size * sizeof(unsigned int));

    // 傳入待統計的數組
    checkCudaErrors(cuMemcpyHtoD(dev_source_data_array, source_data_array, ARRAY_LENGTH * sizeof(unsigned char)));
    checkCudaErrors(cuMemcpyHtoD(dev_result_array, result_array, STATIS_size * sizeof(unsigned int)));

    // CPU 計數
    QueryPerformanceCounter(&large_interger);
    CPU_c3_alloc = large_interger.QuadPart;
#pragma endregion


#pragma region(執行核函數，并在核函數完成時檢查錯誤報告)
    // 定義網格的大小（block_rect）、塊的大小（thread_rect）
    dim3 block_rect(max(, ARRAY_LENGTH / max_thread_per_block), , );
    dim3 thread_rect(min(max_thread_per_block, ARRAY_LENGTH), , );

    cout << "block_rect :\t" << block_rect.x << "\t" << block_rect.y << "\t" << block_rect.z << "\t" << endl;
    cout << "thread_rect :\t" << thread_rect.x << "\t" << thread_rect.y << "\t" << thread_rect.z << "\t" << endl;

    void *arr[] = { (void *)&dev_thread_index_array, (void*)&dev_source_data_array, \
        (void *)&dev_result_array, (void *)&dev_clock_counter_array };

    // 啟動核函數
    checkCudaErrors(cuLaunchKernel(kernel_addr, \
        block_rect.x, block_rect.y, block_rect.z, \
        thread_rect.x, thread_rect.y, thread_rect.z, \
        , , \
        &arr[], ));

    // 同步，檢查核函數計算過程
    checkCudaErrors(cuCtxSynchronize());


    // CPU 計數
    QueryPerformanceCounter(&large_interger);
    CPU_c4_caculate = large_interger.QuadPart;
#pragma endregion


#pragma region(傳出資料，輸出結果)
    checkCudaErrors(cuMemcpyDtoH(thread_index_array, dev_thread_index_array, ARRAY_LENGTH * sizeof(int)));
    checkCudaErrors(cuMemcpyDtoH(result_array, dev_result_array, STATIS_size * sizeof(unsigned int)));
    checkCudaErrors(cuMemcpyDtoH(clock_counter_array, dev_clock_counter_array,  * ARRAY_LENGTH * sizeof(unsigned int)));


    // 輸出結果
    /*
    */
    float temp_float;
    printf("\r\n\tthread index\tresult\t\tclock counter\tcost time(ms)\r\n");
    for (int i = ; i < ARRAY_LENGTH; i++)
    {
        printf("\t%u\t", thread_index_array[i]);
        printf("\t%u\t\t", source_data_array[i]);
        printf("%u\t", \
        ((clock_counter_array[i + ARRAY_LENGTH] - clock_counter_array[i])));

        if (clock_counter_array[i + ARRAY_LENGTH] > clock_counter_array[i])
        {
            temp_float = ((float)((clock_counter_array[i + ARRAY_LENGTH] - clock_counter_array[i])))\
            / (float)(GPU_clock_rate / );
        }
        else
        {
            temp_float = ((float)((MAXUINT32 + clock_counter_array[i + ARRAY_LENGTH] - clock_counter_array[i])))\
            / (float)(GPU_clock_rate / );
        }

        cost_time_array[i] = temp_float;
        printf("\t%f\r\n", cost_time_array[i]);
    }
    cout << endl;
    for (int i = ; i < ; i++)
    {
        printf("%d\t%u\r\n", i, result_array[i]);
    }


    // CPU 計數
    QueryPerformanceCounter(&large_interger);
    CPU_c5_copyout = large_interger.QuadPart;
#pragma endregion


#pragma region(釋放記憶體、顯存)
    // 釋放顯存
    checkCudaErrors(cuMemFree(dev_thread_index_array));
    checkCudaErrors(cuMemFree(dev_source_data_array));
    checkCudaErrors(cuMemFree(dev_result_array));
    checkCudaErrors(cuMemFree(dev_clock_counter_array));
    // 釋放記憶體
    free(thread_index_array);
    free(source_data_array);
    free(result_array);
    free(clock_counter_array);
    free(cost_time_array);

    // CPU 計數
    QueryPerformanceCounter(&large_interger);
    CPU_c6_end = large_interger.QuadPart;
#pragma endregion


#pragma region(CPU結束計時，輸出結果)

    printf("編譯.cu時長：\t%.4lf 毫秒\r\n", (double)(CPU_c2_compile - CPU_c1_start) / CPU_counter_frequency * );
    printf("配置設定顯存時長：\t%.4lf 毫秒\r\n", (double)(CPU_c3_alloc - CPU_c2_compile) / CPU_counter_frequency * );
    printf("GPU計算時長：\t%.4lf 毫秒\r\n", (double)(CPU_c4_caculate - CPU_c3_alloc) / CPU_counter_frequency * );
    printf("傳出資料時長：\t%.4lf 毫秒\r\n", (double)(CPU_c5_copyout - CPU_c4_caculate) / CPU_counter_frequency * );
    printf("釋放記憶體時長：\t%.4lf 毫秒\r\n", (double)(CPU_c6_end - CPU_c5_copyout) / CPU_counter_frequency * );
    printf("CPU計時總時長：\t%.4lf 毫秒\r\n\r\n", (double)(CPU_c6_end - CPU_c1_start) / CPU_counter_frequency * );
#pragma endregion



    return ;
}



// 檢查顯示卡硬體屬性
void check_Cuda_information(int main_argc, char ** main_argv)
{
    // 裝置ID
    int devID;
    // 裝置屬性
    cudaDeviceProp deviceProps;

    //
    cout << "argc = " << main_argc << endl;
    for (int i = ; i < main_argc; i++)
    {
        printf("argv[%d] = %s\r\n", i, main_argv[i]);
    }
    cout << endl;

    // 擷取裝置ID
    devID = findCudaDevice(main_argc, (const char **)main_argv);

    // 擷取GPU資訊
    checkCudaErrors((CUresult)cudaGetDeviceProperties(&deviceProps, devID));
    cout << "devID = " << devID << endl;
    // 顯示卡名稱
    cout << "CUDA device is \t\t\t" << deviceProps.name << endl;
    // 每個 線程塊（Block）中的最大線程數
    cout << "CUDA max Thread per Block is \t" << deviceProps.maxThreadsPerBlock << endl;
    max_thread_per_block = deviceProps.maxThreadsPerBlock;
    // 每個 多處理器組（MultiProcessor）中的最大線程數
    cout << "CUDA max Thread per SM is \t" << deviceProps.maxThreadsPerMultiProcessor << endl;
    // GPU 中 SM 的數量
    cout << "CUDA SM counter\t\t\t" << deviceProps.multiProcessorCount << endl;
    // 線程束大小
    cout << "CUDA Warp size is \t\t" << deviceProps.warpSize << endl;
    // 每個SM中共享記憶體的大小
    cout << "CUDA shared memorize is \t" << deviceProps.sharedMemPerMultiprocessor << "\tbyte" << endl;
    // 每個Block中寄存器的數量
    cout << "CUDA register per block is \t" << deviceProps.regsPerBlock << "\t" << endl;
    // 每個SM中寄存器的數量
    cout << "CUDA register per SM is \t" << deviceProps.regsPerMultiprocessor << "\t" << endl;
    // GPU時鐘頻率
    //  int    clockRate;       /**< Clock frequency in kilohertz */
    cout << "GPU clock frequency is \t\t" << deviceProps.clockRate << "\tkHz" << endl;
    GPU_clock_rate = deviceProps.clockRate * ;


    cout << endl;
}


void print_PTX(char * PTX_ptr, int size_of_PTX)
{
    // 輸出PTX
    printf("\r\nHere is the PTX code:\r\n\r\n");
    for (int i = ; i < size_of_PTX; i++)
    {
        printf("%c", *(PTX_ptr + i));
    }
    printf("\r\nAbove is the PTX code--------------------------\r\n\r\n");

    return;
}

參考：

1.《CUDA并行程式設計》機械工業出版社

2.NVIDIA官方原子操作說明文檔

3.CUDA Toolkit Documation

4.CUDA Toolkit Documation Atomic篇

5.GM107晶片白皮書 NVIDIA GeForce GTX 750 Ti Whitepaper

CUDA學習筆記(5) 原子操作

繼續閱讀

yolov7 tensorrt模型加速部署【實戰】

C#建立線程三種方式

Android中更新UI的三種方式

linxu下CUDA靜态庫-上

四大元件之——server（服務）及Aidl（接口定義語言）

iOS - 解決This application is modifying the autolayout engine from a background thread

This application is modifying the autolayout engine from a background thread 原因

This application is modifying the autolayout engine from a background thread, which can lead

go中channel和goroutine以及select的配合使用

線程通信和程序通信差別（線程程序差別）

程序與線程及通信總結

《cuda并行程式設計》勘誤（3）

《cuda并行程式設計》勘誤（2）

ThreadLocal 線程局部變量ThreadLocal線程局部變量

ImportError: libcublas.so.9.0: cannot open shared object file: No such file or directory（完美解決）

一種解決思路： ImportError: libcublas.so.10.0: cannot open shared object file: No such file