一篇文章理解CUDA架构、编程与进阶使用

本文章详细介绍了CUDA的架构和基础编程方法，并对他的进阶优化方法进行了简单介绍，以便大家对CUDA编程有一个整体的认知。

---

一、CUDA架构

下图为GPU硬件模型：
一块GPU包括3级：GPU、多核处理器、线程处理器

• 一个GPU包含多个多核处理器（SM，图中的Mlultiprocessor），GPU的内存是全局内存global memory（可以被所有线程访问）
• 一个多核处理器包含多个线程处理器，多核处理器的内存是共享内存shared memory（编程时划分好block后，一个block内的所有线程可以访问共享内存）
• 线程处理器，最基本的计算单元，有自己的局部内存和寄存器，只能自己访问

在CUDA编程时，我们经常用到thread, block, grid，其中thread对应硬件上的线程处理器，grid对应一块GPU。而block可以由我们自定义维度，对应到硬件上，其实是由一个多核处理器中的多个线程处理器组合而成，可以将一个多核处理器划分为多个block。

线程束(warp)是最基本的执行单元，一个warp包含32个基本的计算单元-线程thread，也就是说比如我发一个指令，那么线程束中的32个thread将会并行执行该指令。（所以在我们划分blocksize的时候，一般都会设置成32的倍数）

二、CUDA编程基础

CUDA编程并行计算整体流程

1. 在GPU上分配显存，将CPU上的数据拷贝到显存上
2. 利用核函数完成GPU显存中数据的计算
3. 将显存中的计算结果拷贝回CPU内存中

从矩阵加法和矩阵乘法来学习CUDA编程的基本框架（耐心看完下面代码，基本都有注解，可以模仿范式编写自己的代码）

1.矩阵加法

计算矩阵加法：C = A + B，设A B为一维矩阵，长度为n

//核函数（即在GPU中执行的函数/用__global__申明）
__global__
void vecAddKernel(float* A_d, float* B_d, float* C_d, int n)
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;  //计算线程ID
    if (i < n) C_d[i] = A_d[i] + B_d[i];  //筛选ID小于n的线程，即例如线程1计算C_d[1] = A_d[1] + B_d[1]
}

//主函数
int main(int argc, char *argv[]) {

    int n = 10000;

    size_t size = n * sizeof(float);

    // 在CPU上分配内存
    float *a = (float *)malloc(size);
    float *b = (float *)malloc(size);
    float *c = (float *)malloc(size);

	//初始化a b的值（将需要计算的向量放到分配好的内存中）
    for (int i = 0; i < n; i++) {
        float af = rand() / double(RAND_MAX);
        float bf = rand() / double(RAND_MAX);
        a[i] = af;
        b[i] = bf;
    }

	//在GPU上分配显存（格式按照 参考下面代码，size为需要分配的显存大小）
    float *da = NULL;
    float *db = NULL;
    float *dc = NULL;
    cudaMalloc((void **)&da, size);
    cudaMalloc((void **)&db, size);
    cudaMalloc((void **)&dc, size);

	//将CPU上初始化的a b值拷贝到GPU上
    cudaMemcpy(da,a,size,cudaMemcpyHostToDevice);
    cudaMemcpy(db,b,size,cudaMemcpyHostToDevice);
	
	//划分GPU的block和Grid
    int threadPerBlock = 256;  //一个warp大小为32，一般设置为32的倍数
    int blockPerGrid = (n + threadPerBlock - 1)/threadPerBlock;  //根据划分的blocksize计算gridsize
	
	//调用核函数
    vecAddKernel <<< blockPerGrid, threadPerBlock >>> (da, db, dc, n);
   	//将GPU上的计算结果拷贝回CPU
    cudaMemcpy(c,dc,size,cudaMemcpyDeviceToHost);
	
	//释放GPU显存资源
    cudaFree(da);
    cudaFree(db);
    cudaFree(dc);
    
	//释放CPU内存资源
    free(a);
    free(b);
    free(c);
    return 0;
}

2.矩阵乘法

计算矩阵乘法：C = A * B，矩阵A的维度为M*K，矩阵B的维度为K*N

#define M 512
#define K 512
#define N 512

void initial(float *array, int size)
{
	for (int i = 0; i < size; i++)
	{
		array[i] = (float)(rand() % 10 + 1);
	}
}

