《CUDA C编程权威指南》—读书笔记—2.3线程组织

GPU高速运算的很大一个特点就是，利用足够多的线程来隐藏内存操作等各类延迟。选择合适的网格和块来组织线程结构，对内核性能会产生极大的影响。以矩阵加法为例。
在全局内存中，一个矩阵往往以行优先的方式进行存储，即下一行接着上一行以线性方式在物理内存中存储。如下图，为一个8*6的矩阵。

在kernel中，一个线程被分配一个矩阵元素进行运算 。对于一个二维组织的线程来说，我们需要在编程时明确知道一下内容:
1、线程和块索引
2、矩阵中给定点的坐标
3、全局线性内存的偏移量
当我们确定了线程结构后，对于一个二维的矩阵（图像数据也是类似），我们可以将线程索引与矩阵坐标一一映射。

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

再将索引映射到实际的内存中，取得实际的坐标

idx=ix+iy*nx

其实不难理解，就是相当于将一个矩阵（图像）分块，利用块的维度和块内的绝对坐标来映射原来未分块时的坐标。以上的nx是矩阵在x方向上的数量，ny是在y方向上的数量。若是类似opecv读取的图像行列，nx对应cols。
以二维网格二维块的矩阵加法为例子。

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include<algorithm>

int nx = 1 << 12;
int ny = 1 << 12;

void arrayGenerate(float* array, int n) {
	srand((unsigned int)time(0));
	for (int i = 0; i < n; i++) {
		array[i] = (float)(rand() & 0xFF) / 10.0f;
	}
}

void arraySumCpu(float* A, float* B, int n, float* res) {
	for (int i = 0; i < n; i++) {
		res[i] = A[i] + B[i];
	}
}

void comp(float *CPU, float* GPU, int n) {
	double epsilon = 1.0E-8;
	bool match = 1;
	for (int i = 0; i < n; i++) {
		if (abs(CPU[i] - GPU[i]) > epsilon) {
			match = 0;
			printf("Arrays do not match\n");
			printf("host %5.2f gpu %5.2f at currentt %d\n", CPU[i], GPU[i], i);
			break;
		}
	}
	if (match) printf("Array  match\n");
}

__global__ void arraySumGpu(float* d_A, float* d_B, float* d_C, int x,int y) {
	unsigned int ix = threadIdx.x+blockIdx.x*blockDim.x;
	unsigned int iy = threadIdx.y + blockIdx.y*blockDim.y;
	unsigned int idx = ix + iy*x;
	if (ix<x&&iy<y) 
		d_C[idx] = d_A[idx] + d_B[idx];
}

int main(int argc, char **argv) {
	//printf("%s starting...\n", argv[0]);

	//set device
	int dev = 0;
	cudaSetDevice(dev);

	float *A, *B, *gpuRes;
	A = (float*)malloc(sizeof(float)*nx*ny);
	B = (float*)malloc(sizeof(float)*nx*ny);
//	hostRes = (float*)malloc(sizeof(float)*nx*ny);
	gpuRes = (float*)malloc(sizeof(float)*nx*ny);

	//初始化数组
	arrayGenerate(A, nx*ny);
	arrayGenerate(B, nx*ny);


	float *d_A, *d_B, *d_C;
	cudaMalloc((float**)&d_A, sizeof(float)*nx*ny);
	cudaMalloc((float**)&d_B, sizeof(float)*nx*ny);
	cudaMalloc((float**)&d_C, sizeof(float)*nx*ny);

	cudaMemcpy(d_A, A, sizeof(float)*nx*ny, cudaMemcpyHostToDevice);
	cudaMemcpy(d_B, B, sizeof(float)*nx*ny, cudaMemcpyHostToHost);

	int block_x = 32;
	int block_y = 32;
	dim3 block(block_x,block_y);
	dim3 grid((nx+block.x-1)/block.x, (ny+block.y-1)/block.y);
	
	cudaEvent_t start, stop;
	cudaEventCreate(&start);
	cudaEventCreate(&stop);
	float elapsedTime;
	cudaEventRecord(start, 0);
	arraySumGpu << <grid, block >> > (d_A, d_B, d_C, nx, ny);
	cudaEventRecord(stop, 0);
	cudaEventSynchronize(stop);
	cudaEventElapsedTime(&elapsedTime, start, stop);
	printf("time：%fms\n", elapsedTime);

	cudaMemcpy(gpuRes, d_C, sizeof(float)*nx*ny, cudaMemcpyDeviceToHost);

	//comp(hostRes, gpuRes, nx*ny);

	cudaFree(d_A);
	cudaFree(d_B);
	cudaFree(d_C);
	free(A);
	free(B);
	free(gpuRes);

	system("pause");
	return 0;
}

由于上述代码运行的平台为Fermi架构的GT 820M，性能是在够呛。书上使用的是1638416384的数据，在我这台老爷机上，哪怕将block的thread数开启到最大的3232=1024，二维grid的维度也超过了我的硬件限制。因此，缩减数据量为40964096。不同线程结构运行时间如下。1616的内核配置会让块数又超过限制，这里就不让老爷机献丑啦哈哈哈。

内核配置	内核运行时间
（32，32）	25.843519ms
（32，16）	17.604671ms