CUDA编程练习（三）向量内积

Simple_Bright&Man

已于 2024-06-04 20:55:12 修改

阅读量332

点赞数 10

分类专栏： CUDA编程文章标签：算法 python c++

于 2024-05-28 22:28:55 首次发布

本文链接：https://blog.csdn.net/SimpleBright_Man/article/details/139279382

版权

CUDA编程专栏收录该内容

14 篇文章 2 订阅

订阅专栏

1 想法

相乘没什么好说的，取值是平行的，a取i，b也取i

不行，脑子开始内讧了。
有以下规约想法：

元素太多，每个线程在数组上跨grid规模取值求和，一口气就能把元素数量；
线程足够多，直接block内折半向下归约，公共内存里面会出现block数量的元素（或者放在共享内存开辟的第一个值的位置）；
原子加法，或者压到单block内接着归约；
考虑线程束，线程束规约一次性压缩32倍（？？？有用么这，不还是SM一波一波的一次性压缩一半，记得复习）

代码写满了注释，要理解应该是只考验耐心:>

2 代码

#include <iostream>
#include <assert.h>
#include <cstdlib>
#include <cmath>

#define VECTOR_LENGTH 100000000
#define THREAD_PER_BLOCK 256
#define BLOCK_NUM 64
#define MAX_ERR 1e-4
#define CUDA_CALL(x) do { if((x) != cudaSuccess) { printf("Error at %s:%d\n",__FILE__,__LINE__); return EXIT_FAILURE;}} while(0)

void cuda_error_check(const char* prefix, const char * postfix) {
    if (cudaPeekAtLastError()!=cudaSuccess) {
        printf("\n%s%s%s", prefix, cudaGetErrorString(cudaGetLastError()), postfix);
        cudaDeviceReset();
    }
}

__global__ void dummy_vectorDotProduct(float *out, float *a, float *b, int n) 
{
    float sum=0;  // 寄存器
    for(int i = 0; i < n; i++)
    {
        sum = sum +  a[i] * b[i];
    }
    *out = sum;
}

__global__ void sharedAtomic_vectorDotProduct(float *out, float *a, float *b, int n) // 考虑到共享内存，一个block一个结果
{
    int leap= gridDim.x*blockDim.x;
    int bid = blockIdx.x*blockDim.x+threadIdx.x;
    int tid = threadIdx.x;
    
    __shared__ float ssTmp[THREAD_PER_BLOCK];  // 共存
    ssTmp[tid] = 0;
    // 同步修改
    __syncthreads();
    
    for(int i = bid; i < n; i+=leap)
    {
        ssTmp[tid] += a[i] * b[i];  // 全部加到一个grid里面，
        // __syncthreads();  // 猜猜这个东西有多重要？加了就会导致GPU原地爆炸？BUT WHY？？
    }
    __syncthreads();
    
    // 对共享内存进行归约，限制于block内
    for(int stride = blockDim.x/2; stride>0; stride/=2) 
    {
        if (tid<stride) {
            ssTmp[tid] += ssTmp[tid+stride];
        }
        __syncthreads();
    }
    if (tid==0) atomicAdd(out, ssTmp[0]);
    
}

__global__ void sharedVec_vectorDotProduct(float *out_SsVec, float *a, float *b, int n) // 考虑到共享内存，先存到公共内存
{
    int leap=gridDim.x*blockDim.x;
    int bid = blockIdx.x*blockDim.x+threadIdx.x;
    int tid = threadIdx.x;
    
    __shared__ float ssTmp[THREAD_PER_BLOCK];  // 共存
    ssTmp[tid] = 0;
    // 同步修改
    __syncthreads();
    
    for(int i = bid; i < n; i+=leap)
    {
        ssTmp[tid] += a[i] * b[i];  // 全部加到一个grid里面，
        // __syncthreads();  // 猜猜这个东西有多重要？加了就会导致GPU原地爆炸？BUT WHY？？
    }
    __syncthreads();
    
    // 对共享内存进行归约
    for(int stride = blockDim.x/2; stride>0; stride/=2) 
    {
        if (tid<stride) {
            ssTmp[tid] += ssTmp[tid+stride];
        }
        __syncthreads();
    }
    if (tid==0) out_SsVec[blockIdx.x]=ssTmp[0];
}


