CUDA笔记
2006年，NVIDIA公司发布了CUDA，CUDA是建立在NVIDIA的GPUs上的一个通用并行计算平台和编程模型，基于CUDA编程可以利用GPUs的并行计算引擎来更加高效地解决比较复杂的计算难题。近年来，GPU最成功的一个应用就是深度学习领域，基于GPU的并行计算已经成为训练深度学习模型的标配。此外CUDA也支持其它编程语言，如C/C++，Python，Fortran等语言。
GPU架构特点
串行计算和并行计算
高性能计算的关键是利用多核处理器进行并行计算，当我们求解一个计算机程序任务时，很自然的想法就是将任务分解成一系列小任务然后一一完成。
在串行计算中，就会让处理器每次只处理一个计算任务，之后再处理下一个，直到结束
但是串行计算缺点明显：如果我们有多核处理器，就可以利用多核处理器同时处理多个任务，而且小任务间并没有关联，那此时串行编程就失去意义了。
为了进一步加快大任务的计算速度，我们可以把独立的模块分配到不同的处理器上进行同时计算（并行），最后再将这些结果整合，完成一次计算任务
所以一个程序可不可以进行并行计算，关键就在于要分析出程序所含执行模块，其中那些是独立的，哪些是强依赖强耦合的。独立的模块使用并行计算加快速度，强依赖强耦合的模块就使用串行编程，即串行+并行完成一次高性能计算。
CPU&GPU
这里我们谈谈CPU和GPU之间的区别和各自的特点，我们在谈论并行和串行计算时多次谈到多核这个概念。所以我们从“核”这个角度来切入，首先CPU是由专为顺序串行处理而优化的几个核心组成的，GPU则是由数以千计的更小、更高效的核心组成，这些核心专为同时处理多任务设计。
也就是，CPU虽然每个核心自身能力极强，处理任务上非常强悍，无奈他核心少，在并行计算上表现不佳；反观GPU，虽然他的每个核心的计算能力不算强，但他胜在核心非常多，可以同时处理多个计算任务，在并行计算的支持上做得很好。
GPU和CPU的不同硬件特点决定了他们的应用场景，CPU是计算机的运算和控制的核心，GPU主要用作图形图像处理。图像在计算机呈现的形式就是矩阵，我们对图像的处理其实就是操作各种矩阵进行计算，而很多矩阵的运算其实可以做并行化，这使得图像处理可以做得很快，因此GPU在图形图像领域也有了大展拳脚的机会。
下图表示的就是一个多GPU计算机硬件系统，可以看出，一个GPU内存就有很多个SP和各类内存，这些硬件都是GPU进行高效并行计算的基础。
我们再从数据角度切入，谈谈CPU和GPU各自的特点：
CPU需要很强的通用性来处理各种不同的数据类型，比如整型、浮点数等，同时它又必须擅长处理逻辑判断所导致的大量分支跳转和中断处理，所以CPU其实就是一个能力很强的伙计，他能把很多事处理得妥妥当当，当然啦我们需要给他很多资源供他使用（各种硬件），这也导致了CPU不可能有太多核心（核心总数不超过16）。
GPU面对的则是类型高度统一的、相互无依赖的大规模数据和不需要被打断的纯净的计算环境，GPU有非常多核心（费米架构就有512核），虽然其核心的能力远没有CPU的核心强，但是胜在多，在处理简单计算任务时呈现出“人多力量大”的优势，这就是并行计算的魅力。
整理一下两者特点就是：
- CPU：擅长流程控制和逻辑处理，不规则数据结构，不可预测存储结构，单线程程序，分支密集型算法
- GPU：擅长数据并行计算，规则数据结构，可预测存储模式
在现在计算机体系架构中，要完成CUDA并行计算，单靠GPU是不可能的，GPU并不是一个独立运行的计算平台，而需要与CPU协同工作，可以看成是CPU的协处理器，因此当我们在说GPU并行计算时，其实是指的基于CPU+GPU的异构计算架构。在异构计算架构中，GPU与CPU通过PCIe总线连接在一起来协同工作，CPU所在位置称为为主机端（host），而GPU所在位置称为设备端（device），如下图所示。
可以看到GPU包括更多的运算核心，其特别适合数据并行的计算密集型任务，如大型矩阵运算，而CPU的运算核心较少，但是其可以实现复杂的逻辑运算，因此其适合控制密集型任务。
另外，CPU上的线程是重量级的，上下文切换开销大，但是GPU由于存在很多核心，其线程是轻量级的。因此，基于CPU+GPU的异构计算平台可以优势互补，CPU负责处理逻辑复杂的串行程序，而GPU重点处理数据密集型的并行计算程序，从而发挥最大功效。
具体一点，异构计算的意思就是不同体系结构的处理器相互协作完成计算任务。CPU负责总体的程序流程，而GPU负责具体的计算任务，当GPU各个线程完成计算任务后，将GPU那边计算得到的结果拷贝到CPU端，完成一次计算任务。
应用程序利用GPU实现加速的总体分工就是：密集计算代码（约占5%的代码量）由GPU负责完成，剩余串行代码由CPU负责执行
CUDA编程基础
CUDA线程模型
要深刻理解kernel，必须要对kernel的线程层次结构有一个清晰的认识。首先GPU上很多并行化的轻量级线程。kernel在device上执行时实际上是启动很多线程，一个kernel所启动的所有线程称为一个网格（grid），同一个网格上的线程共享相同的全局内存空间，grid是线程结构的第一层次，而网格又可以分为很多线程块（block），一个线程块里面包含很多线程，这是第二个层次。
即(Grid>Block>Thread)
线程两层组织结构如下图所示，这是一个gird和block均为2-dim的线程组织
grid和block都是定义为dim3类型的变量，dim3可以看成是包含三个无符号整数（x，y，z）成员的结构体变量，在定义时，缺省值初始化为1。因此grid和block可以灵活地定义为1-dim，2-dim以及3-dim结构，对于图中结构（主要水平方向为x轴），定义的grid和block如下所示，kernel在调用时也必须通过执行配置<<<grid, block>>>来指定kernel所使用的线程数及结构。
dim3 grid(3, 2);
dim3 block(5, 3);
kernel_fun<<< grid, block >>>(prams...);
<<<grid,block>>>本质就是一个多级索引的方法：
- 第一级索引是(grid.xIdx, grid.yIdy)，对应上图例子就是(1, 1)，通过它我们就能找到了这个线程块的位置
- 启动二级索引(block.xIdx, block.yIdx, block.zIdx)来定位到指定的线程
所以，一个线程需要两个内置的坐标变量（blockIdx，threadIdx）来唯一标识，它们都是dim3类型变量，其中blockIdx指明线程所在grid中的位置，而threaIdx指明线程所在block中的位置，如图中的Thread (1,1)满足：
threadIdx.x = 1
threadIdx.y = 1
blockIdx.x = 1
blockIdx.y = 1
SP&SM
- 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有非常严格的限制，也就限制了并行能力。
上面说到了kernel的线程组织层次，那么一个kernel实际上会启动很多线程，这些线程是逻辑上并行的，但是在物理层却并不一定。这其实和CPU的多线程有类似之处，多线程如果没有多核支持，在物理层也是无法实现并行的。但是好在GPU存在很多CUDA核心，充分利用CUDA核心可以充分发挥GPU的并行计算能力。GPU硬件的一个核心组件是SM，前面已经说过，SM是英文名是 Streaming Multiprocessor，翻译过来就是流式多处理器。SM的核心组件包括CUDA核心，共享内存，寄存器等，SM可以并发地执行数百个线程，并发能力就取决于SM所拥有的资源数。
当一个kernel被执行时，它的gird中的线程块被分配到SM上，一个线程块只能在一个SM上被调度。SM一般可以调度多个线程块，这要看SM本身的能力。那么有可能一个kernel的各个线程块被分配多个SM，所以grid只是逻辑层，而SM才是执行的物理层。
SM采用的是SIMT (Single-Instruction, Multiple-Thread，单指令多线程)架构，基本的执行单元是线程束（warps)，线程束包含32个线程，这些线程同时执行相同的指令，但是每个线程都包含自己的指令地址计数器和寄存器状态，也有自己独立的执行路径。所以尽管线程束中的线程同时从同一程序地址执行，但是可能具有不同的行为，比如遇到了分支结构，一些线程可能进入这个分支，但是另外一些有可能不执行，它们只能死等，因为GPU规定线程束中所有线程在同一周期执行相同的指令，线程束分化会导致性能下降。当线程块被划分到某个SM上时，它将进一步划分为多个线程束，因为这才是SM的基本执行单元，但是一个SM同时并发的线程束数是有限的。这是因为资源限制，SM要为每个线程块分配共享内存，而也要为每个线程束中的线程分配独立的寄存器。所以SM的配置会影响其所支持的线程块和线程束并发数量。
总之，就是grid和thread block只是逻辑划分，一个kernel的所有线程其实在物理层是不一定同时并发的。所以kernel的grid和block的配置不同，性能会出现差异，这点是要特别注意的。还有，由于SM的基本执行单元是包含32个线程的线程束，所以block大小一般要设置为32的倍数。
block组织结构
有时候，我们要知道一个线程在blcok中的全局ID，此时就必须还要知道block的组织结构，这是通过线程的内置变量blockDim来获得。它获取线程块各个维度的大小。
对于一个2-dim的block((D_x,D_y))，线程((x,y))的ID值为((x+y,*D_x))，如果是3-dim的block((D_x,D_y,D_z))，线程((x,y,z))的ID值为((x+y*D_x+z*D_x*D_y))。另外线程还有内置变量gridDim，用于获得网格块各个维度的大小。
block是软件概念，一个block只会由一个sm调度，程序员在开发时，通过设定block的属性，告诉GPU硬件，我有多少个线程，线程怎么组织。而具体怎么调度由sm的warps scheduler负责，block一旦被分配好SM，该block就会一直驻留在该SM中，直到执行结束。一个SM可以同时拥有多个blocks，但需要序列执行。下图显示了GPU内部的硬件架构：
kernel的这种线程组织结构天然适合vector,matrix等运算，如我们将利用上图2-dim结构实现两个矩阵的加法，每个线程负责处理每个位置的两个元素相加，代码如下所示。线程块大小为(16, 16)，然后将N*N大小的矩阵均分为不同的线程块来执行加法运算。
// Kernel定义
__global__ void MatAdd(float A[N][N], float B[N][N], float C[N][N])
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
if (i < N && j < N)
C[i][j] = A[i][j] + B[i][j];
}
int main()
{
...
// Kernel 线程配置
dim3 threadsPerBlock(16, 16);
dim3 numBlocks(N / threadsPerBlock.x+1, N / threadsPerBlock.y+1);
// kernel调用
MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
...
}
CUDA内存模型
- 每个线程都用自己的registers（寄存器）
- 每个线程都有自己的local memory（局部内存）
- 每个线程块内都有自己的shared memory（共享内存），所有线程块内的所有线程共享这段内存资源，且其生命周期与线程块一致
- 每个grid都有自己的global memory（全局内存），不同线程块的线程都可使用
- 每个grid都有自己的只读内存块：constant memory（常量内存）和texture memory（纹理内存），不同线程块的线程都可使用
线程访问这几类存储器的速度是
在进行CUDA编程前，可以先检查一下自己的GPU的硬件配置，这样才可以有的放矢，可以通过下面的程序获得GPU的配置属性：
vim test.cu
//test.cu
#include "device_launch_parameters.h"
#include <iostream>
int main()
{
int deviceCount;
cudaGetDeviceCount(&deviceCount);
for(int i=0;i<deviceCount;i++)
{
cudaDeviceProp devProp;
cudaGetDeviceProperties(&devProp, i);
std::cout << "使用GPU device " << i << ": " << devProp.name << std::endl;
std::cout << "设备全局内存总量： " << devProp.totalGlobalMem / 1024 / 1024 << "MB" << std::endl;
std::cout << "SM的数量：" << devProp.multiProcessorCount << std::endl;
std::cout << "每个线程块的共享内存大小：" << devProp.sharedMemPerBlock / 1024.0 << " KB" << std::endl;
std::cout << "每个线程块的最大线程数：" << devProp.maxThreadsPerBlock << std::endl;
std::cout << "设备上一个线程块（Block）种可用的32位寄存器数量： " << devProp.regsPerBlock << std::endl;
std::cout << "每个EM的最大线程数：" << devProp.maxThreadsPerMultiProcessor << std::endl;
std::cout << "每个EM的最大线程束数：" << devProp.maxThreadsPerMultiProcessor / 32 << std::endl;
std::cout << "设备上多处理器的数量： " << devProp.multiProcessorCount << std::endl;
std::cout << "======================================================" << std::endl;
}
return 0;
}
//使用nvcc编译
nvcc test.cu -o test
//运行
./test
// 输出如下
use GPU device0:NVIDIA GeForce RTX 2080 Ti
device global memory total: 11016MB
num of SM:68
maxThreadsPerBlock:1024
regsPerBlock:65536
maxThreadsPerMultiProcessor:1024
maxThreadsPerMultiProcessor/32:32
multiProcessorCount:68
======================================================
use GPU device1:NVIDIA GeForce RTX 2080 Ti
device global memory total: 11019MB
num of SM:68
maxThreadsPerBlock:1024
regsPerBlock:65536
maxThreadsPerMultiProcessor:1024
maxThreadsPerMultiProcessor/32:32
multiProcessorCount:68
======================================================
CUDA编程模型
如何在GPU上运行程序或者函数？
在CUDA中，host和device是两个重要的概念，我们用host指代CPU及其内存，而用device指代GPU及其内存。CUDA程序中既包含host程序，又包含device程序，它们分别在CPU和GPU上运行。同时，host与device之间可以进行通信，这样它们之间可以进行数据拷贝。典型的CUDA程序的执行流程如下：
- 分配host内存，并进行数据初始化；
- 分配device内存，并从host将数据拷贝到device上；
- 调用CUDA的核函数在device上完成指定的运算；
- 将device上的运算结果拷贝到host上；
- 释放device和host上分配的内存。
上面流程中最重要的一个过程是调用CUDA的核函数来执行并行计算，kernel是CUDA中一个重要的概念，kernel是在device上线程中并行执行的函数，核函数用__global__符号声明，在调用时需要用<<<grid, block>>>来指定kernel要执行的线程数量，在CUDA中，每一个线程都要执行核函数，并且每个线程会分配一个唯一的线程号thread ID，这个ID值可以通过核函数的内置变量threadIdx来获得。
由于GPU实际上是异构模型，所以需要区分host和device上的代码，在CUDA中是通过函数类型限定词开区别host和device上的函数，主要的三个函数类型限定词如下：
__global__：在device上执行，从host中调用（一些特定的GPU也可以从device上调用），返回类型必须是void，不支持可变参数参数，不能成为类成员函数。注意用__global__定义的kernel是异步的，这意味着host不会等待kernel执行完就执行下一步。__device__：在device上执行，单仅可以从device中调用，不可以和__global__同时用。__host__：在host上执行，仅可以从host上调用，一般省略不写，不可以和__global__同时用，但可和__device__，此时函数会在device和host都编译。
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）
怎么用代码表示线程组织模型？
可以用dim3类来表示网格和线程块的组织方式，网格grid可以表示为一维和二维格式，线程块block可以表示为一维、二维和三维的数据格式。
dim3 DimGrid(100, 50); //5000个线程块，维度是100*50
dim3 DimBlock(4, 8, 8); //每个线层块内包含256个线程，线程块内的维度是4*8*8
怎么计算线程ID？
1.使用N个线程块，每一个线程块只有一个线程
dim3 dimGrid(N);
dim3 dimBlock(1);
此时的线程号的计算方式就是
threadId = blockIdx.x;
其中threadId的取值范围为0到N-1。对于这种情况，我们可以将其看作是一个列向量，列向量中的每一行对应一个线程块。列向量中每一行只有1个元素，对应一个线程。
2.使用M×N个线程块，每个线程块1个线程
由于线程块是2维的，故可以看做是一个M*N的2维矩阵，其线程号有两个维度，即：
dim3 dimGrid(M,N);
dim3 dimBlock(1);
其中
blockIdx.x 取值0到M-1
blcokIdx.y 取值0到N-1
这种情况一般用于处理2维数据结构，比如2维图像。每一个像素用一个线程来处理，此时需要线程号来映射图像像素的对应位置，如
pos = blockIdx.y * blcokDim.x + blockIdx.x; //其中gridDim.x等于M
3.使用一个线程块，该线程具有N个线程
dim3 dimGrid(1);
dim3 dimBlock(N);
此时线程号的计算方式为
threadId = threadIdx.x;
其中threadId的范围是0到N-1，对于这种情况，可以看做是一个行向量，行向量中的每一个元素的每一个元素对应着一个线程。
4.使用M个线程块，每个线程块内含有N个线程
dim3 dimGrid(M);
dim3 dimBlock(N);
这种情况，可以把它想象成二维矩阵，矩阵的行与线程块对应，矩阵的列与线程编号对应，那线程号的计算方式为
threadId = threadIdx.x + blcokIdx*blockDim.x;
上面其实就是把二维的索引空间转换为一维索引空间的过程。
5.使用M×N的二维线程块，每一个线程块具有P×Q个线程
dim3 dimGrid(M, N);
dim3 dimBlock(P, Q);
这种情况其实是我们遇到的最多情况，特别适用于处理具有二维数据结构的算法，比如图像处理领域。
其索引有两个维度
threadId.x = blockIdx.x*blockDim.x+threadIdx.x;
threadId.y = blockIdx.y*blockDim.y+threadIdx.y;
上述公式就是把线程和线程块的索引映射为图像像素坐标的计算方法。
编程实例
向量相加
知道了CUDA编程基础，我们就来个简单的实战，利用CUDA编程实现两个1024*1024的float向量的加法
CPU实现
我们首先来看一下CPU如何串行实现这个任务
vim add_cpu.cu
#include <iostream>
#include <stdlib.h>
#include <sys/time.h>
#include <math.h>
using namespace std;
int main()
{
struct timeval start, end;
gettimeofday( &start, NULL );
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;
gettimeofday( &end, NULL );
int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
cout << "total time is " << timeuse/1000 << "ms" <<endl;
return 0;
}
nvcc add_cpu.cu -o add_cpu
./add_cpu
# output
max_error is 0
total time is 16ms
GPU实现
如果使用GPU采用并行计算执行此任务，又会有什么表现呢？
在实现之前，先简单介绍一下CUDA编程中内存管理API。首先是在device上分配内存的cudaMalloc函数：
cudaError_t cudaMalloc(void** devPtr, size_t size);
这个函数和C语言中的malloc类似，但是在device上申请一定字节大小的显存，其中devPtr是指向所分配内存的指针。同时要释放分配的内存使用cudaFree函数，这和C语言中的free函数对应。另外一个重要的函数是负责host和device之间数据通信的cudaMemcpy函数：
cudaError_t cudaMemcpy(void* dst, const void* src, size_t count, cudaMemcpyKind kind)
其中src指向数据源，而dst是目标区域，count是复制的字节数，其中kind控制复制的方向：cudaMemcpyHostToHost, cudaMemcpyHostToDevice, cudaMemcpyDeviceToHost及cudaMemcpyDeviceToDevice，如cudaMemcpyHostToDevice将host上数据拷贝到device上。
这里grid和block都设计为1-dim，首先定义kernel如下：
// 两个向量加法kernel，grid和block均为一维
__global__ void add(float A[], float B[], float C[], int n)
{
// 获取全局索引
int i = blockDim.x * blockIdx.x + threadIdx.x;
C[i] = A[i] + B[i];
}
其中stride是整个grid的线程数，有时候向量的元素数很多，这时候可以将在每个线程实现多个元素（元素总数/线程总数）的加法，相当于使用了多个grid来处理，这是一种grid-stride loop方式，不过下面的例子一个线程只处理一个元素，所以kernel里面的循环是不执行的。下面我们具体实现向量加法：
编程要点：
- 每个Block中的Thread数最大不超过512；
- 为了充分利用SM，Block数尽可能多，至少>100
vim add_gpu.cu
#include "cuda_runtime.h"
#include <stdlib.h>
#include <iostream>
#include <sys/time.h>
using namespace std;
__global__ void Plus(float A[], float B[], float C[], int n)
{
int i = blockDim.x * blockIdx.x + threadIdx.x;
C[i] = A[i] + B[i];
}
int main()
{
struct timeval start, end;
gettimeofday( &start, NULL );
int N = 1024 * 1024;
int nBytes = N * sizeof(float);
// 申请host内存
float *x, *y, *z;
x = (float*)malloc(nBytes);
y = (float*)malloc(nBytes);
z = (float*)malloc(nBytes);
// 初始化数据
for (int i = 0; i < N; ++i)
{
x[i] = 90.0;
y[i] = 10.0;
}
// 申请device内存
float *d_x, *d_y, *d_z;
cudaMalloc((void**)&d_x, nBytes);
cudaMalloc((void**)&d_y, nBytes);
cudaMalloc((void**)&d_z, nBytes);
// 将host数据拷贝到device
cudaMemcpy((void*)d_x, (void*)x, nBytes, cudaMemcpyHostToDevice);
cudaMemcpy((void*)d_y, (void*)y, nBytes, cudaMemcpyHostToDevice);
// 定义kernel的执行配置
dim3 blockSize(256);
dim3 gridSize((N + blockSize.x - 1) / blockSize.x);
// 执行kernel
Plus<< < gridSize, blockSize >> >(d_x, d_y, d_z, N);
// 将device得到的结果拷贝到host
cudaMemcpy((void*)z, (void*)d_z, nBytes, cudaMemcpyDeviceToHost);
// 检查执行结果
float maxError = 0.0;
for (int i = 0; i < N; i++)
maxError += fabs(100.0 - z[i]);
cout << "max error is: " << maxError << endl;
// 释放device内存
cudaFree(d_x);
cudaFree(d_y);
cudaFree(d_z);
// 释放host内存
free(x);
free(y);
free(z);
gettimeofday( &end, NULL );
int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
cout << "total time is " << timeuse/1000 << "ms" <<endl;
return 0;
}
nvcc add_gpu.cu -o add_gpu
./add_gpu
# output
max error is: 0
total time is 125ms
由上面的例子看出，使用CUDA编程时我们看不到for循环了，因为CPU编程的循环已经被分散到各个thread上做了，所以我们也就看到不到for一类的语句。从结果上看，CPU的循环计算的速度比GPU计算快多了，原因就在于CUDA中有大量的内存拷贝操作（数据传输花费了大量时间，而计算时间却非常少），如果计算量比较小的话，CPU计算会更合适一些。
这里我们的向量大小为1024，而block大小为256，那么grid大小是4096，kernel的线程层级结构如下图所示：
在上面的实现中，我们需要单独在host和device上进行内存分配，并且要进行数据拷贝，这是很容易出错的。好在CUDA 6.0引入统一内存（Unified Memory）来避免这种麻烦，简单来说就是统一内存使用一个托管内存来共同管理host和device中的内存，并且自动在host和device中进行数据传输。CUDA中使用cudaMallocManaged函数分配托管内存：
cudaError_t cudaMallocManaged(void **devPtr, size_t size, unsigned int flag=0);
利用统一内存，可以将上面的程序简化如下：
vim add_gpu_memory.cu
#include "cuda_runtime.h"
#include <stdlib.h>
#include <iostream>
#include <sys/time.h>
using namespace std;
__global__ void Plus(float A[], float B[], float C[], int n)
{
int i = blockDim.x * blockIdx.x + threadIdx.x;
C[i] = A[i] + B[i];
}
int main()
{
struct timeval start, end;
gettimeofday( &start, NULL );
int N = 1024 * 1024;
int nBytes = N * sizeof(float);
// 申请托管内存
float *x, *y, *z;
cudaMallocManaged((void**)&x, nBytes);
cudaMallocManaged((void**)&y, nBytes);
cudaMallocManaged((void**)&z, nBytes);
// 初始化数据
for (int i = 0; i < N; ++i)
{
x[i] = 90.0;
y[i] = 10.0;
}
// 定义kernel的执行配置
dim3 blockSize(256);
dim3 gridSize((N + blockSize.x - 1) / blockSize.x);
// 执行kernel
add << < gridSize, blockSize >> >(x, y, z, N);
// 同步device 保证结果能正确访问
cudaDeviceSynchronize();
// 检查执行结果
float maxError = 0.0;
for (int i = 0; i < N; i++)
maxError += fabs(100.0 - z[i]);
cout << "max error is: " << maxError << endl;
// 释放内存
cudaFree(x);
cudaFree(y);
cudaFree(z);
gettimeofday( &end, NULL );
int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
cout << "total time is " << timeuse/1000 << "ms" <<endl;
return 0;
}
nvcc add_gpu_memory.cu -o add_gpu_memory
./add_gpu_memory
# output
max error is: 0
total time is 145ms
相比之前的代码，使用统一内存更简洁了，值得注意的是kernel执行是与host异步的，由于托管内存自动进行数据传输，这里要用cudaDeviceSynchronize()函数保证device和host同步，这样后面才可以正确访问kernel计算的结果。
矩阵相加
即对两个矩阵对应坐标的元素相加后的结果存储在第三个的对应位置的元素上
CPU实现
vim madd_cpu.cu
#include <stdlib.h>
#include <iostream>
#include <sys/time.h>
#include <math.h>
#define ROWS 1024
#define COLS 1024
using namespace std;
int main()
{
struct timeval start, end;
gettimeofday( &start, NULL );
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;
gettimeofday( &end, NULL );
int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
cout << "total time is " << timeuse/1000 << "ms" <<endl;
return 0;
}
nvcc madd_cpu.cu -o madd_cpu
./madd_cpu
# output
max_error is 0
total time is 14ms
GUP实现
vim madd_gpu.cu
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <sys/time.h>
#include <stdio.h>
#include <math.h>
#define Row 1024
#define Col 1024
__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];
}
}
int main()
{
struct timeval start, end;
gettimeofday( &start, NULL );
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 memory
cudaMalloc((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 value
for (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(16, 16);
dim3 blockNumber( (Col + threadPerBlock.x - 1)/ threadPerBlock.x, (Row + threadPerBlock.y - 1) / threadPerBlock.y );
printf("Block(%d,%d) Grid(%d,%d).
", 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
", 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
", max_error);
gettimeofday( &end, NULL );
int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
printf("total time is %d ms
", timeuse/1000);
return 0;
}
nvcc madd_gpu.cu -o madd_gpu
./madd_gpu
# output
Block(16,16) Grid(64,64).
max_error is 0
total time is 190 ms
从结果看出，CPU计算时间还是比GPU的计算时间短。这里需要指出的是，这种二维数组的程序写法的效率并不高（虽然比较符合我们的思维方式），因为我们做了两次访存操作。所以一般而言，做高性能计算一般不会采取这种编程方式。
矩阵乘法
两矩阵相乘，左矩阵第一行乘以右矩阵第一列（分别相乘，第一个数乘第一个数），乘完之后相加，即为结果的第一行第一列的数，依次往下算，直到计算完所有矩阵元素
设输入矩阵为(A)和(B)，要得到(C=A imes B)。实现思路是每个线程计算(C)的一个元素值(C_{i,j})
CPU实现
vim mtime_cpu.cu
#include <iostream>
#include <stdlib.h>
#include <sys/time.h>
#define ROWS 1024
#define COLS 1024
using namespace std;
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;
}
}
int main()
{
struct timeval start, end;
gettimeofday( &start, NULL );
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);
gettimeofday( &end, NULL );
int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
cout << "total time is " << timeuse/1000 << "ms" <<endl;
return 0;
}
nvcc mtime_cpu.cu -o mtime_cpu
./mtime_cpu
# output
total time is 5005ms
GPU实现
对于矩阵运算，应该选用grid和block为2-D的。首先定义矩阵的结构体：
// 矩阵类型，行优先，M(row, col) = *(M.elements + row * M.width + col)
struct Matrix
{
int width;
int height;
float *elements;
};
然后实现矩阵乘法的核函数，这里我们定义了两个辅助的__device__函数分别用于获取矩阵的元素值和为矩阵元素赋值，具体代码如下：
// 获取矩阵A的(row, col)元素
__device__ float getElement(Matrix *A, int row, int col)
{
return A->elements[row * A->width + col];
}
// 为矩阵A的(row, col)元素赋值
__device__ void setElement(Matrix *A, int row, int col, float value)
{
A->elements[row * A->width + col] = value;
}
// 矩阵相乘kernel，2-D，每个线程计算一个元素
__global__ void matMulKernel(Matrix *A, Matrix *B, Matrix *C)
{
float Cvalue = 0.0;
int row = threadIdx.y + blockIdx.y * blockDim.y;
int col = threadIdx.x + blockIdx.x * blockDim.x;
for (int i = 0; i < A->width; ++i)
{
Cvalue += getElement(A, row, i) * getElement(B, i, col);
}
setElement(C, row, col, Cvalue);
}
最后我们采用统一内存编写矩阵相乘的测试实例：
vim mtime_gpu.cu
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <sys/time.h>
#include <iostream>
#include <math.h>
using namespace std;
// 矩阵类型，行优先，M(row, col) = *(M.elements + row * M.width + col)
struct Matrix
{
int width;
int height;
float *elements;
};
// 获取矩阵A的(row, col)元素
__device__ float getElement(Matrix *A, int row, int col)
{
return A->elements[row * A->width + col];
}
// 为矩阵A的(row, col)元素赋值
__device__ void setElement(Matrix *A, int row, int col, float value)
{
A->elements[row * A->width + col] = value;
}
// 矩阵相乘kernel，2-D，每个线程计算一个元素
__global__ void matMulKernel(Matrix *A, Matrix *B, Matrix *C)
{
float Cvalue = 0.0;
int row = threadIdx.y + blockIdx.y * blockDim.y;
int col = threadIdx.x + blockIdx.x * blockDim.x;
for (int i = 0; i < A->width; ++i)
{
Cvalue += getElement(A, row, i) * getElement(B, i, col);
}
setElement(C, row, col, Cvalue);
}
int main()
{
struct timeval start, end;
gettimeofday( &start, NULL );
int width = 1 << 10;
int height = 1 << 10;
Matrix *A, *B, *C;
// 申请托管内存
cudaMallocManaged((void**)&A, sizeof(Matrix));
cudaMallocManaged((void**)&B, sizeof(Matrix));
cudaMallocManaged((void**)&C, sizeof(Matrix));
int nBytes = width * height * sizeof(float);
cudaMallocManaged((void**)&A->elements, nBytes);
cudaMallocManaged((void**)&B->elements, nBytes);
cudaMallocManaged((void**)&C->elements, nBytes);
// 初始化数据
A->height = height;
A->width = width;
B->height = height;
B->width = width;
C->height = height;
C->width = width;
for (int i = 0; i < width * height; ++i)
{
A->elements[i] = 1.0;
B->elements[i] = 2.0;
}
// 定义kernel的执行配置
dim3 blockSize(32, 32);
dim3 gridSize((width + blockSize.x - 1) / blockSize.x,
(height + blockSize.y - 1) / blockSize.y);
// 执行kernel
matMulKernel << < gridSize, blockSize >> >(A, B, C);
// 同步device 保证结果能正确访问
cudaDeviceSynchronize();
// 检查执行结果
float maxError = 0.0;
for (int i = 0; i < width * height; ++i)
maxError = fmax(maxError, fabs(C->elements[i] - 2 * width));
cout << "最大误差: " << maxError << endl;
gettimeofday( &end, NULL );
int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
cout << "total time is " << timeuse/1000 << "ms" <<endl;
return 0;
}
nvcc mtime_gpu.cu -o mtime_gpu
./mtime_gpu
# output
最大误差: 0
total time is 152ms
从这个矩阵乘法任务可以看出，我们通过GPU进行并行计算的方式仅花费了152ms，但是CPU串行计算方式却花费了5005ms，计算速度提升了五十多倍，可见并行计算的威力！
另一个版本的GPU实现
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <sys/time.h>
#include <stdio.h>
#include <math.h>
#define Row 1024
#define Col 1024
__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;
}
int main()
{
struct timeval start, end;
gettimeofday( &start, NULL );
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 memory
int *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 value
for (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(16, 16);
dim3 blockNumber((Col+threadPerBlock.x-1)/ threadPerBlock.x, (Row+threadPerBlock.y-1)/ threadPerBlock.y );
printf("Block(%d,%d) Grid(%d,%d).
", 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);
gettimeofday( &end, NULL );
int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
printf("total time is %d ms
", timeuse/1000);
return 0;
}
TAHT‘S ALL.TANKS FOR READING.