cuda c和c++的简单介绍

原文地址，纯翻译

https://developer.nvidia.com/blog/easy-introduction-cuda-c-and-c/

　　这是cuda并行计算平台 c和c++接口系列的第一篇文章。学习前要求熟练掌握c，针对cuda fortran编程的帖子也会同步更新。这两个系列将涵盖cuda平台上并行计算基本概念。从这里开始，除非有特别说明，否则我将使用属于"cuda c"作为cuda c和c++的简写。cuda c本质上是带有一些拓展的c/c++，允许使用多个并行线程在gpu上执行的函数。

CUDA编程模型基础

　　在接触cuda c代码之前，那些刚接触cuda的人最好先了解cuda编程模型基本描述和其中的一些术语。

　　cuda编程模型是一种同时使用cpu和gpu的异构模型。在cuda中，host指cpu及其内存，device指gpu及其内存，host上运行的代码可以管理host和device上的内存、启动kernel(核函数)，这些kernel是device上运行的函数，他们由gpu上的许多线程并发运行。

　　鉴于cuda编程模型的异构性质，cuda c程序的普遍操作顺序是：

　　1.声明和分配host端，device端内存

　　2.初始化host端数据

　　3.将数据从host端传送到device端

　　4.执行一个或多个核函数

　　5.将结果从device端传回host端

　　记住这个操作流程，让我们来看一个cuda c的例子

第一个cuda c程序

　　上一篇文章中，我介绍了六种SAXPY(Scalar Alpha X Plus Y)的方法,其中就包括了cuda c版本，SAXPY表示单精度A*X+Y，对于并行计算来说是一个很好的hello world程序。在这篇文章中，我将展示一个cuda SAXPY的更完整版本，详细说明做了什么以及为什么这样做，完整的SAXPY代码如下：

#include <stdio.h>

__global__ void saxpy(int n,float a,float *x,float *y)
{
    int i=blockIdx.x*blockDim.x+threadIdx.x;
    if(i<n)
    {
        y[i]=a*x[i]+y[i];
    }
}

int main(void)
{
    int N=1<<20;    // 1左移20
    float *x,*y,*d_x,*d_y;
    x=(float*)malloc(N*sizeof(float));
    y=(float*)malloc(N*sizeof(float));

    cudaMalloc(&d_x,N*sizeof(float));
    cudaMalloc(&d_y,N*sizeof(float));

    for(int i=0;i<N;++i)
    {
        x[i]=1.0f;
        y[i]=2.0f;
    }

    cudaMemcpy(d_x,x,N*sizeof(float),cudaMemcpyHostToDevice);
    cudaMemcpy(d_y,y,N*sizeof(float),cudaMemcpyHostToDevice);

    // perform SAXPY in 1M elements
    int threads=256;
    int blocks=(N+threads-1)/threads;   //  上取整
    saxpy<<<blocks,threads>>>(N,2.0f,d_x,d_y);
    
    cudaMemcpy(y,d_y,N*sizeof(float),cudaMemcpyDeviceToHost);

    float maxError=0.0f;
    for(int i=0;i<N;++i)
    {
        maxError=max(maxError,abs(y[i]-4.0f));
    }
    printf("Max error: %f\n",maxError);

    cudaFree(d_x);
    cudaFree(d_y);
    free(x);
    free(y);

    return 0;
}

　　函数saxpy就是在gpu上运行的核函数，main函数是host端代码，让我们先从host端代码开始讨论这个程序

Host端代码

　　主函数声明了两对数组

float *x,*y,*d_x,*d_y;
x=(float*)malloc(N*sizeof(float));
y=(float*)malloc(N*sizeof(float));

cudaMalloc(&d_x,N*sizeof(float));
cudaMalloc(&d_y,N*sizeof(float));

　　指针x和y指向host端数组，使用malloc函数申请内存空间，d_x,d_y数组指向device端数组，使用cuda runtime api中的cudaMalloc函数申请空间。cuda编程中host端和device端具有独立的内存空间，两者都可以通过host端代码进行管理（cuda c内核还可以在支持他的设备上分配设备内存）。

　　host端代码随后初始化host端数组。在此，我们将x设为元素全部为1的数组，y设为元素全为2的数组

for(int i=0;i<N;++i)
{
    x[i]=1.0f;
    y[i]=2.0f;
}

　　要初始化device端数组，我们只需使用cudaMemcpy函数将数据从x，y拷贝到device端对应数组d_x,d_y上，这个过程就像c语言中的memcpy函数，唯一的区别就是cudaMemcpy需要地四个参数来指定数据拷贝到方向（host端到device端还是device端到host端），在此例中，我们使用cudaMemcpyHostToDevice表示数据从host拷贝到device端。

cudaMemcpy(d_x,x,N*sizeof(float),cudaMemcpyHostToDevice);
cudaMemcpy(d_y,y,N*sizeof(float),cudaMemcpyHostToDevice);

