第三章 简单CUDA程序的基本框架

        第二章中介绍了Hello World程序CUDA的线程组织，但printf函数只有在调试的时候才用到，本章将通过数组相加的计算详解CUDA程序的基本框架。

3.1 数组相加

        考虑一个同样长度的一维数组的对应元素之和，其C++程序代码如下所示:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

const double epsilon=1e-15;
const double a=1.23;
const double b=2.34;
const double c=3.57;

void add(double *a,double *b,double *c,const int number)
{
    for(int i=0;i<number;i++)
        {
            c[i]=a[i]+b[i];
        }
}

void check(const double *z,const int Number)
    {
        for(int i=0;i<Number;i++)
            {
                if( abs(z[i]-c)>epsilon)
                    {
                        throw i;
                    }
            }
    }
int main()
{
    const int Number=1e10;
    const int ByteNumber=Number*sizeof(int);
    double *x=(double*)malloc(ByteNumber);
    double *y=(double*)malloc(ByteNumber);
    double *z=(double*)malloc(bytenumber);
    for(int i=0;i<Number;i++)
        {
            x[i]=a;
            y[i]=b;
        }
    add(x,y,z,Number);
    try
    {
        check(z,Number);
    }
    catch(int error)
    {
        printf("数据在%d位置处发生了错误",error);
        return 0;
    }
    printf("no error");
    free(x);
    free(y);
    free(z);
    return 0;

    
}

使用g++编译运行该程序，将在屏幕上输出“no error”,表示计算结果正确，现在对该程序进行如下解释:
        (1)在主函数中32-34行首先定义了3个双精度浮点数类型的指针变量，然后将其指向由函数malloc()(<stdlib.h>进行了定义)分配的内存，从而得到3个长度为1e8的一维数组，将需要2.4GB主机内存，后面的CUDA也会用到这么大的设备内存，若主机内存不够多，可以适当调整数组的大小。
        (2)在35-39行将x，y数组每个元素分别初始化为1.23和2.34.
        (3)在第40行调用自定义函数add()计算数组x和数组y的和，并将结果保存在z中。
        (4)在第43行用一个自定义函数check()检查数组z中的每个元素是否都是正确值3.57.且判断两个浮点数特别是双精度浮点数不能使用‘=’，而要将两者的差的绝对值与一个很小的数进行判断。
        (5)51-53用以释放内存。

3.2 CUDA程序的基本框架

        一个典型的CUDA程序框架见以下代码：
 

头文件包含
常量定义(或宏定义)
C++自定义函数和CUDA核函数的声明(原型)
int main()
{
    分配主机和设备内存
    初始化主机中的数据
    将某些数据从主机复制到设备
    调用核函数在设备中进行计算
    将某些数据从设备复制到主机
    释放主机和设备内存
}
C++自定义函数和CUDA核函数的定义（实现）

在上述CUDA程序的基本框架中，有很多内容还没有介绍。但是，先把CUDA求数组和的全部源代码列出来，之后在逐步讲解：
 

#include<stdio.h>
#include<math.h>
const double epsilon=1.0e-15;
const double a=1.23;
const double b=2.34;
const double c=3.57;

