cuda 编程入门

1. 简介

首先， 大家都知道在做深度学习以及科学计算时利用GPU算法，可以大规模的提高运算效率。我们首先来了解一下CPU与GPU的区别。

上图为CPU 与GPU 架构上的区别，绿色部分可以理解为逻辑计算部分，黄色为控制区域， 而红色为存储区域，最明显的区别就是CPU拥有更多的存储，尤其是缓存，而GPU拥有更多的计算单元，简单来说单个GPU拥有更多的计算资源。而相对于比较简单的大规模运算，单个GPU可以媲美与集群。由于CPU丰富的缓存，CPU比较适合来计算相对比较复杂的计算，而且一般来说CPU的计算精度更加的高。

在大规模科学计算中，利用GPU的优势的确非常大（主要是能耗以及资源利用效率上面），于此同时利用计算过程中需要各个CPU之间的通信是比较耗时的，但是GPU每个计算单元的相互之间的通信是非常困难的。

本文默认读者已经配置好cuda编译环境。

cuda 下载地址：https://developer.nvidia.com/cuda-downloads

安装指南：https://docs.nvidia.com/cuda/

2. Cuda 第一个程序Helloworld

在介绍第一个程序中我们需要了解一下几个概念。

kernal 代表任务的个数，一个block中拥有多个thread，多个block组成一个grid。每个thread都可以进行逻辑计算。目前每个block中可以建立的thread数量是有限的目前最好的GPU为2048个，而一般的最多只有1024个。这里只是简单介绍，如果先了解更多可以参考博文：https://blog.csdn.net/junparadox/article/details/50540602

#include<stdio.h>

//function define for GPU __global__ is the key word for cuda
__global__ void helloFromGpu(void){
    int index;
    index = blockIdx.x *blockDim.x +threadIdx.x;
    /* 
        blockIdx blockDim threadIdx is variation in system
    */ 
    //printf("%d\n", i);
    printf("Hello World From GPU Thread at %d!\n", index);
}

int main(void){
    //hello from GPU
    printf("Hello World from CPU!!\n");
    // Using the function of GPU <<<block num, thread num>>
    helloFromGpu <<<10, 10>>>();
    cudaDeviceSynchronize();
    return 0;
}

然后利用nvcc编译即可， 

nvcc hello.cu -o hello

接下来我们来解释一下每一行代码的含义：

__global__ void helloFromGpu(void)

函数的定义与C语言基本一样，唯一的的区别多了一个关键字，__global__的含义是对于GPU和CPU都可用，与之对于的关键字还有__device__ 和 __host__ 只有设备端可用（GPU）, 和 只有主机端可用（CPU）对于初学者来说只定义__global__就可以了。

index = blockIdx.x *blockDim.x +threadIdx.x;

这边我们只介绍一维的情况，threadIdx.x 为每个block中thread的ID从[0, 1, 2 .......], blockIdx.x 是block的ID， blockDim.x 是每个block的大小，在一维情况下等同于thread的数量。index为GPU的索引。

helloFromGpu <<<10, 10>>>();

函数的调用与C语言的基本类似，唯一多的是<<<block, thread>>>。

cudaDeviceSynchronize();

这一条命令一定不可以省略，其作用为等待GPU完成任务在进行下一步，否则程序在调用完GPU函数后直接进行下一步，这样就不会有GPU的输出。

3. 数组点乘

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>

#define N (1<<24)
#define blocksize 1024
#define blocknumb (N/blocksize)

// This is a function to check the copy data is right
#define checkCudaAPIErrors(F) if((F) != cudaSuccess) \
{printf("Error at line %d in file %s: %s\n", __LINE__, __FILE__, cudaGetErrorString(cudaGetLastError())); exit(-1);}

// The function compute dot mutiply in CPU
void vecDot(double *a, double *b, double *host_sum){
    int i;
    double sum = 0.0;
    for (i = 0; i < N; i++){
        sum += a[i]*b[i];
    }
    *host_sum = sum;
}

