CUDA笔记

Posted on | In | 本文总阅读量 次

CUDA笔记

CUDA函数前缀

​ CUDA使用cu作为文件类型后缀，在CUDA中，有三种常见的前缀__device__、__global__、__host__，分别代表CUDA的三种运行场景，如下表所示：

限定符 执行 调用 备注
global 设备端执行 可以从主机调用也可以从计算能力3以上的设备调用 必须有一个void的返回类型
device 设备端执行 设备端调用
host 主机端执行 主机调用 可以省略

注意，一个函数可以同时被多个前缀所修饰，如CUDA 10浮点数的转换：

1
__host__ __device__ __half__ __float2half( const float a ) throw()

上面这个函数可以在host端被调用，也可以在device端被调用，我们可以通过一个函数前缀判断这个函数的运行环境。

CUDA shared memory

cuda内存模型

​ CUDA内存模型如上，最底层的DRAM代表global memory。在global memory部分，数据对齐仍然是很重要的话题，当使用L1（L1和Shared Memory是共享的）的时候，对齐问题可以忽略，但是非连续内存仍然会降低性能。在某些情况下，非连续的访问不可避免，使用shared memory 可以提高系统的性能。

GPU上的内存有两种：

  • On-board memory(板载显存)：主要包括全局内存（global memory）、本地内存（local memory）、常量内存（constant memory）、纹理内存（texture memory）等
  • On-chip memory（片上内存）：主要包括寄存器（register）和共享内存（shared memory）

片上内存往往比板载内存快很多，因此我们可以使用shared来进行编程，其主要作用有：

  • 线程间交流通道
  • 可编程的cache
  • 转换数据的临时存储器，以减少全局内存访问

shared memory(SMEM)是GPU的重要组成之一，物理上，每个SM包含一个当前正在执行的block中所有thread共享的低延迟的内存池。SMEM使得同一个block中的thread可以相互合作，重用on-chip数据，并且能够减少kernel需要的global memory带宽。

​ 由于shared memory和L1要比L2和global memory更接近SM，shared memory的延迟要比global memory低20到30倍，带宽大约高10倍。

​ 当一个block开始执行时，GPU会分配其一定数量的shared memory，这个shared memory的地址空间会由block的所有thread共享。所以，使用越多的shared memory，能够并行的active就越少。

Shared Memory Allocation（SMEM分配）

​ 我们可以动态或者静态地分配shared memory，其声明可以在kernel内部，也可以作为全局变量。其标识符为__shared__。下面这句声明了一个2D的浮点型数组：

1
__shared__ float tile[size_x][size_y];

如果在kernel中声明的话，其作用域就是在kenel内，否则就是对所有kernel有效。如果shared memory的大小在编译期未知的话，可以使用extern关键字修饰，如下面声明一个未知大小的1D数组：

1
extern __shared__ int tile[];

然后在每个kernel调用时需要使用kernel<<<grid, size*sizeof(int)>>>(...)，而且，只有一维数组才能这样使用。

Shared Memory Banks and Access Mode（共享内存的带宽和访问模式）

pinned memory(固定内存)

在CUDA编程中，内存拷贝是一个非常费时的一个动作。

CUDA-copy-mem

从上图我们可以看出：

  1. CPU和GPU之间的总线是PCIe，双向传输
  2. CPU和GPU之间的数据拷贝是使用DMA机制来实现

我们使用cudaMalloc为GPU分配内存，malloc为CPU分配内存，除此之外，CUDA还提供了自己独有的机制来分配host内存：cudaHostAlloc()，它与malloc不同的是，malloc分配的是可分页的主机内存，而cudaHostAlloc分配的是页锁定的主机内存，也称作固定内存（pinned memory）、不可分页内存。它的一个重要特点是操作系统不会对这块内存进行分页并交换到磁盘上，从而保证了内存始终驻留在物理内存中。因此，操作系统能够安全地使某个应用程序访问该内存的物理地址，因为其不会被破坏或重新定位。

由于GPU知道内存的物理地址，因此就可以使用DMA技术来在GPU和CPU之间复制数据。当使用可分页的内存进行复制时（使用malloc），CUDA驱动程序仍会通过dram把数据传给GPU，但这时数据会执行两遍：

  1. 从可分页内存复制一块到临时的页锁定内存
  2. 从这个页锁定内存复制到GPU上

