第二章

1.坐标映射

从前面的博文我们可以知道，global memory是可以划分成网格（一个程序对应一个网格），网格由块组成，块由线程组成。一个块内的线程可以相互访问，相互等待。

通过对前面并行向量加法的分析，我们知道，网格和块的大小会影响核函数的性能，这一篇博文将随着书本来探究如何组织网格和块从而获得更高效的性能。

还是以矩阵加法为例，在矩阵加法中，传统的是使用二维网格和二维块的布局来分配线程的归属，但这种分配方法不能获得最好的性能，我们可以自然地提出以下几种不同的分配方法：
1.由二维线程块构成的二维网格（就是矩阵里面套矩阵）
2.由一维线程块构成的一维网格（就是长条套长条）
3.由一维线程块构成的二维网格（就是矩阵里面套长条）

那么首先，我们要了解如何在各种奇奇怪怪的布局下得到线程的索引，该索引对于任意一个global memory分割方式来说，都是确定的。所以我更愿意将其称之为线程的坐标，但为了统一性，下文还是称作为索引。

对于一个矩阵来说，假定分割方式是二维线程块组成的二维网格，我们把线程当作矩阵中的一个元素。设某个线程在矩阵的坐标映射为(ix,iy)那么我们有

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

有点抽象，我们来具体地分析一下
这个式子其实写成
ix=blockIdx.x*blockDim.x+threadIdx.x
更好

首先，我们从大的单位：块，找起。线程位于第n个块的第i个线程，每个块有k个线程，那么它的坐标应该表示为（从1开始计数）：
x=n*k+i
这个还是很抽象，我们再具体一点

我们假设一个学校的一个年级有18个班，班号从1排到18，又假设每个班都有50人，学号从1排到50。现在要给全年级的学生分配一个唯一的年级识别号，那么16班32号的年级识别号就为：
16*50+32=832
也就是它在年级中的索引为832。
这样就好理解了吧？
y轴方向是一样的情况，就不再赘述了。

我们再往上看，刚才分析的是矩阵，现在我要将这个矩阵拍扁（就像将数字图像转化为一个unsigned char *一样），那么就把矩阵的每一行拿出来（从第2行开始），将它的头和上一行的尾相接，形成一个新的行。
又抽象了，还是继续拿学校做例子。

我们假设这个学校有三个年级，每个年级的分布情况都和上述一样，那么我现在要给全校的学生分配一个唯一的校识别号。那么假设有一名学生是二年级的16班的32号，那么它的校识别号应该为：
2*16*50+32=1632

所以索引的概念就被建立起来了。
有了这些储备知识之后，我们继续看矩阵中的线程到global memory的唯一映射的公式：

idx=ix*ny+iy;

（或者是idx=iy*nx+ix;）
夭寿了哟，这个又是什么鬼啊？
莫急，对比一下刚刚我们所说的概念，或者直接就把公式带进去，可以得到：
idx=(threadIdx.x+blockIdx.x*blockDim.x)*nx+threadIdx.y+blockIdx.y*blockDim.y
那么nx其实就是在x方向上的线程的最多个数。

首先我们要注意一件事，所有的坐标都是以0为起始的，这一点与刚才我们所说的学校的例子不相同。x轴的方向水平向右，y轴的方向竖直向下。

我们假设一个分割方式为2*3的网格，每个矩阵的分割方式为4*1，如下图：

那么我们可以得到21这个点的ix为5，iy为2，也就是(5,2)这个点，而nx=8，所以
idx=2*8+5=21
这就可以很直观地理解了吧？
那么怎么生成这样的分割方式呢？

int nx=8;
int ny=4;
dim3 block(4,2); //生成x方向上有4个线程，y方向上有2个的块
dim3 grid((nx+block.x-1)/block.x,(ny+block.y-1)/block.y); 

最后一句的意思是这样的：
生成x方向上可以装下8个线程的最小的x方向上的块数，每块装4个，所以x方向上就是2块。y方向上可以装下4个线程的最小y方向上的块数，每块2个，所以y方向上就是2块。所以生成了一个2*2的网格。

2.二维网格和二维块矩阵求和

我们可以尝试用二维网格和二维块来对矩阵求和（也就是矩阵+矩阵）。也就是说，每个二维平面上的线程点处理对应坐标的矩阵元素。
Example 2-7

#include <cuda_runtime.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>
#include <device_launch_parameters.h>
#define abs(a) a>0 ? a : -a
#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)); \
        exit(-10*error); \
	} \
}
#define calCoor(a) threadIdx.##a+blockIdx.##a*blockDim.##a
double cpuSecond()
{
	return clock();
}
void sumMatrixOnHost(float *A, float *B, float *C, const int nx, const int ny) {
	float *ia = A;
	float *ib = B;
	float *ic = C;

	for (int iy = 0; iy < ny; iy++) {
		for (int ix = 0; ix < nx; ix++) {
			ic[ix] = ia[ix] + ib[ix];
		}
		ia += nx; ib += nx; ic += nx;
	}
}
void initialData(float *ip, int size) {	
	//generate different seed for random number	
	time_t t;
	srand((unsigned)time(&t));
	for (int i = 0; i < size; i++) {
		ip[i] = (float)(rand() & 0xFF) / 10.0f;
	}
}

