CUDA_01组织并行线程

CUDA编程中最重要的一环就是如何让数据并行计算，要想实现较为高效的并行计算就需要在CUDA编程时手动为每个数据分配线程，这就涉及到如何在不同的grid 和 block 中计算出合适线程的id号。

CUDA组织并行线程

文章参考

概念

CUDA编程中最重要的一环就是如何让数据并行计算，要想实现较为高效的并行计算就需要在CUDA编程时手动为每个数据分配线程，这就涉及到如何在不同的grid 和 block 中计算出合适线程的id号。

回顾GPU中一个 kernel 的大致框架：

所在一般在编程时需要事先声明grid, block, thread的大小：

dim3 block(4,2);  // 声明一个 block 其中包含 4 行 2 列 = 8 个线程
dim3 grid(2,3);  // 声明一个 grid 其中包含 2 行 3 列 = 6 个block
// 所以总共可分配的线程数为： 8 * 6 = 48 个

所以在得到一个block中一个线程的索引号（threadIdx）后，需要还需要根据 gridDim. gridIdx, blockDim, blockIdx 计算获得其全局的索引号。详细: grid、block、thread的关系及thread索引的计算

这边以简单的二维block为例介绍：

计算方式就如图中所示。

因为CUDA 每一个线程执行相同的代码，就是异构计算中的多线程单指令， 如果每个不同的线程执行同样的代码，又处理同一组数据，CUDA常用的做法是让不同的线程对应不同的数据，也就是用线程的全局全局标号对应不同组的数据。因为无论时设备还是主机内存都是线性存在的，比如二维矩阵 $(8 \times 6)$， 存储在内存中可能就是如下结构：

我们要做的管理就是：

• 矩阵和块索引（计算线程的全局索引）；

• 计算矩阵中给定点的坐标 (ix, iy);

• (ix, iy) 对应的线性内存的位置；

线性位置的计算方式：

$$idx = ix + iy \times nx$$

这样的分配方式就将数组矩阵中 (ix, iy) 位置处的数据分配到线程标号为 ix,iy 处去计算。

例子

随机生成一个 6 行 8列的数组，总共有 48 个数字，如下：

实际内存是一维线性的，如下：

$$[29.6 , ... , 22.16 , 26.49 , ... , 23.45 , 30.92 , ... , ...]$$

block 和 grid 的大小为：

dim3 blcok(4,2);
dim3 grid(2,3);

这里创建了一个二维block和 二维线程，即使用 2 x 3 个block ，每个 block中再开辟 4x 2 个线程， 则 2 x 3 x 4 x 2 = 48，然后将这48个数字分配进去，每个数字一个线程,就相当于每个线程进行一个数字的操作，进行的计算的并行化。

有如下形式：

按照下面的方式打印索引：

__global__ void printThreadIndex(float *A, const int nx, const int ny)
{
	int ix = threadIdx.x + blockIdx.x*blockDim.x;  
	int iy = threadIdx.y + blockIdx.y*blockDim.y;
	unsigned int idx = iy * nx + ix; 
//nx = 8, ny= 6 在进一步计算全局线程标号，是为了将这48个数据分别匹配到不同的线程中去
	printf("thread_id(%d,%d) block_id(%d,%d) coordinate(%d,%d)"
		"global index %2d ival %2d\n", threadIdx.x, threadIdx.y,
		blockIdx.x, blockIdx.y, ix, iy, idx, A[idx]);
}

//dim3 blcok(4,2);
//dim3 grid(2,3);
//int nx = 8;
//int ny = 6;
printThreadIndex << <grid, block >> > (A_dev, nx, ny);

对应的计算示意图如下：

程序打印输出结果如下：

设置不同的block或者 grid大小，计算的性能也会不一样的，如下：

可以看到在GPU上计算明显比CPU上计算快很多。

附录

时间计算的比较：

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "freshman.h"
#include <stdio.h>

void sumMatrix2D_CPU(float * MatA, float * MatB, float * MatC, int nx, int ny)
{
	float * a = MatA;
	float * b = MatB;
	float * c = MatC;
	for (int j = 0; j < ny; j++)
	{
		for (int i = 0; i < nx; i++)
		{
			c[i] = a[i] + b[i];
		}
		c += nx;
		b += nx;
		a += nx;
	}
}
__global__ void sumMatrix(float * MatA, float * MatB, float * MatC, int nx, int ny)
{
	int ix = threadIdx.x + blockDim.x*blockIdx.x;
	int iy = threadIdx.y + blockDim.y*blockIdx.y;
	int idx = ix + iy * ny;
	if (ix < nx && iy < ny)
	{
		MatC[idx] = MatA[idx] + MatB[idx];
	}
}

