CUDA多线程编程

一、头文件

1、cuda_runtime.h

CUDA 目前有两种不同的 API：Runtime API 和 Driver API，两种 API 各有其适用的范围。由于 runtime API 较容易使用，一开始我们会以 runetime API 为主。
cuda_runtime.h头文件一般就是Runtime API中，运行时的API和其参数的定义。（如果使用Driver API，驱动API则头文件使用cuda.h）。

2、device_launch_parameters.h

要使用threadIdx、blockIdx、blockDim等内置变量时要在头文件里导入device_launch_parameters.h

二、对设备函数的一些操作

1、thrust容器

（1）声明与使用
Trust 提供了两个vector容器：host_vector 与 device_vector。按照命名规则，host_vector位于主机端，device_vector位于GPU设备端。Trust的vector容器与STL中的容器类似，是通用的容器，可以存储任何数据类型，可以动态调整大小。
代码示例：

thrust::device_vector<float>& ZXY

其中，thrust::告诉编译器在thrust命名空间中查找函数与类；device_vector表示使用的容器位于GPU设备端；定义的容器类型是float；容器中存放的是ZXY的地址。
应用实例：

thrust::device_vector<float> d_xds(xds, xds + DNU);

初始化一段内存，d_xds表示初始化的容器名称，其中的xds表示起始指针，xds + DNU表示结束指针。
（2）清除
clear()
代码示例：

thrust::device_vector<float> vol(hvol, hvol + siz);
vol.clear();

其中vol为一个容器，第二行代码clear()表示清楚其中数据，但容器整体结构并未改变。

三、CUDA流

介绍：
1、流可以看成是在设备上work的一个队列，host端将work加入队列，然后继续添加。设备在资源free时，开始调度streams里面的work。
2、CUDA的操作也是在流里面，比如Kernel的启动，内存的拷贝。
3、在同一个流里面的操作是有序的（FIFO），不可以重叠了。
4、在不同的流里面的操作是无序的，可以重叠。
用法：

cudaStream_t stream; //声明一个stream
cudaStreamCreate(&stream); //分配stream
cudaStreamDestroy(stream); //取消分配的stream，在stream中的work完成后同步host端。

除非是特别指定了stream，所有的调用都放在默认流里。这个默认流通常指的是“Stream 0”。
参考网页

四、核函数

官方的文档称为函数执行环境标识符Function execution space specifiers，也就是他指明了这段函数是在哪里被调用的。

__global__

这个前缀修饰的函数是核函数，这些函数在CPU端调用，在GPU上执行。
注意：
（1）修饰的函数必须采用void返回值，并且需要在调用时制定运行的参数 (也就是<<<>>>里的block数和thread)；
（2）任何对__global__函数的调用都必须指定该调用的执行配置。执行配置定义将用于在该设备上执行函数的网格和块的维度，以及相关的流。
（3）函数是异步的，这也代表着函数没被执行完就返回了控制权，所以测量核函数的时间需要同步操作才能获得准确的结果。
代码示例：

 naive_copyToTwoVolumes << <gid, blk >> >(函数的输入参数)

其中，尖括号作用为线程配置，gid类型若为dim3，指定网格的维度和大小，gid.x * gid.y * gid.z 等于所启动的块数量；blk的类型若为 dim3，指定各块的维度和大小，Db.x * Db.y * Db.z 等于各块的线程数量；这些参数并不是传递给设备代码的参数，而是告诉运行时如何启动设备代码。传递给设备代码本身的参数是放在圆括号中传递的。

五、cuda kernal中的block、grid、thread

threadIdx、blockIdx、blockDim和gridDim

threadIdx  //uint3类型，表示一个线程的索引。
blockIdx  //uint3类型，表示一个线程块的索引，一个线程块中通常有多个线程。
blockDim  //dim3类型，表示线程块的大小。
gridDim  //dim3类型，表示网格的大小，一个网格中通常有多个线程块。

以上涉及到线程、线程块、线程格的知识，下面是它们之间关系的示意图。grid表示线程格，一个线程格内包含许多block线程块，而一个线程块内又包含许多线程thread，其中线程块和线程都可以是一维、二维或三维的。