__global__ void sumMatrixOnGPU2D(float *MatA, float *MatB, float *MatC, int nx, int ny) {
	unsigned int ix = calCoor(x);
	unsigned int iy = calCoor(y);
	unsigned int idx = iy * nx + ix;

	if (ix < nx&&iy < ny)
		MatC[idx] = MatA[idx] + MatB[idx];
}
void checkResult(float *hostRef, float *gpuRef, const int N) {
	double epslion = 1.0E-8;
	bool match = 1;
	for (int i = 0; i < N; i++) {
		if (abs(hostRef[i] - gpuRef[i]) > epslion) {
			match = 0;
			printf("Arrays do not match!\n");
			printf("host %5.2f gpu %5.2f at current %d\n", hostRef[i], gpuRef[i], i);
			break;
		}
	}
	if (match) printf("Arrays match.\n\n");
}
int main(int argc, char **argv) {
	printf("%s Starting...\n", argv[0]);

	//get device information
	int dev = 0;
	cudaDeviceProp deviceProp;
	CHECK(cudaGetDeviceProperties(&deviceProp, 0));
	printf("Using Device %d: %s\n", dev, deviceProp.name);
	CHECK(cudaSetDevice(dev));

	//set up data size of matrix
	int nx = 1 << 14;
	int ny = 1 << 14;

	int nxy = nx * ny;
	int nBytes = nxy * sizeof(float);
	printf("Matrix size: nx %d ny %d\n", nx, ny);

	//malloc host memory
	float *h_A, *h_B, *hostRef, *gpuRef;
	h_A = (float *)malloc(nBytes);
	h_B = (float *)malloc(nBytes);
	hostRef = (float *)malloc(nBytes);
	gpuRef = (float *)malloc(nBytes);

	double iStart, iElaps;
	//initialize data at host side	
	iStart = cpuSecond();
	initialData(h_A, nxy);
	initialData(h_B, nxy);
	iElaps = cpuSecond() - iStart;
	memset(hostRef, 0, nBytes);
	memset(gpuRef, 0, nBytes);
	//add matrix at host side for result checks	
	iStart = cpuSecond();
	sumMatrixOnHost(h_A, h_B, hostRef,nx,ny);
	iElaps = cpuSecond() - iStart;

	//malloc device global memory	
	float *d_MatA, *d_MatB, *d_MatC;
	cudaMalloc((void **)&d_MatA, nBytes);
	cudaMalloc((void **)&d_MatB, nBytes);
	cudaMalloc((void **)&d_MatC, nBytes);
	//transfer data from host to device	
	cudaMemcpy(d_MatA, h_A, nBytes, cudaMemcpyHostToDevice);
	cudaMemcpy(d_MatB, h_B, nBytes, cudaMemcpyHostToDevice);
	
	//invoke kernel at host side
	int dimx = 32;
	int dimy = 32;
	dim3 block(dimx, dimy);
	dim3 grid((nx + block.x - 1) / block.x, (ny + block.y - 1) / block.y);

	iStart = cpuSecond();
	sumMatrixOnGPU2D <<<grid, block >>> (d_MatA, d_MatB, d_MatC, nx,ny);
	cudaThreadSynchronize();
	iElaps = cpuSecond() - iStart;
	printf("sumMatrixOnGPU2D <<<(%d,%d),(%d,%d)>>> Time elapsed %f sec\n", grid.x,grid.y,block.x,block.y,iElaps);
	//copy kernel result back to host side	
	cudaMemcpy(gpuRef, d_MatC, nBytes, cudaMemcpyDeviceToHost);
	//check device results	
	checkResult(hostRef, gpuRef, nxy);
	//free device global memory	
	cudaFree(d_MatA);
	cudaFree(d_MatB);
	cudaFree(d_MatC);
	//free host memory	
	free(h_A);
	free(h_B);
	free(hostRef);
	free(gpuRef);

	//reset device
	cudaDeviceReset();
	return(0);
}

我们来分析一下这个代码，nx=1<<14，ny=1<<14，也就是2^13=16384。所以会生成一个x方向有16384个线程，y方向上有16384个线程的网格。接下来就是划分块数，我们可以看到，x方向上每块有32个线程，y方向上每块有32个线程，所以x方向上要划16384/32=512块，y方向上要划16384/32=512块。我们来通过程序的结果来验证我们的想法对不对：

可知我们的推测是正确的。

再来看核函数部分，首先通过宏Calcoor计算出x,y方向上的ix,iy，然后通过传进去的nx,ny的任意一个，求出当前调用核函数线程的idx。然后计算对应矩阵idx的元素。请记住，一定要加上边界判断条件。因为你不确定当前的线程坐标有没有越出你要计算的矩阵的坐标（就是慎防大网格算小矩阵越界的情况）。