目录

1. 使用一个二维网格和二维块的矩阵加法

1.1 关键代码

1.2 完整代码

1.3 运行时间

2. 使用一维网格和一维块的矩阵加法

2.1 关键代码

2.2 完整代码

2.3 运行时间

3. 使用二维网格和一维块的矩阵矩阵加法

3.1 关键代码

3.2 完整代码

3.3 运行时间

---

1. 使用一个二维网格和二维块的矩阵加法

这里建立二维网格(2,3)+二维块(4,2)为例，使用其块和线程索引映射矩阵索引。

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

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

 比如要获取矩阵元素(col,row) = (2,4) ,其全局索引是34，映射到矩阵坐标上，

ix = 2 + 0*3=2； iy = 0 + 2*2=4. 然后再映射到全局内存idx = 4*8 + 2 = 34.

1.1 关键代码

(1) 先固定二维线程块尺寸，再利用矩阵尺寸结合被分块尺寸，推导出二维网格尺寸。

// config
int dimx = 32;
int dimy = 32;
dim3 block(dimx, dimy);  // 二维线程块(x,y)=(4,2)
dim3 grid((nx + block.x - 1) / block.x, (ny + block.y - 1) / block.y); // 二维网格(2,3)
// 直接nx/block.x = 8/4=2. (8+4-1)/4=2. (9+4-1)/4=3

(2) 前面建立好了二维网格和二维线程块，根据公式去求矩阵索引。

// 去掉了循环
// 利用公式，使用二维网格、二维线程块映射矩阵索引。
__global__ void sumMatrixOnDevice2D(float* d_a, float* d_b, float* d_c, const int nx, const int ny)
{
	// 二维网格和二维块，映射到矩阵坐标
	unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
	unsigned int iy = threadIdx.y + blockIdx.y * blockDim.y;
	// 由矩阵坐标, 映射到全局坐标(都是线性存储的)
	unsigned int idx = iy * nx + ix;  // 坐标(ix, iy),前面由iy行，每行有nx个元素
	// 相加
	if (ix < nx && iy < ny)  // 配置线程的可能过多，这里防止越界。
	{
		d_c[idx] = d_a[idx] + d_b[idx];
	}
}

1.2 完整代码

#include "cuda_runtime.h"
#include "device_launch_parameters.h"  // threadIdx

#include <stdio.h>    // io
#include <time.h>     // time_t clock_t
#include <stdlib.h>  // rand
#include <memory.h>  //memset


#define CHECK(call)                                   \
{                                                     \
    const cudaError_t error_code = call;              \
    if (error_code != cudaSuccess)                    \
    {                                                 \
        printf("CUDA Error:\n");                      \
        printf("    File:       %s\n", __FILE__);     \
        printf("    Line:       %d\n", __LINE__);     \
        printf("    Error code: %d\n", error_code);   \
        printf("    Error text: %s\n",                \
            cudaGetErrorString(error_code));          \
        exit(1);                                      \
    }                                                 \
}

/// <summary>
/// 矩阵相加，线性存储的二维矩阵
/// </summary>
/// <param name="h_a"></param>
/// <param name="h_b"></param>
/// <param name="h_c"></param>
/// <param name="nx"></param>
/// <param name="ny"></param>
void sumMatrixOnHost(float* h_a, float* h_b, float* h_c, const int nx, const int ny)
{
	float* ia = h_a;
	float* ib = h_b;
	float* ic = h_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;  // 移动到下一行，ia下一行的第一个索引变成了0.
	}
}

// 去掉循环
__global__ void sumMatrixOnDevice2D(float* d_a, float* d_b, float* d_c, const int nx, const int ny)
{
	// 二维网格和二维块，映射到矩阵坐标
	unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
	unsigned int iy = threadIdx.y + blockIdx.y * blockDim.y;
	// 由矩阵坐标, 映射到全局坐标(都是线性存储的)
	unsigned int idx = iy * nx + ix;  // 坐标(ix, iy),前面由iy行，每行有nx个元素
	// 相加
	if (ix < nx && iy < ny)  // 配置线程的可能过多，这里防止越界。
	{
		d_c[idx] = d_a[idx] + d_b[idx];
	}
}

