CUDA程序设计(一)

为什么需要GPU

几年前我启动并主导了一个项目，当时还在谷歌，这个项目叫谷歌大脑。该项目利用谷歌的计算基础设施来构建神经网络。

规模大概比之前的神经网络扩大了一百倍，我们的方法是用约一千台电脑。这确实使深度学习取得了相当大的进展。用到相当多的

计算机。不久之后我发现，之前我并没意识到，用一千台电脑是一项非常昂贵的技术。因此，我和我的朋友，意识到，利用一种

不同的技术，仅用三台电脑，而非一千台，就可以做到这点，而秘诀就是利用GPU技术。 

                                                              ---Andrew Ng [The Big Talk：深度学习与在线教育]

GPGPU是最近几年才流行起来的Project，并且在科学计算上使用频繁。

深度学习在06年立项，当时的数据集小的可怜，只需要CPU仿真一下。

09年，Stanford的CV实验室公布ImageNet后，深度学习开始拥抱大数据。

浅层结构+大数据尚且依赖云计算，深度结构+大数据的需求则是云计算提供不了的。

12年，Hinton组的Alex开源了cuda_convnet，正式大规模使用GPU。衍生出Theano、Caffe等项目。

在今天，如果你要修改这些项目代码，包括自己的创新改良，不会CUDA是不行的。

这也是深度学习带来的一个全新领域，它要求研究者不仅要理论强，建模强，程序设计能力也要过硬，不能纸上谈兵。

单从这点来看，深度学习依然和传统机器学习学派分道而行（多使用CPU多线程、云计算）。

大规模分布式CPU云应该尽量避免，尽管面子大，唬人效果好，但毕竟计算效率低，平台组建麻烦。

高效的程序设计，应当从单机中，榨取更多的计算力，进而考虑多机联合。

 

CUDA物理计算体系

1.1 CUDA特性

CUDA可以说是介于NVIDIA Driver和DirectX等库的中间产物。

传统游戏开发者，无须考虑太多GPU硬件架构，毕竟A卡、N卡水火不容。

通过设备商提供的Driver，DirectX或者OpenCL可以将不同架构的硬件底层封装，提供给开发者共有的API。

而CUDA，则可以看成是Driver、DX之流的上级。

为了泛型编程(C、C++、Fortran多语言）、以及榨取更多的计算力，NVIDIA对OpenCL进行的改装，

贴合自己的GPU硬件架构，量身定做出CUDA。

较游戏程序员不同，CUDA程序员主要工作，就是把握硬件架构，在算法理论时间复杂度下，

将算法串行执行体系，改组为并行执行体系。违背硬件的并行设计，只会让你的程序跑的百倍之慢。

1.2 NVIDIA GPU三大经典架构（Fermi、Kepler、Maxwell)

与海贼王里CP9道力值 (道力可以反映一个人的战斗力，但不能完全决定一个人的实力)设定相似：

NVIDIA为它在不同成长阶段卖出的产品，规定了计算能力体系：

计算能力1.0是跑CUDA的最低条件，这一时期的代表作是8800GT家族。

Fermi架构的计算能力是2.0, Kepler是3.0，Maxwell是4.0。

1.2.1 年轻有为的Fermi（费米）

Fermi架构，也就是GTX 4XX系卡，这算是CUDA的一个重要转折点——它把流处理器（SP）改名为"CUDA核心"。

NVIDIA的[Fermi架构图]：

NVIDIA的物理计算体系以分为SM、CUDA核心两大部分。

官方称SM为Stream Multiprocessors（流多处理器），民间多译为SM计算阵列。

物理计算体系对于程序员是透明的，需要调度NVIDIA提供的逻辑计算体系（线程Grid、线程Block、线程Thread)

对于一个指定好逻辑计算体系的任务，GPU会按照负载均衡的原则，将任务平均分至各个SM阵列，

对于每个SM阵列，就调度它手下那一伙CUDA核心干活。

由于每个SM阵列的CUDA核心有限，NVIDIA规定，Fermi架构，每个SM最多并行执行1024个线程。

当然，实际任务中，每个SM会分到几百万个线程，这时候，就只能小部分并行，然后再串行了。

Fermi 1.0架构，官方设计是16组SM，512SP，然而旗舰GTX480最后只弄出了15组，480SP，顺次阉割出GTX470、GTX460。

Fermi 2.0架构，旗舰GTX580，总算达到设计图要求，达到了16组，512SP，顺次阉割出了GTX570，GTX560。

1.2.2 老当益壮的Kepler（开普勒）

Kepler架构最大变化在于， 对每个SM阵列，将SP数量扩大到6倍，达到192SP。谓之曰SMX阵列。

