GPU加速原理浅析及代码实现

一、CUDA简介
二、GPU架构特点
三、CUDA线程模型
四、CUDA内存模型
五、CUDA编程规范
- **第一个要掌握的编程要点**：**我们怎么写一个能在GPU跑的程序或函数呢？**
- **第二个要掌握的编程要点**：**CPU和GPU间的数据传输怎么写？**
- **第三个编程要点**：**grid、block、thread的关系及thread索引的计算**
六、CUDA计算加速测试
- **1. 首先，编写程序查看GPU的硬件配置信息**：
- **2. 第一个计算任务**：将两个元素数目为1024*1024的float数组相加。
- **3. 第二个计算任务**：矩阵加法。
- **4. 第三个计算任务**：矩阵乘法。
- **5. 测试程序汇总**

一、CUDA简介

CUDA(Compute Unified Device Architecture)是一种由NVIDIA推出的通用并行计算架构，该架构使GPU(Graphics Processing Unit)能够对复杂的计算问题做性能速度优化。

二、GPU架构特点

首先我们先来了解一下串行计算和并行计算。众所周知，高性能计算的关键是利用多核处理器进行并行计算。
当我们求解一个计算机程序任务时，我们很自然的想法就是将该任务分解成一系列小任务，把这些小任务一一完成。在串行计算时，我们的想法就是让我们的处理器按次序处理一个又一个的计算任务，直至完成最后一个小任务，也就意味着这个大的程序任务完成了。如下图所示，就是我们用串行编程思想求解问题的步骤。

图1 串行计算流程示意图

串行计算的缺点不言自明，如果我们拥有多核处理器，我们就可以利用多核处理器同时处理多个子任务(前提是这些子任务不需要相互依赖，比如线程1的计算任务不需要用到线程2的计算结果)。为了加快大任务的计算速度，我们可以把一些独立的模块分配到不同的处理器上进行同时计算(这就是并行)，最后再将这些结果进行整合，完成一次任务计算。下图就是将一个大的计算任务分解为多个小任务，然后将独立的小任务分配到不同处理器进行并行计算，最后再通过串行程序把结果汇总完成这次的总的计算任务。

图2 串并行混合计算流程示意图

所以，一个程序可不可以进行并行计算，关键就在于我们要分析出该程序可以拆分出哪几个执行模块，这些执行模块哪些是独立的，哪些又是强依赖强耦合的，独立的模块我们可以试着设计并行计算，充分利用多核处理器的优势进一步加速我们的计算任务，强耦合模块我们就使用串行编程，利用串行+并行的编程思路完成一次高性能计算。
接下来我们分析CPU和GPU的区别和各自的性能特点，我们在谈并行、串行计算时多次谈到“多核”的概念，现在我们先从“核”的角度开始这个话题。首先CPU是专为顺序串行处理而优化的几个核心组成。而GPU则由数以千计的更小、更高效的核心组成，这些核心专门为同时处理多任务而设计，可高效地处理并行任务。

CPU虽然每个核心自身能力极强，处理任务上非常强悍，无奈他核心少，在并行计算上表现不佳；反观GPU，虽然他的每个核心的计算能力不算强，但他胜在核心非常多，可以同时处理多个计算任务，在并行计算的支持上做得很好。

GPU和CPU的不同硬件特点决定了他们的应用场景，CPU是计算机的运算和控制的核心，GPU主要用作图形图像处理。图像在计算机呈现的形式就是矩阵，我们对图像的处理其实就是操作各种矩阵进行计算，而很多矩阵的运算其实可以做并行化，这使得图像处理可以做得很快，因此GPU在图形图像领域也有了大展拳脚的机会。

图3 CPU和GPU架构示意图

现在再从数据处理的角度来对比CPU和GPU的特点。CPU需要很强的通用性来处理各种不同的数据类型，比如整型、浮点数等，同时它又必须擅长处理逻辑判断所导致的大量分支跳转和中断处理，所以CPU其实就是一个能力很强的伙计，他能把很多事处理得妥妥当当，当然啦我们需要给他很多资源供他使用（各种硬件），这也导致了CPU不可能有太多核心（核心总数不超过16）。而GPU面对的则是类型高度统一的、相互无依赖的大规模数据和不需要被打断的纯净的计算环境，GPU有非常多核心（费米架构就有512核），虽然其核心的能力远没有CPU的核心强，但是胜在多，在处理简单计算任务时呈现出“人多力量大”的优势，这就是并行计算的魅力。

CPU：擅长流程控制和逻辑处理，不规则数据结构，不可预测存储结构，单线程程序，分支密集型算法。
GPU：擅长数据并行计算，规则数据结构，可预测存储模式。

现在的计算机体系架构中，要完成CUDA并行计算，单靠GPU一人之力是不能完成计算任务的，必须借助CPU来协同配合完成一次高性能的并行计算任务。

一般而言，并行部分在GPU上运行，串行部分在CPU运行，这就是异构计算。具体一点，异构计算的意思就是不同体系结构的处理器相互协作完成计算任务。CPU负责总体的程序流程，而GPU负责具体的计算任务，当GPU各个线程完成计算任务后，我们就将GPU那边计算得到的结果拷贝到CPU端，完成一次计算任务。

图4 How GPU acceleration works?

如上图所示，应用程序利用GPU实现加速的总体分工就是：密集计算代码（约占5%的代码量）由GPU负责完成，剩余串行代码由CPU负责执行。

三、CUDA线程模型

下面我们介绍CUDA的线程组织结构。首先我们都知道，线程是程序执行的最基本单元，CUDA的并行计算就是通过成千上万个线程的并行执行来实现的。下面的结构简图说明了GPU的不同层次的结构。

图5 CUDA的线程组织结构简图

CUDA的线程模型从小到大来总结为：

Thread：线程，并行的基本单位

Thread Block：线程块，互相合作的线程组，线程块有如下几个特点：
允许彼此同步

可以通过共享内存快速交换数据

以1维、2维或者三维组织

Grid：一组线程块
以1维、2维组织

共享全局内存

Kernel：在GPU上执行的核心程序，这个kernel函数是运行在某个Grid上的。
One kernel 对应 One Grid

每一个block和每一个thread都有自己的ID，我们可以通过相应的索引找到相应的线程块和线程。
threadIdx, blockIdx

Block ID：1D, 2D or 3D

Thread ID：1D, 2D or 3D