void initialData(float* p, const int N)
{
	//generate different seed from random number
	time_t t;
	srand((unsigned int)time(&t));  // 生成种子

	for (int i = 0; i < N; i++)
	{
		p[i] = (float)(rand() & 0xFF) / 10.0f;  // 随机数
	}
}


void checkResult(float* hostRef, float* deviceRef, const int N)
{
	double eps = 1.0E-8;
	int match = 1;
	for (int i = 0; i < N; i++)
	{
		if (hostRef[i] - deviceRef[i] > eps)
		{
			match = 0;
			printf("\nArrays do not match\n");
			printf("host %5.2f gpu %5.2f at current %d\n", hostRef[i], deviceRef[i], i);
			break;
		}
	}
	if (match)
		printf("\nArrays match!\n");
}


int main(void)
{
	// get device info
	int device = 0;
	cudaDeviceProp deviceProp;
	CHECK(cudaGetDeviceProperties(&deviceProp, device));
	printf("Using device: %d %s", device, deviceProp.name);  // 卡号0的显卡名称。
	CHECK(cudaSetDevice(device));  // 设置显卡号

	// set matrix dimension. 2^14 = 16384行列数
	int nx = 1 << 13, ny = 1 << 13, nxy = nx * ny;
	int nBytes = nxy * sizeof(float);

	// 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);  // 设备端拷回的数据
	// init data
	initialData(h_a, nxy);
	initialData(h_b, nxy);
	memset(hostRef, 0, nBytes);
	memset(gpuRef, 0, nBytes);
	clock_t begin = clock();
	// add matrix on host side for result checks.
	sumMatrixOnHost(h_a, h_b, hostRef, nx, ny);
	printf("\ncpu: %f s\n", (double)(clock() - begin) / CLOCKS_PER_SEC);

	// malloc device memory
	float* d_mat_a, * d_mat_b, * d_mat_c;
	cudaMalloc((void**)&d_mat_a, nBytes);
	cudaMalloc((void**)&d_mat_b, nBytes);
	cudaMalloc((void**)&d_mat_c, nBytes);

	// transfer data from host to device
	cudaMemcpy(d_mat_a, h_a, nBytes, cudaMemcpyHostToDevice);
	cudaMemcpy(d_mat_b, h_b, nBytes, cudaMemcpyHostToDevice);

	// config
	int dimx = 32;
	int dimy = 32;
	dim3 block(dimx, dimy);  // 二维线程块(x,y)=(4,2)
	dim3 grid((nx + block.x - 1) / block.x, (ny + block.y - 1) / block.y); // 二维网格(2,3)
	// 直接nx/block.x = 8/4=2. (8+4-1)/4=2.

	// invoke kernel
	begin = clock();
	sumMatrixOnDevice2D << <grid, block >> > (d_mat_a, d_mat_b, d_mat_c, nx, ny);
	CHECK(cudaDeviceSynchronize());
	printf("\ngpu: %f s\n", (double)(clock() - begin) / CLOCKS_PER_SEC);

	// check kernel error
	CHECK(cudaGetLastError());

	// copy kernel result back to host side
	cudaMemcpy(gpuRef, d_mat_c, nBytes, cudaMemcpyDeviceToHost);

	// check result
	checkResult(hostRef, gpuRef, nxy);

	// free memory
	cudaFree(d_mat_a);
	cudaFree(d_mat_b);
	cudaFree(d_mat_c);
	free(h_a);
	free(h_b);
	free(hostRef);
	free(gpuRef);

	// reset device
	cudaDeviceReset();

	return 0;
}

1.3 运行时间

加法运行速度提高了8倍。数据量越大，提升越明显。

2. 使用一维网格和一维块的矩阵加法

 为了使用一维网格和一维块，需要写一个新的核函数，其中每个线程处理ny个数据元素。如上图，一维网格是水平方向的一维，一维块是水平方向的一维。