__global__ void vec_sumReduction(float *out_SsVec, float *out, int n) // 对公共内存的求和捏，目前考虑总数小于共享内存的量的情况，单block
{
    int tid = threadIdx.x;
    // 读取到共存，然后进行压缩
    __shared__ float ssTmp[THREAD_PER_BLOCK]; // 这里的blockDim一定是大于32的2的幂次方
    if (tid<n) {
        ssTmp[tid] = out_SsVec[tid];
    }
    else {
        ssTmp[tid] = 0;
    }
    // 同步修改
    __syncthreads();
    
    for(int stride = blockDim.x/2; stride>0; stride/=2)  // 规约到线程束上了，呀~
    {
        if (tid<stride) {
            ssTmp[tid] += ssTmp[tid+stride];
        }
        __syncthreads();
    }
    if (tid==0) *out = ssTmp[0];
}

double cpuGpu_gap(float *a, float *b, int n) {
    //double tmp=0;
    //for (int i=0; i<n; i++) {
    //    tmp += abs(a[i]-b[i]);
    //}
    
    return *a-*b;
}

int main() {
    // 设置使用什么卡
    cudaSetDevice(3);
    // 声明变量，结果都用double表示，double不了一点，不是，说好支持atomicAdd的double运算呢？？
    float *a, *b;
    float *cpuout, *device_out, *device_out1, *device_out2;
    float *d_a, *d_b;
    float *d_out, *d_middle;
    // 初始化cpu端的空间
    a = new float[VECTOR_LENGTH];
    b = new float[VECTOR_LENGTH];
    cpuout = new float; *cpuout=0;
    device_out = new float; device_out1 = new float; device_out2 = new float;
    
    // 初始化
    for(int i = 0; i < VECTOR_LENGTH; i++)
    {
        a[i] = std::rand()/double(RAND_MAX)/100;
        b[i] = std::rand()/double(RAND_MAX)/100;  // 防止溢出最蠢的办法，不过，也能看看精确度上的差异
    }
    // cpu端计算
    for (int i=0; i<VECTOR_LENGTH; i++) {
        *cpuout += a[i]*b[i];
    }
    
    printf("The result of cpu:\t%f\n", *cpuout);
           
           
    // 创建gpu端空间，复制到gpu端
    cudaMalloc(&d_a, sizeof(float) * VECTOR_LENGTH);
    cudaMalloc(&d_b, sizeof(float) * VECTOR_LENGTH);
    cudaMalloc(&d_middle, sizeof(float) * BLOCK_NUM);
    cudaMalloc(&d_out, sizeof(float));  // 这个倒是不用多大的内存
    cudaMemcpy(d_a, a, sizeof(float) * VECTOR_LENGTH, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, b, sizeof(float) * VECTOR_LENGTH, cudaMemcpyHostToDevice);
    
    //记录时间
    float time_gpu, time_gpu1, time_gpu2;
    cudaEvent_t start, stop;
    
    // -----------------------------------------
    //创建Event
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    //记录当前时间
    cudaEventRecord(start);
    
    dummy_vectorDotProduct<<<1,1>>>(d_out, d_a, d_b, VECTOR_LENGTH);
    CUDA_CALL(cudaGetLastError());
    
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);    //等待事件完成。记录之前的任务
    cudaEventElapsedTime(&time_gpu, start, stop);    //计算时间差
    printf("The time for dummy_vectorAdd:\t%f(ms)\n", time_gpu);
    // 转移回CPU计算差值
    cudaMemcpy(device_out, d_out, sizeof(float), cudaMemcpyDeviceToHost);
    printf("The gap for dummy_vectorAdd:\t%f\n", cpuGpu_gap(cpuout, device_out, VECTOR_LENGTH));
    cudaEventDestroy(start);    //消除Event
    cudaEventDestroy(stop);
    
    //  ------------------------------------------
    *device_out=0;
    cudaMemcpy(d_out,device_out, sizeof(float), cudaMemcpyHostToDevice);
    //创建Event
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    //记录当前时间
    cudaEventRecord(start);
    
    sharedAtomic_vectorDotProduct<<<BLOCK_NUM,THREAD_PER_BLOCK>>>(d_out, d_a, d_b, VECTOR_LENGTH);  // SM 72 线程最大1024，来，一个SM负载6个！
    
    CUDA_CALL(cudaGetLastError());
    cuda_error_check("Error ", " returned from literal startup kernel");
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);    //等待事件完成。记录之前的任务
    cudaEventElapsedTime(&time_gpu1, start, stop);    //计算时间差
    printf("The time for sharedAtomic:\t%f(ms)\n", time_gpu1);
    // 转移回CPU计算差值
    cudaMemcpy(device_out1, d_out, sizeof(float), cudaMemcpyDeviceToHost);
    printf("The gap for sharedAtomic:\t%f\n", cpuGpu_gap(cpuout, device_out1, VECTOR_LENGTH));
    cudaEventDestroy(start);    //消除Event
    cudaEventDestroy(stop);
    
    //  ------------------------------------------
    //创建Event
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    //记录当前时间
    cudaEventRecord(start);
    
    sharedVec_vectorDotProduct<<<BLOCK_NUM,THREAD_PER_BLOCK>>>(d_middle, d_a, d_b, VECTOR_LENGTH);  // SM 72 线程最大1024，来，一个SM负载6个！
    vec_sumReduction<<<1,512>>>(d_middle, d_out, BLOCK_NUM);
    
    CUDA_CALL(cudaGetLastError());
    
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);    //等待事件完成。记录之前的任务
    cudaEventElapsedTime(&time_gpu2, start, stop);    //计算时间差
    printf("The time for sharedReduction:\t%f(ms)\n", time_gpu2);
    // 转移回CPU计算差值
    cudaMemcpy(device_out2, d_out, sizeof(float), cudaMemcpyDeviceToHost);
    printf("The gap for sharedReduction:\t%f\n", cpuGpu_gap(cpuout, device_out2, VECTOR_LENGTH));
    cudaEventDestroy(start);    //消除Event
    cudaEventDestroy(stop);
    
    
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_out);
    cudaFree(d_middle);
    
    std::cout<<"以上就是并行与不并行的速度差距，还有原子操作和归约到底的差距"<<std::endl;
}

