并行计算(二)——CUDA

一、简介

CUDA是NVIDIA提供的一种通用的并行计算平台和编程模型，使用CUDA可以像在CPU上一样使用GPU进行编程。

CUDA要介绍的话东西实在太多了，而且GPU的工作原理和CPU尽管是有些相似的，但是实际使用的思路和CPU却可能完全不同，这里也只能简单讲一点。CUDA C编程和普通C语言也没有什么太多的不同，由于CPU和GPU使用的二进制指令不同，因此使用CUDA C编程时需要指定编译器哪些代码编译成CPU的代码，哪些代码编译成GPU的代码。通常我们通过在函数前添加__device__和__global__来告知编译器该函数是在GPU上执行的，而具体的编译则由CUDA的编译器来配合C编译器共同完成。当在函数前添加__device__时则表明该函数在GPU上调用GPU上执行，而__global__则是告知编译器该函数在CPU上调用GPU上执行，另外CUDA还有一个__host__来声明函数是CPU上的函数，但是通常函数声明时只要不在前面加__device__和__global__则默认函数为CPU上的函数。

由于GPU是大规模并发执行的，因此在编写GPU函数时和上一篇的OpenMP很像，我们都必须确定每个线程的索引，并对相应索引的数据进行操作。CUDA的线程分为线程块(block)和线程(thread)两层，一个线程块由多个线程组成，通过线程块和线程编号即可获得需要操作的数据的下标。CUDA本身还提供了多维的线程块和线程索引，具体的用法以及执行模式后面的例子会有体现。

二、示例

1、矢量相加

#include "cuda_runtime.h"
#include <iostream>
#include "device_launch_parameters.h"
#define FL float
#define N 1024
__global__ void add(FL* x,FL* y,FL* z)
{
    unsigned i=blockIdx.x*blockDim.x+threadIdx.x;
    z[i]=x[i]+y[i];
    //blockIdx.x是执行当前代码的线程块编号
    //blockDim.x是线程块内的线程数量
    //threadIdx.x是当前线程号
    //i是通过线程块编号和线程编号获取要计算加和的数据的位置
    //此时只要保证线程块数量与每个线程块中线程总数的乘积是数组大小即可
    //通常也许数组大小并不能被线程块内线程数整除，这时加一个边界检查即可
}
int main()
{
    FL *x,*y,*z,*gpu_x,*gpu_y,*gpu_z;
    x=(FL*)malloc(sizeof(FL)*N);//在主存中申请内存
    y=(FL*)malloc(sizeof(FL)*N);
    z=(FL*)malloc(sizeof(FL)*N);
    x[10]=y[10]=100;
    cudaMalloc(&gpu_x,sizeof(FL)*N);申请显存
    cudaMalloc(&gpu_y,sizeof(FL)*N);
    cudaMalloc(&gpu_z,sizeof(FL)*N);
    cudaMemcpy(gpu_x,x,sizeof(FL)*N,cudaMemcpyHostToDevice);//将数据从主存拷贝至显存
    cudaMemcpy(gpu_y,y,sizeof(FL)*N,cudaMemcpyHostToDevice);

    add<<<8,128>>>(gpu_x,gpu_y,gpu_z);//调用GPU函数，在GPU中计算矢量和，尖括号中是线程块数量和每个线程块中的线程数

    cudaMemcpy(z,gpu_z,sizeof(FL)*N,cudaMemcpyDeviceToHost);//将计算结果从显存拷贝至主存
    std::cout<<z[10]<<'\n';
    free(x);
    free(y);
    free(z);
    cudaFree(gpu_x);
    cudaFree(gpu_y);
    cudaFree(gpu_z);
    return 0;
}

在Linux下使用nvcc进行编译，然后执行即可。这里多说一句，我是用笔记本测试的，笔记本Linux用独显比较麻烦，如果碰巧你和我一样用的是不够国际化的不知名发行版，那可能你会和我一样遇到很多坑爹的情况。

笔记本上Linux带N卡的机器基本上就是三个方案：

第一是不使用独显，纯靠集显输出，这种情况下CUDA程序是不能正常运行的
第二种就是使用大黄蜂，这个东西就是平时用集显，然后你可以通过在运行的命令前加optirun来使用独显执行程序
第三种就是nvidia-primer这个是使用独显计算集显输出(貌似是这样的)这种方案表现上就是所有时候都是用独显在运行程序。

