第2章 CUDA编程模型

2.1 CUDA编程模型概述

CUDA编程模型提供了一个计算机架构抽象作为应用程序和硬件之间的桥梁。

通信抽象是程序与编程模型实现之间的分界线，它通过专业的硬件原语和操作系统的编译器或库来实现。

利用编程模型所编写的程序指定了程序的各组成部分是如何共享信息及相互协作的。

编程模型从逻辑上提供了一个特定的计算机架构，通常它体现在编程语言或编程环境中。

CUDA另外利用GPU架构的计算能力提供了以下几个特有功能：
        1.一种通过层次结构在GPU中组织线程的方法；
        2.一种通过层次结构在GPU中访问存储器的方法；

以程序员的角度可从以下几个不同的角度来看待并行计算：
        1.专业领域层；
        2.框架逻辑层；
        3.硬件层；

在编程和算法设计的过程中，你最关心的应是在领域层如何解析数据和函数，以便在并行运行环境中能正确、高效地解决问题。

当进入编程阶段，你的关注点落在如何组织并发线程。在这个阶段，你需要从逻辑层面来思考，以确保你的线程和计算能正确地解决问题。

在硬件层，通过理解线程是如何映射到核心可以帮助提供其性能。

2.1.1 CUDA编程结构

在一个异构环境中包含多个CPU与GPU，每个CPU与GPU的内存都有一条PCI-E总线分隔开。因此，需要注意区分以下内容：
        1.主机：CPU及其内存
        2.设备：GPU及其内存

从CUDA6.0开始，Nvidia提出名为统一寻址(Unified Memory)的编程模型的改进，它连接了主机内存和设备内存空间，可使用单个指针访问CPU和GPU内存，无须彼此之间手动拷贝数据。本章节中，重要的是应学会如何为主机和设备分配内存空间以及如何在CPU和GPU之间拷贝共享数据。

内核(Kernel)是CUDA编程模型的一个重要组成部分，运行在GPU上运行。在此背景下，CUDA调度管理核函数。在主机上，基于应用程序数据以及GPU的性能定义如何让设备实现算法功能。这样做的目的使你专注于算法的逻辑，且在创建和管理大量的GPU线程时不必拘于细节。

多数情况下，主机可以独立地对设备进行操作。内核一旦启动，管理权立刻返回给主机，释放CPU来执行由设备运行的并行代码实现的额外的任务。CUDA编程模型主要是异步的，因此在GPU上进行的运算可以与主机-设备通信重叠。一个典型的CUDA程序包括由并行代码互补的串行代码。

一个典型的CUDA程序遵循以下模式：
        1.将数据从CPU拷贝到GPU内存；
        2.调用核函数对存储在GPU内存中的数据进行拷贝；
        3.将数据从GPU内存传送回CPU内存；

2.1.2 内存管理

cudaMalloc()函数负责向设备分配一定字节的线性内存，并以devPtr(void**)的形式返回指向所分配的内存指针。

cudaMemcpy()函数从src指向的源存储区复制一定数量的字节到dst指向的目标存储区。复制方向由kind指定。此函数是以同步方式执行，因为在cudaMemcpy函数返回以及传输操作完成之前主机应用程序是阻塞的。

除了内核启动之外的CUDA调用都会返回一个错误的枚举类型cudaError_t。可以用cudaGetErrorString()函数将错误代码转换为可读的错误信息。

一个简化的GPU内存模型，主要包括两部分：全局内存和共享内存。

CUDA编程模型最显著的一个特点就是揭示了内存层次结构。每个GPU设备都有用于不同用途的存储结构。

GPU内存层次中，最主要的两种内存是全局内存和共享内存。全局内存类似于CPU的系统内存，而共享内存类似于CPU的缓存。然而，GPU的共享内存可以由CUDA C的内核直接控制。

nvcc封装了几种内部编译工具，CUDA编译器命令行选项在不同阶段启动不同的工具完成编译工作。-Xcomplier用于指定命令行选项是指向C编译器还是预处理编译器。 -std=c99是指按照C99标准编写的。

cudaFree()函数释放GPU的内存。

数组求和案例：

#include <cstdlib>
#include <cmath>
#include <stdio.h>
#include <time.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <cuda_runtime_api.h>

__constant__  int d_nElem = 1024;


#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)); \
	}\
}

void sumArraysOnHost(float* A, float* B, float* C, const int N)
{
	for (int idx = 0; idx < N; ++idx)
		C[idx] = A[idx] + B[idx];
	printf("OK!\n");
}