2.1 关键代码

(1) 建立垂直方向的一维块，再是水平方向一维网格

// config
int dimx = 32;
int dimy = 1;
dim3 block(dimx, dimy);  // 一维线程块(x,y)=(32,1)，其中每一个线程都处理ny个元素。
// 一维网格((nx+block.x-1)/block.x,1)
dim3 grid((nx + block.x - 1) / block.x, 1); 

 (2) 前面建立好了一维网格和一维线程块，根据公式去求矩阵索引。

 假设blockDim = (32,1). gridDim = (1024,1). 比如当前的threadIdx.x = 10, 由于前面有一个块，当前的位置映射矩阵索引

ix = threadIdx.x + blockIdx.x * blockDim.x = 10 + 1 * 32 = 42

由于一个线程处理ny个元素（所以这里有一个循环，ix不变，iy从第一行开始到ny行）。

__global__ void sumMatrixOnDevice1D(float* d_a, float* d_b, float* d_c, const int nx, const int ny)
{
	unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
	// 一个线程处理ny个数据
	if (ix < nx)
	{
		for (int iy = 0; iy < ny; iy++)  // 处理ix这一列元素。
		{
			int idx = iy * nx + ix;  // 第iy行前面有iy*nx个元素，再加上当前行第ix个
			d_c[idx] = d_a[idx] + d_b[idx];
		}
	}
}

2.2 完整代码

#include "cuda_runtime.h"
#include "device_launch_parameters.h"  // threadIdx

#include <stdio.h>    // io
#include <time.h>     // time_t clock_t
#include <stdlib.h>  // rand
#include <memory.h>  //memset


#define CHECK(call)                                   \
{                                                     \
    const cudaError_t error_code = call;              \
    if (error_code != cudaSuccess)                    \
    {                                                 \
        printf("CUDA Error:\n");                      \
        printf("    File:       %s\n", __FILE__);     \
        printf("    Line:       %d\n", __LINE__);     \
        printf("    Error code: %d\n", error_code);   \
        printf("    Error text: %s\n",                \
            cudaGetErrorString(error_code));          \
        exit(1);                                      \
    }                                                 \
}

/// <summary>
/// 矩阵相加，线性存储的二维矩阵
/// </summary>
/// <param name="h_a"></param>
/// <param name="h_b"></param>
/// <param name="h_c"></param>
/// <param name="nx"></param>
/// <param name="ny"></param>
void sumMatrixOnHost(float* h_a, float* h_b, float* h_c, const int nx, const int ny)
{
	float* ia = h_a;
	float* ib = h_b;
	float* ic = h_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;  // 移动到下一行，ia下一行的第一个索引变成了0.
	}
}

// 去掉了循环
__global__ void sumMatrixOnDevice2D(float* d_a, float* d_b, float* d_c, const int nx, const int ny)
{
	// 二维网格和二维块，映射到矩阵坐标
	unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
	unsigned int iy = threadIdx.y + blockIdx.y * blockDim.y;
	// 由矩阵坐标, 映射到全局坐标(都是线性存储的)
	unsigned int idx = iy * nx + ix;  // 坐标(ix, iy),前面由iy行，每行有nx个元素
	// 相加
	if (ix < nx && iy < ny)  // 配置线程的可能过多，这里防止越界。
	{
		d_c[idx] = d_a[idx] + d_b[idx];
	}
}

__global__ void sumMatrixOnDevice1D(float* d_a, float* d_b, float* d_c, const int nx, const int ny)
{
	unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
	// 一个线程处理ny个数据
	if (ix < nx)
	{
		for (int iy = 0; iy < ny; iy++)
		{
			int idx = iy * nx + ix;  // 第iy行前面有iy*nx个元素，再加上当前行第ix个
			d_c[idx] = d_a[idx] + d_b[idx];
		}
	}
}

void initialData(float* p, const int N)
{
	//generate different seed from random number
	time_t t;
	srand((unsigned int)time(&t));  // 生成种子

	for (int i = 0; i < N; i++)
	{
		p[i] = (float)(rand() & 0xFF) / 10.0f;  // 随机数
	}
}