而页锁定内存只需执行后面这一步，因而速度提高了一倍。

当我们在调用cudaMemcpy(dst, src, …)时，程序会自动检测dest或者src是否为Pinned Memory，若不是，则会先将内容拷进一不可见的Pinned Memory中，然后再进行传输。可以手动指定Pinned Memory，对应的API为cudaHostAlloc(address, size, option)分配地址，cudaFreeHost(pointer)释放地址。注意，分配的内存都是在Host端。

不过，把所有的malloc都替换为cudaHostAlloc()是不对的。固定内存是一把双刃剑，当使用固定内存时，虚拟内存的功能就会失去，尤其是在应用程序中都使用固定内存时，会导致系统内存很快被耗尽，影响自身以及其他应用程序的执行。

所以，建议针对cudaMemcpy()调用中的源内存或者目标内存，才使用页锁定内存，并且在不再使用他们的时候立即释放，而不是关闭应用程序的时候才释放。

零拷贝内存

通常来说，主机不能直接访问设备变量，同时设备变量也不能直接访问主机变量。但是有一个例外：零拷贝内存，主机和设备都可以访问零拷贝内存。它允许你将主机内存直接映射到GPU内存空间上。因此，当你对GPU上的内存解引用时，如果它是基于GPU的，那么你就将获得了全局内存的高速带宽（200GB/s），如果是零拷贝内存的，它会提交一个PCI-E读取事务，很长时间之后，主机会通过PCI-E总线返回数据，带宽降低到16GB/s。

当程序是计算密集型时，零拷贝可能是一项非常有用的技术，它节省了设备显式传输的时间，事实上，它可以将计算和数据传输流水化，而且无需执行显式的内存管理。

当使用零拷贝内存来共享主机和设备间的数据时，必须同步主机和设备间的内存访问，同时更改主机和设备的零拷贝内存中的数据将导致不可预料的后果。

零拷贝内存是固定（不可分页）内存，该内存映射到设备地址空间中，使用以下函数创建一个到固定内存的映射（和锁页内存类似）：

1
cudaError_t cudaHostAlloc(void **pHost, size_t count, unsigned int flags);

这个函数分配了count字节的主机内存，该内存是页面锁定的内存且设备可访问的，用这个函数分配的内存必须用cudaFreeHost函数释放。flags参数可以对以分配的特殊属性进一步进行配置：

  • cudaHostAllocDefault：使cudaHostAlloc函数的行为与cudaMallocHost函数一致。
  • cudaHostAllocPortable：可以返回被CUDA上下文所使用的固定内存，而不仅仅是执行内存分配。
  • cudaHostAllocWriteCombined：返回写结合内存，该内存可以在某些系统配置上通过PCIe总线更快地传输，但是它在大多数主机上不能被有效地读取。因此，写结合内存对缓冲区来说是一个很好的选择，该内存通过设备使用映射的固定内存或主机到设备的传输。
  • cudaHostAllocMapped：零拷贝内存，该标志返回可以实现主机写入和设备读取被映射到设备地址空间的主机内存。

分配好零拷贝内存之后，就可以使用下列函数获取映射到固定内存的设备指针了：

1
cudaError_t cudaHostGetDevicePointer(void **pDevice, void *pHost, unsigned int flags);

该函数返回一个在pDevice中的设备指针，该指针可以在设备上被引用以访问映射得到的固定主机内存。如果设备不支持映射得到的固定内存，则该函数失效。flag保留，始终为0.

在频繁进行读写操作时，使用零拷贝内存作为设备内存的补充将显著降低性能。因为每一次映射到零拷贝内存的传输必须经过PCIe总线，与全局内存相比，延迟也显著增加。

在使用零拷贝内存时，需要检查设备是否支持固定内存映射，cudaDeviceProp的canMapHostMemory成员是一个bool类型值，true表示支持固定内存映射。

一个使用例子如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#include <numeric>
#include <stdio.h>
#include <stdlib.h>
void checkCUDAError(const char *msg) {
    cudaError_t err = cudaGetLastError();
    if (cudaSuccess != err) {
        fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString(err));
        exit(EXIT_FAILURE);
    }
}

