CUDA-4

CUDA-4

这次关注的是Shared MemorySynchronization. 这主要是在线程之间有交流的时候会用到。 这次是以内积为例来说明,其两个向量先作对应的×操作,然后再把各自的和加起来。

先看一下普通的c语言如何写。

#include <stdio.h>

#define N 100

int dot(int *a, int *b){
    int temp=0;
    for(int i=0;i<N;i++){
        temp += (a[i]*b[i]);

    }   

    return temp;
}


int main(){
    int a[N],b[N];
    //int c[N];
    for(int i=0;i<N;i++){
        a[i] = i*i;
        b[i] = i*2;
    
    }   

    
    int ret;
    ret = dot(a, b); 
    printf("the result is %d",ret);

    return 0;
}

可以运行,但是感觉写得太难看了。

可以看到,如果向量的长度很大的话,运行时间将线性增长。下面看看用cuda写的,下面是根据书上改的。

第一种方式。


#include <stdio.h>

#define imin(a,b) (a<b?a:b)


const int N = 33*1024;
const int threads_per_block=256;    //就是blockDim.x的值
const int blocks_per_grid = imin(32, (N-1+threads_per_block)/threads_per_block);   //即看看是否需要'tile'操作,或者理解成
要用一块砖还是多块砖


//kernel function


__global__ void dot(float *a, float *b, float *c){

    __shared__ float cache[threads_per_block];

    int tid = threadIdx.x + blockIdx.x*blockDim.x;
    int cacheIndex = threadIdx.x;

    float temp = 0;
    while(tid<N){
        temp += a[tid]*b[tid];
        tid += blockDim.x*gridDim.x;   //如果只用一块砖,就不用这个操作了。
    }

    cache[cacheIndex] = temp;

    __syncthreads();  // 同步

    // 算一个block中的结果
    int M = blockDim.x;
    while(M>=1){
        if(cacheIndex<M && cacheIndex >=1)
            cache[0] += cache[cacheIndex];
        __syncthreads();
        M -= 1; 
    }

    //c来存当前block的结果
    if(cacheIndex == 0){
        c[blockIdx.x] = cache[0];
    }

}


int main(){
    float *a, *b, c, *partial_c;
    float *dev_a, *dev_b, *dev_partial_c;

    // allocate memory on the cpu
    a = (float*)malloc(N*sizeof(float));
    b = (float*)malloc(N*sizeof(float));
    partial_c = (float*)malloc(blocks_per_grid*sizeof(float));

    //allocate the memory on the gpu
    cudaMalloc((void**)&dev_a, N*sizeof(float));
    cudaMalloc((void**)&dev_b, N*sizeof(float));
    cudaMalloc((void**)&dev_partial_c, blocks_per_grid*sizeof(float));

    // fill in the host memory with data
    for(int i=0;i<N;i++){
        a[i] = i;
        b[i] = i*2;
    }
    //copy the arrays to gpu    

    cudaMemcpy(dev_a, a, N*sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(dev_b, b, N*sizeof(float), cudaMemcpyHostToDevice);
    //注意是目的地在前。

    //run on the gpu
    dot<<<blocks_per_grid, threads_per_block>>>(dev_a, dev_b, dev_partial_c);

    // copy the reult to cpu
    cudaMemcpy(partial_c, dev_partial_c, blocks_per_grid*sizeof(float), cudaMemcpyDeviceToHost);

    // cal the final result on cpu
    c = 0;
    for(int i=0;i<blocks_per_grid;i++){
        c+=partial_c[i];

    }

    // compare the results
    #define sum_squares(x) (x*(x+1)*(2*x+1)/6)

    printf("GPU: %.6g, CPU: %.6g\n", c, 2 * sum_squares((float)(N-1)));

    //free memory

    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_partial_c);

    free(a);
    free(b);
    free(partial_c);

    return 0;
                                                                                            
}


因为在同一个block中的thread是都能够看得到cache的,所以可以只在第0个线程那里做处理就可以了。