void initialData(float* ip, int size)
{
	time_t time(0);
	srand(time);
	for (int i = 0; i < size; ++i)
	{
		ip[i] = static_cast<float>(rand()) / RAND_MAX;
	}
}

void checkResult(float* hostRef, float* gpuRef, const int N)
{
	double epsilon = 1.0e-6;
	int match = 1;
	for (int i = 0; i < N; ++i)
	{
		if (abs(hostRef[i] - gpuRef[i]) > epsilon)
		{ 
			printf("Error: hostRef[%d]: %f, gpuRef[%d]: %f\n", i, hostRef[i], i, gpuRef[i]);
			match = 0;
			break;
		}
	}

	if (match)
		printf("Arrays match!\n\n");

}

__global__ void sumArraysOnDevice(float* A, float* B, float* C)
{
	unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;

	if (x < d_nElem)
	{
		C[x] = A[x] + B[x];
	}
}

int main(int argc, char** argv)
{
	// set Device
	int nGPUs;
	cudaGetDeviceCount(&nGPUs);

	// declare variables
	int nElem = 1 << 18;
	int nElem_local = (nElem + nGPUs - 1) / nGPUs;


	size_t nBytes = nElem * sizeof(float);

	float* h_A, * h_B, *h_C, *gpu_ref;
	float* d_A, * d_B, * d_C;

	// memeory allocation on host
	h_A = static_cast<float*>(malloc(nBytes));
	h_B = static_cast<float*>(malloc(nBytes));
	h_C = static_cast<float*>(malloc(nBytes));
	gpu_ref = static_cast<float*>(malloc(nBytes));

	// memory allocation on device
	cudaMalloc((void**)&d_A, nBytes);
	cudaMalloc((void**)&d_B, nBytes);
	cudaMalloc((void**)&d_C, nBytes);

	// data initialization
	initialData(h_A, nElem);
	initialData(h_B, nElem);

	// data transfer from host to device
	cudaMemcpy(d_A, h_A, nBytes, cudaMemcpyHostToDevice);
	cudaMemcpy(d_B, h_B, nBytes, cudaMemcpyHostToDevice);
	cudaMemcpyToSymbol(&d_nElem, &nElem, 1, sizeof(int));

	sumArraysOnHost(h_A, h_B, h_C, nElem);

	// set kernel function configuration on device
	dim3 block(128);
	dim3 grid((nElem + block.x - 1) / block.x);

	// run kernel
	clock_t startTime, elapseTime;
	startTime = clock();
	for (int i = 0; i < nGPUs; ++i)
	{
		unsigned int offset = i * nElem_local;
		sumArraysOnDevice << <grid, block >> > (d_A + offset, d_B + offset, d_C + offset);
		cudaDeviceSynchronize();
	}
	elapseTime = (clock() - startTime) / CLOCKS_PER_SEC;
	printf("%9.8f\n", static_cast<float>(elapseTime));

	// data transfer from device to host
	cudaMemcpy(gpu_ref, d_C, nBytes, cudaMemcpyDeviceToHost);

	// result check
	checkResult(h_C, gpu_ref, nElem);

	// free memory
	free(h_A);
	free(h_B);
	free(h_C);
	free(gpu_ref);
	cudaFree(d_A);
	cudaFree(d_B);
	cudaFree(d_C);

	cudaDeviceReset();

	return 0;
}

不同的存储空间：使用CUDA进行编程最常犯错就是对不同内存空间的不恰当的引用。对于在GPU上被分配的内存来说，设备指针在主机代码可能并没有被引用。比如你执行了gpuRef = d_C。但在cuda6.0提出了统一寻址后，则可以使用一个指针来同时访问CPU与GPU的内存。

2.1.3 线程管理

了解如何组织线程是CUDA编程的一个关键部分。CUDA编程模型是一个三层的线程层次结构，由线程、线程块和线程块网格构成。

由一个内核启动的所有线程统称为一个网格。同一个网格中所有线程共享相同的全局内存空间。一个网格由多个线程块构成，一个线程块包含一组线程，同一线程块内的线程协作可以通过以下方式来实现：
        1.同步；
        2.共享内存；

注意不同块内的线程不能直接协作。

线程依靠以下两个坐标来区分彼此：
        1. blockIdx;(线程块在线程网格内索引)；
        2. threadIdx;(块内的线程索引)；

基于上述坐标，你可以将部分数据分配给不同的线程。

