CUDA 编程
GPU 与 CPU
- CPU：擅长流程控制和逻辑处理，不规则数据结构，不可预测存储结构，单线程程序，分支密集型算法
- GPU：擅长数据并行计算，规则数据结构，可预测存储模式
异构编程
现在的计算机体系架构中，要完成CUDA并行计算，单靠GPU一人之力是不能完成计算任务的，必须借助CPU来协同配合完成一次高性能的并行计算任务。
一般而言，并行部分在GPU上运行，串行部分在CPU运行，这就是异构计算。具体一点，异构计算的意思就是不同体系结构的处理器相互协作完成计算任务。CPU负责总体的程序流程，而GPU负责具体的计算任务，当GPU各个线程完成计算任务后，我们就将GPU那边计算得到的结果拷贝到CPU端，完成一次计算任务。
CUDA线程模型
CUDA的线程模型从小往大来总结就是：
- Thread：线程，并行的基本单位
- Thread Block：线程块，互相合作的线程组，线程块有如下几个特点：
- 允许彼此同步
- 可以通过共享内存快速交换数据
- 以1维、2维或3维组织
- Grid：一组线程块
- 以1维、2维组织
- 共享全局内存
Kernel：在GPU上执行的核心程序，这个kernel函数是运行在某个Grid上的。
理解kernel，必须要对kernel的线程层次结构有一个清晰的认识。首先GPU上很多并行化的轻量级线程。kernel在device上执行时实际上是启动很多线程，一个kernel所启动的所有线程称为一个网格（grid），同一个网格上的线程共享相同的全局内存空间，grid是线程结构的第一层次，而网格又可以分为很多线程块（block），一个线程块里面包含很多线程，这是第二个层次。线程两层组织结构如上图所示，这是一个gird和block均为2-dim的线程组织。grid和block都是定义为dim3类型的变量，dim3可以看成是包含三个无符号整数（x，y，z）成员的结构体变量，在定义时，缺省值初始化为1。因此grid和block可以灵活地定义为1-dim，2-dim以及3-dim结构，kernel调用时也必须通过执行配置<<<grid, block>>>来指定kernel所使用的网格维度和线程块维度。
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有非常严格的限制，也就限制了并行能力。
需要指出，每个SM包含的SP数量依据GPU架构而不同，Fermi架构GF100是32个，GF10X是48个，Kepler架构都是192个，Maxwell都是128个。
简而言之，SP是线程执行的硬件单位，SM中包含多个SP，一个GPU可以有多个SM（比如16个），最终一个GPU可能包含有上千个SP。这么多核心“同时运行”，速度可想而知，这个引号只是想表明实际上，软件逻辑上是所有SP是并行的，但是物理上并不是所有SP都能同时执行计算（比如我们只有8个SM却有1024个线程块需要调度处理），因为有些会处于挂起，就绪等其他状态，这有关GPU的线程调度。
每个线程由每个线程处理器（SP）执行
线程块由多核处理器（SM）执行
一个kernel其实由一个grid来执行，一个kernel一次只能在一个GPU上执行
block是软件概念，一个block只会由一个sm调度，程序员在开发时，通过设定block的属性，告诉GPU硬件，我有多少个线程，线程怎么组织。而具体怎么调度由sm的warps scheduler负责，block一旦被分配好SM，该block就会一直驻留在该SM中，直到执行结束。一个SM可以同时拥有多个blocks，但需要序列执行。
CUDA内存模型
CUDA中的内存模型分为以下几个层次：
- 每个线程都用自己的registers（寄存器）。
- 每个线程都有自己的local memory（局部内存）。
- 每个线程块内都有自己的shared memory（共享内存），所有线程块内的所有线程共享这段内存资源。
- 每个grid都有自己的global memory（全局内存），不同线程块的线程都可使用。
- 每个grid都有自己的constant memory（常量内存）和texture memory（纹理内存）），不同线程块的线程都可使用。
| 存储器 | 位置 | 拥有缓存 | 访问权限 | 变量生存周期 |
|---|---|---|---|---|
| register | GPU片内 | N/A | Device可读/写 | 与thread相同 |
| Local memory | 板载显存 | 无 | Device可读/写 | 与thread相同 |
| Shared memory | GPU片内 | N/A | Device可读/写 | 与block相同 |
| Constant memory | 板载显存 | 有 | Device可读，host要读写 | 可在程序中保持 |
| Texture memory | 板载显存 | 有 | Device可读，host要读写 | 可在程序中保持 |
| Global memory | 板载显存 | 无 | Device可读/写, host可读/写 | 可在程序中保持 |
| Host memory | Host内存 | 无 | host可读/写 | 可在程序中保持 |
| Pinned memory | Host内存 | 无 | host可读/写 | 可在程序中保持 |
线程访问这几类存储器的速度是register > local memory >shared memory > global memory。
CUDA编程模型
一个CUDA程序的可以分为两个部分: 在CPU上运行的Host程序；在GPU上运行的Device程序。两者拥有各自的存储器。GPU上运行的函数又被叫做kernel函数，通过__global__关键字声名。
__ device__设备端执行，设备端调用__global__设备端执行，主机端调用__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）
计算线程号
使用N个线程块，每一个线程块只有一个线程，即
dim3 dimGrid(N);
dim3 dimBlock(1);
此时的线程号的计算方式就是
threadId = blockIdx.x;
其中threadId的取值范围为0到N-1。对于这种情况，我们可以将其看作是一个列向量，列向量中的每一行对应一个线程块。列向量中每一行只有1个元素，对应一个线程。
使用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
使用一个线程块，该线程具有N个线程
dim3 dimGrid(1);
dim3 dimBlock(N);
此时线程号的计算方式为：
threadId = threadIdx.x;
其中threadId的范围是0到N-1，对于这种情况，可以看做是一个行向量，行向量中的每一个元素的每一个元素对应着一个线程。
使用M个线程块，每个线程块内含有N个线程
dim3 dimGrid(M);
dim3 dimBlock(N);
这种情况，可以把它想象成二维矩阵，矩阵的行与线程块对应，矩阵的列与线程编号对应，那线程号的计算方式为：
threadId = threadIdx.x + blcokIdx*blockDim.x;
上面其实就是把二维的索引空间转换为一维索引空间的过程。
使用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;
上述公式就是把线程和线程块的索引映射为图像像素坐标的计算方法。