nclude <stdio.h>

#define imin(a,b) (a<b?a:b)


const int N = 33*1024;
const int threads_per_block=256;    //就是blockDim.x的值
const int blocks_per_grid = imin(32, (N-1+threads_per_block)/threads_per_block);   //即看看是否需要'tile'操作,或者>理解成要用一块砖还是多块砖


//kernel function


__global__ void dot(float *a, float *b, float *c){

    __shared__ float cache[threads_per_block];

    int tid = threadIdx.x + blockIdx.x*blockDim.x;
    int cacheIndex = threadIdx.x;

    float temp = 0;
    while(tid<N){
        temp += a[tid]*b[tid];
        tid += blockDim.x*gridDim.x;   //如果只用一块砖,就不用这个操作了。
    }

    cache[cacheIndex] = temp;

    __syncthreads();  // 同步

   //c来存当前block的结果
    if(cacheIndex == 0){
        for(int i=1;i<blockDim.x;i++){
            cache[0] += cache[i];
        }
        __syncthreads();
        c[blockIdx.x] = cache[0];
    }

}


int main(){
    float *a, *b, c, *partial_c;
    float *dev_a, *dev_b, *dev_partial_c;

    // allocate memory on the cpu
    a = (float*)malloc(N*sizeof(float));
    b = (float*)malloc(N*sizeof(float));
    partial_c = (float*)malloc(blocks_per_grid*sizeof(float));

    //allocate the memory on the gpu
    cudaMalloc((void**)&dev_a, N*sizeof(float));
    cudaMalloc((void**)&dev_b, N*sizeof(float));
    cudaMalloc((void**)&dev_partial_c, blocks_per_grid*sizeof(float));

    // fill in the host memory with data
    for(int i=0;i<N;i++){
        a[i] = i;
        b[i] = i*2;
    }
    //copy the arrays to gpu    

    cudaMemcpy(dev_a, a, N*sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(dev_b, b, N*sizeof(float), cudaMemcpyHostToDevice);
    //注意是目的地在前。

    //run on the gpu
    dot<<<blocks_per_grid, threads_per_block>>>(dev_a, dev_b, dev_partial_c);

    // copy the reult to cpu
    cudaMemcpy(partial_c, dev_partial_c, blocks_per_grid*sizeof(float), cudaMemcpyDeviceToHost);

    // cal the final result on cpu
    c = 0;
    for(int i=0;i<blocks_per_grid;i++){
        c+=partial_c[i];

    }

    // compare the results
    #define sum_squares(x) (x*(x+1)*(2*x+1)/6)

    printf("GPU: %.6g, CPU: %.6g\n", c, 2 * sum_squares((float)(N-1)));

    //free memory

    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_partial_c);

    free(a);
    free(b);
    free(partial_c);

    return 0;
                                             
}


第三种是书上的思想,更好一些,其把每个block中的256个,做了一个处理。加的次数变成log级别的了。

即用下面的代码来处理当前block中的结果。

int i = blockDim.x/2;
    while(i != 0){
        if(cacheIndex<i){
            cache[cacheIndex] += cache[cacheIndex+i];
        }
        __syncthreads();
        i/=2;

    }
    


这个地方还是


 //c来存当前block的结果
    if(cacheIndex == 0){
        c[blockIdx.x] = cache[0];
    }


回看整个过程

相当于把代码复制了32(gridDim.x)份, 每一份是一个block,然后在每一个block中有一个shared 的cache,这个cache在该block为所有的线程所共享。但是注意不同的block,cache是不同的。每一个block中的结果算完之后做一个小的汇总, 存在dev_partial_c中,最后再把每个block中的算的小的汇总再来一次大的汇总。这个大的汇总可以在cpu上来做,因为数量级已经很小了。

打赏,谢谢~~

取消

感谢您的支持,我会继续努力的!

扫码支持
扫码打赏,多谢支持~

打开微信扫一扫,即可进行扫码打赏哦