void checkResult(float* hostRef, float* deviceRef, const int N)
{
	double eps = 1.0E-8;
	int match = 1;
	for (int i = 0; i < N; i++)
	{
		if (hostRef[i] - deviceRef[i] > eps)
		{
			match = 0;
			printf("\nArrays do not match\n");
			printf("host %5.2f gpu %5.2f at current %d\n", hostRef[i], deviceRef[i], i);
			break;
		}
	}
	if (match)
		printf("\nArrays match!\n");
}


int main(void)
{
	// get device info
	int device = 0;
	cudaDeviceProp deviceProp;
	CHECK(cudaGetDeviceProperties(&deviceProp, device));
	printf("Using device: %d %s", device, deviceProp.name);  // 卡号0的显卡名称。
	CHECK(cudaSetDevice(device));  // 设置显卡号

	// set matrix dimension. 2^14 = 16384行列数
	int nx = 1 << 13, ny = 1 << 13, nxy = nx * ny;
	int nBytes = nxy * sizeof(float);

	// 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);  // 设备端拷回的数据
	// init data
	initialData(h_a, nxy);
	initialData(h_b, nxy);
	memset(hostRef, 0, nBytes);
	memset(gpuRef, 0, nBytes);
	clock_t begin = clock();
	// add matrix on host side for result checks.
	sumMatrixOnHost(h_a, h_b, hostRef, nx, ny);
	printf("\ncpu: %f s\n", (double)(clock() - begin) / CLOCKS_PER_SEC);

	// malloc device memory
	float* d_mat_a, * d_mat_b, * d_mat_c;
	cudaMalloc((void**)&d_mat_a, nBytes);
	cudaMalloc((void**)&d_mat_b, nBytes);
	cudaMalloc((void**)&d_mat_c, nBytes);

	// transfer data from host to device
	cudaMemcpy(d_mat_a, h_a, nBytes, cudaMemcpyHostToDevice);
	cudaMemcpy(d_mat_b, h_b, nBytes, cudaMemcpyHostToDevice);

	// config
	int dimx = 32;
	int dimy = 1;
	dim3 block(dimx, dimy);  // 一维线程块(x,y)=(32,1)，其中每一个线程都处理ny个元素。
	dim3 grid((nx + block.x - 1) / block.x, 1); // 一维网格((nx+block.x-1)/block.x,1)

	// invoke kernel
	begin = clock();
	//sumMatrixOnDevice2D << <grid, block >> > (d_mat_a, d_mat_b, d_mat_c, nx, ny);
	sumMatrixOnDevice1D << <grid, block >> > (d_mat_a, d_mat_b, d_mat_c, nx, ny);
	CHECK(cudaDeviceSynchronize());
	printf("\ngpu: %f s\n", (double)(clock() - begin) / CLOCKS_PER_SEC);

	// check kernel error
	CHECK(cudaGetLastError());

	// copy kernel result back to host side
	cudaMemcpy(gpuRef, d_mat_c, nBytes, cudaMemcpyDeviceToHost);

	// check result
	checkResult(hostRef, gpuRef, nxy);

	// free memory
	cudaFree(d_mat_a);
	cudaFree(d_mat_b);
	cudaFree(d_mat_c);
	free(h_a);
	free(h_b);
	free(hostRef);
	free(gpuRef);

	// reset device
	cudaDeviceReset();

	return 0;
}

2.3 运行时间

3. 使用二维网格和一维块的矩阵矩阵加法

当使用一个包含一维块的二维网格时，每个线程都只关注一个数据元素并且网格的第二个维数等于ny。比如块的维度是(32,1)，网格的维度是((nx+32-1)/32, (ny+1-1)/1) = ((nx+32-1)/32, ny).  

这可以看作是含有一个二维块的二维网格的特殊情况，其中块的第二个维数是1. 

利用块和网格索引，映射矩阵索引，如上图，blockIdx.y等于矩阵索引iy： 

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