Kernel在device上执行时实际上是启动很多线程，一个kernel所启动的所有线程称为一个网格(grid)，同一个网格上的线程共享相同的全局内存空间。grid是线程结构的第一层次，而网格又可以分为很多线程块(block)，一个线程块里面包含很多线程，这是第二个层次。线程的两层组织结构如图5所示，这是一个grid和block均为2-dim的线程组织。grid和block都是定义为dim3类型的变量，dim3可以看成是包含三个无符号整数(x,y,z)成员的结构体变量，在定义时，缺省值初始化为1。因此，gird和block可以灵活地定义为1-dim，2-dim和3-dim结构，kernel调用时也必须通过执行配置<<<grid, block>>>来指定kernel所使用的网格维度和线程块维度。以上图为例，如何通过<<<grid, block>>>标记方式索引到我们想要的那个线程呢？CUDA的这种<<<grid, block>>>其实就是一个多级索引的方法，第一级索引是(grid.xIdx, grid.yIdy)，对应上图的例子就是(1, 1)，通过它我们就能找到这个线程块的位置，然后我们启动二级索引(block.xIdx, block.yIdy, block.zIdz)来定位到指定的线程。这就是CUDA的线程组织结构。

SP：最基本的处理单元，streaming processor，也成为CUDA core。最后具体的指令和任务都是在SP上处理的。GPU进行并行计算，也就是很多个SP同时做处理。
SM：多个SP加上其他的一些资源组成一个streaming multiprocessor，也叫GPU大核，其他资源如：warp scheduler，register，shared memory等。SM可以看作GPU的心脏(对比CPU核心)，register和shared memory是SM的稀缺资源。CUDA将这些资源分配给所有驻留在SM中的threads。因此，这些有限的资源就使每个SM中active warps有非常严格的限制，也就限制了并行能力。

简而言之，SP是线程执行的硬件单位，SM中包含多个SP，一个GPU可以有多个SM（比如16个），最终一个GPU可能包含有上千个SP。这么多核心“同时运行”，速度可想而知。这个“同时运行”只是想表明：实际上，软件逻辑上是所有SP是并行运行的，但是物理上并不是所有SP都能够同时执行计算（比如我们只有8个SM却有1024个线程块需要调度），因为有些会处于挂起，就绪等其他状态，这有关GPU的线程调度。
下面这个图将从硬件角度和软件角度解释CUDA的线程模型。

图6 CUDA线程模型解析图

每个线程由每个线程处理器(SP)执行

线程块由多核处理器(SM)执行

一个Kernel其实由一个grid来执行，一个kernel一次只能在一个GPU上执行

block是软件概念，一个block只会由一个SM调度，程序员在开发时，通过设定block的属性，告诉GPU硬件，我有多少个线程，线程该怎么组织。而具体怎么调度由SM的warp scheduler负责，block一旦被分配好SM，该block就会一直驻留在该SM中，直到执行结束。一个SM可以同时拥有多个blocks，但需要序列执行。下图显示了GPU内部的硬件架构：

图7 GPU内部硬件架构图

四、CUDA内存模型

CUDA中的内存模型分为以下几个层次：

每个线程都用自己的register（寄存器）
每个线程都有自己的local memory（局部内存）
每个线程块内都有自己的shared memory（共享内存），所有线程块内的所有线程共享这段内存资源
每个grid都有自己的global memory（全局内存），不同线程块的线程都可以使用
每个grid都有自己的constant memory（常量内存）和texture memory（纹理内存），不同线程块的线程都可以使用

线程访问这几类存储器的速度是 register>local memory>shared memory>global memory
下图表示上述内存在计算机架构中的所在层次。

图8 内存架构层次图

五、CUDA编程规范

那么，CUDA如何编程呢？？？
我们先捋一捋常见的CUDA术语：

第一个要掌握的编程要点：我们怎么写一个能在GPU跑的程序或函数呢？

通过关键字就可以表示某个程序在CPU上跑还是在GPU上跑！如下表所示，比如我们用__global__定义一个kernel函数，就是CPU上调用，GPU上执行，注意__global__函数的返回值必须设置为void。

第二个要掌握的编程要点：CPU和GPU间的数据传输怎么写？

首先介绍在GPU内存分配回收内存的函数接口：

cudaMalloc()：在设备端分配 global memory
cudaFree()：释放存储空间

CPU的数据和GPU端数据做数据传输的函数接口是一样的，他们通过传递的函数实参（枚举类型）来表示传输方向：
cudaMemcpy(void *dst, void *src, size_t nbytes, enum cudaMemcpyKind direction)
enum cudaMemcpyKind:

cudaMemcpyHostToDevice（CPU到GPU）
cudaMemcpyDeviceToHost（GPU到CPU）
cudaMemcpyDeviceToDevice（GPU到GPU）

第三个编程要点：grid、block、thread的关系及thread索引的计算

1. 首先，列出CUDA计算线程索引的一般公式：

// CUDA thread index:
int blockId = blockIdx.z*(gridDim.x*gridDim.y) + blockIdx.y*gridDim.x + blockIdx.x;
int threadId = blockId*(blockDim.x*blockDim.y*blockDim.z) + threadIdx.z*(blockDim.x*blockDim.y) + threadIdx.y*blockDim.x + threadIdx.x;

从上到下分别是：
block的3，2，1维；thread的3，2，1维。
如果没有对应维度，删除对应的计算部分即可。
维度的最小值是1，但是索引的值最小是0。
2. 接下来，我们来看一段简单代码，其功能是对两个数组求和，并保存到另一个数组。

#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <iostream>using namespace std;// 二：线程执行代码
__global__ void vector_add(float* vec1, float* vec2, float* vecres, int length)
{int tid = threadIdx.x;if (tid < length) {vecres[tid] = vec1[tid] + vec2[tid];}
}int main()
{const int length = 16; // 数组长度为16float a[length], b[length], c[length]; // host中的数组for (int i = 0; i < length; i++) // 初始赋值{ a[i] = b[i] = i;}float* a_device, *b_device, *c_device; // device中的数组cudaMalloc((void**)&a_device, length * sizeof(float)); // 分配内存cudaMalloc((void**)&b_device, length * sizeof(float));cudaMalloc((void**)&c_device, length * sizeof(float));cudaMemcpy(a_device, a, length * sizeof(float), cudaMemcpyHostToDevice); // 将host数组的值拷贝给device数组cudaMemcpy(b_device, b, length * sizeof(float), cudaMemcpyHostToDevice);// 一：参数配置dim3 grid(1, 1, 1), block(length, 1, 1); // 设置参数vector_add<<<grid,block>>>(a_device, b_device, c_device, length); // 启动kernelcudaMemcpy(c, c_device, length * sizeof(float), cudaMemcpyDeviceToHost); // 将结果拷贝到hostfor (int i = 0; i < length; i++) // 打印出来方便观察{                          cout << c[i] << " ";}cout << endl;// system("pause");return 0;
}

