CUDA编程模型-以数组加法为例
通于一而万事毕，无心得而鬼神服。
本文以数组加法为例，简述CUDA的编程模型
基本概念
确保了解《CUDA的线程结构》中的内容再进行下去，因为对于线程结构本章会省略一部分描述
编程模型
首先，编程模型是写程序时可被控制的部分，通过编程模型可以控制计算设备的行为；其次，编程模型对于编译器和库函数来说是一层通信抽象；
再次，对于GPU来说，编程模型可分为
- 核函数
- 内存管理
- 线程管理
- 流
最后，对于Cuda编程，可抽象出几个过程
- 领域层（业务主导要实现的目标）
- 逻辑层（基于逻辑组织线程代码）
- 硬件层（思考如何将线程映射到机器上，以及思考如何提升性能）
接下来，会从以下方面说明Cuda的编程结构
- 内存
- 线程
核函数
- 启动
- 编写
- 验证
- 错误处理
内存管理
首先应该理解异构环境下内存的样貌：通过PCIe总线分割和通信的两种内存：主机内存和设备内存
一个异构环境，通常有多个CPU多个GPU，他们都通过PCIe总线相互通信，也是通过PCIe总线分隔开的。所以我们要区分一下两种设备的内存：
- 主机：CPU及其内存
- 设备：GPU及其内存
注意这两个内存从硬件到软件都是隔离的（CUDA6.0 以后支持统一寻址），我们目前先不研究统一寻址，我们现在还是用内存来回拷贝的方法来编写调试程序，以巩固大家对两个内存隔离这个事实的理解。
CUDA的内存管理是C风格的，列表如下：
| 标准C函数 | CUDA C 函数 | 说明 |
|---|---|---|
| malloc | cudaMalloc | 内存分配 |
| memcpy | cudaMemcpy | 内存复制 |
| memset | cudaMemset | 内存设置 |
| free | cudaFree | 释放内存 |
另外，CUDA中设备内存也是分层次的，分为单网格中公用的全局内存和单线程块中公用的共享内存，大致如图：
线程管理
在《CUDA的线程结构》中，曾经用草图展示过cuda线程的组织形式，这里用更好的图展示一下：
另外，不同线程块的线程物理隔离无法相互影响，同一个线程块中的线程可以实现以下两种操作：
- 同步
- 共享内存
最后再重申一下《CUDA的线程结构》中提过的，如何标识一个线程：
- blockIdx（线程块在线程网格内的位置索引）
- threadIdx（线程在线程块内的位置索引）
这两个内置结构体基于 uint3 定义，包含三个无符号整数xyz的结构
还有blockIdx中三个字段的范围和threadIdx中三个字段的范围：
- gridDim
- blockDim
dim3类型(基于uint3定义的数据结构)的变量，也包含三个字段xyz
另外，线程块在3个轴上的个数和一个线程块的线程数是有限制的：
- 线程块在xyz三轴上分别最多：
2^31-1,65535,65535 - 一个线程块线程总数最多1024，在三维上最多:
1024,1024,64
核函数
之前hello world时已经用到了核函数，但是没有给一个明确定义，这里给出：
核函数就是在CUDA模型上诸多线程中运行的那段串行代码，这段代码在设备上运行，用NVCC编译，产生的机器码是GPU的机器码，所以我们写CUDA程序就是写核函数
- 第一步我们要确保核函数能正确的运行产生正切的结果
- 第二步优化CUDA程序的部分，无论是优化算法，还是调整内存结构，线程结构都是要调整核函数内的代码，来完成这些优化的。
接下来将以数组加法为例来实现和使用核函数，以及应用一些内存管理的函数
单机版数组加法
无需多言，直接上代码，明确我们要干什么：
/**
* 单机版:两数组相加
*
*/
#include<iostream>
const double a=1.23;
const double b=3.45;
const double c=4.68;
const int N = 100000000;
const double EPSILON = 1.0e-15;
// 函数1: 数组相加(数组默认是动态分配的)
void add(const double *x,const double *y,double *z,const int N) {
for (int i=0;i<N;i++) {
z[i] = x[i]+y[i];
}
}
// 函数2: 检查结果
void check(const double *z,const int N) {
bool has_error = false;
for (int i=0;i<N;i++) {
if (fabs(z[i]-c) > EPSILON) {
has_error = true;
break;
}
}
printf("%s\n",has_error?"Has errors":"No errors");
}
int main() {
//创建数组xyz
double* x = new double[N];
double* y = new double[N];
double* z = new double[N];
//初始化数组
for (int i=0;i<N;i++) {
x[i] = a;
y[i] = b;
z[i] = 0;
}
//计算数组
add(x,y,z,N);
//检查结果
check(z,N);
//删除数组指针
delete[] x;
delete[] y;
delete[] z;
return 0 ;
}CUDA程序和内存管理
先看代码，比起hello world，这是一个完整的CUDA程序，它包含了下图所示的结构：
/**
* cuda版:两数组相加
*
*/
#include<iostream>
const double a=1.23;
const double b=3.45;
const double c=4.68;
const int N = 100000000;
const int M = sizeof(double)*N;
const double EPSILON = 1.0e-15;
// 函数1: 数组相加(数组默认是动态分配的)
__global__ void add(const double *x,const double *y,double *z) {
const int n = blockDim.x*blockIdx.x + threadIdx.x;
z[n] = x[n] + y[n];
}
// 函数2: 检查结果
void check(const double *z,const int N) {
bool has_error = false;
for (int i=0;i<N;i++) {
if (fabs(z[i]-c) > EPSILON) {
has_error = true;
break;
}
}
printf("%s\n",has_error?"Has errors":"No errors");
}
int main() {
//主机: 创建数组xyz
double* x = new double[N];
double* y = new double[N];
double* z = new double[N];
//初始化数组
for (int i=0;i<N;i++) {
x[i] = a;
y[i] = b;
z[i] = 0;
}
//主机: 创建设备数组xd,yd
double* xd;
double* yd;
cudaMalloc((void **)&xd,M);
cudaMalloc((void **)&yd,M);
//主机: 复制主机数据到设备
cudaMemcpy(xd,x,M,cudaMemcpyHostToDevice);
cudaMemcpy(yd,y,M,cudaMemcpyHostToDevice);
//主机: 指定设备运行参数
//const dim3 block_size(128,1,1);
const int block_size=128;
const dim3 grid_size(N/block_size,1,1);
//设备: 计算数组
add<<<grid_size,block_size>>>(xd,yd,yd);
//主机: 复制设备结果到主机
cudaMemcpy(z,yd,M,cudaMemcpyDeviceToHost);
//检查结果
check(z,N);
//删除数组指针
cudaFree(xd);
cudaFree(yd);
delete[] x;
delete[] y;
delete[] z;
return 0 ;
}先把核函数的部分放在一边，代码里有以下几点可以注意：
- 隐式设备初始化：了解过opengl或d3d程序的可以明显感觉，CUDA是没有显式初始化设备(即GPU)的函数。在第一次调用一个和设备管理及版本查询功能无关的运行时 API函数时，设备将自动初始化。
内存管理一节已经大致了解了内存管理的相关函数，这里就是应用的地方，需要注意的是cudaMalloc和cudaMemcpy
- cudaMalloc的第一个参数是指针的指针
- cudaMemcpy的最后一个参数规定了数据传输的方式（从名称就可以看出）
CUDA程序和核函数
然后来看核函数，可分为三个方面：
- 启动核函数
- 编写核函数
- 验证核函数
首先，CPU是一个控制者，运行核函数，要从CPU发起，第一，核函数是异步的，第二，如《CUDA的线程结构》所说，三个尖括号<<<grid,block>>>内是对设备代码执行的线程结构的配置：
- 网格的布局（线程块）
- 线程块的布局（线程）
其次，同步函数有两种
cudaDeviceSynchronize显式同步- 涉及到设备端不执行完，主机没办法进行，比如内存拷贝函数，隐式同步
然后看到核函数的写法，cuda中不仅有核函数，还有一些其他类型，通过不同的标识表明，如核函数是__global__：
| 限定符 | 执行 | 调用 | 备注 |
|---|---|---|---|
| global | 设备端执行 | 可以从主机调用也可以从计算能力3以上的设备调用 | 必须有一个void的返回类型 |
| device | 设备端执行 | 设备端调用 | |
| host | 主机端执行 | 主机调用 | 可以省略 |
有些函数可以同时定义为 device 和 host ，这种函数可以同时被设备和主机端的代码调用，主机端代码调用函数很正常，设备端调用函数与C语言一致，但是要声明成设备端代码，告诉nvcc编译成设备机器码，同时声明主机端设备端函数，那么就要告诉编译器，生成两份不同设备的机器码。
其他方面核函数基本和普通函数一致，但有一些限制：
- 除了统一寻址，只能访问设备内存
- 必须有void返回类型
- 不支持可变数量的参数
- 不支持静态变量
- 显示异步行为
- 不可作为成员，但可以包装
- 在动态并行的前提下，可以递归调用
另外设备函数可用__noinline__ __forceinline__来作为内联函数和非内联函数
CUDA小技巧，当我们进行调试的时候可以把核函数配置成单线程的：
下面是一些例子：
- 函数的使用
- 斐波那契数列，本程序意义不大，旨在展示
__device__ __host__的共用 - 斐波那契数列，本程序意义不大，旨在展示递归的使用
/**
* 函数的使用
*/
#include<iostream>
const double a=1.23;
const double b=3.45;
const double c=4.68;
const int N = 100000000;
const int M = sizeof(double)*N;
const double EPSILON = 1.0e-15;
// 函数1: 数组相加(数组默认是动态分配的)
//TODO: 注意到三种写法效率是一样的
__device__ double addfunc1(double x,double y) {
return x+y;
}
__device__ double addfunc2(double x,double y,double *z) {
*z = x+y;
}
__device__ double addfunc3(double x,double y,double &z) {
z= x+y;
}
__global__ void add(const double *x,const double *y,double *z) {
const int n = blockDim.x*blockIdx.x + threadIdx.x;
if (n<N) {
z[n] = x[n] + y[n];
//z[n] = addfunc1(x[n],y[n]);
//addfunc2(x[n],z[n],&z[n]);
//addfunc3(x[n],y[n],z[n]);
}
}/**
* 斐波那契数列，本程序意义不大，旨在展示__device__ __host__的共用
*
*/
#include<iostream>
const int N = 10;
const int M = sizeof(int)*N;
// 函数1: 斐波那契数列
__device__ __host__ void fibN(int* A, const int N) {
if (N<0) return;
if (A[N]==-1) {
if (N>1) {
fibN(A,N-1);
fibN(A,N-2);
A[N] = A[N-1] + A[N-2];
}
if (N<=1) {
A[N]=N;
}
}
}
__global__ void fibNN(int* A, const int N) {
fibN(A, N);
}
// 函数2: 检查结果
//0、1、1、2、3、5、8、13、21、34
void check(const int *A,const int N) {
for (int i=0;i<N;i++) {
std::cout<<A[i]<<" ";
}
std::cout<<std::endl;
}
int main() {
//主机: 创建数组xyz
int* x = new int[N];
//初始化数组
for (int i=0;i<N;i++) {
x[i] = -1;
}
//主机: 创建设备数组xd,yd
int* xd;
cudaMalloc((void **)&xd,M);
//主机: 复制主机数据到设备
cudaMemcpy(xd,x,M,cudaMemcpyHostToDevice);
//主机: 指定设备运行参数
const dim3 block_size(128,1,1);
const dim3 grid_size(128,1,1);
//设备: 计算数组
fibNN<<<grid_size,block_size>>>(xd,N-1);
//主机: 复制设备结果到主机
cudaMemcpy(x,xd,M,cudaMemcpyDeviceToHost);
//检查结果
check(x,N);
for (int i=0;i<N;i++) {
x[i] = -1;
}
fibN(x,N-1);
check(x,N);
//删除数组指针
cudaFree(xd);
delete[] x;
return 0 ;
}/**
* cuda版:两数组相加
*
*/
#include<iostream>
const double a=1.23;
const double b=3.45;
const double c=4.68;
const int N = 100000000;
const int M = sizeof(double)*N;
const double EPSILON = 1.0e-15;
// 函数1: 数组相加(数组默认是动态分配的)
__global__ void add(const double *x,const double *y,double *z) {
const int n = blockDim.x*blockIdx.x + threadIdx.x;
if (n<N) {
z[n] = x[n] + y[n];
}
}
// 函数2: 检查结果
void check(const double *z,const int N) {
bool has_error = false;
for (int i=0;i<N;i++) {
if (fabs(z[i]-c) > EPSILON) {
has_error = true;
break;
}
}
printf("%s\n",has_error?"Has errors":"No errors");
}
int main() {
//主机: 创建数组xyz
double* x = new double[N];
double* y = new double[N];
double* z = new double[N];
//初始化数组
for (int i=0;i<N;i++) {
x[i] = a;
y[i] = b;
z[i] = 0;
}
//主机: 创建设备数组xd,yd
double* xd;
double* yd;
cudaMalloc((void **)&xd,M);
cudaMalloc((void **)&yd,M);
//主机: 复制主机数据到设备
cudaMemcpy(xd,x,M,cudaMemcpyHostToDevice);
cudaMemcpy(yd,y,M,cudaMemcpyHostToDevice);
//主机: 指定设备运行参数
//const dim3 block_size(128,1,1);
const int block_size=128;
const dim3 grid_size(N/block_size,1,1);
std::cout << "Size: " << N/block_size <<"=="<<block_size <<std::endl;
//设备: 计算数组
//TODO: 注意到这里使用了会超过N的执行参数
add<<<800000,block_size>>>(xd,yd,yd);
//主机: 复制设备结果到主机
cudaMemcpy(z,yd,M,cudaMemcpyDeviceToHost);
//检查结果
check(z,N);
//删除数组指针
cudaFree(xd);
cudaFree(yd);
delete[] x;
delete[] y;
delete[] z;
return 0 ;
}/**
* 斐波那契数列，本程序意义不大，旨在展示递归的使用
*
*/
#include<iostream>
const int N = 10;
const int M = sizeof(int)*N;
// 函数1: 斐波那契数列
__global__ void fibN(int* A, const int N) {
if (N<0) return;
if (A[N]==-1) {
if (N>1) {
fibN<<<1,1>>>(A,N-1);
fibN<<<1,1>>>(A,N-2);
A[N] = A[N-1] + A[N-2];
}
if (N<=1) {
A[N]=N;
}
}
}
// 函数2: 检查结果
//0、1、1、2、3、5、8、13、21、34
void check(const int *A,const int N) {
for (int i=0;i<N;i++) {
std::cout<<A[i]<<" ";
}
std::cout<<std::endl;
}
int main() {
//主机: 创建数组xyz
int* x = new int[N];
//初始化数组
for (int i=0;i<N;i++) {
x[i] = -1;
}
//主机: 创建设备数组xd,yd
int* xd;
cudaMalloc((void **)&xd,M);
//主机: 复制主机数据到设备
cudaMemcpy(xd,x,M,cudaMemcpyHostToDevice);
//主机: 指定设备运行参数
const dim3 block_size(128,1,1);
const dim3 grid_size(128,1,1);
//设备: 计算数组
fibN<<<grid_size,block_size>>>(xd,N-1);
//主机: 复制设备结果到主机
cudaMemcpy(x,xd,M,cudaMemcpyDeviceToHost);
//检查结果
check(x,N);
//删除数组指针
cudaFree(xd);
delete[] x;
return 0 ;
}
//0 1 1 0 -2 -2 -2 -4 -2 -6 结果错误 但是可见可以这么用验证和错误检测
所有编程都需要对错误进行处理，早期的编码错误，编译器会帮我们搞定，内存错误也能观察出来，但是有些逻辑错误很难发现，甚至到了上线运行时才会被发现，而且有些厉害的bug复现会很难，不总出现，但是很致命，而且CUDA基本都是异步执行的，当错误出现的时候，不一定是哪一条指令触发的，这一点非常头疼；这时候我们就需要对错误进行防御性处理了
错误检测宏如下：
#pragma once
#include <cstdio>
#define CHECK(call) \
do \
{ \
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); \
} \
} while (0)
cudaGetErrorString可把error_code转为相应的提示信息
用这个宏函数包装大部分CUDA运行时API函数。有一个例外是cudaEventQuery()函数，因为它很有可能返回cudaErrorNotReady,但又不代表程序出错了
#include "error.cuh"
#include <cmath>
#include <cstdio>
const double EPSILON = 1.0e-15;
const double a = 1.23;
const double b = 2.34;
const double c = 3.57;
void __global__ add(const double *x, const double *y, double *z, const int N);
void check(const double *z, const int N);
int main(void)
{
const int N = 100000000;
const int M = sizeof(double) * N;
double *h_x = (double*) malloc(M);
double *h_y = (double*) malloc(M);
double *h_z = (double*) malloc(M);
for (int n = 0; n < N; ++n)
{
h_x[n] = a;
h_y[n] = b;
}
double *d_x, *d_y, *d_z;
CHECK(cudaMalloc((void **)&d_x, M));
CHECK(cudaMalloc((void **)&d_y, M));
CHECK(cudaMalloc((void **)&d_z, M));
//TODO:此处错误复制方向搞反了
CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyDeviceToHost));
CHECK(cudaMemcpy(d_y, h_y, M, cudaMemcpyDeviceToHost));
const int block_size = 128;
const int grid_size = (N + block_size - 1) / block_size;
add<<<grid_size, block_size>>>(d_x, d_y, d_z, N);
CHECK(cudaMemcpy(h_z, d_z, M, cudaMemcpyDeviceToHost));
check(h_z, N);
free(h_x);
free(h_y);
free(h_z);
CHECK(cudaFree(d_x));
CHECK(cudaFree(d_y));
CHECK(cudaFree(d_z));
return 0;
}
void __global__ add(const double *x, const double *y, double *z, const int N)
{
const int n = blockDim.x * blockIdx.x + threadIdx.x;
if (n < N)
{
z[n] = x[n] + y[n];
}
}
void check(const double *z, const int N)
{
bool has_error = false;
for (int n = 0; n < N; ++n)
{
if (fabs(z[n] - c) > EPSILON)
{
has_error = true;
}
}
printf("%s\n", has_error ? "Has errors" : "No errors");
}用上述方法不能捕捉调用核函数的相关错误，因为核函数不返回任何值(回顾一下，核函数必须用void修饰)。有一个方法可以捕捉调用核函数可能发生的错误
#include "error.cuh"
#include <cmath>
#include <cstdio>
const double EPSILON = 1.0e-15;
const double a = 1.23;
const double b = 2.34;
const double c = 3.57;
void __global__ add(const double *x, const double *y, double *z, const int N);
void check(const double *z, const int N);
int main(void)
{
const int N = 100000000;
const int M = sizeof(double) * N;
double *h_x = (double*) malloc(M);
double *h_y = (double*) malloc(M);
double *h_z = (double*) malloc(M);
for (int n = 0; n < N; ++n)
{
h_x[n] = a;
h_y[n] = b;
}
double *d_x, *d_y, *d_z;
CHECK(cudaMalloc((void **)&d_x, M));
CHECK(cudaMalloc((void **)&d_y, M));
CHECK(cudaMalloc((void **)&d_z, M));
CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice));
CHECK(cudaMemcpy(d_y, h_y, M, cudaMemcpyHostToDevice));
//TODO:错误在线程块的大小太大
const int block_size = 1280;
const int grid_size = (N + block_size - 1) / block_size;
add<<<grid_size, block_size>>>(d_x, d_y, d_z, N);
CHECK(cudaGetLastError());
CHECK(cudaDeviceSynchronize());
CHECK(cudaMemcpy(h_z, d_z, M, cudaMemcpyDeviceToHost));
check(h_z, N);
dim3 a;
free(h_x);
free(h_y);
free(h_z);
CHECK(cudaFree(d_x));
CHECK(cudaFree(d_y));
CHECK(cudaFree(d_z));
return 0;
}
void __global__ add(const double *x, const double *y, double *z, const int N)
{
const int n = blockDim.x * blockIdx.x + threadIdx.x;
if (n < N)
{
z[n] = x[n] + y[n];
}
}
void check(const double *z, const int N)
{
bool has_error = false;
for (int n = 0; n < N; ++n)
{
if (fabs(z[n] - c) > EPSILON)
{
has_error = true;
}
}
printf("%s\n", has_error ? "Has errors" : "No errors");
}CHECK(cudaGetLastError());
CHECK(cudaDeviceSynchronize());其中，第一条语句的作用是捕捉第二个语句之前的最后一个错误，第二条语句的作用是同步主机与设备。之所以要同步主机与设备，是因为核函数的调用是异步的，即主机发出调用核函数的命令后会立即执行后面的语句，不会等待核函数执行完毕。