iy = threadIdx.y + blockIdx.y * blockDim.y = threadIdx.y + 0 * 1 = threadIdx.y 

再利用矩阵索引，映射全局内存索引:

idx = iy * nx + ix;   // 当前行iy，块外的前面有iy * nx个元素，块内索引是ix

3.1 关键代码

(1) 建立一维线程块和二维网格

// config
int dimx = 32;
int dimy = 1;
dim3 block(dimx, dimy);  // 一维线程块(x,y)=(32,1)，其中每一个线程都处理一个元素。
// ((nx+32-1)/32, (ny+1-1)/1) = ((nx+32-1)/32, ny)
dim3 grid((nx + block.x - 1) / block.x, ny);

（2）利用块和网格索引，映射矩阵索引，再映射全局内存索引

__global__ void sumMatrixOnDevice2DGrid1DBlock(float* d_a, float* d_b, float* d_c, const int nx, const int ny)
{
	unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
	unsigned int iy = blockIdx.y;  // threadIdx.y + blockIdx.y * blockDim.y = threadIdx.y + 0 * 1
	unsigned int idx = iy * nx + ix;
	if (ix < nx && iy < ny)
	{
		d_c[idx] = d_a[idx] + d_b[idx];
	}
}

3.2 完整代码

 

#include "cuda_runtime.h"
#include "device_launch_parameters.h"  // threadIdx

#include <stdio.h>    // io
#include <time.h>     // time_t clock_t
#include <stdlib.h>  // rand
#include <memory.h>  //memset


#define CHECK(call)                                   \
{                                                     \
    const cudaError_t error_code = call;              \
    if (error_code != cudaSuccess)                    \
    {                                                 \
        printf("CUDA Error:\n");                      \
        printf("    File:       %s\n", __FILE__);     \
        printf("    Line:       %d\n", __LINE__);     \
        printf("    Error code: %d\n", error_code);   \
        printf("    Error text: %s\n",                \
            cudaGetErrorString(error_code));          \
        exit(1);                                      \
    }                                                 \
}

/// <summary>
/// 矩阵相加，线性存储的二维矩阵
/// </summary>
/// <param name="h_a"></param>
/// <param name="h_b"></param>
/// <param name="h_c"></param>
/// <param name="nx"></param>
/// <param name="ny"></param>
void sumMatrixOnHost(float* h_a, float* h_b, float* h_c, const int nx, const int ny)
{
	float* ia = h_a;
	float* ib = h_b;
	float* ic = h_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;  // 移动到下一行，ia下一行的第一个索引变成了0.
	}
}

// 去掉了循环
__global__ void sumMatrixOnDevice2D(float* d_a, float* d_b, float* d_c, const int nx, const int ny)
{
	// 二维网格和二维块，映射到矩阵坐标
	unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
	unsigned int iy = threadIdx.y + blockIdx.y * blockDim.y;
	// 由矩阵坐标, 映射到全局坐标(都是线性存储的)
	unsigned int idx = iy * nx + ix;  // 坐标(ix, iy),前面由iy行，每行有nx个元素
	// 相加
	if (ix < nx && iy < ny)  // 配置线程的可能过多，这里防止越界。
	{
		d_c[idx] = d_a[idx] + d_b[idx];
	}
}

__global__ void sumMatrixOnDevice1D(float* d_a, float* d_b, float* d_c, const int nx, const int ny)
{
	unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
	// 一个线程处理ny个数据
	if (ix < nx)
	{
		for (int iy = 0; iy < ny; iy++)
		{
			int idx = iy * nx + ix;  // 第iy行前面有iy*nx个元素，再加上当前行第ix个
			d_c[idx] = d_a[idx] + d_b[idx];
		}
	}
}

__global__ void sumMatrixOnDevice2DGrid1DBlock(float* d_a, float* d_b, float* d_c, const int nx, const int ny)
{
	unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
	unsigned int iy = blockIdx.y;  // threadIdx.y + blockIdx.y * blockDim.y = threadIdx.y + 0 * 1
	unsigned int idx = iy * nx + ix;
	if (ix < nx && iy < ny)
	{
		d_c[idx] = d_a[idx] + d_b[idx];
	}
}