每个SMX阵列，包含192个CUDA核心，单次并行吞吐量是2048个线程。

Kepler 1.0架构，官方设计是15组SM，2880SP，然而旗舰GTX580最后只弄出了8组，1536SP，顺次阉割出GTX570、GTX560。

Kepler 2.0架构，旗舰GTX680，总算达到设计图要求，达到了15组，2880SP，顺次阉割出了GTX670，GTX660。

特别版，GTX Titan Z，直接把两块GK110并在一起，合出了30组，5760SP，同时支持双精度浮点计算。

其阉割掉双精度之后，就是GTX690。

值得一提的是，GTX游戏卡直接把双精度阉割掉了，因为只有Tesla做科学计算的时候，才会用双精度浮点运算。

整体来看，Kepler并没有提升SM阵列数量（在产量情况下，GTX 5XX居然只造出了8组）。但是每个阵列的吞吐量扩大到2倍（1024-2048）。

速度也较之提升（CUDA核心扩大到6倍）。

在GTX 6XX系中，15组SM阵列，总算达到了Fermi的两倍吞吐量。

1.2.3 "长者"风范的Maxwell（麦克斯韦）

Maxwell架构是老黄的无奈之举。因为台积电把20nm工艺让给了ARM系（Apple和高通）。

还是基于28nm的Maxwell，继续在SM上大刀阔斧闹改革，将192SP降低为128SP，谓之曰SMM阵列。

Maxwell 最新架构，官方设计是16组SM，2048SP，为旗舰GTX980，顺次阉割出GTX970，GTX960。

特别版，GTX TitanX，24组SM，3072SP，较之TitanZ，阉掉了双精度浮点数支持。

TitanX是老黄在GTC 2015向DL界主推的一块民用卡，因为DL无需高精度浮点，用Tesla太奢侈。

(仔细分析老黄那自信满满进攻DL领域的表情，拿TitanX打游戏实在有点.....)

降低了SP数量之后，Maxwell得到了功耗和性能之间的权衡调整。

1.3 设备信息

编译CUDA Examples下的deviceQuery，你可以看到GPU的基本信息。

红字标出的，是设计并行程序时需要参考的几个重要设备参数。

 

CUDA逻辑计算体系

2.1 线程网格(Grid)、线程块（Block)、线程(Thread)、线程束（Warp)

2.1.1 内核函数

内核函数是并行计算中最基本的单元函数，其特点是：

统一的处理逻辑代码，分布并行掌控不同区域的数据，以此达到多区域数据联动并行执行。

NVIDIA为了CPU在逻辑上，能调度GPU计算的函数，规定了统一的格式。

以__global__限定符为始，声明：__global__ void helloworld()。

__global__意思为，GPU执行，CPU调用

调用时，需要分配  <<<线程块，块内线程数>>>。

如执行helloword，使用1个线程块，块内使用256个线程，则

helloworld<<<1,256>>>

2.1.2 线程网格(Grid)

线程网格在编程时并不存在，它只是抽象上的并行网格体系。

不同种类的内核函数，每种内核函数调度数个的线程块，这数个线程块逻辑上被判为一个Grid。

2.1.3 线程块(Block)

线程块是一个3D结构，强调3D坐标系时，需要以dim3类型声明三维大小。

dim3是个结构体, 成员x、y、z，代表方向轴尺度。

如helloworld<<<dim3(1,1,1),256>>>。

当然，大部分操作基本使用的是1D坐标系，线程块默认全部扩展到X轴上。

一般写成helloworld<<<1,256>>>。

通常在内核函数内，需要获取线程块编号，以便对数据集的不同区域处理，四大重要属性：

☻dim3 gridDim(不是指有多少Grid，而是指一个Grid有多少Block)

☻dim3 blockDim(不是指有多少Block，而是指一个Block有多少Thread)

☻dim3 blockIdx

☻dim3 threadIdx

对于1D坐标系，有int tid= (blockDim.x*blockIdx.x) + threadIdx.x;

tid指明当前线程的编号，是内核函数里最基本的控制变量。

int step=(blockDim.x * gridDim.x);