若线程格和线程块都是三维矩阵。这里假设线程格是一个3×4×5的三维矩阵， 线程块是一个4×5×6的三维矩阵，则有：
（1）gridDim
gridDim.x、gridDim.y、gridDim.z分别表示线程格各个维度的大小：

gridDim.x=3
gridDim.y=4
gridDim.z=5

（2）blockDim
blockDim.x、blockDim.y、blockDim.z分别表示线程块中各个维度的大小：

blockDim.x=4
blockDim.y=5
blockDim.z=6

（3）blockIdx
blockIdx.x、blockIdx.y、blockIdx.z分别表示当前线程块所处的线程格的坐标位置
（4）threadIdx
threadIdx.x、threadIdx.y、threadIdx.z分别表示当前线程所处的线程块的坐标位置
通过 blockIdx.x、blockIdx.y、blockIdx.z、threadIdx.x、threadIdx.y、threadIdx.z就可以完全定位一个线程的坐标位置了。
一般的线程索引方式：

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

假如现在我们想得到第6个block中第五个线程，则根据定义式有idz=5+6*8，threadIdx.x线程id为5（0开始），blockIdx.x块id为6，blockDim.x表示块的维度，即一共有8个线程块。

若定义的线程是三维的，则继续定义其他方向有：

int idx = threadIdx.y + blockIdx.y * blockDim.y;
int idy = threadIdx.z + blockIdx.z * blockDim.z;

如果不好理解可以看下图三维grig、block的示意图

                  三维grig、block的示意图

六、cuda拷贝、申请内存

1、cudaMalloc

在GPU内分配内存，这个函数与CPU中的malloc相似，可以先理解malloc的作用，malloc用法为：
int *a = (int )malloc(nsizeof(int))，返回的是一个int型指针，指向大小为n个int型数据的连续内存地址的首地址，可以理解为a是这个数组的首地址。
cudaMalloc也是十分相似的作用，例如cudaMalloc(float(**)&addr,n*sizeof(float))，意思是在GPU内申明一段n个大小的float型数组，addr这个变量中存的就是用户在GPU中声明的float型数组的首地址。

2、cudaMallocPitch

主要解决所需分配的矩阵行大小不是128bytes倍数的问题。函数主要作用与cudamalloc相同，不同的是cudaMallocPitch每行分配256bytes。

cudaMallocPitch(void** devPtr, size_t* pitch, size_t widthInBytes, size_t height)

devPtr：开辟矩阵的数据的头指针
pitch：分配存储器的宽度，以字节为单位(cuda的返回值)
pitch是指设备端分配内存时每行的字节数，其计算方式为将每行元素数乘以每个元素的大小，然后向上对齐到内存对齐的大小。例如：数据宽度是33，数据类型是float，则pitch=334（一个float类型数据占4bytes大小）+内存对齐的大小+1。内存对齐的大小也叫做多分配的字节，一般是32/64，由计算机硬件决定。
width：分配矩阵的列数
height：分配矩阵的行数
在设备上分配widthInBytes * height字节的线性内存，并返回分配内存的指针devPtr。
PS：cudaMallocPitch()就是为了解决每行首地址是否是global memory对齐段的问题，如果用cudaMallocPitch()来分配N=33(即列数为33)的矩阵时（需要33*4=132个字节），那么申请的矩阵每一行大小会变成256个Bytes(0-131为我们需要使用的空间，132-255未使用)，而不是cudaMalloc中的132个Bytes，这样分配以后，每行的首地址将会是与globla memory分段地址对齐的（都是128的整数倍），warp在访问的时候就可以对齐了！因此函数将确保在任何给出的行中对应的指针是连续的。
原因

3、cudaMemcpy

用于在主机（Host）和设备（Device）之间往返的传递数据，使用方式：

主机到设备：cudaMemcpy(d_A,h_A,nBytes,cudaMemcpyHostToDevice)
设备到主机：cudaMemcpy(h_A,d_A,nBytes,cudaMemcpyDeviceToHost)

改行代码表示从h_A存储区中复制nBytes个字节到d_A中。
可以根据最后一个参数看出，源指针和目标指针分别是设备指针-主机指针、主机指针-设备指针。
注意：该函数是同步执行函数，在未完成数据的转移操作之前会锁死并一直占有CPU进程的控制权