一些知名发行版上nvidia-primer会适配的比较好，如果你恰巧是用的知名发行版，那么恭喜你只需要装好驱动、装好CUDA、装好nvidia-primer就可以了，其他的什么都不用管了。如果你和我一样用的发行版恰好不支持nvidia-primer那么你可能就需要想我一样使用大黄蜂了，不过大黄蜂用着也很简单，你用nvcc编译好二进制文件之后optirun ./a.out就可以了。(如果你的二进制文件名不是a.out就改成对应的名字就好)

CUDA的基本计算就是这样，需要注意的是GPU上的计算数据都来自显存，当你要使用CPU中申请的内存中数据时需要在显存中申请相同大小的内存，并将数据拷贝到显存，同样计算结果也要从显存拷贝到主存才能使用cout、printf等函数输出。另外主存和显存的数据拷贝是一个非常耗时的过程，因此CUDA程序尽量要避免不必要的主存与显存间的数据交互。如果算法本身主存和显存间频繁的数据交互不可避免时，就需要你去评估GPU加速的效果是否能够掩盖显存和主存间数据交互的开销了。

二、矩阵转置

这个测试基本能展示GPU的工作和访存方式了。

首先我们用一维数组来存储矩阵，以行优先的形式进行存储，即同一行的数据是连续的，矩阵一行一行的在一维数组中顺序排列。

#include "cuda_runtime.h"
#include <iostream>
#include "device_launch_parameters.h"
#include <ctime>
using namespace std;
#define FL float
#define N 8192
#define I 4133
#define J 3729
__global__ void kernel_trans(FL* x,FL* y);
void trans(FL* x,FL* y);
int main()
{
    FL *x,*y,*gpu_x,*gpu_y;
    x=(FL*)malloc(sizeof(FL)*N*N);
    y=(FL*)malloc(sizeof(FL)*N*N);
    x[I*N+J]=11323;
    cudaMalloc(&gpu_x,sizeof(FL)*N*N);
    cudaMalloc(&gpu_y,sizeof(FL)*N*N);
    cudaMemcpy(gpu_x,x,sizeof(FL)*N*N,cudaMemcpyHostToDevice);
    cudaMemcpy(gpu_y,y,sizeof(FL)*N*N,cudaMemcpyHostToDevice);

    for(int i=0;i<10;i++) trans(gpu_x,gpu_y);

    cudaMemcpy(x,gpu_x,sizeof(FL)*N*N,cudaMemcpyDeviceToHost);
    cudaMemcpy(y,gpu_y,sizeof(FL)*N*N,cudaMemcpyDeviceToHost);
    cout<<x[I*N+J]<<'\t'<<y[J*N+I]<<'\n';
    free(x);
    free(y);
    cudaFree(gpu_x);
    cudaFree(gpu_y);
    return 0;
}

这是测试函数的主体，根据trans_kernel和trans两个函数的不同可以选择不同的方式进行转置。

1、简单行优先

首先我们测试一下按行读取按列写入的转置

__global__ void kernel_trans(FL* x,FL* y)
{
    unsigned i=blockDim.x*blockIdx.x+threadIdx.x,
             j=blockDim.y*blockIdx.y+threadIdx.y,
             s_pos,d_pos;
    s_pos=i*N+j;
    d_pos=j*N+i;
    y[d_pos]=x[s_pos];
}
void trans(FL* x,FL* y)
{
    kernel_trans<<<dim3(256,256),dim3(32,32)>>>(x,y);
}

这里由于矩阵是8192x8192的，而CUDA中线程块的最大线程数是1024所以我们就按照所有线程块内线程拉满来进行分块，也就是使用256x256的二维线程块，而每个线程块内则是32x32的二维线程。编译完成后在执行命令前nvprof可以记录你的函数执行时间，因为这里测试的是GPU上的性能，因此只需要考察kernel_trans一个函数的用时即可。

使用optirun nvprof ./a.out即可得到kernel_trans的用时，结果如下

Time(%)      Time     Calls       Avg       Min       Max  Name
 49.88%  665.86ms        10  66.586ms  64.950ms  78.212ms  kernel_trans(float*, float*)

我们可以看到这个kernel_trans的平均用时是

66.586ms

2、简单列优先

接下来就是换核函数了，我们试一下列优先读取，按列读按行写

__global__ void kernel_trans(FL* x,FL* y)
{
    unsigned i=blockDim.x*blockIdx.x+threadIdx.x,
             j=blockDim.y*blockIdx.y+threadIdx.y,
             s_pos,d_pos;
    s_pos=i*N+j;
    d_pos=j*N+i;
    y[s_pos]=x[d_pos];
}
void trans(FL* x,FL* y)
{
    kernel_trans<<<dim3(256,256),dim3(32,32)>>>(x,y);
}

