第二章 CUDA编程模型

• 写一个CUDA程序
• 执行一个核函数
• 用网格和线程块组织线程
• GPU性能测试

2.1 CUDA编程模型概述

CUDA编程是一种单指令多线程的变成范式，何为单指令多线程？我的个人理解是某一行代码（或者是写了某个函数），但是会有同时多个“对象”来并行执行这个函数，是并行执行而不是C++多线程里面的并发执行，这个“对象”其实就是CUDA的每一个执行单元（线程束、Core单元）,至于什么是线程束和Core单元，我们后面会讲到，先插个眼，至于C++多线程，我的理解是让CPU来高效的切换上下文执行不同线程体内的方法，属于并发，总之他们最大的区别就是平台不同（GPU vs CPU）。免责声明，我也是刚刚入门的小白，理解可能有误，如有偏差或者错误，也烦请指正。

2.1.1 CUDA编程结构

CUDA是一种C语言扩展的语言，所以整个编程环境包含GPU和CPU，整体步骤分为三个部分：

1. CPU初始化GPU资源，传输数据
2. GPU计算
3. GPU给CPU返回结果

因此我们使用以下进行区分

1. 主机：CPU及其内存
2. 设备：GPU及其内存

在CUDA6.0之前，CPU（主机）和GPU（设备）之间的内存是相互独立的，即一个CPU分配的指针无法访问GPU的内存数据，如果需要将主机的数据拷贝到GPU，需要使用特定的指令进行拷贝，而CUDA6.0以后提出了“统一寻址”的内存模型，将CPU和GPU的内存空间连接起来，可以使用单个指针访问CPU和GPU的内存空间，无序彼此进行手动拷贝，第四章会进行介绍，但是我们还是应该学习一下如何规范地在GPU和CPU之间进行内存分配和数据拷贝。

举个例子：

1、在主机上有一块内存host_a,在GPU上也有一块内存dev_a;
2、你想要吧host_a 的数据拷贝到dev_a上,进行GPU的运算
3、你需要执行 GPU 内存拷贝命令（dev_a，host_a,host_a的长度）
而不是直接另dev_a = host_a,否则会出现系统崩溃。

下面进入正题：
简述一下CUDA编程的流程：
图2-1所示为CUDA程序的执行流程，首先是由CPU初始化GPU资源，然后往GPU传输数据，然后GPU的内核启动，开始执行核函数进行计算（此时程序控制权同时返回给CPU，详见第六章），然后将计算结果返回到CPU。

2.1.2 内存管理

CUDA编程模型假设系统是由一个主机和一个设备组成，各自拥有独立的内存，表2-1列出了标准C函数以及CUDA C函数的内存操作指令

C语言	CUDA C
malloc	cudaMalloc
memcpy	cudaMemcpy
memset	cudaMemset
free	cudaFree

1、CUDA执行内存分配的函数原型为：

cudaError_t cudaMalloc (void** devPtr,size_t size)

该函数负责向设备分配一定字节的线性内存，并以devPtr的形式返回指向所分配内存的指针。cudaMalloc与标准C语言中的malloc函数几乎一样。

2、cudaMemcpy函数负责主机和设备之间的数据传输，其函数原型为：

cudaError_t cudaMemcpy(void* dst,const void* src,size_t count, cudaMemcpyKind kind)

这个函数比较特殊的地方就在于最后一个参数cudaMemcpyKind kind了，它代表数据拷贝的方向：

• cudaMemcpyHostToHost 主机到主机
• cudaMemcpyHostToDevice 主机带设备
• cudaMemcpyDeviceToHost 设备到主机
• cudaMemcpyDeviceToDevice 设备到设备

此外，这个函数时同步的，可以理解，毕竟是拷贝数据，在cudaMemcpy返回以前，主机的程序应该是阻塞的，如果拷贝成功则返回cudaSuccess,失败返回cudaErrorMemoryAllocation。
我们可以用：

char* cudaGetErrorString(cudaError_t error)

这个方法来将错误码转化为可读的错误信息，与strerror类似。

下图所示为简单的CUDA模型：
主机（CPU）进行设备内存分配、数据拷贝到设备，将设备数据拷贝回主机以及设备的内存释放。

代码demo : 使用CPU和GPU进行两个数组的相加：
CPU：sumArraysOnhost.c

c代码我就不手打了，书里P22-23页,也合并到下面GPU程序中了

GPU：sumArraysOnhost.cu