4、cudaMemcpy2D

函数调用方式：

cudaMemcpy2D(void* dst, size_t dpitch, const void* src, size_t spitch,
size_t width, size_t height, enum cudaMemcpyKind kind);

dst: 拷贝到的目标矩阵内存头指针
dpitch: dst指向的2D数组中的内存宽度，以字节为单位，是cuda为了读取方便，对齐过的内存宽度，可能大于一行元素占据的实际内存。实际值大小为dpitchsizeof(数据类型)。
**src ：**源矩阵内存头指针。
spitch: src指向的2D数组中的内存宽度，以字节为单位，spitchsizeof(数据类型)。
width: src指向的2D数组中一行元素占据的实际宽度。以字节为单位，等于
width*sizeof(type)
**height: ** src指向的2D数组的行数。
**kind： **拷贝数据的方向,从src指向的内存区域拷贝数据到dst指向的内存区域。
kind表示拷贝方向：
cudaMemcpyHostToHost；
cudaMemcpyDeviceToHost；
cudaMemcpyHostToDevice；
cudaMemcpyDeviceToDevice。

5、cudaMemset2D

函数调用方法：

cudaError_t cudaMemset2D(void * devPtr, size_t pitch, int value, size_t width, size_t height)	

devPtr - 指向2D设备内存的指针
pitch - 2D设备内存的字节数=pitch*sizeof（type）
value - 为每个指定内存设置的值
width - 设置的矩阵宽度（列数），以字节为单位，所以是Width of matrix set (columns in bytes)*sizeof（type）
height - 设置的矩阵高度（行数）

6、cudaMemcpyToSymbol

将数据从host拷贝到global，拷贝方式的不同是由目的内存申请的方式决定的。
申请的是device内存，cudaMemcpyToSymbol拷贝就是从host拷贝到global memory。
申请的是constant内存，cudaMemcpyToSymbol拷贝就是从host拷贝到constant memory。
调用：

cudaMemcpyToSymbol	(const char * symbol,
const void * src, size_t count, size_t offset = 0,
enum cudaMemcpyKind kind = cudaMemcpyHostToDevice	 
)	

Parameters:
symbol - 设备上拷贝目标
src - 原内存地址
count - 拷贝的数据比特数
offset - 起始符号偏移量（bytes）
kind - 拷贝传输方式
可选： cudaMemcpyHostToDevice / cudaMemcpyDeviceToDevice.

类似的函数：连接

七、cuda纹理内存

1、纹理对象的创建

一个纹理对象是用cudaCreateTextureObject()产生的。cudaCreateTextureObject()有4个参数，常用的前三个是
cudaTextureObject_t *texObj：需要生产的纹理对象；
cudaResourceDesc *resDesc：资源描述符，用来获取述纹理数据；
cudaTextureDesc *texDesc：纹理描述符，用来描述纹理参数；

//示例代码，产生纹理对象
cudaTextureObject_t texObj1;

2、 纹理内存的寻址方式

cudaAddressModeClamp：超出范围就用边界值代替，示意： AA | ABCDE | EE
cudaAddressModeBorder：超出范围就用零代替，示意： 00 | ABCDE | 00
cudaAddressModeWrap：重叠模式（循环），示意： DE | ABCDE || AB
cudaAddressModeMirror：镜像模式，示意： BA | ABCDE | ED

3、滤波方式

filtering mode：滤波模式，定义了fetch返回结果的计算方式。有两种模式：cudaFilterModePoint or cudaFilterModeLinear。
cudaFilterModePoint：点模式，返回最接近的一个点，即最近邻插值。
cudaFilterModeLinear：线性模式，即线性插值，对于一维纹理，两点插值；对于二维纹理，四点插值；对于三维纹理，八点插值。线性模式只有在fetch返回浮点类型数据（注意并非指read mode的归一化浮点模式）下才有效

4、纹理内存的访问模式

访问模式有 cudaReadModeNormalizedFloat 和 cudaReadModeElementType 两种。 cudaReadModeNormalizedFloat读取 4 字节整数时会除以 0x8fff（有符号整数）或 0xffff（无符号整数），从而把值线性映射到 [-1.0, 1.0] 区间（有符号整数）或 [0, 1] 区间（无符号整数），读取 2 字节整数时也会发生类似变换，除以 0x8f 或 0xff 。
cudaReadModeElementType则不会发生这种转换，直接读取原数据。