运行结果如下：

结果是正确的，也是我们所能预料到的。上面的代码计算的是两个一维向量的和，由于数组大小是16，所以我们使用了16个线程来计算。

dim3 grid(1, 1, 1), block(length, 1, 1); // 设置参数

图9 线程结构示意图

如上图所示，就是一个一维的block（此处只有x维度上存在16个线程）。所以，内建变量只有一个在起作用，就是threadIdx.x，它的范围是[0, 15]。因此，我们在计算线程索引时，只用内建变量就行了（其他的为0，写了也不起作用）。

// 二：线程执行代码
__global__ void vector_add(float* vec1, float* vec2, float* vecres, int length)
{int blockId = blockIdx.z*(gridDim.x*gridDim.y) + blockIdx.y*gridDim.x + blockIdx.x;int threadId = blockId*(blockDim.x*blockDim.y*blockDim.z) + threadIdx.z*(blockDim.x*blockDim.y) + threadIdx.y*blockDim.x + threadIdx.x;// 代入线程结构的12个参数，上式可以简化为如下形式// int threadId = threadIdx.x;if (threadId< length) {vecres[threadId] = vec1[threadId] + vec2[threadId];}
}

3.下面来看另一个例子：

dim3 grid(1, 1, 1), block(8, 2, 1); // 设置参数

图10 线程结构示意图

由上图可知，这个线程格只有一个一维的线程块，该线程块内的线程是二维的，x的维度为8，y的维度为2，共有8*2=16个线程，如果要用这16个线程来计算数组的累加，当然是可以的，但是这里需要改动一下线程执行代码中的索引计算方式了。

// 二：线程执行代码
__global__ void vector_add(float* vec1, float* vec2, float* vecres, int length)
{int blockId = blockIdx.z*(gridDim.x*gridDim.y) + blockIdx.y*gridDim.x + blockIdx.x;int threadId = blockId*(blockDim.x*blockDim.y*blockDim.z) + threadIdx.z*(blockDim.x*blockDim.y) + threadIdx.y*blockDim.x + threadIdx.x;// 代入线程结构的12个参数，上式可以简化为如下形式// int threadId = threadIdx.y * blockDim.x + threadIdx.x; if (threadId < length) {vecres[threadId] = vec1[threadId] + vec2[threadId];}
}

我们一定要有并行思想，这里有16个线程，kernel启动后，每个线程都有自己的索引号，比如某个线程位于grid中哪个维度的block（即blockIdx.x, blockIdx.y, blockIdx.z），又位于该block的哪个维度的线程（即threadIdx.x, threadIdx.y, threadIdx.z），利用这些线程索引号映射到对应的数组下标，我们要做的工作就是保证这些下标不重复（如果重复的话，那就惨了），最初的那种一维的计算方式就不行了。因此，通过使用threadIdx和blockDim来进行映射（偏移）。如上，blockDim.x = 8, blockDim.y = 2。
4. 下面再来看一个例子：

dim3 grid(1, 1, 1), block(4, 4, 1); // 设置参数

图11 线程结构示意图

由上图可知，这个线程格只有一个一维的线程块，该线程块内的线程是二维的，x的维度为4，y的维度为4，共有4*4=16个线程，如果要用这16个线程来计算数组的累加，当然也是可以的，其线程执行代码中的索引计算方式与上一个例子相同。

// 二：线程执行代码
__global__ void vector_add(float* vec1, float* vec2, float* vecres, int length)
{int blockId = blockIdx.z*(gridDim.x*gridDim.y) + blockIdx.y*gridDim.x + blockIdx.x;int threadId = blockId*(blockDim.x*blockDim.y*blockDim.z) + threadIdx.z*(blockDim.x*blockDim.y) + threadIdx.y*blockDim.x + threadIdx.x;// 代入线程结构的12个参数，上式可以简化为如下形式// int threadId = threadIdx.y * blockDim.x + threadIdx.x; if (threadId < length) {vecres[threadId] = vec1[threadId] + vec2[threadId];}
}

5. Another example：

dim3 grid(16, 1, 1), block(1, 1, 1); // 设置参数

先来个线程结构示意图如下，绿色的区域表示grid，蓝色的区域表示block，图中有一个grid和16个block，每个block都是一维，而且x维度上只有一个线程。

图12 线程结构示意图

显然，线程索引代码应该如下：

// 二：线程执行代码
__global__ void vector_add(float* vec1, float* vec2, float* vecres, int length)
{int blockId = blockIdx.z*(gridDim.x*gridDim.y) + blockIdx.y*gridDim.x + blockIdx.x;int threadId = blockId*(blockDim.x*blockDim.y*blockDim.z) + threadIdx.z*(blockDim.x*blockDim.y) + threadIdx.y*blockDim.x + threadIdx.x;// 代入线程结构的12个参数，上式可以简化为如下形式// int threadId = blockIdx.x; if (threadId < length) {vecres[threadId] = vec1[threadId] + vec2[threadId];}
}

6. Another example：

dim3 grid(4, 1, 1), block(4, 1, 1); // 设置参数

图13 线程结构示意图

此时，线程索引代码应该如下：

// 二：线程执行代码
__global__ void vector_add(float* vec1, float* vec2, float* vecres, int length)
{int blockId = blockIdx.z*(gridDim.x*gridDim.y) + blockIdx.y*gridDim.x + blockIdx.x;int threadId = blockId*(blockDim.x*blockDim.y*blockDim.z) + threadIdx.z*(blockDim.x*blockDim.y) + threadIdx.y*blockDim.x + threadIdx.x;// 代入线程结构的12个参数，上式可以简化为如下形式// int threadId = blockIdx.x * blockDim.x + threadIdx.x; if (threadId < length) {vecres[threadId] = vec1[threadId] + vec2[threadId];}
}