//核函数（传入显存ABC以及维度信息MNK）
__global__ void multiplicateMatrix(float *array_A, float *array_B, float *array_C, int M_p, int K_p, int N_p)
{
	//这里我们划分的lblock和grid是二维的，分别计算线程的二维索引（x方向和y方向的索引）
	int ix = threadIdx.x + blockDim.x*blockIdx.x;//row number，
	int iy = threadIdx.y + blockDim.y*blockIdx.y;//col number

	if (ix < N_p && iy < M_p)  //筛选线程，每个线程计算C中的一个元素，线程的xy索引与C的元素位置索引对应
	{
		float sum = 0;
		for (int k = 0; k < K_p; k++) //C中的某个元素为A中对应行和B中对应列向量的乘积和。
		{
			sum += array_A[iy*K_p + k] * array_B[k*N_p + ix];
		}
		array_C[iy*N_p + ix] = sum;
	}
}

//主函数
int main(int argc, char **argv)
{
	int Axy = M * K;
	int Bxy = K * N;
	int Cxy = M * N;

	float *h_A, *h_B, *hostRef, *deviceRef;
	
	//在CPU上分配内存
	h_A = (float*)malloc(Axy * sizeof(float)); 
	h_B = (float*)malloc(Bxy * sizeof(float));
	h_C = (float*)malloc(Cxy * sizeof(float));

	initial(h_A, Axy);
	initial(h_B, Bxy);
	
	//在GPU上分配显存
	float *d_A, *d_B, *d_C;
	cudaMalloc((void**)&d_A, Axy * sizeof(float));
	cudaMalloc((void**)&d_B, Bxy * sizeof(float));
	cudaMalloc((void**)&d_C, Cxy * sizeof(float));
	
	//将CPU上初始化的a b值拷贝到GPU上
	cudaMemcpy(d_A, h_A, Axy * sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(d_B, h_B, Bxy * sizeof(float), cudaMemcpyHostToDevice);

	//划分GPU的block和Grid
    int dimx = 2;
    int dimy = 2;
    dim3 block(dimx, dimy);
    dim3 grid((M + block.x - 1) / block.x, (N + block.y - 1) / block.y);
    
    //调用核函数
    multiplicateMatrix<<<grid,block>>> (d_A, d_B, d_C, M, K, N);
    
	//将GPU上计算结果拷贝回CPU
    cudaMemcpy(h_C, d_C, Cxy * sizeof(float), cudaMemcpyDeviceToHost);

	//释放GPU显存资源
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    
	//释放CPU内存资源
    free(h_A);
    free(h_B);
    free(h_C);

return (0);
}

三、CUDA进阶 I——利用共享内存加速访存

一般我们将数据发送到GPU后，默认保存到全局内存，而全局内存的读写速度特别慢，这个时候我们将数据从全局内存放到线程块的共享内存中，计算过程中，读取访问速度更快的共享内存，将会大大减少数据访问耗时，提高程序速度。

1.CUDA内存读写速度比较

下列几种内存的架构参见下图：

• 线程寄存器（~1周期）
• Block共享内存（~5周期）
• Grid全局内存（~500周期）
• Grid常量内存（~5周期）
  

2.申请共享内存

前面对比了共享内存和全局内存的访问速度，为了进一步提高访存速度，可以把全局内存一部分数据拷贝到共享内存中（由于共享内存的大小有限，大概只有几十K，所以只能分多次拷贝数据）

申请共享内存的方式分为静态申请和动态申请
申请共享内存关键字：__ shared __
块内共享内存同步：__syncthreads()函数（块内不同线程之间同步）

• 静态申请

__global__ void staticFun(int* d, int n)
{
	__shared__ int s[64];  //静态申请，需要指定申请内存的大小
	int t = treadIdx.x；
	s[t] = d[t];  //将全局内存数据拷贝到申请的共享内存中，之后利用共享内存中的数据参与运算将会比调
	//用全局内存数据参与运算快(由于共享内存有限，不能全部拷贝到共享内存，这其中就涉及到分批拷贝问题了）
	__syncthreads();//需要等所有线程块都拷贝完成后再进行计算
}

staticFun<<1,n>>(d, n);

• 动态申请

__global__ void dynamicFun(int *d, int n)
{
	extern __shared__ int s[]; //动态申请，不需要指定大小，需要加上extern关键字
	int t = threadIdx.x;
	s[t] = d[t];
	__syncthreads();
}

dynamicFun<<1, n, n*sizeof(int)>>(d, n); //动态申请需要在外部指定共享内存大小

上面内容只是让大家对共享内存如何加速运算有一个初步的认识，详细使用方法可以参考我的另外一篇文章：CUDA加速计算矩阵乘法&进阶玩法(共享内存)

后面的内容有待补充。。。有用的话记得点赞搜藏o

四、CUDA进阶 II——利用stream加速大批量文件IO读写耗时

1. 认识CUDA stream

CUDA的stream流，类似我们经常使用CPU时开多线程。

• 当我们使用GPU进行计算时，如果我们没有主动开启stream流，GPU会自动创建默认流来执行核函数，默认流和CPU端的计算是同步的。（也即在CPU执行任务过程中，必须等GPU执行完核函数后，才能继续往下执行）
• 当我们使用GPU进行计算时，我们可以主动开启多个stream流，类似CPU开启多线程。我们可以将大批量文件读写分给多个流去执行，或者用不同的流分别计算不同的核函数。开启的多个流之间是异步的，流与CPU端的计算也是异步的。所以我们需要注意加上同步操作。
  值得注意的是，受PCIe总线带宽的限制，当一个流在进行读写操作时，另外一个流不能同时进行读写操作，但是其他流可以进行数值计算任务。这个有点类似与CPU中的流水线机制。
  

2. CUDA stream API介绍

• 创建一个stream
cudaStream_t stream;
cudaStreamCreate(&stream);
• 将host数据拷贝到device
cudaMemcpyAsync(dst, src, size, type, stream)
• kernel在流中执行
kernel_name<<<grid, block, stream>>>(praments);
• 同步和查询
cudaError_t cudaStreamSynchronize(cudaStream_t stream)
cudaError_t cudaStreamQuery(cudaStream_t stream);
• 销毁流
cudaError_t cudaStreamDestroy(cudaStream_t stream)

上面仅对CUDA stream有一个简单的介绍和认知，CUDA stream使用示例在这篇文章进行了介绍：CUDA优化方案—stream的使用

官方参考文档：https://developer.nvidia.com/blog/gpu-pro-tip-cuda-7-streams-simplify-concurrency/

五、CUDA进阶 III——调用cuBLAS库API进行矩阵计算

cuBLAS是一个BLAS的实现，允许用户使用NVIDIA的GPU的计算资源。使用cuBLAS 的时候，应用程序应该分配矩阵或向量所需的GPU内存空间，并加载数据，调用所需的cuBLAS函数，然后从GPU的内存空间上传计算结果至主机，cuBLAS API也提供一些帮助函数来写或者读取数据从GPU中。
• 列优先的数组，索引以1为基准
• 头文件 include "cublas_v2.h“
• 三类函数（向量标量、向量矩阵、矩阵矩阵）

cuBlas使用范例

int main(int argc, char **argv)
{
		......
        cublasStatus_t status;
        cublasHandle_t handle;
        cublasCreate(&handle);

        float a = 1, b = 0;
        cublasSgemm(
          handle,
          CUBLAS_OP_T,   //矩阵A的属性参数，转置，按行优先
          CUBLAS_OP_T,   //矩阵B的属性参数，转置，按行优先
          M,          //矩阵A、C的行数
          N,          //矩阵B、C的列数
          K,          //A的列数，B的行数，此处也可为B_ROW,一样的
          &a,             //alpha的值
          d_A,            //左矩阵，为A
          K,          //A的leading dimension，此时选择转置，按行优先，则leading dimension为A的列数
          d_B,            //右矩阵，为B
          N,          //B的leading dimension，此时选择转置，按行优先，则leading dimension为B的列数
          &b,             //beta的值
          d_C,            //结果矩阵C
          M           //C的leading dimension，C矩阵一定按列优先，则leading dimension为C的行数
        );
        cudaMemcpy(deviceRef, d_C, Cxy * sizeof(float), cudaMemcpyDeviceToHost);
        cudaDeviceSynchronize();
        ......
}

详细的各种API介绍，推荐参考官方文档：https://docs.nvidia.com/cuda/cublas/index.html