__global__ void sumNum( int *data) {

    int i = threadIdx.x + blockIdx.x * blockDim.x;
    if(i<1000000000){
     data[i]=10;
    }
}
int main(void) {
    size_t size = 1*1000000000 * sizeof(int);//4G
    //1.启用零复制
    cudaSetDeviceFlags (cudaDeviceMapHost);
    int* data;
    //2.分配主机内存
    cudaHostAlloc((void**) &data, size,
            cudaHostAllocWriteCombined | cudaHostAllocMapped);
     checkCUDAError("cudaHostAlloc data");
    
     memset(data, 0, 1*1000000000 * sizeof(int));
    int *gpudata;
    //3.将常规的主机指针转换成指向设备内存空间的指针
    cudaHostGetDevicePointer(&gpudata, data, 0);
    checkCUDAError("cudaHostGetDevicePointer");
    sumNum<<<1000000000/1024+1023, 1024>>>(gpudata);
    //注意！！因为下面要打印出来测试，所以要先同步数据，这个函数可以保证cpu等待gpu的kernel函数结束才往下运行。如果数据暂时用不到，可以在整体结束以后再加这句话。明显等待kernel函数结束会占用程序进行的时间。
    cudaDeviceSynchronize();
    for (int i = 99999999; i < 1000000000; i=i+100000000) {
        printf("%d \n", data[i]);
    }
    //记得零拷贝的free是这个函数
    cudaFreeHost(data);
    return 0;
}

注意零拷贝内存还有另一种实现方式，就是直接对malloc申请的内存使用cudaHostRegister进行标志位设置（设置成固定内存），然后使用cudaHostGetDevicePointer获取其设备指针。代码例子如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#include "arrayTools.h"
#include "cudaCode.h"
#include "MyTimer.h"

template<class T>
__global__ void warmup(T *A, T *B, T *C, const int n)
{
    int tid = blockDim.x*blockIdx.x + threadIdx.x;
    if (tid >= n)return;
    C[tid] = A[tid] + B[tid];
}

template<class T>
__global__ void sumVecOnDeviceZeroCopy(T *A, T *B, T *C, const int n)
{
    int tid = blockDim.x*blockIdx.x + threadIdx.x;
    if (tid >= n)return;
    C[tid] = A[tid] + B[tid];
}

template<class T>
__global__ void sumVecOnDeviceZeroCopyRegister(T *A, T *B, T *C, const int n)
{
    int tid = blockDim.x*blockIdx.x + threadIdx.x;
    if (tid >= n)return;
    C[tid] = A[tid] + B[tid];
}
int main()
{
    int nElem=1<<20;
    int nBytes=nElem*sizeof(float);
    float * d_A, * d_B, * d_C;
    float * h_A, * h_B, * gpuRes;
    h_A=(float*)malloc(nBytes);
    h_B=(float*)malloc(nBytes);
    initialInt(h_A, nElem);
    initialInt(h_B, nElem);
    gpuRes=(float*)malloc(nBytes);
    
    dim3 block(512, 1);
    dim3 grid((nElem + block.x - 1) / block.x);

    //将h_A、h_B、h_C三块内存锁定并获取对应的设备地址，进行计算，这样避免了从一开始就得cudaHostAlloc申请零拷贝内存的不便
    CHECK(cudaHostRegister(h_A, nBytes, cudaHostRegisterMapped));
    CHECK(cudaHostGetDevicePointer((void **)&d_A, h_A, 0));
    CHECK(cudaHostRegister(h_B, nBytes, cudaHostRegisterMapped));
    CHECK(cudaHostGetDevicePointer((void **)&d_B, h_B, 0));
    CHECK(cudaHostRegister(gpuRes, nBytes, cudaHostRegisterMapped));
    CHECK(cudaHostGetDevicePointer((void **)&d_C, gpuRes, 0));
    memset(gpuRes, 0, nBytes);
    sumVecOnDeviceZeroCopyRegister<float> << <grid, block >> > (d_A, d_B, d_C, nElem);
    CHECK(cudaDeviceSynchronize());
    
    compareVec(gpuRes, cpuRes, 0, nElem);
    return 0;
}

集成架构和离散架构

这两种架构是常见的异构计算系统架构。

集成架构：cpu和gpu集成在一个芯片上，并且在物理内存上共享内存，在这种架构中，由于无需在PCIe总线上备份，所以零拷贝内存在性能和可编程性方面更好一些。