7. Another example：

dim3 grid(2, 2, 1), block(2, 2, 1); // 设置参数

图14 线程结构示意图

此时，线程索引代码应该如下：

// 二：线程执行代码
__global__ void vector_add(float* vec1, float* vec2, float* vecres, int length)
{int blockId = blockIdx.z*(gridDim.x*gridDim.y) + blockIdx.y*gridDim.x + blockIdx.x;int threadId = blockId*(blockDim.x*blockDim.y*blockDim.z) + threadIdx.z*(blockDim.x*blockDim.y) + threadIdx.y*blockDim.x + threadIdx.x;// 代入线程结构的12个参数，上式可以简化为如下形式// int threadId = (blockIdx.y * gridDim.x + blockIdx.x) * (blockDim.x * blockDim.y) + threadIdx.y * blockDim.x + threadIdx.x;if (threadId < length) {vecres[threadId] = vec1[threadId] + vec2[threadId];}
}

8. Another example：

dim3 grid(8, 4, 1), block(8, 2, 1); // 设置参数

图15 线程结构示意图

共有848*2 = 512个线程，此时，线程一维索引方式应该如下：

// 二：线程执行代码
__global__ void vector_add(float* vec1, float* vec2, float* vecres, int length)
{int blockId = blockIdx.z*(gridDim.x*gridDim.y) + blockIdx.y*gridDim.x + blockIdx.x;int threadId = blockId*(blockDim.x*blockDim.y*blockDim.z) + threadIdx.z*(blockDim.x*blockDim.y) + threadIdx.y*blockDim.x + threadIdx.x;// 代入线程结构的12个参数，上式可以简化为如下形式// int threadId = (blockIdx.y * gridDim.x + blockIdx.x) * (blockDim.x * blockDim.y) + threadIdx.y * blockDim.x + threadIdx.x;if (threadId < length) {vecres[threadId] = vec1[threadId] + vec2[threadId];}
}

其二维索引方式应该如下：

__global__ void vector_add(float** mat1, float** mat2, float** matres, int width)
{int x = blockIdx.x * blockDim.x + threadIdx.x;int y = blockIdx.y * blockDim.y + threadIdx.y;if (x < width && y < width) {matres[x][y] = mat1[x][y] + mat2[x][y];}
}

上述公式就是把线程和线程块的索引映射为图形像素坐标的计算方法。

六、CUDA计算加速测试

1. 首先，编写程序查看GPU的硬件配置信息：

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#include <stdio.h>
using namespace std;void INFO_GPU()
{int deviceCount;cudaGetDeviceCount(&deviceCount);for (int i = 0; i < deviceCount; i++){cudaDeviceProp devProp;cudaGetDeviceProperties(&devProp, i);cout << "使用GPU device " << i << ": " << devProp.name << endl;cout << "设备全局内存总量： " << devProp.totalGlobalMem / 1024 / 1024 << "MB" << endl;cout << "SM的数量：" << devProp.multiProcessorCount << endl;cout << "每个线程块的共享内存大小：" << devProp.sharedMemPerBlock / 1024.0 << " KB" << endl;cout << "每个线程块的最大线程数：" << devProp.maxThreadsPerBlock << endl;cout << "设备上一个线程块（Block）种可用的32位寄存器数量： " << devProp.regsPerBlock << endl;cout << "每个EM的最大线程数：" << devProp.maxThreadsPerMultiProcessor << endl;cout << "每个EM的最大线程束数：" << devProp.maxThreadsPerMultiProcessor / 32 << endl;cout << "设备上多处理器的数量： " << devProp.multiProcessorCount << endl;cout << "======================================================" << endl;}
}int main()
{INFO_GPU();system("pause");return 0;
}

输出结果如下：

使用GPU device 0: NVIDIA A100-PCIE-40GB
设备全局内存总量： 40683MB
SM的数量：108
每个线程块的共享内存大小：48 KB
每个线程块的最大线程数：1024
设备上一个线程块（Block）种可用的32位寄存器数量： 65536
每个EM的最大线程数：2048
每个EM的最大线程束数：64
设备上多处理器的数量： 108
======================================================

2. 第一个计算任务：将两个元素数目为1024*1024的float数组相加。

首先，我们尝试使用CPU串行来实现此计算任务。

// add_cpu
void add_cpu()
{clock_t start, end;double duration;start = clock();float* A, * B, * C;int n = 1024 * 1024;int size = n * sizeof(float);A = (float*)malloc(size);B = (float*)malloc(size);C = (float*)malloc(size);for (int i = 0; i < n; i++){A[i] = 90.0;B[i] = 10.0;}for (int i = 0; i < n; i++){C[i] = A[i] + B[i];}float max_error = 0.0;for (int i = 0; i < n; i++){max_error += fabs(100.0 - C[i]);}cout << "max_error is " << max_error << endl;end = clock();duration = (double)(end - start) / CLOCKS_PER_SEC;cout << "total time is " << duration * 1000 << "ms" << endl;// 21 mscout << "======================================================" << endl;
}

输出结果如下：

max_error is 0
total time is 21ms
======================================================

max_error = 0意味着计算正确，程序总耗时 21 ms。
接下来我们使用GPU并行来实现此计算任务。

