为什么需要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~40GBs，也就是说，每秒能访存30~40GB数据。
☻笔记本GDDR5的显存，带宽大概是80GBs。
☻GTX960的带宽是148GBs。
☻淘宝杂牌显卡，带宽大概是40~50GBs。
☻GPU的Cache带宽大概是1.5TBs，CPU的寄存器带宽大概是15TBs。
☻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 "); 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 "); 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) != ' ') *des++ = *src++; *des = ' '; } /****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 ", 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} ", c[0], c[1], c[2], c[3], c[4]); }