离散架构：cpu和gpu是分离的，物理内存上也是分离的。数据需要通过PCIe总线进行交互，数据拷贝的耗时和延迟的代价通常也比较大。因此在这种架构下，零拷贝只在特殊情况下有优势。

由于通过映射的固定内存在主机和设备之间是共享的，所以必须同步内存访问来避免任何潜在的数据冲突，这种数据冲突一般是有多线程异步访问相同的内存而引起的。

CUDA编程流程

CUDA编程模型基础

CUDA编程模型是一个异构模型，需要CPU和GPU协同工作。在CUDA中，host和device是两个重要的概念。host指代CPU及其内存，device指代GPU及其内存。CUDA程序既包括host程序，又包括device程序，它们分别在CPU和GPU上运行。同时，host和device之间可以进行通信，这样他们之间可以进行数据拷贝。典型的CUD程序执行流程如下：

  1. 分配host内存，并进行数据初始化
  2. 分配device内存，并从host拷贝数据到device上
  3. 调用CUDA的核函数在device上完成指定的运算
  4. 将device上的运算结果拷贝到host上
  5. 释放host和device上分配的内存

上面最重要的一个过程是调用CUDA的核函数来执行并行计算，kernel是CUDA中一个重要的概念，kernel是在device上线程中并行执行的函数。核函数使用__global__来声明，在调用时需要用<<<grid, block>>>来指定kernel要执行的线程数量，并且每个线程会分配一个唯一的thread ID，这个ID值可以通过核函数内置的变量threadidx来获得。

在CUDA中使用不同的限定词来区分host和device上的函数，三个前缀在前面已经有了说明。

kernel在device上执行时实际上是启动很多线程，一个kernel所启动的所有线程称为一个网格（grid），同一网格上的线程共享相同的全局内存空间，grid是线程的第一层次。同一个网格上的线程又可以分为很多线程块（block），一个线程块里包含很多线程，这是第二层次。两层结构如下图所示：

CUDA线程结构

1
2
3
4
// 对应的代码
dim3 grid(3, 2);
dim3 block(5, 3);
kernel_fun<<<grid, block>>>(params...);

这是一个grid和block均为2-dim的线程组织。grid和block都定义为dim3类型的变量，dim3可以看成是包含三个无符号整数（x, y, z）成员的结构体变量。在定义时，缺省值初始化为1。因此grid和block可以灵活地定义为1-dim，2-dim以及3-dim结构，对于图中结构（水平方向为x轴），展示如上。kernel在调用时也必须通过执行配置<<<grid, block>>>来指定kernel所使用的线程数及结构。

所以，一个线程需要两个内置的坐标变量（blockidx， trheadidx）来唯一标识，它们都是dim3类型变量，其中blockidx指明线程所在的grid中的位置，而threadidx指明线程所在block中的位置。如图中的Thread（1, 1）满足：

1
2
3
4
threadIdx.x = 1;
threadIdx.y = 1;
blockIdx.x = 1;
blockIdx.y = 1;

​ 一个线程块的线程是放在同一个流式多处理器（SM）上的，但是单个SM的资源有限，这导致线程块中的线程数是有限制的，现代GPUs的线程块可支持的线程数达1024个。有时候，我们要知道一个线程在block中的全局ID，就必须知道block的组织结构，这是通过blockDim来获得的。它获取线程块各个维度的大小。对于一个2-dim的block(Dx, Dy)，线程（x，y）的ID值为(x + y * Dx)，此外，如果是3-dim的block(Dx, Dy, Dz)，线程(x, y, z)的ID值为（x + y * Dx + z * Dx * Dy）。另外还有内置变量gridDim，用于获取网格块各个维度的大小。

​ kernel的这种线程组织结构天然适合vector，matrix等运算，如下我们将利用上图2-dim结构实现两个矩阵的加法，每个线程负责处理两个位置的元素相加。线程块大小为(16, 16)，然后将N×N大小的矩阵均分为不同的线程块来执行加法运算。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
//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 threadPerBlock(16, 16);
    dim3 numBlocks(N / threadPerBlock.x, N / threadPerBlock.y);
    //Kernel调用
    MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
}