// add_gpu
__global__ void Plus(float A[], float B[], float C[], int n)
{// CUDA thread index:int blockId = blockIdx.z * (gridDim.x * gridDim.y) + blockIdx.y * gridDim.x + blockIdx.x;int threadId = blockId * (blockDim.x * blockDim.y * blockDim.z) + threadIdx.z * (blockDim.x * blockDim.y) + threadIdx.y * blockDim.x + threadIdx.x;//int threadId = blockDim.x * blockIdx.x + threadIdx.x;C[threadId] = A[threadId] + B[threadId];
}void add_gpu()
{clock_t start, end;double duration;start = clock();float* A, * Ad, * B, * Bd, * C, * Cd;int n = 1024 * 1024;int size = n * sizeof(float);// CPU端分配内存A = (float*)malloc(size);B = (float*)malloc(size);C = (float*)malloc(size);// 初始化数组for (int i = 0; i < n; i++){A[i] = 90.0;B[i] = 10.0;}// GPU端分配内存cudaMalloc((void**)&Ad, size);cudaMalloc((void**)&Bd, size);cudaMalloc((void**)&Cd, size);// CPU的数据拷贝到GPU端cudaMemcpy(Ad, A, size, cudaMemcpyHostToDevice);cudaMemcpy(Bd, B, size, cudaMemcpyHostToDevice);cudaMemcpy(Bd, B, size, cudaMemcpyHostToDevice);// 定义kernel执行配置，（1024*1024/1024）个block，每个block里面有1024个线程dim3 dimBlock(1024);dim3 dimGrid(n / 1024);// 执行kernelPlus << <dimGrid, dimBlock >> > (Ad, Bd, Cd, n);// 将在GPU端计算好的结果拷贝回CPU端cudaMemcpy(C, Cd, size, cudaMemcpyDeviceToHost);// 校验误差float max_error = 0.0;for (int i = 0; i < n; i++){max_error += fabs(100.0 - C[i]);}cout << "max error is " << max_error << endl;// 释放CPU端、GPU端的内存free(A);free(B);free(C);cudaFree(Ad);cudaFree(Bd);cudaFree(Cd);end = clock();duration = (double)(end - start) / CLOCKS_PER_SEC;cout << "total time is " << duration * 1000 << "ms" << endl;// 130mscout << "======================================================" << endl;
}

输出结果如下：

max error is 0
total time is 130ms
======================================================

max_error = 0意味着计算正确，程序总耗时 130 ms。由上面的例子看出，使用CUDA编程时我们看不到for循环了，因为CPU编程的循环已经被分散到各个thread上做了，所以我们也就看到不到for一类的语句。从结果上看，CPU的循环计算的速度比GPU的多线程计算快多了，原因就在于CUDA中有大量的内存拷贝操作（数据传输花费了大量时间，而计算时间却非常少），如果计算量比较小的话，CPU计算会更合适一些。
另外：

每个Block中的Thread数最大不超过1024；
为了充分利用SM，Block数尽可能多，>100。

3. 第二个计算任务：矩阵加法。

下面计算一个稍微复杂的例子，矩阵加法，即对两个矩阵对应坐标的元素相加后的结果存储在第三个的对应位置的元素上。值得注意的是，这个计算任务采用了二维数组的计算方式，注意一下二维数组在CUDA编程中的写法。

// add_cpu_array
void add_cpu_array()
{clock_t start, end;double duration;start = clock();int* A, ** A_ptr, * B, ** B_ptr, * C, ** C_ptr;int total_size = ROWS * COLS * sizeof(int);A = (int*)malloc(total_size);B = (int*)malloc(total_size);C = (int*)malloc(total_size);A_ptr = (int**)malloc(ROWS * sizeof(int*));B_ptr = (int**)malloc(ROWS * sizeof(int*));C_ptr = (int**)malloc(ROWS * sizeof(int*));//CPU一维数组初始化for (int i = 0; i < ROWS * COLS; i++){A[i] = 80;B[i] = 20;}for (int i = 0; i < ROWS; i++){A_ptr[i] = A + COLS * i;B_ptr[i] = B + COLS * i;C_ptr[i] = C + COLS * i;}for (int i = 0; i < ROWS; i++)for (int j = 0; j < COLS; j++){C_ptr[i][j] = A_ptr[i][j] + B_ptr[i][j];}//检查结果int max_error = 0;for (int i = 0; i < ROWS * COLS; i++){//cout << C[i] << endl;max_error += abs(100 - C[i]);}cout << "max_error is " << max_error << endl;end = clock();duration = (double)(end - start) / CLOCKS_PER_SEC;cout << "total time is " << duration * 1000 << "ms" << endl;// 20mscout << "======================================================" << endl;
}

输出结果如下：

max_error is 0
total time is 20ms
======================================================

max_error = 0意味着计算正确，程序总耗时 20 ms。
接下来我们使用GPU并行来实现此计算任务。

// add_gpu_array
__global__ void addKernel(int** C, int** A, int** B)
{int idx = threadIdx.x + blockDim.x * blockIdx.x;int idy = threadIdx.y + blockDim.y * blockIdx.y;if (idx < Col && idy < Row) {C[idy][idx] = A[idy][idx] + B[idy][idx];}
}void add_gpu_array()
{clock_t start, end;double duration;start = clock();int** A = (int**)malloc(sizeof(int*) * Row);int** B = (int**)malloc(sizeof(int*) * Row);int** C = (int**)malloc(sizeof(int*) * Row);int* dataA = (int*)malloc(sizeof(int) * Row * Col);int* dataB = (int*)malloc(sizeof(int) * Row * Col);int* dataC = (int*)malloc(sizeof(int) * Row * Col);int** d_A;int** d_B;int** d_C;int* d_dataA;int* d_dataB;int* d_dataC;//malloc device memorycudaMalloc((void**)&d_A, sizeof(int**) * Row);cudaMalloc((void**)&d_B, sizeof(int**) * Row);cudaMalloc((void**)&d_C, sizeof(int**) * Row);cudaMalloc((void**)&d_dataA, sizeof(int) * Row * Col);cudaMalloc((void**)&d_dataB, sizeof(int) * Row * Col);cudaMalloc((void**)&d_dataC, sizeof(int) * Row * Col);//set valuefor (int i = 0; i < Row * Col; i++) {dataA[i] = 90;dataB[i] = 10;}//将主机指针A指向设备数据位置，目的是让设备二级指针能够指向设备数据一级指针//A 和  dataA 都传到了设备上，但是二者还没有建立对应关系for (int i = 0; i < Row; i++) {A[i] = d_dataA + Col * i;B[i] = d_dataB + Col * i;C[i] = d_dataC + Col * i;}cudaMemcpy(d_A, A, sizeof(int*) * Row, cudaMemcpyHostToDevice);cudaMemcpy(d_B, B, sizeof(int*) * Row, cudaMemcpyHostToDevice);cudaMemcpy(d_C, C, sizeof(int*) * Row, cudaMemcpyHostToDevice);cudaMemcpy(d_dataA, dataA, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);cudaMemcpy(d_dataB, dataB, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);dim3 threadPerBlock(128, 8);// 不大于1024dim3 blockNumber((Col + threadPerBlock.x - 1) / threadPerBlock.x, (Row + threadPerBlock.y - 1) / threadPerBlock.y);printf("Block(%d,%d)   Grid(%d,%d).\n", threadPerBlock.x, threadPerBlock.y, blockNumber.x, blockNumber.y);addKernel << <blockNumber, threadPerBlock >> > (d_C, d_A, d_B);//拷贝计算数据-一级数据指针cudaMemcpy(dataC, d_dataC, sizeof(int) * Row * Col, cudaMemcpyDeviceToHost);int max_error = 0;for (int i = 0; i < Row * Col; i++){//printf("%d\n", dataC[i]);max_error += abs(100 - dataC[i]);}//释放内存free(A);free(B);free(C);free(dataA);free(dataB);free(dataC);cudaFree(d_A);cudaFree(d_B);cudaFree(d_C);cudaFree(d_dataA);cudaFree(d_dataB);cudaFree(d_dataC);printf("max_error is %d\n", max_error);end = clock();duration = (double)(end - start) / CLOCKS_PER_SEC;cout << "total time is " << duration * 1000 << "ms" << endl;// 24 mscout << "======================================================" << endl;
}