void initialData(float* p, const int N)
{
	//generate different seed from random number
	time_t t;
	srand((unsigned int)time(&t));  // 生成种子

	for (int i = 0; i < N; i++)
	{
		p[i] = (float)(rand() & 0xFF) / 10.0f;  // 随机数
	}
}


void checkResult(float* hostRef, float* deviceRef, const int N)
{
	double eps = 1.0E-8;
	int match = 1;
	for (int i = 0; i < N; i++)
	{
		if (hostRef[i] - deviceRef[i] > eps)
		{
			match = 0;
			printf("\nArrays do not match\n");
			printf("host %5.2f gpu %5.2f at current %d\n", hostRef[i], deviceRef[i], i);
			break;
		}
	}
	if (match)
		printf("\nArrays match!\n");
}


int main(void)
{
	// get device info
	int device = 0;
	cudaDeviceProp deviceProp;
	CHECK(cudaGetDeviceProperties(&deviceProp, device));
	printf("Using device: %d %s", device, deviceProp.name);  // 卡号0的显卡名称。
	CHECK(cudaSetDevice(device));  // 设置显卡号

	// set matrix dimension. 2^14 = 16384行列数
	int nx = 1 << 13, ny = 1 << 13, nxy = nx * ny;
	int nBytes = nxy * sizeof(float);

	// 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);  // 设备端拷回的数据
	// init data
	initialData(h_a, nxy);
	initialData(h_b, nxy);
	memset(hostRef, 0, nBytes);
	memset(gpuRef, 0, nBytes);
	clock_t begin = clock();
	// add matrix on host side for result checks.
	sumMatrixOnHost(h_a, h_b, hostRef, nx, ny);
	printf("\ncpu: %f s\n", (double)(clock() - begin) / CLOCKS_PER_SEC);

	// malloc device memory
	float* d_mat_a, * d_mat_b, * d_mat_c;
	cudaMalloc((void**)&d_mat_a, nBytes);
	cudaMalloc((void**)&d_mat_b, nBytes);
	cudaMalloc((void**)&d_mat_c, nBytes);

	// transfer data from host to device
	cudaMemcpy(d_mat_a, h_a, nBytes, cudaMemcpyHostToDevice);
	cudaMemcpy(d_mat_b, h_b, nBytes, cudaMemcpyHostToDevice);

	// config
	int dimx = 32;
	int dimy = 1;
	dim3 block(dimx, dimy);  // 一维线程块(x,y)=(32,1)，其中每一个线程都处理一个元素。
	// ((nx+32-1)/32, (ny+1-1)/1) = ((nx+32-1)/32, ny)
	dim3 grid((nx + block.x - 1) / block.x, ny);

	// invoke kernel
	begin = clock();
	//sumMatrixOnDevice2D << <grid, block >> > (d_mat_a, d_mat_b, d_mat_c, nx, ny);  
	//sumMatrixOnDevice1D << <grid, block >> > (d_mat_a, d_mat_b, d_mat_c, nx, ny);
	sumMatrixOnDevice2DGrid1DBlock << <grid, block >> > (d_mat_a, d_mat_b, d_mat_c, nx, ny);
	CHECK(cudaDeviceSynchronize());
	printf("\ngpu: %f s\n", (double)(clock() - begin) / CLOCKS_PER_SEC);

	// check kernel error
	CHECK(cudaGetLastError());

	// copy kernel result back to host side
	cudaMemcpy(gpuRef, d_mat_c, nBytes, cudaMemcpyDeviceToHost);

	// check result
	checkResult(hostRef, gpuRef, nxy);

	// free memory
	cudaFree(d_mat_a);
	cudaFree(d_mat_b);
	cudaFree(d_mat_c);
	free(h_a);
	free(h_b);
	free(hostRef);
	free(gpuRef);

	// reset device
	cudaDeviceReset();

	return 0;
}

3.3 运行时间