​ 再次强调CUDA内存模型如下，每个线程都有自己的私有本地内存（Local Memory），而每个线程块都有包含共享内存（Shared Memory），可以被县城快中所有的线程共享，其生命周期和线程块一致。此外，所有的线程都可以访问全局内存（Global Memory），以及一些只读内存块：常量内存（Constant Memory）和纹理内存（Texture Memory）。

CUDA内存模型2

​ GPU有很多CUDA核心，可以充分发挥GPU的并行能力。GPU硬件的一个核心组件是SM（Streaming Multiprocess，流式多处理器）。SM的核心组件包括CUDA核心，共享内存，寄存器等，SM可以并发地执行数百个线程，并发能力就取决于SM所拥有的资源数。当一个kernel被执行时，它的grid中的线程块被分配到SM上。一个线程块只能在一个SM上被调度。SM一般可以调度多个线程块。一个kernel上的各个线程可能被分配给多个SM，所以grid只是逻辑层，SM才是物理层。SM采用的是SIMT（Single-Instruction，Multiple-Thread，单指令多线程）架构，基本的执行单元是线程束（wraps），线程束包含32个线程，这些线程同时执行相同的指令，但是每个线程都包含自己的指令地址计数器和寄存器状态，也有自己独立的执行路径。所以尽管线程束中的线程同时从统一程序地址执行，但是可能有不同的行为，比如遇到了分支结构，一些线程可能进入这个分支，但是另外一些有可能不执行，它们只能死等。因为GPU规定线程束中所有线程在同一周期执行相同的指令，线程束分化会导致性能下降。

​ 当线程块被划分到某个SM上时，它将进一步划分为多个线程束，因为这才是SM的基本执行单元，但是一个SM同时并发的线程束数是有限的，这是因为资源限制，SM要为每个线程块分配共享内存，而也要为每个线程分配独立的寄存器。所以SM的配置会影响其所支持的线程块和线程束的并发数量。总之，就是网格和线程块只是逻辑划分，一个kernel的所有线程其实在物理层不一定同时并发。所以kernel的grid和block的配置的不同，性能会出现差异。由于SM的基本执行单元是包含32个线程的线程束，所以block一般要设置为32的倍数。

​ 在进行CUDA编程前，可以先检查自己的GPU的硬件配置，通过下面的程序可以获得GPU的配置属性：

1
2
3
4
5
6
7
8
9
int dev = 0;
cudaDeviceProp devProp;
CHECK(cudaGetDeviceProperties(&devProp, dev));
std::cout << "使用GPU device " << dev << ": " << devProp.name << 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 << "每个EM的最大线程数：" << devProp.maxThreadsPerMultiProcessor << std::endl;
std::cout << "每个SM的最大线程束数：" << devProp.maxThreadsPerMultiProcessor / 32 << std::endl;

向量加法实例

下面介绍一下CUDA编程中的内存管理API，首先是在device上分配内存的cudaMalloc函数:

1
cudaError_t cudaMalloc(void** devPtr, size_t size);

这个函数和C语言中的malloc类似，在device上申请一定字节大小的显存，其中devPtr是指向所分配内存的指针。同时释放内存使用cudaFree函数，这和C语言中的free函数所对应。另外一个重要的函数是负责host和device之间数据通信的cudaMemcpy函数：

1
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如下：

1
2
3
4
5
6
7
8
9
10
// 两个向量加法，grid和block均为一维
__global__ void add(float* x, float* y, float* z, int n){
    // 获取全局索引，避免线程重复处理同一个元素，相当于把二维数组一维化
    int index = threadIdx.x + blockIdx.x * blockDim.x;
    // 步长，每次有多少个线程在处理
    int stride = blockDim.x * gridDim.x;
    for(int i = index; i < n; i += stride){
        z[i] = x[i] + y[i];
    }
}

其中stride是整个grid的线程数，有时候向量的元素数很多，这时可以在一个线程实现多个元素的加法（元素总数/线程总数的加法），相当于使用了多个grid来处理，这是一种grid-strip loop的方式。不过下面的例子一个线程只处理一个元素，所以kernel里面的循环是不执行的。下面实现具体的向量加法：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
int main(){
    int N = 1 << 20;
    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] = 10.0;
        y[i] = 20.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
    add<<<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 = fmax(maxError, fabs(z[i] - 30.0));
    }
    std::cout << "最大误差为：" << maxError << std::endl;
    
    // 释放device内存
    cudaFree(d_x);
    cudaFree(d_y);
    cudaFree(d_z);
    
    // 释放host内存
    free(x);
    free(y);
    free(z);
    
    return 0;
}