输出结果如下：

Block(128,8)   Grid(8,128).
max_error is 0
total time is 24ms
======================================================

max_error = 0意味着计算正确，程序总耗时 24 ms。从结果看出，CPU计算时间还是比GPU的计算时间短。这里需要指出的是，这种二维数组的程序写法的效率并不高（虽然比较符合我们的思维方式），因为我们做了两次访存操作。所以一般而言，做高性能计算一般不会采取这种编程方式。

4. 第三个计算任务：矩阵乘法。

回顾一下矩阵乘法：两矩阵相乘，左矩阵第一行乘以右矩阵第一列（分别相乘，第一个数乘第一个数），乘完之后相加，即为结果的第一行第一列的数，依次往下算，直到计算完所有矩阵元素。

// mul_cpu
void matrix_mul_cpu(float* M, float* N, float* P, int width)
{for (int i = 0; i < width; i++)for (int j = 0; j < width; j++){float sum = 0.0;for (int k = 0; k < width; k++){float a = M[i * width + k];float b = N[k * width + j];sum += a * b;}P[i * width + j] = sum;}
}void mul_cpu()
{clock_t start, end;double duration;start = clock();float* A, * B, * C;int total_size = ROWS * COLS * sizeof(float);A = (float*)malloc(total_size);B = (float*)malloc(total_size);C = (float*)malloc(total_size);//CPU一维数组初始化for (int i = 0; i < ROWS * COLS; i++){A[i] = 80.0;B[i] = 20.0;}matrix_mul_cpu(A, B, C, COLS);end = clock();duration = (double)(end - start) / CLOCKS_PER_SEC;cout << "total time is " << duration * 1000 << "ms" << endl;// 8703 mscout << "======================================================" << endl;
}

输出结果如下：

total time is 8703ms
======================================================

程序总耗时 8703 ms，接下来我们使用GPU并行来实现此计算任务。梳理一下CUDA求解矩阵乘法的思路：因为C=A×B，我们利用每个线程求解C矩阵每个(x, y)的元素，每个线程载入A的一行和B的一列，遍历各自行列元素，对A、B对应的元素做一次乘法和一次加法。

// mul_gpu
__global__ void matrix_mul_gpu(int* M, int* N, int* P, int width)
{int i = threadIdx.x + blockDim.x * blockIdx.x;int j = threadIdx.y + blockDim.y * blockIdx.y;int sum = 0;for (int k = 0; k < width; k++){int a = M[j * width + k];int b = N[k * width + i];sum += a * b;}P[j * width + i] = sum;
}void mul_gpu()
{clock_t start, end;double duration;start = clock();int* A = (int*)malloc(sizeof(int) * Row * Col);int* B = (int*)malloc(sizeof(int) * Row * Col);int* C = (int*)malloc(sizeof(int) * Row * Col);//malloc device memoryint* d_dataA, * d_dataB, * d_dataC;cudaMalloc((void**)&d_dataA, sizeof(int) * Row * Col);cudaMalloc((void**)&d_dataB, sizeof(int) * Row * Col);cudaMalloc((void**)&d_dataC, sizeof(int) * Row * Col);//set valuefor (int i = 0; i < Row * Col; i++) {A[i] = 90;B[i] = 10;}cudaMemcpy(d_dataA, A, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);cudaMemcpy(d_dataB, B, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);dim3 threadPerBlock(128, 8);// 不超过1024dim3 blockNumber((Col + threadPerBlock.x - 1) / threadPerBlock.x, (Row + threadPerBlock.y - 1) / threadPerBlock.y);printf("Block(%d,%d)   Grid(%d,%d).\n", threadPerBlock.x, threadPerBlock.y, blockNumber.x, blockNumber.y);matrix_mul_gpu << <blockNumber, threadPerBlock >> > (d_dataA, d_dataB, d_dataC, Col);//拷贝计算数据-一级数据指针cudaMemcpy(C, d_dataC, sizeof(int) * Row * Col, cudaMemcpyDeviceToHost);//释放内存free(A);free(B);free(C);cudaFree(d_dataA);cudaFree(d_dataB);cudaFree(d_dataC);end = clock();duration = (double)(end - start) / CLOCKS_PER_SEC;cout << "total time is " << duration * 1000 << "ms" << endl;// 26 mscout << "======================================================" << endl;
}

输出结果如下：

Block(128,8)   Grid(8,128).
total time is 26ms
======================================================

从这个矩阵乘法任务可以看出，我们通过GPU进行并行计算的方式仅花费了26 ms，但是CPU串行计算方式却花费了8.7秒，计算速度提升了三百多倍，可见并行计算的威力！

5. 测试程序汇总