int main(int argc, char** argv)
{
	printf("strating...\n");
	initDevice(0);
	int nx = 1 << 12;  //2**12 = 4096
	int ny = 1 << 12;  //2**12 = 4096
	int nxy = nx * ny;  // 4096 * 4096
	int nBytes = nxy * sizeof(float);  // nxy 共占用多少个字节

	//Malloc
	float* A_host = (float*)malloc(nBytes);  // 分配指针空间
	float* B_host = (float*)malloc(nBytes);
	float* C_host = (float*)malloc(nBytes);
	float* C_from_gpu = (float*)malloc(nBytes);
	initialData(A_host, nxy);
	initialData(B_host, nxy);

	//cudaMalloc
	float *A_dev = NULL;
	float *B_dev = NULL;
	float *C_dev = NULL;
	CHECK(cudaMalloc((void**)&A_dev, nBytes));
	CHECK(cudaMalloc((void**)&B_dev, nBytes));
	CHECK(cudaMalloc((void**)&C_dev, nBytes));


	CHECK(cudaMemcpy(A_dev, A_host, nBytes, cudaMemcpyHostToDevice));
	CHECK(cudaMemcpy(B_dev, B_host, nBytes, cudaMemcpyHostToDevice));

	int dimx = 32;
	int dimy = 32;

	// cpu compute
	cudaMemcpy(C_from_gpu, C_dev, nBytes, cudaMemcpyDeviceToHost);
	double iStart = cpuSecond();
	sumMatrix2D_CPU(A_host, B_host, C_host, nx, ny);
	double iElaps = cpuSecond() - iStart;
	printf("CPU Execution Time elapsed %f sec\n", iElaps);

	//==================================== 2d block and 2d grid ================================
	dim3 block_0(dimx, dimy); // blockDim.x = 32， blockDim.y = 32, threadIdx.x(y) in [0,31]
	// (nx - 1)/ block_0.x + 1 相当于 int(nx/block_0.x)
	dim3 grid_0((nx - 1) / block_0.x + 1, (ny - 1) / block_0.y + 1); //128， 128，  gridDim.x(y) = 128, blockIdx.x(y) in [0,128]
	iStart = cpuSecond();
	sumMatrix << <grid_0, block_0 >> > (A_dev, B_dev, C_dev, nx, ny);
	CHECK(cudaDeviceSynchronize());
	iElaps = cpuSecond() - iStart;
	printf("GPU Execution configuration grid and block shape: <<<(%d,%d),(%d,%d)>>> Time elapsed %f sec\n",
		grid_0.x, grid_0.y, block_0.x, block_0.y, iElaps);
	CHECK(cudaMemcpy(C_from_gpu, C_dev, nBytes, cudaMemcpyDeviceToHost));
	checkResult(C_host, C_from_gpu, nxy);
	//============================================================================================

	//=====================================1d block and 1d grid ==================================
	dimx = 32;
	dim3 block_1(dimx);  // 32
	dim3 grid_1((nxy - 1) / block_1.x + 1);  // (4096 * 4096 -1)/ 32 + 1 = 524288
	iStart = cpuSecond();
	sumMatrix << <grid_1, block_1 >> > (A_dev, B_dev, C_dev, nx*ny, 1);
	CHECK(cudaDeviceSynchronize());
	iElaps = cpuSecond() - iStart;
	printf("GPU Execution configuration grid and block shape: <<<(%d,%d),(%d,%d)>>> Time elapsed %f sec\n",
		grid_1.x, grid_1.y, block_1.x, block_1.y, iElaps);
	CHECK(cudaMemcpy(C_from_gpu, C_dev, nBytes, cudaMemcpyDeviceToHost));
	checkResult(C_host, C_from_gpu, nxy);
	//============================================================================================

	//===================================== 2d block and 1d grid =================================
	dimx = 32;
	dim3 block_2(dimx); // 32
	dim3 grid_2((nx - 1) / block_2.x + 1, ny);  // 128, 4096
	iStart = cpuSecond();
	sumMatrix << <grid_2, block_2 >> > (A_dev, B_dev, C_dev, nx, ny);
	CHECK(cudaDeviceSynchronize());
	iElaps = cpuSecond() - iStart;
	printf("GPU Execution configuration grid and block shape: <<<(%d,%d),(%d,%d)>>> Time elapsed %f sec\n",
		grid_2.x, grid_2.y, block_2.x, block_2.y, iElaps);
	CHECK(cudaMemcpy(C_from_gpu, C_dev, nBytes, cudaMemcpyDeviceToHost));
	checkResult(C_host, C_from_gpu, nxy);
	//=============================================================================================
	cudaFree(A_dev);
	cudaFree(B_dev);
	cudaFree(C_dev);
	free(A_host);
	free(B_host);
	free(C_host);
	free(C_from_gpu);
	cudaDeviceReset();
	return 0;
}