5、cudaResourceDesc资源描述符和cudaTextureDesc纹理描述符

//资源描述符（固定格式）
cudaResourceDesc resDesc;
memset(&resDesc, 0, sizeof(resDesc));//初始化
resDesc.resType = cudaResourceTypeArray;//指定对应设备内存的形式为 CUDA数组
resDesc.res.array.array = cuArray;//CUDA数组 对应的赋值形式

//纹理描述符
cudaTextureDesc texDesc;
memset(&texDesc, 0, sizeof(texDesc));
texDesc.addressMode[0] = cudaAddressModeClamp;
texDesc.addressMode[1] = cudaAddressModeClamp;
texDesc.filterMode = cudaFilterModePoint;//最近邻插值法
texDesc.readMode = cudaReadModeElementType;//若选用cudaFilterModeLinear,则readMode=cudaReadModeNormalizedFloat
texDesc.normalizedCoords = 1;//对坐标进行归一化

接着使用cudaCreateTextureObject()函数创建纹理对象

//创建纹理对象
cudaTextureObject_t tex = 0;
cudaCreateTextureObject(&tex, &resDesc, &texDesc, NULL);

6、CUDA array

(1)cuda array分配
CUDA array是专为纹理获取使用的，是不能通过其它方式进行访问的，其分配函数包括cudaMallocArray和cudaMalloc3DArray。

//以cudaMalloc3DArray为例：
cudaMalloc3DArray //CUDA数组的分配。	
(	
	struct cudaArray ** 	array,
	const struct cudaChannelFormatDesc * 	desc,
	struct cudaExtent 	extent,
	unsigned int 	flags = 0	 
)	

参数:
array - 指向设备内存中已分配数组的指针
desc - 请求通道格式
extent - 请求的分配大小(元素中的宽度字段)
flags - 扩展的标志(目前必须为0)

cudaMalloc3DArray(宽度，高度，深度) 能够分配1D、2D或3D数组，例如：
如果高度和深度范围都为零，则分配一个1D数组。对于1D数组，有效的范围是{(1,8192)，0,0}。同理若
如果只有深度范围为零，则分配2D数组。对于2D数组，有效的范围是{(1,65536)，(1,32768)，0}。
3D为：{(1, 2048), (1, 2048), (1, 2048)}。

(2)纹理通道声明
作用是在设备内存中分配CUDA数组。该函数有一个独立的C和c++ API (c++ API被重载)。
官方使用格式：

cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();

尖括号内为数组数据类型，小括号内为分配数组是所需要的参数，channelDesc为分配的数组名称；
cudaCreateChannelDesc函数的内部结构：

  struct cudaChannelFormatDesc {
          int x, y, z, w;
          enum cudaChannelFormatKind f;
        };
__host__​cudaChannelFormatDesc cudaCreateChannelDesc ( int  x, int  y, int  z, int  w, cudaChannelFormatKind f )

​cudaChannelFormatDesc应该是CUDA频道格式设置的意思，cudaCreateChannelDesc是具体格式设置函数，cudaChannelFormatKind f为设置的数组类型，类型有以下几种：

cudaChannelFormatKindSigned = 0
Signed channel format   （有符号型）
cudaChannelFormatKindUnsigned = 1
Unsigned channel format   （无符号型）
cudaChannelFormatKindFloat = 2
Float channel format   （浮点型）
cudaChannelFormatKindNone = 3
No channel format  （无格式）

其他的四个参数x,y,z,w，对于C API函数，这些是每个通道的比特数。这些可以是颜色通道、空间维度或者任何你想用的东西。从cuda文档中，“返回格式为f的通道描述符，以及每个组件x、y、z和w的比特数。”
示例代码：

cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<T>();

表示在CUDA内存中分配数组T
（3）cudaMemcpy3DParms
CUDA 3D内存拷贝参数

//用法
cudaMemcpy3DParms copyParams = { 0 };
prjSize.width = Width;//定义3D尺寸-长宽高
prjSize.height = Height;
prjSize.depth = Depth;