#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
using namespace std;#define COLS 1024
#define ROWS 1024
#define Col  1024
#define Row  1024void INFO_GPU()
{int deviceCount;cudaGetDeviceCount(&deviceCount);for (int i = 0; i < deviceCount; i++){cudaDeviceProp devProp;cudaGetDeviceProperties(&devProp, i);cout << "使用GPU device " << i << ": " << devProp.name << endl;cout << "设备全局内存总量： " << devProp.totalGlobalMem / 1024 / 1024 << "MB" << endl;cout << "SM的数量：" << devProp.multiProcessorCount << endl;cout << "每个线程块的共享内存大小：" << devProp.sharedMemPerBlock / 1024.0 << " KB" << endl;cout << "每个线程块的最大线程数：" << devProp.maxThreadsPerBlock << endl;cout << "设备上一个线程块（Block）种可用的32位寄存器数量： " << devProp.regsPerBlock << endl;cout << "每个EM的最大线程数：" << devProp.maxThreadsPerMultiProcessor << endl;cout << "每个EM的最大线程束数：" << devProp.maxThreadsPerMultiProcessor / 32 << endl;cout << "设备上多处理器的数量： " << devProp.multiProcessorCount << endl;cout << "======================================================" << endl;/*使用GPU device 0: NVIDIA A100 - PCIE - 40GB设备全局内存总量： 40683MBSM的数量：108每个线程块的共享内存大小：48 KB每个线程块的最大线程数：1024设备上一个线程块（Block）种可用的32位寄存器数量： 65536每个EM的最大线程数：2048每个EM的最大线程束数：64设备上多处理器的数量： 108*/}
}void add_cpu()
{clock_t start, end;double duration;start = clock();float* A, * B, * C;int n = 1024 * 1024;int size = n * sizeof(float);A = (float*)malloc(size);B = (float*)malloc(size);C = (float*)malloc(size);for (int i = 0; i < n; i++){A[i] = 90.0;B[i] = 10.0;}for (int i = 0; i < n; i++){C[i] = A[i] + B[i];}float max_error = 0.0;for (int i = 0; i < n; i++){max_error += fabs(100.0 - C[i]);}cout << "max_error is " << max_error << endl;end = clock();duration = (double)(end - start) / CLOCKS_PER_SEC;cout << "total time is " << duration * 1000 << "ms" << endl;// 21 mscout << "======================================================" << endl;
}// add_gpu
__global__ void Plus(float A[], float B[], float C[], int n)
{// CUDA thread index:int blockId = blockIdx.z * (gridDim.x * gridDim.y) + blockIdx.y * gridDim.x + blockIdx.x;int threadId = blockId * (blockDim.x * blockDim.y * blockDim.z) + threadIdx.z * (blockDim.x * blockDim.y) + threadIdx.y * blockDim.x + threadIdx.x;//int threadId = blockDim.x * blockIdx.x + threadIdx.x;C[threadId] = A[threadId] + B[threadId];
}void add_gpu()
{clock_t start, end;double duration;start = clock();float* A, * Ad, * B, * Bd, * C, * Cd;int n = 1024 * 1024;int size = n * sizeof(float);// CPU端分配内存A = (float*)malloc(size);B = (float*)malloc(size);C = (float*)malloc(size);// 初始化数组for (int i = 0; i < n; i++){A[i] = 90.0;B[i] = 10.0;}// GPU端分配内存cudaMalloc((void**)&Ad, size);cudaMalloc((void**)&Bd, size);cudaMalloc((void**)&Cd, size);// CPU的数据拷贝到GPU端cudaMemcpy(Ad, A, size, cudaMemcpyHostToDevice);cudaMemcpy(Bd, B, size, cudaMemcpyHostToDevice);cudaMemcpy(Bd, B, size, cudaMemcpyHostToDevice);// 定义kernel执行配置，（1024*1024/1024）个block，每个block里面有1024个线程dim3 dimBlock(1024);dim3 dimGrid(n / 1024);// 执行kernelPlus << <dimGrid, dimBlock >> > (Ad, Bd, Cd, n);// 将在GPU端计算好的结果拷贝回CPU端cudaMemcpy(C, Cd, size, cudaMemcpyDeviceToHost);// 校验误差float max_error = 0.0;for (int i = 0; i < n; i++){max_error += fabs(100.0 - C[i]);}cout << "max error is " << max_error << endl;// 释放CPU端、GPU端的内存free(A);free(B);free(C);cudaFree(Ad);cudaFree(Bd);cudaFree(Cd);end = clock();duration = (double)(end - start) / CLOCKS_PER_SEC;cout << "total time is " << duration * 1000 << "ms" << endl;// 114mscout << "======================================================" << endl;
}void add_cpu_array()
{clock_t start, end;double duration;start = clock();int* A, ** A_ptr, * B, ** B_ptr, * C, ** C_ptr;int total_size = ROWS * COLS * sizeof(int);A = (int*)malloc(total_size);B = (int*)malloc(total_size);C = (int*)malloc(total_size);A_ptr = (int**)malloc(ROWS * sizeof(int*));B_ptr = (int**)malloc(ROWS * sizeof(int*));C_ptr = (int**)malloc(ROWS * sizeof(int*));//CPU一维数组初始化for (int i = 0; i < ROWS * COLS; i++){A[i] = 80;B[i] = 20;}for (int i = 0; i < ROWS; i++){A_ptr[i] = A + COLS * i;B_ptr[i] = B + COLS * i;C_ptr[i] = C + COLS * i;}for (int i = 0; i < ROWS; i++)for (int j = 0; j < COLS; j++){C_ptr[i][j] = A_ptr[i][j] + B_ptr[i][j];}//检查结果int max_error = 0;for (int i = 0; i < ROWS * COLS; i++){//cout << C[i] << endl;max_error += abs(100 - C[i]);}cout << "max_error is " << max_error << endl;end = clock();duration = (double)(end - start) / CLOCKS_PER_SEC;cout << "total time is " << duration * 1000 << "ms" << endl;// 15mscout << "======================================================" << endl;
}// add_gpu_array
__global__ void addKernel(int** C, int** A, int** B)
{int idx = threadIdx.x + blockDim.x * blockIdx.x;int idy = threadIdx.y + blockDim.y * blockIdx.y;if (idx < Col && idy < Row) {C[idy][idx] = A[idy][idx] + B[idy][idx];}
}void add_gpu_array()
{clock_t start, end;double duration;start = clock();int** A = (int**)malloc(sizeof(int*) * Row);int** B = (int**)malloc(sizeof(int*) * Row);int** C = (int**)malloc(sizeof(int*) * Row);int* dataA = (int*)malloc(sizeof(int) * Row * Col);int* dataB = (int*)malloc(sizeof(int) * Row * Col);int* dataC = (int*)malloc(sizeof(int) * Row * Col);int** d_A;int** d_B;int** d_C;int* d_dataA;int* d_dataB;int* d_dataC;//malloc device memorycudaMalloc((void**)&d_A, sizeof(int**) * Row);cudaMalloc((void**)&d_B, sizeof(int**) * Row);cudaMalloc((void**)&d_C, sizeof(int**) * Row);cudaMalloc((void**)&d_dataA, sizeof(int) * Row * Col);cudaMalloc((void**)&d_dataB, sizeof(int) * Row * Col);cudaMalloc((void**)&d_dataC, sizeof(int) * Row * Col);//set valuefor (int i = 0; i < Row * Col; i++) {dataA[i] = 90;dataB[i] = 10;}//将主机指针A指向设备数据位置，目的是让设备二级指针能够指向设备数据一级指针//A 和  dataA 都传到了设备上，但是二者还没有建立对应关系for (int i = 0; i < Row; i++) {A[i] = d_dataA + Col * i;B[i] = d_dataB + Col * i;C[i] = d_dataC + Col * i;}cudaMemcpy(d_A, A, sizeof(int*) * Row, cudaMemcpyHostToDevice);cudaMemcpy(d_B, B, sizeof(int*) * Row, cudaMemcpyHostToDevice);cudaMemcpy(d_C, C, sizeof(int*) * Row, cudaMemcpyHostToDevice);cudaMemcpy(d_dataA, dataA, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);cudaMemcpy(d_dataB, dataB, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);dim3 threadPerBlock(128, 8);// 不大于1024dim3 blockNumber((Col + threadPerBlock.x - 1) / threadPerBlock.x, (Row + threadPerBlock.y - 1) / threadPerBlock.y);printf("Block(%d,%d)   Grid(%d,%d).\n", threadPerBlock.x, threadPerBlock.y, blockNumber.x, blockNumber.y);addKernel << <blockNumber, threadPerBlock >> > (d_C, d_A, d_B);//拷贝计算数据-一级数据指针cudaMemcpy(dataC, d_dataC, sizeof(int) * Row * Col, cudaMemcpyDeviceToHost);int max_error = 0;for (int i = 0; i < Row * Col; i++){//printf("%d\n", dataC[i]);max_error += abs(100 - dataC[i]);}//释放内存free(A);free(B);free(C);free(dataA);free(dataB);free(dataC);cudaFree(d_A);cudaFree(d_B);cudaFree(d_C);cudaFree(d_dataA);cudaFree(d_dataB);cudaFree(d_dataC);printf("max_error is %d\n", max_error);end = clock();duration = (double)(end - start) / CLOCKS_PER_SEC;cout << "total time is " << duration * 1000 << "ms" << endl;// 20 mscout << "======================================================" << endl;
}// mul_cpu
void matrix_mul_cpu(float* M, float* N, float* P, int width)
{for (int i = 0; i < width; i++)for (int j = 0; j < width; j++){float sum = 0.0;for (int k = 0; k < width; k++){float a = M[i * width + k];float b = N[k * width + j];sum += a * b;}P[i * width + j] = sum;}
}void mul_cpu()
{clock_t start, end;double duration;start = clock();float* A, * B, * C;int total_size = ROWS * COLS * sizeof(float);A = (float*)malloc(total_size);B = (float*)malloc(total_size);C = (float*)malloc(total_size);//CPU一维数组初始化for (int i = 0; i < ROWS * COLS; i++){A[i] = 80.0;B[i] = 20.0;}matrix_mul_cpu(A, B, C, COLS);end = clock();duration = (double)(end - start) / CLOCKS_PER_SEC;cout << "total time is " << duration * 1000 << "ms" << endl;// 8864 mscout << "======================================================" << endl;
}// mul_gpu
__global__ void matrix_mul_gpu(int* M, int* N, int* P, int width)
{int i = threadIdx.x + blockDim.x * blockIdx.x;int j = threadIdx.y + blockDim.y * blockIdx.y;int sum = 0;for (int k = 0; k < width; k++){int a = M[j * width + k];int b = N[k * width + i];sum += a * b;}P[j * width + i] = sum;
}void mul_gpu()
{clock_t start, end;double duration;start = clock();int* A = (int*)malloc(sizeof(int) * Row * Col);int* B = (int*)malloc(sizeof(int) * Row * Col);int* C = (int*)malloc(sizeof(int) * Row * Col);//malloc device memoryint* d_dataA, * d_dataB, * d_dataC;cudaMalloc((void**)&d_dataA, sizeof(int) * Row * Col);cudaMalloc((void**)&d_dataB, sizeof(int) * Row * Col);cudaMalloc((void**)&d_dataC, sizeof(int) * Row * Col);//set valuefor (int i = 0; i < Row * Col; i++) {A[i] = 90;B[i] = 10;}cudaMemcpy(d_dataA, A, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);cudaMemcpy(d_dataB, B, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);dim3 threadPerBlock(128, 8);// 不超过1024dim3 blockNumber((Col + threadPerBlock.x - 1) / threadPerBlock.x, (Row + threadPerBlock.y - 1) / threadPerBlock.y);printf("Block(%d,%d)   Grid(%d,%d).\n", threadPerBlock.x, threadPerBlock.y, blockNumber.x, blockNumber.y);matrix_mul_gpu << <blockNumber, threadPerBlock >> > (d_dataA, d_dataB, d_dataC, Col);//拷贝计算数据-一级数据指针cudaMemcpy(C, d_dataC, sizeof(int) * Row * Col, cudaMemcpyDeviceToHost);//释放内存free(A);free(B);free(C);cudaFree(d_dataA);cudaFree(d_dataB);cudaFree(d_dataC);end = clock();duration = (double)(end - start) / CLOCKS_PER_SEC;cout << "total time is " << duration * 1000 << "ms" << endl;// 27 mscout << "======================================================" << endl;
}int main()
{INFO_GPU();add_cpu();add_gpu();add_cpu_array();add_gpu_array();mul_cpu();mul_gpu();system("pause");return 0;
}

输出结果如下：

使用GPU device 0: NVIDIA A100-PCIE-40GB
设备全局内存总量： 40683MB
SM的数量：108
每个线程块的共享内存大小：48 KB
每个线程块的最大线程数：1024
设备上一个线程块（Block）种可用的32位寄存器数量： 65536
每个EM的最大线程数：2048
每个EM的最大线程束数：64
设备上多处理器的数量： 108
======================================================
max_error is 0
total time is 21ms
======================================================
max error is 0
total time is 122ms
======================================================
max_error is 0
total time is 21ms
======================================================
Block(128,8)   Grid(8,128).
max_error is 0
total time is 23ms
======================================================
total time is 8849ms
======================================================
Block(128,8)   Grid(8,128).
total time is 27ms
======================================================