线程网络坐标变量是基于uint3定义的CUDA内置的向量类型，是一个包含3个无符号整数的结构。

网格和块的维度由下列两个内置变量指定。blockDim和gridDim。

通常，一个线程网格会被组织成为线程块的二维数组形式，一个线程块会被组织成线程的三维数组形式。

线程网格和线程块均使用3个dim3类型的无符号整型字段，而未使用的字段将被初始化为1且忽略不计。

在CUDA程序中有两组不同的网格和块变量：手动定义的dim3数据类型和预定义的uint3数据类型。手动定义的dim3类型的网格和块变量仅在主机端可见，而uint3类型的内置预初始化的网格和块变量仅在设备端可见。

dim3 block(32, 16);
dim3 grid(12,34,23);

区分主机端和设备端的网格和块变量的访问是很重要的。例如，声明一个主机端的块变量，按如下定义它的坐标并对其进行访问: block.x/y/z。

在设备端，已经预定义了内置块变量的大小：blockDim.x/y/z。

对于一个给定的数据大小，确定网格和块尺寸的一般步骤为;
        1. 确定块的大小；
        2. 在已知数据大小和块大小的基础上计算网格维度；

要确定块尺寸，通常需要考虑：
        1.内核的性能特性；
        2.GPU资源的限制；

CUDA的特点之一就是通过编程模型揭示了一个多层的线程层次结构，而且网格和块的维数会影响性能。

网格和块的维度存在几个限制因素，对于块大小的一个主要限制因素就是可利用的计算资源，如寄存器，共享资源等。某些限制可以通过查询GPU设备撤回。

网格和块从逻辑上代表了一个核函数的线程层次结构。

2.1.4 启动一个CUDA核函数

CUDA内核调用时C语言函数调用的延伸，<<<>>>运算符是核函数的执行配置。

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

CUDA编程模型揭示了线程层次结构。利用执行配置可以指定线程在GPU上调度运行的方式。执行配置的第一个值是网格维度，也就是启动块的数目。第二个值是块维度，也就是每个块中线程的数目。通过指定网格核块的维度，可以进行以下配置;
        1. 内核中线程的数目；
        2. 内核中使用的线程布局；

同一个块中的线程之间可以相互协作，不同块内的线程不能直接协作。

由于数据在全局内存中是线性存储的，可以利用线程索引进行以下操作;
        1. 在网格中标识一个唯一的线程；
        2. 建立线程和数据元素之间的映射关系；

核函数的调用与主机线程是异步的，即核函数调用后，控制权立刻返回主机端。

也可以调用以下函数来强制主机端程序等待所有的核函数执行结束，cudaDeviceSynchronize()函数。

一些CUDA运行时API在主机和设备之间是隐式同步的。当使用cudaMemcpy()函数即是如此。

2.1.5 编写核函数

核函数是在设备端执行的代码。在核函数中，需要为一个线程规定要进行的计算以及要进行的数据访问。

__global__ void kernel_name(argument list);

以上是用__global__声明的核函数。核函数必须返回一个void类型。

限定符	执行	调用	备注
__global__	设备端	主机端调用，也可以从计 算能力为3的设备中调用	必须有一个void返回类型
__device__	设备端	设备端
__host__	主机端	主机端	可以省略

__device__和__host__限定符可以一齐使用，这样函数可以同时在主机端和设备端进行编译。

CUDA核函数的限制;
        1. 只能访问设备内存；
        2. 必须返回void类型；
        3. 不支持可变数量的参数；
        4. 不支持静态变量；
        5. 显示异步行为；

2.1.6 验证核函数

需要一个主机函数来验证核函数的结果。

除了许多可用的调试工具外，还有两个非常简单实用的方法可用于验证核函数。首先，使用printf函数。其次，可以将执行参数设置为<<<1, 1>>>，因此强制用一个块和一个线程执行核函数，这模拟了串行执行程序。对于调试和验证结果是否正确是非常有用的。

2.1.7 处理错误

由于CUDA调用是异步的，所以有时可能很难确定某个错误是由那一步程序引起的。定义一个错误处理宏封装所有的CUDA API调用，这简化了错误检查过程。

检查宏仅在调试为目的，因为在和核函数启动后添加这个检查点会阻塞主机端线程，使该检查点成为全局屏障。

2.1.8 编译和执行

nvcc  name.cu -o name

2.2 给核函数记时

最简单的方法是在主机端使用一个CPU或GPU计时器来计算内核的执行时间。

2.2.1 用CPU计时器计时