这里偷个懒，就直接互换了s_pos和d_pos的位置，其实严谨一点应该是互换两个变量的值才对。同样测一下用时，平均用时是

50.46ms

可以看到列优先读取是快于行优先读取的，这是为什么呢？

这就是前面说的GPU的执行和访存模式的问题了。GPU中有多个SM，每个SM内部又有多个SP更多的细节这里也不多讲了。GPU在执行时会把一个线程块内的每32个线程组织成一个线程束(warp)发到一个SM上执行，而一个SM可以同时并发多个线程块的线程束，这样一来多个SM就可以并发更多的线程块中的线程束，从而实现GPU上的大规模并发。对于一个SM来说尽管在任务调度上它是并发多个SM的，但是其中能进行运算的就是SP，SP的数量限制了物理上一个SM的计算能力，对于一个SM来讲它能够驻留最大的线程束的数量是多于它物理上能够并发的线程束数量的，SM通过交替执行其上驻留的不同线程束来实现并发同时隐藏IO。列优先读取时不同线程束都是读取的不同行的数据，而GPU访存时和CPU相似有一个Cache Line的概念，一次取出多个数据放入缓存中。当多个线程束在SM上交替执行时，第一个线程束在执行时就会缓存多行的数据，当轮换到另一个线程束时其所需要的数据可能已经在Cache中了。而行优先存储时由于其一个线程束中的数据都是连续的，而对于转置这种操作每个线程又只执行一次，因此每个线程束执行时都需要从显存中取去数据，性能就会有所下降。

3、优化线程块内线程数

之前的测试我们是拉满每个线程块中的线程数量，这里我们适当减少每个线程块中的线程数，而适当增加线程块的数量，我们将每个线程块内的线程数改为8x32共256个线程，这时线程块数量就变成了1024x256

__global__ void kernel_trans(FL* x,FL* y)
{
    unsigned i=blockDim.x*blockIdx.x+threadIdx.x,
             j=blockDim.y*blockIdx.y+threadIdx.y,
             s_pos,d_pos;
    s_pos=i*N+j;
    d_pos=j*N+i;
    y[s_pos]=x[d_pos];
}
void trans(FL* x,FL* y)
{
    kernel_trans<<<dim3(1024,256),dim3(8,32)>>>(x,y);
}

用时为

49.64ms

我们看到性能再一次得到一点提升，这次又是为什么呢？这里就纯粹是线程束调度的问题了，之前有讲过SM会并发多个线程块的多个线程束，但是对于SM来讲其上能驻留的线程块数量和线程束数量都是有限制的，我的卡比较老了计算能力只有3.5，一个SM可以驻留64个线程束，和16个线程块，这样算下来当一个线程块包含4个线程束时恰好一个SM可以一次并发16个完整的线程块。而4个线程束一个线程束32个线程总共就是128个线程。由于一级缓存的Cache Line是128bit即32个单精度浮点数，因此我们使用4x32这种二维线程网格来计算，但是实际测试中我发现以4x32这种线程网格计算时速度反而比32x32还要慢，具体原因尚不清楚。我的一个猜测是GPU的显存写入可能是与读取相同采用一次32bit来写入的，因此1列8行的数据转置成一行时就恰好可以一次写入，当我使用8x32这种线程网格时，性能就有了些许提升。但是考虑到这种情况缓存利用率并不算高，因此性能提升并不明显。

4、指令流水优化

再进一步就是流水优化了，GPU也存在和CPU类似的指令流水因此一个线程执行多条指令就可以获得更好的流水，进一步提高性能。这里做一个简单的操作，考虑到一个SM中一级缓存的大小，要充分利用一级缓存，按照8x32的线程来组织线程块时，可以一个线程块执行四个线程块的任务，于是我们就可以使用如下的核函数。

__global__ void kernel_trans(FL* x,FL* y)
{
    unsigned i=4*blockDim.x*blockIdx.x+threadIdx.x,
             j=blockDim.y*blockIdx.y+threadIdx.y,
             in=j*N+i,out=i*N+j;
    y[out]=x[in];
    y[out+blockDim.x]=x[in+blockDim.x*N];
    y[out+2*blockDim.x]=x[in+blockDim.x*2*N];
    y[out+3*blockDim.x]=x[in+blockDim.x*3*N];
}
void trans8(FL* x,FL* y)
{
    kernel_trans<<<dim3(256,256),dim3(8,32)>>>(x,y);
}

这时测得的核函数用时为

24.24ms

可以看到当一个线程中进行多个操作时我们又进一步提升了性能。