并行计算(二)——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加速的效果是否能够掩盖显存和主存间数据交互的开销了。
2、矩阵转置
这个测试基本能展示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.586ms2、简单列优先
接下来就是换核函数了，我们试一下列优先读取，按列读按行写
__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来说尽管在任务调度上它是并发多个线程束的，但是其中能进行运算的就是SP，SP的数量限制了物理上一个SM的计算能力，对于一个SM来讲它能够驻留最大的线程束的数量是多于它物理上能够并发的线程束数量的，SM通过交替执行其上驻留的不同线程束来实现并发同时隐藏IO。列优先读取时不同线程束都是读取的不同行的数据，而GPU访存时和CPU相似有一个Cache Line的概念，一次取出多个数据放入缓存中。当多个线程束在SM上交替执行时，第一个线程束在执行时就会缓存多行的数据，当轮换到另一个线程束时其所需要的数据可能已经在Cache中了。而行优先存储时由于其一个线程束中的数据都是连续的，而对于转置这种操作每个线程又只执行一次，因此每个线程束执行时都需要从显存中取去数据，性能就会有所下降。
3、共享内存的使用
CUDA可以通过在核函数中声明局部变量的语句前加__shared__来在一个线程块内申请共享内存，这部分内存本质上是一块可控的L1 Cache，CUDA将每一个SM上的一级缓存分为共享内存和普通一级缓存两部分，因此使用共享内存和使用一级缓存的性能完全相同。我们就利用共享内存来进一步加速转置。
这里要说一下做法，首先我们申请一个共享内存，然后读取显存中数据并写入到共享内存中，当一个线程块完成了全部共享内存的写入后，再统一把共享内存中的数据写入显存中。这样通过一些巧妙的设计可以几乎达到显存的读写都是顺序的，而转置中不可避免的间断读写都可以在L1 Cache中完成，这样速度就能得到很大的提升。
基本转置思路就是共享内存以转置后的形式划分为二维数组(本例依旧以上面测得最高性能的8x32的子块来转置，那么转置后就是32x8的子块，共享内存就申请为[32][8]这样的二维数组)，读取数据时采用按列读取(根据前面的测试我们可以看到这种读取方式更加有优势)然后将读取的数据转置并存入共享内存相应的位置，最后在完成全部读取之后将共享内存中的数据直接按顺序写入显存。具体代码如下
__global__ void kernel_trans(FL* x,FL* y)
{
unsigned i=blockDim.x*blockIdx.x,
j=blockDim.y*blockIdx.y,
s_pos=i*N+j+threadIdx.y/4*N+8*(threadIdx.y%4)+threadIdx.x,
d_pos=j*N+i+threadIdx.y*N+threadIdx.x;
__shared__ FL t[32][8];
t[(threadIdx.y%4)*8+threadIdx.x][threadIdx.y/4]=x[s_pos];
__syncthreads();
y[d_pos]=t[threadIdx.y][threadIdx.x];
}
void trans(FL* x,FL* y)
{
kernel_trans<<<dim3(1024,256),dim3(8,32)>>>(x,y);
}这时我们测试一下性能，用时为
41.09ms可以看到当使用了共享内存之后性能得到了很大提升。
4、利用L1 Cache代替共享内存
前面已经说了共享内存本质上就是L1 Cache，只不过CUDA通过共享内存为我们提供了一个显式控制L1 Cache的途径，因此在转置这种简单操作中其实通过巧妙的设计，我们完全可以不使用共享内存而达到相同的性能。
首先GPU的显存读写都是以32B进行的，以单精度计算就是8个单精度浮点数。GPU执行时是将每个线程块中的线程以32个为一组组织成线程束(Warp)的形式在SM上执行的，一个SM可以执行多个线程块上的多个线程束。我们可以将GPU中的计算与CPU类比，一个线程束(Warp)类似于一个CPU线程，而一个SM相当于一个CPU核心(带有超线程)。那么类比我之前一篇CPU转置的设计，我们同样可以考虑让一个线程束来完成8x8子块转置，从而达到一个较高的性能。但是由于一个线程束只能包含32个线程，因此我们先试着让它完成8x4的子块转置。这样的话我们可以保证一次转置的显存写入都是完整，至于读取的浪费我们就适当的牺牲一下，从上面的测试来看，这一部分对性能的影响其实不大，因为一个线程块在全部执行完之前是常驻，块内不同线程束的交替并不容易导致缓存更新，所以不同线程束的读取可以认为都是从一级缓存中读取的，多读几次也无妨。具体的设计如下。
__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_trans1<<<dim3(1024,256),dim3(8,32)>>>(x,y);
}我们看到这次使用的线程块依旧是8x32的，这里说一下二维线程索引的线程束划分。对于二维线程来说其本质还是一维线程，CUDA以列优先来组织二维线程索引，即对于[threadIdx.x][threadIdx.y]这一线程其相当于[threadIdx.y*blockDim.x+threadIdx.x]这个一维线程，所以当我们使用8x32这种线程划分时，一个线程束完成的是8x4的子块转置，刚好是我们想要的子块大小。
这时我们测一下速度，用时为
43.32ms可以看到这时速度与前面利用共享内存转置相比略有降低，但是与未使用共享内存相比却性能提升了很大，这也从侧面证明了我们的猜想，就是对于转置这种简单操作，我们完全可以利用执行顺序上的设计来通过一级缓存代替使用共享内存得到更高性能的核函数。
接下来我们再试一下一个线程束内以4x8来转置，这时可以认为是一次子块转置的每一行读取都是完整的，但是一次缓存写入的利用率只有一半。
__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_trans1<<<dim3(12048,256),dim3(4,32)>>>(x,y);
}这时测一下速度，用时为
41.35ms可以看到这次速度就与使用共享内存基本相同了，因此可以推测GPU转置过程中提高显存读取带宽的利用率似乎更加有效。