可以使用gettimeofday系统调用来创建一个CPU计时器，以获取系统的时钟时间。程序中需要添加time.h文件。

由于GPU的可扩展性，你需要借助块和线程的索引来计算一个按行优先的数组索引i，并在核函数中添加限定条件来检验索引值是否越界。(if i < N)

2.2.2 用nvprof工具计时

自CUDA5.0以来，NVDIA提供了nvprof的命令行分析工具，可以帮助从应用程序的cpu和gpu活动情况中获取时间线信息，其包括内核执行、内存传输以及CUDA API的调用。其用法如下;

nvprof [nvprof_args] <application> [application_args]

nvprof显示的结果前半部分来自于程序的输出，后半部分来自于nvprof的输出。

nvprof是一个能帮助你理解在执行应用程序时所花费的时间主要用在何处的强大工具。

对于HPC工作负载，理解程序中通信比的计算是非常重要的。如果你的应用程序用于计算的时间大于数据传输的所用的时间，那么或许可以压缩这些操作，并完全隐藏与传输数据有关的延迟。如果你的应用程序用于计算的时间少于数据传输所用的时间，那么需要尽量减少主机和设备之间的传输。也可以用CUDA流和事件来压缩计算量和通信量。

由nvprof得到的计数器可以帮助你获取应用程序的指令和内存吞吐量。然后可以判定你的应用程序的性能受限于算法还是受限于内存带宽。

2.3 组织并行线程

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

通常情况下，一个矩阵用行优先的方法在全局内存进行线性存储。在矩阵加法的核函数中，一个线程通常被分配一个数据元素来处理。首先要完成的任务是使用线程索引从全局内存访问指定的数据。通常情况下，需要管理三种索引;
        1. 线程和块索引；
        2. 矩阵中给定点的坐标；
        3. 全局线性内存中的偏移量；

第一步，可以用以下公式把线程和块索引映射到矩阵坐标上；

ix = threadIdx.x + blockIdx.x * blockDim.x;
iy = threadIdx.y + blockIdx.y * blockDim.y;

第二步：可以用一下公式把矩阵坐标映射到全局内存中的索引上/存储单元上

idx = ix + iy * nx;

// 2D grid and 2D block
__global__ void sumMatrixOnDevice(float* A, float* B, float* C)
{
	unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
	unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
	unsigned int idx = x + y * d_N;

	if (x < d_N && y < d_M)
	{
		C[idx] = A[idx] + B[idx];
	}
}

// 1D grid and 1D block
__global__ void sumMatrixOnDevice(float* A, float* B, float* C)
{
	unsigned int ix = blockDim.x * blockIdx.x + threadIdx.x;

	if (idx < d_N)
	{
        for (int iy = 0; iy < d_M; iy++)
        {
            unsigned int idx = ix + iy * d_N;
		    C[idx] = A[idx] + B[idx];
        }

	}
}

// 2D grid and 1D block
__global__ void sumMatrixOnDevice(float* A, float* B, float* C)
{
	unsigned int ix = blockDim.x * blockIdx.x + threadIdx.x;
    unsigned int iy = blockIdx.y;
    unsigned int idx = ix + iy * d_N;

	if (idx < d_N && iy < d_M)
	{
	    C[idx] = A[idx] + B[idx];
	}
}

注意，2Dgrid and 1D block唯一优点是每个线程减少一个整数乘法和一次整数加法的计算。

2.4 设备管理

在本节中，你将通过以下两种方式查询和管理GPU设备;
        1. CUDA运行时API函数；
        2. NVIDIA系统管理界面命令行实用程序；

2.4.1 使用运行时API查询GPU信息

cudaGetDeviceProperties()函数查询GPU的所有信息；

2.4.2 确定最优GPU

通过查询每个设备的流多处理器的个数

2.4.4 在运行时设置设备

使用环境变量CUDA_VISIBLE_DEVICES，就可以在运行时指定所选的GPU且无需更改应用程序。例如CUDA_VISIBLE_DEVICES = 2，系统将会屏蔽其他GPU。CUDA_VISIBLE_DEVICES = 2，3，系统将会只适用ID为2,3的设备。

2.5 总结

与C语言中的并行编程相比，CUDA程序中的线程层次结构是其独有的结构。通过一个抽象的多级线程层次结构，CUDA能够控制一个大规模并行环境。你也学习到网格和线程块的尺寸对于内核性能有很大影响。

学习如何组织线程是CUDA编程的重点之一，理解网格和线程块的启发式的最好方法就是编写程序。