#include <cuda_runtime.h>  // cuda头文件
#include <stdio.h>
#include <time.h>
#include <memory.h>
#include <stdlib.h>
// c语言中两个一维矩阵相加
void sumArrays(float* a, float* b, float* res, const int size)
{
	for (int i = 0; i < size; i += 4)
	{
		res[i] = a[i] + b[i];
		res[i + 1] = a[i + 1] + b[i + 1];
		res[i + 2] = a[i + 2] + b[i + 2];
		res[i + 3] = a[i + 3] + b[i + 3];
	}
}

// CUDA 核函数，和c语言代码的最大区别是：不需要for循环了！！！
// 这就是CUDA编程模型是单指令多线程的魅力吧
__global__ void sumArraysGPU(float* a, float* b, float* res)
{
	// 每个线程在执行时，都能获取到它们自身的blockIdx、blocKDim
	// threadIdx等内置变量，得到一个唯一的索引，
	//然后去读取它对应的矩阵值 
	int i = blockIdx.x * blockDim.x + threadIdx.x;
	res[i] = a[i] + b[i];
}

// 检查cpu和cuda执行的加法结果是否一致
void checkResult(float* hostRef, float* gpuRef, const int N)
{
	double epsilon = 1.0E-8;
	for (int i = 0; i < N; i++)
	{
		if (abs(hostRef[i] - gpuRef[i]) > epsilon)
		{
			printf("Results don\'t match!\n");
			printf("%f(hostRef[%d] )!= %f(gpuRef[%d])\n", hostRef[i], i, gpuRef[i], i);
			return;
		}
	}
	printf("Check result success!\n");
}
// 矩阵的随机初始化
void initialData(float* ip, int size)
{
	time_t t;
	srand((unsigned)time(&t));
	for (int i = 0; i < size; i++)
	{
		ip[i] = (float)(rand() & 0xffff) / 1000.0f;
	}
}
// 矩阵的随机初始化
void initialData_int(int* ip, int size)
{
	time_t t;
	srand((unsigned)time(&t));
	for (int i = 0; i < size; i++)
	{
		ip[i] = int(rand() & 0xff);
	}
}

int main() {
	// main函数就是主机代码
	// 首先初始化GPU设备
	int dev = 0;
	cudaDeviceProp deviceProp;
	cudaGetDeviceProperties(&deviceProp, dev);
	printf("Using device %d: %s\n", dev, deviceProp.name);

	// 分配cpu和GPU内存
	int nElem = 1 << 14;
	printf("Vector size:%d\n", nElem);
	int nByte = sizeof(float) * nElem;
	float* a_h = (float*)malloc(nByte);
	float* b_h = (float*)malloc(nByte);
	float* res_h = (float*)malloc(nByte);
	float* res_from_gpu_h = (float*)malloc(nByte);
	memset(res_h, 0, nByte);
	memset(res_from_gpu_h, 0, nByte);

	float* a_d, * b_d, * res_d;
	cudaMalloc((float**)&a_d, nByte);
	cudaMalloc((float**)&b_d, nByte);
	cudaMalloc((float**)&res_d, nByte);
	
	//初始化矩阵
	initialData(a_h, nElem);
	initialData(b_h, nElem);
	
	// CPU 数据拷贝到GPU，cudaMemcpyHostToDevice
	cudaMemcpy(a_d, a_h, nByte, cudaMemcpyHostToDevice);
	cudaMemcpy(b_d, b_h, nByte, cudaMemcpyHostToDevice);
	
	// 配置CUDA，每个线程块包含1024个线程
	dim3 block(1024);
	
	// 一共有 nElem / 1024个网格块，
	// 所以总共是nElem个线程，和数组大小一致
	// 因此第i个线程就去读取数组1和数组2的第i个值，然后相加
	// 把结果存储到结果的第i个值中
	// 这样复杂度就是O(1)了，不过好像书里没有提及GPU的复杂度？？？ 
	dim3 grid(nElem / block.x);
	
	// 调用核函数
	sumArraysGPU << <grid, block >> > (a_d, b_d, res_d);
	printf("Execution configuration<<<%d,%d>>>\n", grid.x, block.x);
	
	// 把结果拷贝回主机内存,这里有一个隐式同步
	cudaMemcpy(res_from_gpu_h, res_d, nByte, cudaMemcpyDeviceToHost);
	sumArrays(a_h, b_h, res_h, nElem);
	
	// 比较cpu和gpu的计算结果
	checkResult(res_h, res_from_gpu_h, nElem);
	
	// 释放内存
	cudaFree(a_d);
	cudaFree(b_d);
	cudaFree(res_d);

	free(a_h);
	free(b_h);
	free(res_h);
	free(res_from_gpu_h);
	return 0;
}