　　在运行完核函数后，为了将数据d_y从device端拷贝回host端的y，我们使用cudaMemcpy函数，第四个参数指定为cudaMemcpyDeviceToHost

cudaMemcpy(y,d_y,N*sizeof(float),cudaMemcpyDeviceToHost);

启动核函数

　　saxpy核函数由以下代码启动

saxpy<<<blocks,threads>>>(N,2.0f,d_x,d_y);

　　<<<>>>中间的信息是核函数的执行配置，具体指有多少设备线程并行执行内核。在cuda中，软件有一个线程层次结构，它模仿线程处理器在gpu上的分组方式。在cuda编程模型中，我们通常说通过grid来启动线程block。执行配置中的第一个参数指定grid中的block数目，第二个参数指定一个block中的线程数。

　　block和grid可以通过给dim3传值构造成1维，2维或3维度数据结构（由cuda定义的具有x,y,z策划那个圆的简单结构），但是对于这个简单的例子，我们只需要一个维度，所以我们传递一个整数。在这个例子中，我们设置一个block中包含256个线程，使用取整算法来计算得到处理数组N个元素所需要的block数目（(N+256-1)/256）。

　　对于数组中元素的数量不能被线程块大小整除的情况，内核代码必须检查越界内存访问。（我认为就是检测取整多出来的那部分）

释放内存

　　程序结束时，我们需要释放掉所有申请的内存空间。对于用cudaMalloc()申请的device端内存，使用cudaFree()来释放。对于host端申请单内存，使用free()来释放。

cudaFree(d_x);
cudaFree(d_y);
free(x);
free(y);

device端代码

　　现在我们来看核函数代码

__global__ void saxpy(int n,float a,float *x,float *y)
{
    int i=blockIdx.x*blockDim.x+threadIdx.x;
    if(i<n)
    {
        y[i]=a*x[i]+y[i];
    }
}

　　在cuda中，我们使用__global__声明符来声明核函数。device代码中定义的变量不需要指定为device变量，因为他们被指定驻留在device上。（我的理解是核函数里声明的变量，他的生命周期就是该线程上执行的核函数生命周期）。在本例中，n，a和i变量将由每个线程存储在一个寄存器中，并且指针x和y必须是指向device内存地址空间的指针。因为当我们从host端启动核函数时，我们要将d_x,d_y传递给kernel。前两个变量n和a，我们没有明确的代码将它们从host端传到device端，由于函数参数在c/c++中默认按值传递，因此cuda运行时可以自动处理这些值到device的传输。cuda runtime api这个特性使得gpu上启动核函数非常自然和容易，几乎与调用c函数一样。

　　核函数saxpy的代码只有两行，前面我们提到过，核函数是由多线程并发执行的。如果我们希望每个线程处理数组的一个元素，那么我们需要一种区分和识别每个线程的方法。cuda定义了变量blockDim，blockIdx和threadIdx。这些预定义变量的类型为dim3，类似于主机代码中的执行配置参数。预定义变量blockDim包含在核函数启动的第二个执行配置参数中，指定block中的线程个数（维度）。预定义变量threadIdx和blockIdx分别表示block内的线程索引和grid内地block索引。表达方式：

int i=blockIdx.x*blockDim.x+threadIdx.x;

　　生成用于访问数组元素的全局索引。我们在这个例子中没有使用gridDim，gridDim表示核函数启动的第一个执行配置参数中指定的grid尺寸（一个grid里面有多少个block）

　　在使用索引访问数组元素之前，会根据元素个数n检查索引是否合法，以确保没有内存越界访问。如果数组中的元素不能被block中的线程数整除，（由于向上取整）核函数启动的线程数会大于数组大小，此时需要进行检查。核函数的第二行逐个元素执行SAXPY操作，除了边界检查外，他的结果等于host端通过循环实现的SAXPY

if(i<n) y[i]=a*x[i]+y[i];

编译和运行代码

　　cuda c编译器，nvcc是NVIDIA CUDA Toolkit（http://www.nvidia.com/content/cuda/cuda-toolkit.html）的一部分。为了编译我们的SAXPY例子，我们将代码保存为.cu格式，文件名是saxpy.cu。我们可以使用nvcc编译它

nvcc -o saxpy saxpy.cu

　　我们可以运行代码

./saxpy

总结与结论

　　通过使用cuda c对saxpy简单实现，我们现在了解了cuda编程的基础知识。将c代码移植到cuda c上只需要几个c扩展：核函数的__global__声明符；核函数启动时的执行配置（<<<blocks,threads>>>）；用于识别和区分并行执行核函数的gpu线程的内置设备变量blockDim,blockIdx和threadIdx。

　　异构cuda编程模型的一个优点是，可以增量的将现有代码从c移植到cuda c，一次一个内核。

　　在本系列的下一篇文章中，我么将了解一些性能测量和指标。