__global__ void add(const double *x,const double *y,const double *z)
{
    int tid=blockDim.x*blockIdx.x+threadIdx.x;
    z[tid]=x[tid]+y[tid];
}
void check(double *z,const int n)
{
    bool haserror=false;
    for(int i=0;i<n;i++)
        {
            if(abs(z[i]-c)>epsilon)
                haserror=true;
        }
printf("%s\n",haserror?"ERROR":"NO ERROR");
}
int main()
{
    const int n=1e8;
    const int m=sizeof(double)*n;
    double *h_x=(double*)malloc(m);
    double *h_y=(double*)malloc(m);
    double *h_z=(double*)malloc(m);
    for(int i=0;i<n.i++)
        {
            h_x[i]=a;
            h_y[i]=b;
        }
    double *d_x,*d_y,*d_z;
    cudaMalloc((void**)d_x,m);
    cudaMalloc((void**)d_y,m);
    cudaMalloc((void**)d_z,m);
    cudaMemcpy(d_x,h_x,m,cudaMemcpyHostToDevice);
    cudaMemcpy(d_y,h_y,m,cudaMemcpyHostToDevice);
    
    const int block_size=128;
    const int grid_size=n/128;
    add<<<grid_size,block_size>>>(d_x,d_y,d_z);
    cudaMemcpy(h_z,d_z,m,cudaMemcpyDeviceToHost);
    check(h_z,n);
    free(h_x);
    free(h_y);
    free(h_z);
    cudaFree(d_x);
    cudaFree(d_y);
    cudaFree(d_z);
    return 0;
    
}

用nvcc编译该程序时，并指定与GeForce MX450对应的计算能力(可选取自己所用的GPU的计算能力)：
        nvcc -arch=sm_75 add1.cu
将得到一个可执行文件a.exe.运行该程序得到的输出与之前的C++程序得到的结果一样，说明达到了预期的效果。
        值得注意的是，当使用较大的数据量时，网格大小相对很大。本例中1e8/128=781250.若用CUDA8.0，会使用默认的2.0的计算能力，且网格大小gridDim.x的限值为65535，小于本例中所使用的值，将会导致程序无法正确的运行。

3.2.1 隐形的设备初始化

        在CUDA的运行时API中，没有明显地初始化设备(GPU)的函数。在第一次调用一个和设备管理及版本查询功能无关的运行时API函数时，设备将自动初始化。

3.2.2 设备内存的分配和释放

        在上述程序中，首先在主机中定义了三个数组并进行了初始化。与C++版本一致。然后在35-38行定义了三个数组指针并分配了内存(显存)。若不看后面不知道三个指针指向哪个内存，只有用cudaMalloc()函数才能确定他们指向设备中哪个内存，而不是主机中的内存。该函数涉及一个CUDA运行时API函数。所有CUDA运行时API函数都以cuda开头，本文涉及较少的内容，完整API运行时函数参照：https://docs.nvidia.com/cuda/cuda-runtime-api。
        C++中分配内存可由malloc()(<stdlib.h>中)进行主机内存的分配，同样在CUDA中可由cudaMalloc()进行分配。该函数的原型如下:
        cudaError_t cudaMalloc(void ** address,size_t size);
其中:
        (1)第一个参数address是待分配设备内存的指针。由于内存(地址)本身就是一个指针，所以待分配设备内存的指针就是指针的指针，即双重指针。
        (2)第二个参数就是待分配内存的字节数。
        (3)返回值是一个错误代号。如果调用成功，则返回cudaSuccess,否则返回一个代表错误的代号。
        该函数为某个变量分配size个字节的线性内存(linear memory)。在36-38行忽略了cudaMalloc的返回值，且采用sizeof(double),在一般情况下该值为8。
        调用函数cudaMalloc()时传入的第一个参数(void**)&d_x。可以知道d_x是一个double类型的指针，那么他的地址&d_x就是double类型的双重指针，通过void**进行强转，也可以不写强转从而可以调用如下:cudaMalloc(&d_x,m);
        之所以使用双重指针，是因为函数的功能是改变指针d_x本身的内容，而不是改变d_x所指向内存缓冲区中的变量值，从而必须将d_x的地址&d_x传给函数cudaMalloc()才能达到这种效果。且通过h_前缀标识主机变量，d_标识设备变量。
        在C++中分配和释放内存分别采用malloc()和free(),从而在cuda函数中使用cudaMalloc()和cudaFree()进行内存的分配和释放，该函数的原型是:
        cudaError_t cudaFree(void *address);