我这里的环境是w10,vs2019，显卡是RTX3050Ti,cuda10,和书里的linux环境有点区别，如何配置环境可以百度一下
cuda的代码是以.cu为后缀，然后用nvcc编译器进行编译
我的编译和运行结果如下：

2.1.3 线程管理

我们在主机调用核函数后，它的执行会移动到设备上，也就是我们最开始提到的CUDA的每一个执行单元,有一个内核启动的所有线程统称为一个网格（grid）,一个网格中有许多的线程块（block）,每个线程块组织了多个线程用来执行核函数。同一个线程块内的线程可以通过：同步以及共享内存进行协作。
每一个线程都要一个独立的ID（索引），这个ID是根据两个坐标进行计算得到：blockIdx（线程块在网格内的索引）和threadIdx（线程在线程块内的索引）,这些索引是核函数的内置变量，也就是说在启动核函数时，每个线程都能知道自己的blockIdx和threadIdx。
在CUDA中，网格和线程块最多是可以通过三维来组织，，我们也可通过blickDim 和 gridDim来查询网格以及线程块的维度，但是常见的是一维和二维。具体看一下这个例子

#include <cuda_runtime.h>
#include <stdio.h>
int main(int argc,char ** argv)
{
  int nElem=1024;
  dim3 block(1024);
  dim3 grid((nElem-1)/block.x+1);
  printf("grid.x %d block.x %d\n",grid.x,block.x);

  block.x=512;
  grid.x=(nElem-1)/block.x+1;
  printf("grid.x %d block.x %d\n",grid.x,block.x);

  block.x=256;
  grid.x=(nElem-1)/block.x+1;
  printf("grid.x %d block.x %d\n",grid.x,block.x);

  block.x=128;
  grid.x=(nElem-1)/block.x+1;
  printf("grid.x %d block.x %d\n",grid.x,block.x);

  cudaDeviceReset();
  return 0;
}

2.1.4 启动一个CUDA核函数

kernel_name <<< grid,block >>> (argument list);

这行代码的kernel_name是一个核函数<<< >>>里面的参数是核函数的执行配置，第一个值grid是网格的维度，第二个值block的线程块维度，例如你要配置网格为一维4个网格块，每一个网格块包含8个线程，那么你就可以按照如下填写核函数的配置：

dim3 grid(4)  // dim3是一个uint3的数据结构，可以创建三维的数据，也就是x,y,z 
dim3 block(8)  // eg dim grid(3,4,2) 这就是创建了一个三维的数据，x维度为3，y维度为4，z维度为2
kernel_name <<< grid,block >>> (argument list);

一个网格包含四个红色框（线程块），一个线程块包含8个线程，它们都是按照一维组织的。

还是这张图，他是一个二维的网格，网格的行数为2，列数为3；网格中的线程块，行数为3，列数为5，因此创建的规则如下：

dim3 grid(3,2);  
dim3 block(5,3);
kernel_name <<< grid,block >>> (argument list);

要注意一点，GPU网格、线程块的排列规则和线性代数的矩阵规则有点不同，矩阵的第一维行数，第二维代表列数，例如A[4][3] 是代表一个四行三类的矩阵，但是block(4,3)代表的是一个三行四列的线程块，感觉有点像和矩阵反着来。
然后例如图中的一个线程Thread(4,2),代表这线程在这一个线程块中的行id threadIdx.x = 4, threadIdx.y = 2;总之第一维代表x方向,第二维位代表y方向，如果有第三维就代表z方向。网格和线程块的表示规则相同。

异步行为：CUDA不同于C语言，所有的核函数都是异步的!!!
在主机传输完数据、调用核函数后，内核调度相应的线程块进行异步的计算，这个时候控制权又返回到主机这边，相当于GPU自己去干活了，CPU自己也要干活（主机程序继续往下走了），但是如果主机走完程序那就要退出了（return 0）,而GPU的任务还没干完，这个时候就需要让主机进行等待，用到下面这个API来进行阻塞同步：

cudaError_t cudaDeviceSynchronized(void);

调用这个函数也称为显式同步，而像cudaMemcpy这些函数，本身是有同步机制在的，成为隐式同步。