// The function using GPU
__global__ void vecDotGPU(double *a, double *b, double *sub_sum){
    // The index of GPU
    int gid = blockDim.x * blockIdx.x + threadIdx.x;
    // create the share memory for the same block different thread
    __shared__ double component[blocksize];
    // store some block mutiply in a array (share memory)
    component[threadIdx.x] = a[gid] * b[gid];
    // A very important function wait all the thread complete the work to make sure get right result
    __syncthreads();
    // for different thread result can add get a sum result
    // This may difficult to understand easily 
    for (int i = (blocksize >>1); i>0; i=(i>>1)){
        if(threadIdx.x < i)
            component[threadIdx.x] += component[threadIdx.x + i];
        __syncthreads();
    }
    // get the result for every block
    if (threadIdx.x == 0){
        sub_sum[blockIdx.x] = component[0];
    }
}
int main(void){
    int i;
    double *host_a, *host_b;
    double *host_sum;
    struct timeval start;
    struct timeval end;
    double elapsedTime;
    
    double *device_a, *device_b;
    double gpu_sum = 0;
    double *host_subSum;
    double *device_subSum;

    host_a = (double *)malloc(sizeof(double) * N);
    host_b = (double *)malloc(sizeof(double) * N);

    host_sum = (double *)malloc(sizeof(double));
    // init  host_a and host_b
    for (i = 0; i < N; i++){
        host_a[i] = (double)rand()/RAND_MAX;
        host_b[i] = (double)rand()/RAND_MAX;
    }

    vecDot(host_a, host_b, host_sum);

    // GPU 
    gettimeofday(&start, NULL);

    host_subSum = (double *)malloc(sizeof(double) * blocknumb);

    // create memory in GPU
    cudaMalloc((void**)&device_a, sizeof(double) *N);
    cudaMalloc((void**)&device_b, sizeof(double) *N);

    cudaMalloc((void**)&device_subSum, sizeof(double) *blocknumb);

    // copy data from CUP to GPU (host to device)
    checkCudaAPIErrors(cudaMemcpy(device_a, host_a, sizeof(double) * N, cudaMemcpyHostToDevice));
    checkCudaAPIErrors(cudaMemcpy(device_b, host_b, sizeof(double) * N, cudaMemcpyHostToDevice));

    // Call the GPU function
    vecDotGPU<<<blocknumb, blocksize>>>(device_a, device_b, device_subSum);
    //copy dat from GPU to CPU (device to host)
    checkCudaAPIErrors(cudaMemcpy(host_subSum, device_subSum, sizeof(double) * blocknumb, cudaMemcpyDeviceToHost));
    gettimeofday(&end, NULL);
    // free the memeory of GPU
    cudaFree(device_a);
    cudaFree(device_b);
    cudaFree(device_subSum);
    // Add all block result get the sum
    for (i =0; i < blocknumb; i++){
        gpu_sum += host_subSum[i];
    }

    elapsedTime  = (end.tv_sec - start.tv_sec) * 1000.0;
    elapsedTime += (end.tv_usec - start.tv_usec) / 1000.0;
    printf("GPU running time is %f ms\n", elapsedTime);
    printf("CPU sum result is %lf\n", *host_sum);
    printf("GPU sum result is %lf\n", gpu_sum);
    // free the memory of CPU
    free(host_a);
    free(host_b);
    free(host_sum);
    free(host_subSum);

}

利用以下命令编译:

nvcc vectorDotMul.cu  -o vectorDotMul

__shared__ double component[blocksize];

首先来介绍以下子函数中的语法: 首先我们要知道gpu中不同block数据通信是非常困难的，而相同block 不同thread之间的通信需要借助与share memory，作为数组点乘，我们需要把两个数组相乘之后把所有乘积相加，我们这里做法是利用share memory 将相同block里面的乘积相加之后，再在CPU上把不同的block结果相加，从而得到最终结果。

__syncthreads();

这个函数的作用是等待同一block上所有thread完成之前的操作，主要是由于在求和过程中，不同thread之间的数据会有相互依赖，如果不加这个函数，会出现有的thread还在计算上一步结果，而有的thread已经在计算下面的结果，会出现计算错误。

接下来介绍一下同一block下每个thread相加算法，由于算法不是特别明确，我们只介绍主要过程，具体详细内容读者可以通过代码自行测试理解。

主要是将得到的乘积数组分成两部分，然后将前半部分 加 后半部分， 然后将索引长度减半，直到索引长度为1。

ps: 我们的threadNum 为2的n次方。

cudaMalloc((void**)&device_a, sizeof(double) *N);

在GPU内存中开辟sizeof(double) * N的memory，并且将这块内存赋给device_a。

cudaMemcpy(device_a, host_a, sizeof(double) * N, cudaMemcpyHostToDevice)

将CPU中host_a 复制到 GPU中 device_a

cudaMemcpy(host_subSum, device_subSum, sizeof(double) * blocknumb, cudaMemcpyDeviceToHost)

 将GPU中的数据copy回CPU。

cudaFree(device_subSum);

Free GPU中的内存空间。