这里的参数address就是待释放的设备内存变量(不是双重指针)，返回值是一个错误变量，如果调用成功，返回cudaSuccess。
        主机内存也可以用new 和delete运算符进行分配和释放。在分配和释放内存的时候必须进行匹配，也就是说cudaMalloc()和cudaFree()必须成对存在，若将cudaFree()换成free()则也会输出NO ERROR，但是程序退出之前，会出现所谓的段错误(segmentation fault)。
        从计算能力2.0开始，CUDA还允许在核函数内部使用malloc()和free()动态的分配一定数量的全局内存，但是程序性较差。

3.2.3 主机与设备之间的数据传递

        在分配了设备内存之后，就可以将某些数据从主机传递到设备中去了，第39-40行将主机的h_x和h_y中的数据复制到设备中的相应变量d_x和d_y所指的缓冲区中，这里用到了CUDA的运行时API函数cudaMemcpy().其原型如下:
cudaError_t cudaMemcpy(void *dst,const void *src,size_t count,enum cudaMemcpykind kind);
其中:
        (1)第一个参数dst表示目标地址
        (2)第二个参数是src的源地址
        (3)第三个参数count是复制数据的字节数
        (4)第四个参数kind是一个枚举类型的变量，标志数据传递方向。他只能取以下几个值:
                1)cudaMemcpyHostToHost,表示从主机复制到主机。
                2)cudaMemcpyHostToDevice,表示从主机复制到设备。
                3)cudaMemcpyDeviceToHost,表示从设备复制到主机。
                4)cudaMemcpyDeviceToDevice,表示从设备复制到设备。
                5)cudaMemcpyDefault,表示根据指针dst和src所指向地址自动判断数据的传输方向，这要求系统具有统一虚拟寻址的功能(unified virtual addressing)的功能(要求64位的主机)。CUDA正逐步放弃对32位主机的支持，所以一般情况下数据传输方向是没有问题的。
        (5)返回值是一个错误代号。若调用成功则返回cudaSuccess.
        (6)该函数的作用是将一定字节数的数据从源地址所指向的缓冲区中复制到目标地址所指的缓冲区。
        通过观察第39行，就是将h_x指向的主机内存中的m个字节的数据复制到d_x所指的设备内存中，因为这里的源地址是主机中的内存，目标地址是设备中的地址，多以第四个参数必须是cudaMemcpyDeviceToHost或cudaMemcpyDefault.

3.2.4 核函数中数据与线程的对应

        将数据从主机传递到设备之后，就可以调用核函数在设备中进行计算，42-44行确定了核函数的执行配置:使具有128个线程的一维线程块，一共有1e8/128个线程块。仔细对比第一章中的add.cpp和本次的add1.cu的设备端函数，可以看出，将主机中的函数改为设备中的核函数非常简单:基本上是去掉一层循环。主机函数中需要通过循环的方式对每一个元素进行操作。在设备的核函数中，用“单指令-多线程”的方式编写代码，顾可以去掉循环，只需要将数组元素和线程指标一一对应即可。
        例如，使用了该语句：
        const int tid=blockDim.x*blockIdx.x+threadIdx.x;
来确定对应方式。赋值右边出现内建变量，左边的n是需要用的数组元素指标。第0个线程块的blockDim.x个线程对应于0-blockDim.x-1的编号，第一个线程块的blockDim.x个线程对应blockDim.x-2*blockDim.x-1的编号。线程块指标与数据一一对应之后，就可以对数组元素进行操作了z[n]=x[n]+y[n].
        

3.2.5 核函数的要求

        核函数是CUDA编程中最重要的部分，现在列出编写核函数时要注意的几点。
        (1)返回必然是void，可以使用return 但是不能返回任何值。
        (2)必须有限定符__global__,也可以加上static等限定符，次序可以使随意的。
        (3)函数名无特殊要求，不支持可变数量的参数列表，参数列表必须确定。
        (4)除非采用统一内存机制，不然传递给核函数的数组指针必须指向设备内存。
        (5)核函数不能成为一个类的成员，通常的一个做法是用一个包装函数调用核函数，而将包装函数定义为类的成员。
        (6)在计算能力3.5之前，核函数之间不能相互调用，但在3.5之后引入了动态并行机制，在核函数中可以调用其他核函数甚至是自己。
        (7)核函数的调用都是在设备中执行，且调用的时候必须指定之星配置。