copyParams是定义的结构体的名称，需要自定义。
注意：
创建extent时，在旧版本中宽度Width以字节为单位，即必须乘上sizeof(DTYPE)，新版本已经不用乘以sizeof(DTYPE)。
（4）make_cudaPitchedPtr
pitched pointer工具，也是绑定3D纹理内存的一种工具：
因为在当前的CUDA版本中，3D的线性内存是无法直接绑定到texture memory，一维的可以，因此，需要将数据首先放进一个3D的CUDA array，然后将3D CUDA array绑定到texture memory上。

copyParams.srcPtr = make_cudaPitchedPtr(
		(void*)sourceData, prjSize.width*sizeof(type)
		prjSize.width, prjSize.height);

参数说明：
(void*)sourceData–需要传递的数据
prjSize.width*sizeof(type)–维度参数1，注意这个需要乘以数据类型的字节大小
prjSize.width, prjSize.height–维度参数2、3，3D数据这里就是前两个维度

具体原理为，如果访问数组元素*u[x][y][z]，通过pitched pointer访问则是u_p[x+y*pitch+ z*pitch*height ]。 显然，这里pitch=width，因此当创建pitched pointer时我们需要将width和height作为参数传递给函数make_cudaPitchedPtr()。

在这里尤其要注意的是，pitched pointer指向的array与传统的C语言数组的存储方式不同，C语言访问元素u[x][y]*[z]是通过u[y*width*depth+x*depth+z]。因此为了正确读取所需元素，建议逆序建立pitched pointer，即depth和height：

copyParams.srcPtr = make_cudaPitchedPtr((void*)u, array_depth*sizeof(float), array_depth, array_height);

此时相当于数组u[x][y][z]被转置，在CUDA3D array中对应元素为u[z][y][x]。
PS：有些程序，在数据建立的时候就是建立成ZXY形式，后期再经过处理出一个ZYX形式，这里就不用逆序建立了。按照第一个建立方式即可。
参考博文
（5）cuda释放
cudaFreeArray()释放。

八、设备端函数

1、cudaSetDevice

显卡选择函数，当计算机中有多个GPU时，可以用这个函数选择应用哪个GPU进行并行计算。

cudaSetDevice(0)

2、CUDA_SAFE_CALL

CUDA中一种保护机制，可以用cudaerror代替。

3、cudaDeviceReset

销毁本进程中所有线程在对应GPU设备上的资源分配与状态。

4、cudaThreadSynchronize();

块内通信:通过共享内存进行通信，块内每个线程都能访问共享存储器，不同块的线程不能通信。
__syncthreads()； 当某个线程执行到该函数时，进入等待状态，直到同一线程块（Block）中所有线程都执行到这个函数为止，即一个__syncthreads()相当于一个线程同步点，确保一个Block中所有线程都达到同步，然后线程进入运行状态。
调用 cudaThreadSynchronize()函数,会使cpu处于等待状态,等待所有的线程都执行完毕.但是,cudaThreadSynchronize()函数并不能在kernel中使用。因为CUDA API和host代码是异步的，cudaDeviceSynchronize可以用来停住CPU并等待CUDA中的操作完成。
也就是说cudaThreadSynchronize()是用来同步线程的，而cudaDeviceSynchronize是用来同步整个设备代码的。
参考网页

5、__syncthreads()

__syncthreads()是cuda的内建函数，用于块内线程通信。
参考：https://blog.csdn.net/jyl1999xxxx/article/details/68950846

6、错误返回函数

cudaGetLastError：
cudaGetLastError函数用于返回最新的一个运行时调用错误。
其返回的值有很多种可能，常见的：
cudaSuccess：成功无异常
cudaErrorSetOnActiveProcess：由异步编程引起的错误
这表明用户在调用非设备管理操作(分配内存和启动内核是非设备管理操作的例子)初始化CUDA运行后，调用了cudaSetDevice()、cudaSetValidDevices()、cudaSetDeviceFlags()等等设置操作。
详细返回值及解析：参考索引文章
cudaGetErrorString：
对于任何CUDA错误，都可以通过函数cudaGetErrorString函数来获取错误的详细信息。