CUDA reduce 并行规约求和

规约求和原理如下, 把半的数据加到前一半上, 经过N次迭代, 最后归为一个数

那么为什么可以用CUDA 加速呢?

1. 一个block中的线程可以同步处理, 比如一个block中的32个数:
   第一次迭代:
   data[0] += data[16]
   data[1] += data[17]
   …
   data[15] += data31]
   可以同步进行;
   第二次迭代:
   data[0] += data[8]
   data[1] += data[9]
   …
   data[7] += data[15];
   可以同步进行;
   第N次迭代:
   data[0] += data[1]

2. 将数据分到不同的block中, 先对每一个block单独并行处理, 最后再处理每一个block 的结果, 例如, 1024 * 1024 大小的数据, 每一个block处理1024个数据, 一共1024 个block, 处理完后剩下1024个数据, 接着在处理, 就有加速的效果:

每一个block 块作为 global 核函数的处理单元, 在一个block中能保证一个block 中的数据同步.

代码如下:

__global__ void global_reduce_kernel(float *d_in, float *d_out, int maxsize){
    int my_id = threadIdx.x + blockDim.x * blockIdx.x;
    int tid = threadIdx.x;
    if(my_id > maxsize)return;
    for( int s = blockDim.x / 2; s > 0; s>>=1){//迭代N 次
        if(tid < s){ // 为什么用tid 约束, 因为能保持block同步.
            if((my_id + s) > maxsize){
                continue;
            }else{
                d_in[my_id] += d_in[my_id + s]; // 将后半段数据加到前半段
            }
        }
        __syncthreads();  // 保证block中数据同步
    }
    if(tid == 0){ // 处理每一个block 块
        d_out[blockIdx.x] = d_in[my_id];
    }
}

tid : 每一个block 中的线程id
my_id : 所有 block 中的线程id
程序执行时会执行 blocks * threads 次 global 核函数, 区分开的就是通过
threadIdx.x 和 blockIdx.x

2. 使用 shared 共享内存

将每一个数据拷贝到对应的block 中的 shared 共享内存中. 这样减少全局内存的访问,提升速度.

__global__ void share_reduce_kernel(float *d_in, float *d_out, int maxsize){
    extern __shared__ float sdata[]; // 每一个block中共享一个sdata 
    int tid = threadIdx.x;           // 每一个block中线程id
    int mid = threadIdx.x + blockDim.x * blockIdx.x; //所有block中的访问id
    sdata[tid] = d_in[mid]; //将所有block中的数据复制到每一个block中的共享内存sdata中, 这是共享内存精华
    __syncthreads(); // 等待block中所有数据拷贝完
    for(int i = blockDim.x / 2; i > 0; i>>=1){
        if(tid < i){
            if(mid + i > maxsize){
                continue;
            }else{
                sdata[tid] += sdata[tid + i];  // 已经将每一个block的数据拷贝到sdata, 因此只需要在block中操作.
            }
        }
    }
    __syncthreads(); // 等待block中所有数据处理完

    if(tid == 0){
        d_out[blockIdx.x] = sdata[0]; // 每一个block 中的共享内存中的求和结果赋值给d_out
    }
}

在一个核函数中声明 共享内存 extern shared float sdata[], 其实就是在一个block 块中分配一块共享内存, 如果有N 个block块, 就相当分配了N 个共享内存块, 每一个 block 中共享同一个共享内存.

sdata[tid] = d_in[mid];
将不同block 中的数据放置到不同block中的共享内存中.可同步进行,几乎不耗费时间,但要同步.
往后的处理只在block中进行,减少全局内存访问,即: sdata[tid] += sdata[tid + i];

注意:
规约求和时, 每个block的线程数量为 2 的 N 次方, 否则可能导致对不齐没有归到一个数上结果出错.

附上全代码:

#include<iostream>
using namespace std;

__global__ void share_reduce_kernel(float *d_in, float *d_out, int maxsize){
    extern __shared__ float sdata[]; // 每一个block中共享一个sdata 
    int tid = threadIdx.x;           // 每一个block中线程id
    int mid = threadIdx.x + blockDim.x * blockIdx.x; //所有block中的访问id
    sdata[tid] = d_in[mid]; //将所有block中的数据复制到每一个block中的共享内存sdata中, 这是共享内存精华
    __syncthreads(); // 等待block中所有数据拷贝完
    for(int i = blockDim.x / 2; i > 0; i>>=1){
        if(tid < i){
            if(mid + i > maxsize){
                continue;
            }else{
                sdata[tid] += sdata[tid + i];  // 已经将每一个block的数据拷贝到sdata, 因此只需要在block中操作.
            }
        }
    }
    __syncthreads(); // 等待block中所有数据处理完

    if(tid == 0){
        d_out[blockIdx.x] = sdata[0]; // 每一个block 中的共享内存中的求和结果赋值给d_out
    }
}

