CUDA学习(十一):原子操作实现向量内积

---

博主CUDA学习系列汇总传送门（持续更新）：编程语言|CUDA入门

---

本文章为 《 GPU编程与优化 大众高性能计算》的读书笔记，例子也都取自书中教程。

---

  向量内积运算中，需要用到的原子操作函数是原子加法函数，GPU的向量内积运算包含的两个累加部分都可以采用原子操作替换，分别是block内thread间归约和block间归约，当然也可以只用一次原子操作。

  原子操作又可以分为全局存储原子操作和共享存储原子操作，正好能对应向量内积中的两个累加部分。

---

一、原子操作的概念

摘自 CUDA atomic原子操作

• CUDA的原子操作可以理解为对一个变量进行“读取-修改-写入”这三个操作的一个最小单位的执行过程，这个执行过程不能够再分解为更小的部分，在它执行过程中，不允许其他并行线程对该变量进行读取和写入的操作。基于这个机制，原子操作实现了对在多个线程间共享的变量的互斥保护，确保任何一次对变量的操作的结果的正确性。
• 原子操作确保了在多个并行线程间共享的内存的读写保护，每次只能有一个线程对该变量进行读写操作，一个线程对该变量操作的时候，其他线程如果也要操作该变量，只能等待前一线程执行完成。原子操作确保了安全，代价是牺牲了性能。
• 原子操作的结果是直接写回存储，返回的值是旧值，所以获取新值需要从新从存储中读取数据。通俗讲，每个线程操作一次共享内存，都需要对这个共享内存进行读-写操作。

原子操作函数种类丰富，具体如下：摘自https://blog.csdn.net/q583956932/article/details/78826987

函数名	作用
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

二、一次原子操作替换两次归约

两次归约在《CUDA学习(十)：向量内积的多种方法实现》中已经提过。

这里直接用原子的加法操作对所有线程的乘积结果进行归约。

• 优点：kernel函数相对比较简单，每个线程求得乘积结果tmp后，直接用原子加法计算即可
• 缺点：原子操作对数据访问是串行的，频繁的原子操作会非常耗时

/* 五、两次归约替换一次原子操作 */

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include <iostream>

const int N = 2048;

const int nBlockNum = 16;//开32个block
const int threadnum = 32;//开32个线程

template <typename T>
__global__ void dot_gpu_4(T *a, T *b, T *c, int n) // c是必须是经过初始化后的，否则结果不对
{
    const int nStep = gridDim.x * blockDim.x; // 跳步的步长，几个block*block内部线程的数量=跳步线程量
    double dTemp = 0.0;
    int nTidIdx = blockIdx.x * blockDim.x + threadIdx.x;//当前线程在全局线程中的索引
    while (nTidIdx < n) // 当前线程需要循环多少次计算结果
    {
        dTemp += a[nTidIdx] * b[nTidIdx];
        nTidIdx += nStep;
    }
    atomicAdd(c, dTemp); // 原子操作
}

int main()
{
    float a[N], b[N];
    float c = 0.0;
    for(int i = 0; i < N; ++i) // 为数组a、b赋值
    {
        a[i] = i * 1.0;
        b[i] = 1.0;
    }
    float *d_a = NULL, *d_b = NULL, *d_c = NULL;
    cudaMalloc(&d_a, N *sizeof(float));
    cudaMemcpy(d_a, a, N *sizeof(float), cudaMemcpyHostToDevice);

    cudaMalloc(&d_b, N *sizeof(float));
    cudaMemcpy(d_b, b, N *sizeof(float), cudaMemcpyHostToDevice);


    cudaMalloc(&d_c, sizeof(float));
    cudaMemset(&d_c, 0, sizeof(float)); //因为原子是累加得到最终结果，所以结果必须赋初值0
    dot_gpu_4<<< nBlockNum, threadnum >>>(d_a, d_b, d_c, N);
  
    cudaMemcpy(&c, d_c, sizeof(float), cudaMemcpyDeviceToHost);
 	cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);
    printf("a * b = %f\n", c);
    printf("hello world\n");
    return 0;
}

三、Block内归约，block间原子操作

• 原子操作对数据访问是串行的，为了减少全局存储操作数量，在上面全局原子操作代码中加入block内低线程归约，最后对每个block的归约结果进行原子加法运算

书上实验表明，block内归约，block间原子操作的kernel函数的性能提升明显。

/* 五、两次归约替换一次原子操作 */

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include <iostream>

const int N = 2048;

const int nBlockNum = 16;//开32个block
const int nThreadnum = 32;//开32个线程

template <typename T>
__global__ void dot_gpu_5(T *a, T *b, T *c, int n) // c是必须是经过初始化后的，否则结果不对
{
    __shared__ T tmp[nThreadnum];
    const int nStep = gridDim.x * blockDim.x; // 跳步的步长，几个block*block内部线程的数量=跳步线程量
    double dTemp = 0.0;
    int nTidIdx = blockIdx.x * blockDim.x + threadIdx.x;//当前线程在全局线程中的索引
    while (nTidIdx < n) // 当前线程需要循环多少次计算结果
    {
        dTemp += a[nTidIdx] * b[nTidIdx];
        nTidIdx += nStep;
    }
    tmp[threadIdx.x] = dTemp;//将每个线程的内积放入到共享存储中
    __syncthreads(); // 同步操作，等待上面执行完成

    int i = nThreadnum / 2; // 低线程归约
    while( i != 0)
    {
        if(threadIdx.x < i)
        {
            tmp[threadIdx.x] += tmp[threadIdx.x + i];
        }
        __syncthreads();
        i /= 2;
    }
    // 注意，每个block只能原子操作一次，所以需要==0或者1或者线程内任意一个数字，否则c = 和*numThread
    if(threadIdx.x == 0) 
    {
        atomicAdd(c, tmp[0]);
    }

}

int main()
{
     float a[N], b[N];
    float c = 0.0;
    for(int i = 0; i < N; ++i) // 为数组a、b赋值
    {
        a[i] = i * 1.0;
        b[i] = 1.0;
    }
    float *d_a = NULL, *d_b = NULL, *d_c = NULL;
    cudaMalloc(&d_a, N *sizeof(float));
    cudaMemcpy(d_a, a, N *sizeof(float), cudaMemcpyHostToDevice);

    cudaMalloc(&d_b, N *sizeof(float));
    cudaMemcpy(d_b, b, N *sizeof(float), cudaMemcpyHostToDevice);


    cudaMalloc(&d_c, sizeof(float));
    cudaMemset(&d_c, 0, sizeof(float));//因为原子是累加得到最终结果，所以结果必须赋初值0
    dot_gpu_5<<< nBlockNum, nThreadnum >>>(d_a, d_b, d_c, N);
  
    cudaMemcpy(&c, d_c, sizeof(float), cudaMemcpyDeviceToHost);
 	cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);
    printf("a * b = %f\n", c);
    printf("hello world\n");
    return 0;
}