3 结果

The result of cpu:      1605.164795
The time for dummy_vectorAdd:   1890.833130(ms)
The gap for dummy_vectorAdd:    -0.000122
The time for sharedAtomic:      0.975808(ms)
The gap for sharedAtomic:       -894.899170
The time for sharedReduction:   0.890496(ms)
The gap for sharedReduction:    -894.899170
以上就是并行与不并行的速度差距，还有原子操作和归约到底的差距

4 总结+尚存的问题

并行很好，没什么好说的，那么该怎么用共享内存？怎么用原子加法？怎么用公共内存（寄存器）？这些是核心问题；
拆成两个核函数，反而可能会造成性能下降（10^6规模和一个e的规模结果是不一样的），因为异步性，原子加法反而没有那么的串行，不能过分追求并行化，还得考虑组织成本；
说好计算能力高的显卡支持double的原子加法呢？？？估计是编译的时候选项不对，待老夫换个IDE（换个选项就可以了，确实）
变量名规范点，尤其是这种涉及不太熟悉的变量的情况；
你还没尝试成功线程束展开（估计也就是if分支没写对。。。）；——2024.6.4可能线程束并不能在现在的显卡上随便展开哦
你还没尝试迭代写法（元素个数大于的情况下，继续迭代，小于一个block的时候，用小规模的kernel算）