编写函数计算度为i的结点数目_CUDA程序设计入门（二）—

上期内容：SM和SP

原本这期的标题为“任务分解与执行模式”。因为该标题覆盖的知识面太广，导致坑很难填，所以我改成了“核函数”，本期重点讨论核函数。

这期我们将接触线程（thread）、线程束（threading warp）、线程块（threading block）、线程网格（threading grid）、核函数（kernel）等概念，内容较上期更多。知乎及各大博客上不少同类文章对这些概念进行了详细的介绍，大家也可读一读其他作者的文章。

1. 任务分解

任务分解是并行程序设计的必要步骤，把一个串行的任务划分为多个子任务，使子任务并行执行。任务分解可分为域分解和功能分解。其中域分解将算例使用到的数据分割为多个数据块，由不同的调度单位处理不同的数据块；而功能分解将算例分割为多个步骤，由不同的调度单位处理不同的步骤。在实际应用中，由于域分解较易实现负载均衡，且适用于GPU、Xeon Phi、FPGA等加速卡，多数情况下采用域分解。

但是给定一个算例，你如何对它做域分解呢？这就和调度单位有关了。CUDA的设计师们设计了多层次的调度单位，包括线程、线程束、线程块、线程网格。其中，每个CUDA核函数占用一个线程网格的资源，线程网格下又包括线程块等其他调度单位。

CUDA调度单位的层次关系

2. 核函数示例

下面我们通过实现一个简单的CUDA程序来体会如何使用CUDA核函数。这个程序的功能是获取3维网格中的所有线程的线程号。线程号的计算方法如下：

其中线程块号表示线程块在网格内的相对位置，通过下式计算：

其中，gridDim.x、gridDim.y、gridDim.z分别为当前网格中x轴、y轴、z轴方向的线程块数目（1维线程块和2维线程块可分别看作是gridDim.y=1、gridDim.z=1和gridDim.z=1的情况），blockIdx.x、blockIdx.y、blockIdx.z分别为线程块在当前网格中x轴、y轴、z轴方向上的坐标。

线程块大小就是一个线程块的尺寸，通过下式计算：

其中，blockDim.x、blockDim.y、blockDim.z分别为当前线程块中x轴、y轴、z轴方向的线程数目。

块内线程号是线程相对于其所在分块的位置，起到“偏移地址”的作用，通过下式计算：

其中，threadIdx.x、threadIdx.y、threadIdx.z分别表示线程在所在线程块中x轴、y轴、z轴方向上的坐标。

该程序的任务分解方案是将所有的线程号平均分摊每个线程计算，每个线程负责计算自己的线程号。

2.1 编写核函数

该核函数的功能为计算出当前线程号，并保存计算结果。核函数不能返回任何值（返回类型必须是void），因此我们需要借助参数将计算结果“带出去”。那么对应的代码可写成：

#include <cuda_runtime.h>
#include <device_launch_parameters.h>

/**
 * @brief cuda kernel -- compute id of a thread
 * @param array that stores thread ids
 * @return return is not allowed
 */
__global__ void computeThreadID(unsigned int* threadID);

__global__ void computeThreadID(unsigned int* threadID)
{
    int tid = (blockIdx.z * gridDim.y * gridDim.x +
        blockIdx.y * gridDim.x + blockIdx.x) *
        blockDim.z * blockDim.y * blockDim.x +
        threadIdx.z * blockDim.y * blockDim.x +
        threadIdx.y * blockDim.x + threadIdx.x;
    threadID[tid] = tid;
}

其中，“__global__”前缀告诉NVCC编译器这是一个由CPU和GPU共同执行的核函数，由CPU调用并将函数中的指令发射给GPU。

除“__global__”外，类似的前缀还有：

“__host__”表示函数仅在host端（CPU）执行，相当于不加任何前缀。
“__device__”表示函数仅在device端（GPU）执行，由“__device__”修饰的函数可以被核函数调用。

核函数对应网格内所有的线程共享核函数中的指令，但不同的线程中blockIdx、threadIdx的值是不同的。

2.2 调用核函数

这一步需要对网格和线程块进行设置，由于本例的数据规模可以在GPU可承受的范围内随意设置，所以我们随便把网格大小为2*3*4，线程块大小设为5*6*7。另外，我们还需要通过cudaMalloc、cudaFree函数分别在GPU上分配/释放内存空间，通过cudaMemcpy在host端和device端之间传输数据。最后还要把结果打印出来。那么对应的代码可写成：

/*********************************************************************
 * @file 	check_gpuinfo.cu
 * @brief 	fetch thread ID
 * @author 	Bin Qu
 * @email 	benquickdenn@foxmail.com
 * @date	2019-12-1
 * you can reedit or modify this file
*********************************************************************/

#include <cstdio>
#include <cstdlib>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>

#include "kernels.h"

/**
 * @brief display thread ID
 */
void dispThreadID();

/**
 * @brief main entry
 * @return exit status
 */
int main(int argc, char** argv)
{
	dispThreadID();
	return 0;
}

void dispThreadID()
{
	/* initialize grid */
	dim3 gridSize(2, 3, 4);

	/* initialize block */
	dim3 blockSize(5, 6, 7);

	/* allocate memory on host */
	const unsigned int memSpace = gridSize.z * gridSize.y * gridSize.x *
		blockSize.z * blockSize.y * blockSize.x;
	unsigned int* threadID;
	threadID = (unsigned int*)std::malloc(memSpace * sizeof(unsigned int));
	
	/* allocate memory on device */
	unsigned int* cuThreadId;
	cudaMalloc((void**)&cuThreadId, memSpace * sizeof(unsigned int));

	/* copy data from host to device */
	/* in this application, there is no need to copy data from host to device */
	//cudaMemcpy(cuThreadId, threadID, memSpce * sizeof(unsigned int), cudaMemcpyHostToDevice)

	/* call kernel */
	computeThreadID<<<gridSize, blockSize>>>(cuThreadId);

	/* copy data from device to host */
	cudaMemcpy(threadID, cuThreadId, memSpace * sizeof(unsigned int), cudaMemcpyDeviceToHost);

	/* free memory on device */
	cudaFree(cuThreadId);

	/* display thread id */
	for (int i = 0; i < memSpace; i++)
		std::printf("%d, ", threadID[i]);
	std::printf("rn");

	/* free memory on host */
	std::free(threadID);
}

由于线程的数目过多，这里不展示完整的结果。最大的线程号为5039，它刚好等于

参考资料

[1] 陈国良. 并行计算:结构·算法·编程[M]. 2011.

下期预告：

实现简单的应用：向量加