__global__ void global_reduce_kernel(float *d_in, float *d_out, int maxsize){
    int my_id = threadIdx.x + blockDim.x * blockIdx.x;
    int tid = threadIdx.x;
    if(my_id > maxsize)return;
    printf("maxsize = %d \n", maxsize);
    //printf("blockDim.x = %d \n", blockDim.x);
    for( int s = blockDim.x / 2; s > 0; s>>=1){
        printf(" s = %d \n ", s);
        if(tid < s){
            if((my_id + s) > maxsize){
                //d_in[my_id] = d_in[my_id];
                continue;
            }else{
                d_in[my_id] += d_in[my_id + s];
            }

            printf("tid = %d , my_id = %d , my_id + s = %d , d_in[my_id] = %d \n", tid, my_id, my_id+s, d_in[my_id]);
            //d_in[my_id] += d_in[my_id + s];
            
        }
        __syncthreads();
    }
    if(tid == 0){
        printf(" blockIdx.x = %d , my_id = %d \n ", blockIdx.x, my_id);
        d_out[blockIdx.x] = d_in[my_id];
    }
}


void test_cuda_device(){
    int deviceCount;
    cudaGetDeviceCount(&deviceCount);
    if (deviceCount == 0) {
        fprintf(stderr, "error: no devices supporting CUDA.\n");
        exit(EXIT_FAILURE);
    }
    int dev = 0;
    cudaSetDevice(dev);

    cudaDeviceProp devProps;
    if (cudaGetDeviceProperties(&devProps, dev) == 0)
    {
        printf("Using device %d:\n", dev);
        printf("%s; global mem: %dB; compute v%d.%d; clock: %d kHz\n",
               devProps.name, (int)devProps.totalGlobalMem, 
               (int)devProps.major, (int)devProps.minor, 
               (int)devProps.clockRate);
    }
}

void reduce(float *d_out, float* d_intermediate, float *d_in, int size){
    int max_threadPerBlock = 4;
    int threads = max_threadPerBlock;
    int blocks = (size + max_threadPerBlock - 1) / max_threadPerBlock;
    //global_reduce_kernel<<<blocks, threads>>>(d_in, d_intermediate, size);
    share_reduce_kernel<<<blocks, threads>>>(d_in, d_intermediate, size);

    int remain_size = blocks;
    threads = (blocks + 1) / 2 * 2;
    blocks = 1;

    cudaDeviceSynchronize();

    float temp[20];
    cudaMemcpy(temp, d_intermediate, 10 * sizeof(float), cudaMemcpyDeviceToHost);
    for(int k = 0; k < 20; k++){
         std::cout << " k = " << k << "   " <<temp[k]<< std::endl;
    }

    std::cout << " blocks = " << blocks << "  threads = " <<threads<< std::endl;
    //global_reduce_kernel<<<blocks, 8>>>(d_intermediate, d_out, remain_size);
    share_reduce_kernel<<<blocks, 8>>>(d_intermediate, d_out, remain_size);
}
int main(){

    
    test_cuda_device();
    int d_size = 18;
    float *hdata = new float[d_size];
    float sum = 0;
    for(int i = 0; i < d_size; i++){
        hdata[i] = i;
        sum += hdata[i];
    }
    std::cout << " sum = " << sum << std::endl;

    float *ddata_in, *ddata_intermediate, *ddata_out;
    cudaMalloc((void**)&ddata_in, sizeof(float) * d_size);
    cudaMalloc((void**)&ddata_intermediate, sizeof(float) * d_size);
    cudaMalloc((void**)&ddata_out, sizeof(float) * 1);

    cudaMemcpy(ddata_in, hdata, d_size * sizeof(float), cudaMemcpyHostToDevice);

    reduce(ddata_out, ddata_intermediate, ddata_in, d_size);

    float  hdata_out;
    cudaMemcpy(&hdata_out, ddata_out, 1 * sizeof(float), cudaMemcpyDeviceToHost);
    std::cout << " hdata_out = " << hdata_out << std::endl;
    return 0;
}