2.1.5 编写核函数

什么是核函数：
只能由设备（CUDA）调用的函数成为核函数，其特点是必须用 __global__关键字声明，而且必须返回 void
核函数的特点：

1. 仅能范围设备的内存
2. 返回void
3. 不支持可变参数类型
4. 不支持静态变量
5. 异步

下表总结了CUDA C程序的函数类型限定符

限定符	执行	调用
global	设备端	可以从主机端调用，也可以从计算能力为3的设备调用
device	设备端	仅从设备端调用，或者核函数（global）中调用
host	主机端	仅从主机端调用，可以省略这个关键字

此外__device__ 和 __host__关键字可以一起使用，声明一个可以在主机、设备端同时调用的函数。

2.1.6 验证核函数

2.1.7 处理错误

本小节是使用了一个宏定义来进行错误码判断，要注意宏定义的写法：只能写在一行，因此每一行结尾都要带一个换行符 \ 表示是同一行

#define CHECK(call)\
{\
  const cudaError_t error=call;\
  if(error!=cudaSuccess)\
  {\
      printf("ERROR: %s:%d,",__FILE__,__LINE__);\
      printf("code:%d,reason:%s\n",error,cudaGetErrorString(error));\
      exit(1);\
  }\
}

2.2 给核函数计时

2.2.1 cpu给核函数计时

在c语言中，我们常用gettimeofday来记录耗时

#include <sys/time.h>  // linux
$include <time.h>     // windows

struct timeval start,end;
gettimeofday(&start,NULL);
dosomething();
gettimeofday(&end,NULL);

int elaps_time = (end.tv_sec-start.tv_sec)*100000+\
(end.tv_usec-start.tv_usec); // 单位是微秒

2.2.2 nvprof给核函数计时

由于笔者电脑的GPU和CUDA版本较高，已经无法使用nvprof这个工具了，所以相关的内容只能从书上了解，自己没有实践。

现在取而代之的好像叫ncu和nsyc,后面再进行补充，插个眼。

2.3 组织并行线程

GPU的优势就在于它可以组织多个线程进行真正并行的计算，所以如果选择合适的网格和线程块大小，就能够对内核性能产生很大的影响，本节通过一个矩阵加法的例子来进一步说明这一点。对于矩阵的运算，传统的想法是使用一个包含二维网格与二维线程块的布局来组织线程，但是这种方法并非能获取到最佳性能，因此本节将采用一下布局进行比较：

1. 二维网格+二维线程块
2. 一维网格+一维线程块
3. 二维网格+一维线程块

2.3.1 使用块和线程建立矩阵索引

对于一个矩阵加法，往往是通过给一个线程分配指定的数据元素进行计算，因此首先需要确定的是某个线程的索引，以及它要访问的矩阵（存在内存中）的地址。

对于具体的某个线程，它自己知道的是它的threadIdx信息，blockIdx信息和blockDim信息，所以该线程在当前线程块中的索引为：
block_idx = threadIdx.y*blockDim.x + threadIdx

而如果它想要知道自己在当前CUDA全局下的索引，他需要先计算：
ix = threadIdx.x + blockIdx.xblockDim.x
iy = threadIdx.y + blockdix.yblockDim.y

ifx = iynx + ix 其中nx为矩阵的列数（blockDim.x * grid.x）
我的理解是，加入要计算矩阵A[m][n]+B[m][n],那么会创建一个总线程数量为mn的CUDA设备，分别让每个线程去计算对应所以的矩阵值。

2.4 设备管理

英伟达提供了几个API用于查询和管理GPU设备：

2.4.1 使用cudaDeviceProp结构体保存GPU信息

cudaErrot_t cudaGetDeviceProperties(cudaDeviceProp* prop, int device);

其中cudaDeviceProp结构体返回了GPU设备的属性.

2.4.2 确定最优的GPU

通过cudaGetDeviceProperties方法找到多个GPU中处理器最大的那个，代码就不放了

2.4.3 nivdia-smi

具体的查询信息就看书吧

总结

本章是CUDA的基础，介绍了CUDA的编程模型，线程的组织方式以及内存管理和信息查询的API，这个是最基础的内容了。后面三四五六章是CUDA关于内存、缓冲的介绍，难度直线提升，加油！各位看官如果觉得有用，麻烦动动发财小手点个赞再走。

最后感谢一下谭升大佬的博客: link，受益匪浅，代码也是参考大佬的github,如有兴趣可以去star一下他！