这里我们向量的大小是1<<20，而block的大小为256,那么grid大小为（2^20 / 256）= 4096，kernel的线程层级如下图所示：

kenel线程层级

同时，block不是越大越好，要适当选择。

在上面的实现中，我们需要单独在host和device上进行内存分配，并且进行数据拷贝，这是很容易出错的，在CUDA 6.0引入了统一内存(Unified Memory)来避免这种麻烦，简单地说就是统一内存使用一个托管内存来共同管理host和device的内存，并且自动在host和device中进行数据传输。CUDA中使用cudaMallocManaged函数分配托管内存：

1
cudaError_t cudaMallocManaged(void **devPtr, size_t size, unsigned int flag = 0);

利用统一内存，可以将上面的程序简化如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
int main(){
    int N = 1 << 20;
    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] = 10.0;
        y[i] = 20.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 = fmax(maxError, fabs(z[i] - 30.0));
    }
    std::cout << "最大误差为：" << maxError << std::endl;
    
    // 释放内存
    cudaFree(x);
    cudaFree(y);
    cudaFree(z);
    
    return 0;
}

需要注意的是kernel执行和host是异步的，由于托管内存自动进行数据传输，这里要用cudaDeviceSynchronize()函数保证device和host同步，这样后面才可以正确访问kernel计算的结果。

矩阵乘法实例

设输入矩阵为A和B，得到的矩阵为C = A×B。实现的思路就是每个线程计算C的一个元素的值Ci,j，对于矩阵运算，应该选用grid和block为2-D的，首先定义矩阵的结构体：

1
2
3
4
5
6
// 矩阵类型，行优先，M(row, col) = *(M.elements + row * M.width + col)
struct Matrix {
    int width;
    int height;
    float *elements;
}

矩阵乘法的实现图如下：

矩阵乘法实现图

然后实现矩阵的核函数，这里我们定义了两个辅助的__device__函数分别用于获取矩阵的元素值和为矩阵元素赋值，具体代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
// 获取矩阵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);
}

最后我们使用统一内存来编写矩阵相乘的测试实例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
int main(){
    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));
    }
    std::cout << "最大误差为：" << maxError << std::endl;
    
    return 0;
}

这里设计的线程block为（32, 32），grid大小为（32, 32），最终测试结果如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
nvprof cuda9.exe
==16304== NVPROF is profiling process 16304, command: cuda9.exe
最大误差: 0
==16304== Profiling application: cuda9.exe
==16304== Profiling result:
            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
 GPU activities:  100.00%  1.32752s         1  1.32752s  1.32752s  1.32752s  matMulKernel(Matrix*, Matrix*, Matrix*)
      API calls:   83.11%  1.32762s         1  1.32762s  1.32762s  1.32762s  cudaDeviceSynchronize
                   13.99%  223.40ms         6  37.233ms  37.341us  217.66ms  cudaMallocManaged
                    2.81%  44.810ms         1  44.810ms  44.810ms  44.810ms  cudaLaunch
                    0.08%  1.3300ms        94  14.149us       0ns  884.64us  cuDeviceGetAttribute
                    0.01%  199.03us         1  199.03us  199.03us  199.03us  cuDeviceGetName
                    0.00%  10.009us         1  10.009us  10.009us  10.009us  cuDeviceTotalMem
                    0.00%  6.5440us         1  6.5440us  6.5440us  6.5440us  cudaConfigureCall
                    0.00%  3.0800us         3  1.0260us     385ns  1.5400us  cudaSetupArgument
                    0.00%  2.6940us         3     898ns     385ns  1.5390us  cuDeviceGetCount
                    0.00%  1.9250us         2     962ns     385ns  1.5400us  cuDeviceGet

==16304== Unified Memory profiling result:
Device "GeForce GT 730 (0)"
   Count  Avg Size  Min Size  Max Size  Total Size  Total Time  Name
    2051  4.0000KB  4.0000KB  4.0000KB  8.011719MB  21.20721ms  Host To Device
     270  45.570KB  4.0000KB  1.0000MB  12.01563MB  7.032508ms  Device To Host