由于CUDA限制每个Block的线程数(2.0以上通常使用1024，以下通常使用512）

所以在常规元素分解模型中，通常把每个Block的线程数设置为常量(固定不动)

这时，有两个策略：

① 其一，不固定Block：

这种方法最为常用，由于CUDA对每个任务而言，对Block数量的限制很松，

如图：

这时候，可以采取为每个线程分配一个元素的方法，用

$BLOCKS=(N + THREADS - 1) / THREADS$

算出一个动态的Block数量的需求，这时候，for(i=tid;i<N;i+=step)等效于for(i=tid;i<N;i+=1)

因为这个循环根本不会执行第二次。

① 其二，固定Block数量：

这时，这时候，为了跑完全部的N个元素，有些Thread会启动人工循环。

i+=step会将元素坐标继续跳转，因为N必然大于step，你不能用+1来取剩余的元素吧？

这两种方法本质上是等效的，由于在物理执行时，同时并行线程最多大概是3072，

几百万、甚至几千万的Block会被CUDA扔到等待队列里，由CUDA自己安排自动循环，

有时甚至比你的人工循环更高效，所以，通常用①的方法，为了保持写法一致，step也会作为默认跳转量。

2.1.4 线程(Thread)

CUDA逻辑体系里最基本的执行单位，等效于CPU的线程。

内核函数一旦被指明了线程块大小，线程大小后，每个线程就分配到了一个内核函数的副本。

区别这些的线程的唯一方法就是线程编号tid，通过tid，让不同线程窥视数据集的不同部分。

用相同的逻辑代码，执行数据空间的不同子集。

2.1.5 线程束(Warp)

线程束对用户透明，它是NVIDIA强行规定的。目前显卡都固定为32。

逻辑上，线程束将32个线程编为一组。

一般微机系统，如8086，它的访存模式是串行的。每一个总线周期，吞一个字节进来。

NVIDIA的GPU在一个总线周期内，能够最大吞32*4=128字节。

前提是当个线程束内的线程，逐序访问显存，这特别需要设计数据存储形式。

使用线程束的目的是掩盖单个总线周期过长的问题，通常要跑500~600个T周期。

 

简单样例: HelloWorld、向量加法

3.1 CUDA程序基本结构

CUDA的编译器NVCC只能编译.cu文件。cu文件可以看作是h&cpp文件的混合体。

所有代码可以写在一个cu文件里，但是毫无可读性。

CUDA最大的亮点在于提供C(C++)接口，所以设计模式应该保持传统C程序风格。

☻风格一：将不同函数写进不同cu文件里，模块化。

看起来挺好的，实际操作起来比较麻烦。

NVCC编译器在生成程序时，不支持跨文件link，也就是说，#include一个cu文件，会导致重定义。

VS中会这么报错：

error LNK2005: "void __cdecl cudaStrcpy(char *,char *)" (?cudaStrcpy@@YAXPAD0@Z) 已经在 gpu_helloworld.cu.obj 中定义

所以，NVCC编译完cu文件后，要在项目属性里把冲突的cu文件从生成中排除，然后link。

☻风格二：h、cpp、cu混用

上图中指明了，h&cpp是由C编译器来编译的。C编译器里不允许#include一个cu文件（不支持)

否则VS会报错

main.obj : error LNK2001: 无法解析的外部符号 _threadIdx

解决方案：

若要引用cu里的函数，在main.cpp里外部extern声明一下，让VS转为NVCC编译器处理。

3.2 GPU版HelloWorld

3.2.1 内存与显存之争

——老板，我要4G显存的显卡。

    ——亲，这是我们最新款显卡GT810，4G显存哦，只要399。您看那个GTX960

          才2G显存，跑游戏卡死了。

——好，我就买这个了。

买显卡，看显存大小，似乎已经成了电脑小白必带的标签。

显存是什么？显卡中最慢的外部存储体，焊在PCB板上，与GPU独立，想来几斤就有几斤。

☻DDRIII的内存，传输带宽是30~40GB\s，也就是说，每秒能访存30~40GB数据。

☻笔记本GDDR5的显存，带宽大概是80GB\s。

☻GTX960的带宽是148GB\s。

☻淘宝杂牌显卡，带宽大概是40~50GB\s。

☻GPU的Cache带宽大概是1.5TB\s，CPU的寄存器带宽大概是15TB\s。

☻GPU大概有50%的T周期浪费在访问显存上。作为GPU上最慢最庞大的存储体，显存带宽极为重要。

★★★结论：根据木桶原理，访存速度远远远远远小于并行处理速度，因而I/O瓶颈是GPU的最大问题。

★★★CPU不可以直接访显存，仅__device__、__global__函数有权访问，所以需要内存显存数据相互转移。

3.2.2 内存与显存大挪移

一般来说，一个CUDA程序必然少不了以下三步：

☻cudaMalloc：创建新的动态显存堆

☻cudaMemcpy：将主机(Host)内存复制到设备(Device)显存

☻显存处理完之后，cudaMemcpy：设备(Device)显存复制回主机(Host)内存，释放显存cudaFree

其中第三步最容易遗忘。要知道，CPU最后是无法使用显存中的数据的。

3.2.3 GPU设备选择

CUDA并行有三大境界

☻多台机器间并行

☻多路显卡间并行

☻单卡多线程并行

一般单机中需要操控多路显卡交火，比如AlexNet就是典型的双路交火结构。

不像孤岛危机之类的游戏那样，CUDA需要手动检测、切换设备，这是一切计算的开头准备工作。

下面的InitCUDA函数演示了单卡计算的设备选择：

#include "cuda_runtime.h"
#include "cstdio"
bool InitCUDA()
{
    int cnt, i;
    cudaGetDeviceCount(&cnt);
    if (!cnt)
    {
        printf("No GPU Device\n");
        return false;
    }
    for (i = 0; i < cnt; i++)
    {
        cudaDeviceProp device;
        if (cudaGetDeviceProperties(&device, i)==cudaSuccess)
            if (device.major >= 1) break;
    }
    if (i == cnt)
    {
        printf("GPU's Computing Capability must higher than 1.0\n");
        return false;
    }
    cudaSetDevice(i);
    return true;
}

只有计算能力>=1.0(8800GT以上)的N卡才能做CUDA计算，默认选择符合条件的第一块卡。

3.2.4 HelloWorld

GPU版HelloWorld，主要目的是演示CUDA基本程序框架：

☻将HelloWorld复制进显存

☻让GPU完成strcpy函数

☻将显存中的HelloWorld转回内存，并且打印

整体代码:

/****kernel.cu****/
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
__global__ void cudaStrcpy(char *des, char *src)  /*内核函数*/
{
    while ((*src) != '\0') *des++ = *src++;
    *des = '\0';
}


/****gpu_helloworld.cu****/
#include "device.cu"
#include "kernel.cu"
#include "cstring"
void helloworld(char *str1, char *str2)
{
    InitCUDA();
    char *dev_str1=0, *dev_str2=0;
    int size = strlen(str1) + 1;
    cudaMalloc((void**)&dev_str1, size);   /*cuda系函数必须放在cu文件里*/
    cudaMalloc((void**)&dev_str2, size);
    cudaMemcpy(dev_str1, str1, size,cudaMemcpyHostToDevice);
    cudaStrcpy<<<1,1>>>(dev_str2, dev_str1); /*单线程块、单线程*/
    cudaMemcpy(str2, dev_str1, size, cudaMemcpyDeviceToHost);
}


/****main.cpp****/
#include "cstdio"
#include "cstring"
extern void helloworld(char *str1, char *str2);
int main()
{
    char src[] = "HelloWorld with CUDA";
    char *des = new char[strlen(src)+1];
    helloworld(src, des);
    printf("%s\n", des);
}

 3.3 GPU版向量加法

向量加法是CUDA 7.0在VS中提供的样例模板，演示了并行算法的经典trick：循环消除。

利用单个线程块中，多个线程并发执行，来消除循环。

时间复杂度估计，不能简单从O(n)迁移到O(1)，因为GPU同时并行量存在限制。

即便是Kepler架构中拥有192SP的SM阵列，理论同时并行量也不过是2048。

/****kernel.cu****/
__global__ void kernel_plus(int *a, int *b, int *c)
{
    int x = threadIdx.x;
    c[x] = a[x] + b[x];
}

/****gpu_vectoradd.cu****/
void vectorAdd(int *a, int *b, int *c,int size)
{
    if (!InitCUDA()) return;
    int *dev_a = 0, *dev_b = 0, *dev_c = 0;
    cudaMalloc((void**)&dev_a, size*sizeof(int));
    cudaMalloc((void**)&dev_b, size*sizeof(int));
    cudaMalloc((void**)&dev_c, size*sizeof(int));
    cudaMemcpy(dev_a, a, size*sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(dev_b, b, size*sizeof(int), cudaMemcpyHostToDevice);
    kernel_plus << <1, size >> >(dev_a, dev_b, dev_c);
    cudaMemcpy(c, dev_c, size*sizeof(int), cudaMemcpyDeviceToHost);
}

/****main.cpp****/
extern void vectorAdd(int *a, int *b, int *c, int size);
int main()
{
    int a[5] = { 1, 2, 3, 4, 5 }, b[5] = { 10, 20, 30, 40, 50 }, c[5] = { 0 };
    vectorAdd(a, b, c,5);
    printf("{1,2,3,4,5} + {10,20,30,40,50} = {%d,%d,%d,%d,%d}\n",
        c[0], c[1], c[2], c[3], c[4]);
}