3.2.6 核函数中if语句的必要性

        如果N是blockDim.x的整数倍(即block_size),不会引发任何问题，因为核函数中的线程数目刚好等于数组元素的个数，当N不是blockDim.x的整数倍的时候，会引发错误。
        若将N改为1e8+1,依然取block_size=128.此时N/128之后，商为781250，余数是1。可以将grid_size定义为781251，使得定义的线程数为1e8+128。可以通过条件语句规避不需要的线程操作，从而设定grid_size=:
        (N-1)/block_size+1;或者(N+block_size-1)/block_size;
以上两个语句等效于：
        int grid_size=(N%block_size==0)?(N/block_size):((N-1+block_size)/block_size);
因为此时线程数(1e8+128)大于总的线程数(1e8+1),如取消if将会导致非法的设备内存操作。可以在核函数中使用return终止语句的进行。也可以将上述代码写成以下形式:
        const int n=blockIdx.x*blockDim.x+threadIdx.x;
        if(n>N) return;
        z[n]=x[n]+y[n];

void __global__ add(const double *a,const double *b,const double*c,const int N)
{
    const int tid=blockDim.x*blockIdx.x+threadIdx.x;
    if(tid>N) return;
    else
        {
    c[tid]=a[tid]+b[tid];
        }
}

3.3 自定义设备函数

        核函数可以调用不带执行配置的自定义函数，这种自定义函数称为设备函数(device function)
.其是在设备中被调用的。与之相比核函数是在设备中执行，但在主机端被调用。

3.3.1 函数执行空间标识符

        在CUDA程序中，以下标识符确定一个函数在哪里被调用，以及哪里被执行:
        (1)__global__修饰函数称为核函数，一般由主机调用，在设备中执行，若采用动态并行，也可以在核函数中调用自己或其他核函数。
        (2)用__device__修饰的函数称为设备函数，只能被核函数或其他设备函数调用，在设备中执行。
        (3)用__host__修饰的函数是主机端普通C++函数，可忽略，有时可用__host__和__device__标识一个函数既是一个C++中的普通函数，也是一个设备函数。可以在主机和核函数中进行调用。
        (4)不能用__global__和__device__同时修饰一个函数。
        (5)不能使用__host__和__global__进行锈蚀。__global__不能与任意host或device修饰词连用。
        (6)编译器可以将__device__修饰的函数定义为内联或非内联元素，可以使用__noinline__和__forceinline__进行改变。
 

3.3.2 为数组相加核函数添加一个设备函数

        

//版本一有返回值
double __device__ add1_device(const double x,const double y)
{
    return x+y;
}

void __global__ add1(const double *x,const double *y,const double *z,const int n)
{
    int tid=threadIdx.x+blockDim.x*blockIdx.x;
    if(tid<n)
    {
        z[tid]=add1_device(x[tid],y[tid]);
    }
}

//版本二通过指针进行
double __device__ add2_device(const double x,const double y,double *z)
{
    *z =x+y;
}

void __global__ add1(const double *x,const double *y,const double *z,const int n)
{
    int tid=threadIdx.x+blockDim.x*blockIdx.x;
    if(tid<n)
    {
        add2_device(x[tid],y[tid],&z[tid]);
    }
}
//版本三引用传值
double __device__ add1_device(const double x,const double y，double & z)
{
    z= x+y;
}

void __global__ add1(const double *x,const double *y,const double *z,const int n)
{
    int tid=threadIdx.x+blockDim.x*blockIdx.x;
    if(tid<n)
    {
        add1_device(x[tid],y[tid],z[tid]);
    }
}