CUDA编程-part1

CUDA启动！

并行计算

理论上可以先学一下openMP或者pthread，毕竟是C/C++编程而且涉及到并行技术。

不同的数据划分会影响计算的性能：

对于即将学习的GPU，其中比较核心的部分可以分为：

核函数
内存管理
线程管理
流

在进行逻辑层面对数据、架构更深理解之前，首先需要在领域层了解如何使用CUDA解决简单的问题。

相较于常规的CPU编程，我们只需要考虑内存和CPU的交互即可（实际上这个也不需要），但是对于GPU（同时包含各种并行结构，包括多线程、MIMD、SIMD、指令级并行，因此可以称为SIMT单指令多线程）编程，一般异构的下甚至会有多CPU和多GPU，而资源调度又是需要经过PCIE总线（先不说统一寻址），因此需要将内存从宿主调用到设备上。

对于GPU而言，GPU容量的指标有两个：

CUDA核心数量
内存大小

对于GPU性能，指标同样有两个：

峰值计算性能
内存带宽

相对于openmp和pthread操作，CUDA编程做了相当复杂的抽象从而是得对架构有更深的理解。

对于CUDA编程，由于其实现了抽象，因此需要了解的地方主要有两个地方：

通过层次结构在GPU中组织线程的方法
通过层次结构在GPU中访问内存的方法

这里有个主从关系需要知道：对于宿主的一个核函数1-1grid，grid1-n块，块1-n线程：

可以理解成block是个进程，也因此每个block之间是同步而且是共享内存的。

一般来说一个网格被分为二维的块而块分为三维的线程。

CUDA模型编程

一个标准的CUDA编程分为三步：

数据从CPU内存拷贝到GPU内存
调用核函数对存储在GPU内存中的数据进行操作
数据从GPU传回到CPU内存中

依次说明。首先是内存管理的过程：

内存管理

cudaMalloc负责分配一定的线性内存，并以devPtr的形式返回指向分配内存的指针。整体的过程如图所示：

从示例代码可以看出了进行cudaMalloc时需要使用二维指针，对于线性的内存这里确实存疑，但是代码是这样写的（关键是地址）：

template<class T>
static __inline__ __host__ cudaError_t cudaMalloc(
  T      **devPtr,
  size_t   size
)
{
  return ::cudaMalloc((void**)(void*)devPtr, size);
}
// part of code
  memset(res_h,0,nByte);
  memset(res_from_gpu_h,0,nByte);

  float *a_d,*b_d,*res_d;
  CHECK(cudaMalloc((float**)&a_d,nByte));
  CHECK(cudaMalloc((float**)&b_d,nByte));
  CHECK(cudaMalloc((float**)&res_d,nByte));

由于设备内存和主机内存并不同，因此完全不能简单的进行：gpuRef=d_C这样的操作。

线程管理

首先大致说明一下核函数是如何运行的：核函数长得很像正常的函数：

kernel_name«<grid,block»>(argument list);

function_name(argument list);

但是与正常不同的是grid和block的定义，也就是接下来需要了解的内容：

CUDA将结构划分为两层，由线程块和线程块网格构成。在调用核函数时可以指定对应的grid和block，此外，grid和block并不是如图所示的二维结构，在代码进行定义的时候使用的是dim3.x\.y\.z每个均为uint。

在调用核函数时，除了使用简单的常量表示，更常见的是用dim3形式来创建并使用：

dim3 block(3);
dim3 grid((nElem+block.x-1)/block.x)

这里在定义时，grid是作为block的倍数的，至于为什么设计，后面再说。

另外一个有趣的点是使用两层的原因，其实也很简单，迎合GPU的架构以及更加方便的调用资源。

最后，要注意的是核函数的运行是异步的，一旦启动会立即执行。因此在有些场合，需要使用一些同步手段：

cudaError_t cudaDeviceSynchronize(void);

这是一个显式的用于同步的函数，当然在实际使用中，资源抢占也会导致不得不同步：

cudaError_t cudaMemcpy(void* dst,const void * src,
  size_t count,cudaMemcpyKind kind);

了解核函数的相关概念之后就可以尝试编写核函数了，与正常的函数类似，核函数多加了一个__global__（限定词），具体如下：

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

除此之外还有以下限制：

只能访问设备内存
必须返回void
不支持可变数量参数
不支持静态变量
显示异步行为

一个简单的例子：

void sumArraysOnHost(float *A, float *B, float *C, const int N) {
  for (int i = 0; i < N; i++)
    C[i] = A[i] + B[i];
}

__global__ void sumArraysOnGPU(float *A, float *B, float *C) {
  int i = threadIdx.x;
  C[i] = A[i] + B[i];
}

上述是一个串行和并行的代码，倒确实不难理解。上面都是关于核函数编写的一些技术型问题，实际上还需要进行一些技巧上的了解：

验证正确性

核函数的验证很简单，与非并行的比较即可，在编写中添加这些内容可以很方便的实现找到问题所在。详细代码查看代码3:

/*
* https://github.com/Tony-Tan/CUDA_Freshman
* 3_sum_arrays
*/
#include <cuda_runtime.h>
#include <stdio.h>
#include "tool.h"


void sumArrays(float * a,float * b,float * res,const int size)
{
  for(int i=0;i<size;i+=4)
  {
    res[i]=a[i]+b[i];
    res[i+1]=a[i+1]+b[i+1];
    res[i+2]=a[i+2]+b[i+2];
    res[i+3]=a[i+3]+b[i+3];
  }
}
__global__ void sumArraysGPU(float*a,float*b,float*res)
{
  int i=threadIdx.x;
  res[i]=a[i]+b[i];
}
int main(int argc,char **argv)
{
  int dev = 0;
  cudaSetDevice(dev);

  int nElem=1024;
  printf("Vector size:%d\n",nElem);
  int nByte=sizeof(float)*nElem;
  float *a_h=(float*)malloc(nByte);
  float *b_h=(float*)malloc(nByte);
  float *res_h=(float*)malloc(nByte);
  float *res_from_gpu_h=(float*)malloc(nByte);
  memset(res_h,0,nByte);
  memset(res_from_gpu_h,0,nByte);

  float *a_d,*b_d,*res_d;
  CHECK(cudaMalloc((float**)&a_d,nByte));
  CHECK(cudaMalloc((float**)&b_d,nByte));
  CHECK(cudaMalloc((float**)&res_d,nByte));

  initialData(a_h,nElem);
  initialData(b_h,nElem);

  CHECK(cudaMemcpy(a_d,a_h,nByte,cudaMemcpyHostToDevice));
  CHECK(cudaMemcpy(b_d,b_h,nByte,cudaMemcpyHostToDevice));

  dim3 block(nElem);
  dim3 grid(nElem/block.x);
  sumArraysGPU<<<grid,block>>>(a_d,b_d,res_d);
  printf("Execution configuration<<<%d,%d>>>\n",block.x,grid.x);

  CHECK(cudaMemcpy(res_from_gpu_h,res_d,nByte,cudaMemcpyDeviceToHost));
  sumArrays(a_h,b_h,res_h,nElem);

  checkResult(res_h,res_from_gpu_h,nElem);
  cudaFree(a_d);
  cudaFree(b_d);
  cudaFree(res_d);

  free(a_h);
  free(b_h);
  free(res_h);
  free(res_from_gpu_h);

  return 0;
}

#define CHECK(call)\
{\
  const cudaError_t error=call;\
  if(error!=cudaSuccess)\
  {\
      printf("ERROR: %s:%d,",__FILE__,__LINE__);\
      printf("code:%d,reason:%s\n",error,cudaGetErrorString(error));\
      exit(1);\
  }\
}

此外考虑到代码优化过程和最终结果，知道最后优化到什么程度也是非常重要的。因此需要使用计时器来获取代码的计算时间等。

在博客中使用了nvprof来获取计算时间（其中也包含了各个函数的运行时间，因此较为详细，可惜不支持8.0以上计算能力的设备）

组织并行线程

上面讲解了基本的CUDA并行编程模型，接下来就是如何组织线程实现并行任务。

首先是需要知道如何实现线程和块的索引：

上图相对来说比较详细和清楚，但是可能存在一些细节上的问题。对于行列索引，有全局索引和局部索引之分。

对于计算机的架构，显然多线程单指令是没有意义的，因此实际上进行CUDA并行计算的时候是将数据进行并行的。要做的就是使用上面的索引来索引数据从而实现数据并行的个功能。

一个矩阵加法的例子，尽管是矩阵加法，但是数据实际上是线性存储的，因此为了实现更加高效的计算，最好是将一行的数据放在一个block上，利用空间连续性提高性能，下面的代码提供了三种调用方式（二维线程块组成二维网格，一维线程块组成一维网格，一维线程块组成二维网格）：

我们首先不具体实现三种调用方式，在进行索引的时候，对于一个二维的示例而言，需要管理三种索引：

线程和块索引
矩阵中给定点的坐标
全局线性内存中的偏移量

遵循上面的原则可以使得编程时减少错误。首先是映射：

对于全局索引：

$$ix=threadIdx.x + blockIdx.x \times blockDim.x$$

$$iy=threadIdx.y + blockIdx.y \times blockDim.y$$

完成块、线程到矩阵的索引之后，接着需要把矩阵坐标映射到全局内存的索引/存储单元上，这个就简单一点了：

$$idx = iy * nx + ix$$


__global__ void sumMatrix(float * MatA,float * MatB,float * MatC,int nx,int ny)
{
    int ix=threadIdx.x+blockDim.x*blockIdx.x;
    int iy=threadIdx.y+blockDim.y*blockIdx.y;
    int idx=ix+iy*ny;
    if (ix<nx && iy<ny)
    {
      MatC[idx]=MatA[idx]+MatB[idx];
    }
}

int main(int argc,char** argv)
{
  printf("strating...\n");
  initDevice(0);
  int nx=1<<12;
  int ny=1<<12;
  int nxy=nx*ny;
  int nBytes=nxy*sizeof(float);

  //Malloc
  float* A_host=(float*)malloc(nBytes);
  float* B_host=(float*)malloc(nBytes);
  float* C_host=(float*)malloc(nBytes);
  float* C_from_gpu=(float*)malloc(nBytes);
  initialData(A_host,nxy);
  initialData(B_host,nxy);

  //cudaMalloc
  float *A_dev=NULL;
  float *B_dev=NULL;
  float *C_dev=NULL;
  CHECK(cudaMalloc((void**)&A_dev,nBytes));
  CHECK(cudaMalloc((void**)&B_dev,nBytes));
  CHECK(cudaMalloc((void**)&C_dev,nBytes));


  CHECK(cudaMemcpy(A_dev,A_host,nBytes,cudaMemcpyHostToDevice));
  CHECK(cudaMemcpy(B_dev,B_host,nBytes,cudaMemcpyHostToDevice));

  int dimx=32;
  int dimy=32;

  // cpu compute
  cudaMemcpy(C_from_gpu,C_dev,nBytes,cudaMemcpyDeviceToHost);
  double iStart=cpuSecond();
  sumMatrix2D_CPU(A_host,B_host,C_host,nx,ny);
  double iElaps=cpuSecond()-iStart;
  printf("CPU Execution Time elapsed %f sec\n",iElaps);

  // 2d block and 2d grid
  dim3 block_0(dimx,dimy);
  dim3 grid_0((nx-1)/block_0.x+1,(ny-1)/block_0.y+1);
  iStart=cpuSecond();
  sumMatrix<<<grid_0,block_0>>>(A_dev,B_dev,C_dev,nx,ny);
  CHECK(cudaDeviceSynchronize());
  iElaps=cpuSecond()-iStart;
  printf("GPU Execution configuration<<<(%d,%d),(%d,%d)>>> Time elapsed %f sec\n",
        grid_0.x,grid_0.y,block_0.x,block_0.y,iElaps);
  CHECK(cudaMemcpy(C_from_gpu,C_dev,nBytes,cudaMemcpyDeviceToHost));
  checkResult(C_host,C_from_gpu,nxy);
  // 1d block and 1d grid
  dimx=32;
  dim3 block_1(dimx);
  dim3 grid_1((nxy-1)/block_1.x+1);
  iStart=cpuSecond();
  sumMatrix<<<grid_1,block_1>>>(A_dev,B_dev,C_dev,nx*ny ,1);
  CHECK(cudaDeviceSynchronize());
  iElaps=cpuSecond()-iStart;
  printf("GPU Execution configuration<<<(%d,%d),(%d,%d)>>> Time elapsed %f sec\n",
        grid_1.x,grid_1.y,block_1.x,block_1.y,iElaps);
  CHECK(cudaMemcpy(C_from_gpu,C_dev,nBytes,cudaMemcpyDeviceToHost));
  checkResult(C_host,C_from_gpu,nxy);
  // 2d block and 1d grid
  dimx=32;
  dim3 block_2(dimx);
  dim3 grid_2((nx-1)/block_2.x+1,ny);
  iStart=cpuSecond();
  sumMatrix<<<grid_2,block_2>>>(A_dev,B_dev,C_dev,nx,ny);
  CHECK(cudaDeviceSynchronize());
  iElaps=cpuSecond()-iStart;
  printf("GPU Execution configuration<<<(%d,%d),(%d,%d)>>> Time elapsed %f sec\n",
        grid_2.x,grid_2.y,block_2.x,block_2.y,iElaps);
  CHECK(cudaMemcpy(C_from_gpu,C_dev,nBytes,cudaMemcpyDeviceToHost));
  checkResult(C_host,C_from_gpu,nxy);


  cudaFree(A_dev);
  cudaFree(B_dev);
  cudaFree(C_dev);
  free(A_host);
  free(B_host);
  free(C_host);
  free(C_from_gpu);
  cudaDeviceReset();
  return 0;
}

另一种验证的手段就是将核函数调用的block和grid设置为1，1从而达到串行的效果。

上述是目前阶段的一整个流程，随着不断复杂，整体也会不断变化。

GPU设备信息

以往的设备信息都是通过nvidia-smi获得的，考虑到应用场景，往往是够的，但是在CUDA编程下需要知道更多的信息：

nvidia-smi -q -i 0

nvidia-smi -q -i 0 -d MEMORY

但是感觉对于CUDA编程，还是使用教程中的程序更好一点。

CUDA执行模型基础

CUDA整体可以理解成是多个SM（流多处理器）组成的，每个SM支持上百个线程并发执行。这是第一个概念。但是在执行的时候，并不是将其直接全部并行使用，而是将其分为32个线程一组的线程束，而且这些线程是同步的（相较于SIMD对多数据进行相同的指令，GPU的SIMT允许同一线程束的多个线程独立执行，从相同的指令地址开始执行但是可能会有不同的行为（后面提到的分支）），这32个线程只能跑一个程序而且提前跑完或者不跑都只能等其他线程完成。简而言之整体的架构：

上图展示除了逻辑和物理设备的对应关系，其中一点要注意的是块内是可以同步的，即一个SM内数据是同步的，不会出现数据冲突问题，但是对于不同块之间是异步的，因此可能会导致冲突问题。完成上面对基本架构的理解之后就是整体架构的详细了解了。

由于一个SM上可以同时具有多个线程束，因此就回到了CPU调度线程时候可能遇到的问题，当有些线程束因为一些理由闲置的时候，SM可以立即从同一SM的其他线程束调度过来并运行（并发的线程束间并没有切换开销，因为硬件资源已经分配到了SM上的所有线程和块中，所以最新调度的线程束状态已经存储在SM上了）

Fermi架构

最早的架构，简单但是够用了：

这里整体的架构就不多说明了，我们同样要关注的是每个SM中的内容：

CUDA核自然是再熟悉不过了，但是这次的核心在于橙色的东西，这个橙色的主要有四个，分别是两组束调度器和指令分配单元，也因此可以知道，实际上对于线程束上还需要分为两组进行调度。

除此之外，由于一个线程块内的切换代价基本是没有的（数据共享），因此可以实现较快的调度和并发从而充分利用GPU，如图所示的灵活任务调度：

Kepler架构

接下来是Fermi架构的升级版，相较于上一代架构，Kepler架构强化了SM并使用了动态并行以及Hyper-Q技术（此外还增加了技术参数，堆料）

每个MX变成SMX，此外，使得内核可以启动内核（可以进行递归操作）

Hyper-Q指的是GPU和CPU间的同步硬件链接，在GPU执行的时候CPU可以做更多的事情，Fermi和Kelper的区别如下：

除了CUDA core的增加，之前单通道调度整体，现在可以多个硬件队列调度grid了。

然而上面仅仅是一些简单的GPU架构理解，后面就可以自己看一些GPU的架构和相关技术参数了。由于不同的架构、参数会带来运行结果的差异，需要根据硬件的相关信息进行参数的调整以及优化（也是一个核心内容）。举个例子，在4核8线程的CPU上，make -j8的速度比make -j4快。

对CUDA编程的理解了解之前可能认为就是使得线程、线程块数量尽可能匹配结构拥有的数量，但是实际上在复杂场景下，还得考虑由于计算差异导致的性能变化，举个例子，一个运算中的sin、cos计算较多，但是考虑到不同架构之间的SFU（特殊功能单元）数量不同，运行的效果和结果也会存在差异。

因此在优化时需要考虑如下内容：

代码的时间、空间复杂度
特殊指令的使用情况
函数调用的频率和持续时间

在实现高性能并行的时候也需要注意程序的健壮性（并行是这样的），而优化的核心在于充分理解硬件架构、测试软件等从而实现最佳的性能。

具体的调试和测试就需要通过nvvp、nvprof等手段进行测试了。

线程束运行本质

前面已经知道了线程束是由dim3指定的，因此可以使用三维坐标指定对应的线程束，然而实际上，线程在硬件层面上还是一维的架构。因此对于块中的每个线程，都可以使用内置变量threadIdx和blockDim来计算：

$$threadIdx.y * blockDim.x + threadIdx.x$$

$$threadIdx.z * blockDim.y * blockDim.x + threadIdx.y * blockDim.x + threadIdx.x$$

由于前面提到了线程会被分成32个一组，也因此如果无法组成32个一组，就会导致多余出来的线程成组并占用独立资源。（举个例子：80个软件线程要占用96个硬件线程）

通常情况下，线程束在执行时会被分配到相同的指令从而处理自己的私有数据，但是在执行指令时，我们有时候会使用到条件语句，即一个线程束中部分运行if而另一部分运行else。

对于CPU而言，由于其是流水线结构，在运行时采用分支预测可以很好的解决if的判断等待问题，但是问题在于，GPU的设计并没有那么复杂，无法实现这种分支预测功能。在核函数中同一个块不同线程运行时面对不同的if所带来的分支运行就是线程束的分化

但是前面提到了，一个指令周期，一个线程束执行的都是相同的指令，但是上面的分化似乎是相悖的。为了解决这个问题，使用了非常简单的方案：等。在别的指令执行if时，别的等待，同时在执行else时另一组等待，因此条件分支越多，性能下降越多。（注意，这里研究的是一个线程束内部的问题，即一个SM要处理的内容，不同线程束之间是不会遇到问题）

显然，这种问题的解决方式也比较简单，就是把线程束的分化变成线程束间的分化，下面是一个例子：

__global__ void mathKernel1(float *c)
{
	int tid = blockIdx.x* blockDim.x + threadIdx.x;

	float a = 0.0;
	float b = 0.0;
	if (tid % 2 == 0)
	{
		a = 100.0f;
	}
	else
	{
		b = 200.0f;
	}
	c[tid] = a + b;
}

__global__ void mathKernel2(float *c)
{
	int tid = blockIdx.x* blockDim.x + threadIdx.x;
	float a = 0.0;
	float b = 0.0;
	if ((tid/warpSize) % 2 == 0)
	{
		a = 100.0f;
	}
	else
	{
		b = 200.0f;
	}
	c[tid] = a + b;
}

详细代码参考一下8_divergence.cu，除了上面采用的手段来减少分支的代价以外，还使用了warmup来使得其之后运行速度正常。而且nvcc实际上默认带有分支预测功能，因此在一定程度上还是做了优化的。我们使用

sudo nvprof –metrics branch_efficiency ./divergence

来检测分支性能，比较幽默地是看来我的设备是没有分支优化的。不过后面应该还会有更加详细的说明。不过从分支效率不是50%可以看出来设备还是做了优化的。

然而，前面提到了部分线程会不运行而另外一些会运行，这时候我们把观察的目标放到线程束内部，除了从计算资源的角度理解线程调度，在学习计算机组成原理的时候还考虑了资源分配相关的线程调度，对于并行计算更加重要，接下来就是对线程束资源分配的相关内容。

前面知道了通过分支优化或者一些手段可以尽可能减少分支等待，但是在有些情况下，无论如何最后都会遇到不得不等待的分支，即存在一些线程束无法执行。对于这些无法执行的线程束，往往是两种可能：就绪态、等待资源态。这里有点像计算机组成原理中的资源调度问题，而线程束能否激活并使用，取决于：程序计数器、寄存器、共享内存因素。简而言之，核函数要求的资源越少-内存越大，同时容纳的线程束就越多。

也因此会出现当核函数较大时，一个SM的资源无法处理一个完整块从而无法启动程序。这时候就需要对架构和代码更深刻的理解和优化了。

一个早期的参数列表，实际上现在随便一个设备的计算能力都爆炒了，写这篇文章时候我的笔记本MX250的计算能力已经是6.0了。如果说资源满足计算进行的条件，接下来就是等待运行了，这里将其称之为活跃的块，活跃的块所包含的线程束就是活跃的线程束。但是一个块内的线程束可以进一步被分为3种类型：

选定的线程束
阻塞的线程束
符合条件的线程束

一个SM上的线程束在每个周期都会选择活跃的线程束然后调度到执行单元上运行。但是有点好玩的地方，一个块上既然数据资源调度好了，为什么还会出现阻塞，实际上主要存在两个原因：

32个CUDA核心可以用于执行
当前指令中所有的参数都就绪

尽管线程束的分配是由设备实现并且由于上下文调度对资源要求较低，但是资源分配往往限制线程束的数量，因此提高性能的另一个核心是尽可能地使得资源满足尽可能多地线程束使其活跃。

其次是利用率的问题，前面提到了需要使得尽可能多地线程束运行。这里引入一个新概念：指令发出和完成之间的时钟周期被称为指令延迟。相较于CPU上执行C语言，CUDA对于指令延迟更加关注，因为CPU核心是为了同时最小化延迟1~2个线程而设计的，但是GPU上存在大量的线程，因此GPU的指令延迟可以被线程束计算隐藏，从而使得计算效率更高。

这里引入一个队列理论的定理：利特尔法则：

$$所需线程束数量=延迟\times 吞吐量$$

一个简单的例子，指令平均延迟周期为5，为了每个周期跑6个线程束，至少需要准备30个未完成的线程束。

然而，单单讨论延迟和吞吐量还是有点以偏概全了，实际上也有复杂的因素在里面：

延迟不必多说，取决于指令的计算性能，吞吐量则取决于内存频率和带宽（计算方式也很简单：吞吐量=带宽*频率）得到带宽之后，带宽*指令延迟即可得到并行的io量。

因此，举个例子，如果一个线程需要4B的数据，那么为了延迟隐藏，就需要准备：

$$74KB/4B=18500 线程数$$

根据上述的例子可以知道，延迟隐藏取决于活跃线程束的数量，这一数量受资源的约束。因此要想实现性能最高的并行，需要控制好两者之间的关系。

接下来是GPU的占用率，类似于CPU的多核处理占用率，GPU的占用率为：

$$占用率=\frac{活跃线程束数量}{最大线程束数量}$$

获得SM中最大的线程束数量可以参考官方的simpleDeviceQuery.cu实现：

#include <cuda_runtime.h>
#include <stdio.h>

int main(int argc,char** argv)
{
    printf("%s Starting ...\n",argv[0]);
    int deviceCount = 0;
    cudaError_t error_id = cudaGetDeviceCount(&deviceCount);
    if(error_id!=cudaSuccess)
    {
        printf("cudaGetDeviceCount returned %d\n ->%s\n",
              (int)error_id,cudaGetErrorString(error_id));
        printf("Result = FAIL\n");
        exit(EXIT_FAILURE);
    }
    if(deviceCount==0)
    {
        printf("There are no available device(s) that support CUDA\n");
    }
    else
    {
        printf("Detected %d CUDA Capable device(s)\n",deviceCount);
    }
    int dev=0,driverVersion=0,runtimeVersion=0;
    cudaSetDevice(dev);
    cudaDeviceProp deviceProp;
    cudaGetDeviceProperties(&deviceProp,dev);
    printf("Device %d:\"%s\"\n",dev,deviceProp.name);
    cudaDriverGetVersion(&driverVersion);
    cudaRuntimeGetVersion(&runtimeVersion);
    printf("  CUDA Driver Version / Runtime Version         %d.%d  /  %d.%d\n",
        driverVersion/1000,(driverVersion%100)/10,
        runtimeVersion/1000,(runtimeVersion%100)/10);
    printf("  CUDA Capability Major/Minor version number:   %d.%d\n",
        deviceProp.major,deviceProp.minor);
    printf("  Total amount of global memory:                %.2f GBytes (%llu bytes)\n",
            (float)deviceProp.totalGlobalMem/pow(1024.0,3),deviceProp.totalGlobalMem);
    printf("  GPU Clock rate:                               %.0f MHz (%0.2f GHz)\n",
            deviceProp.clockRate*1e-3f,deviceProp.clockRate*1e-6f);
    printf("  Memory Bus width:                             %d-bits\n",
            deviceProp.memoryBusWidth);
    if (deviceProp.l2CacheSize)
    {
        printf("  L2 Cache Size:                            	%d bytes\n",
                deviceProp.l2CacheSize);
    }
    printf("  Max Texture Dimension Size (x,y,z)            1D=(%d),2D=(%d,%d),3D=(%d,%d,%d)\n",
            deviceProp.maxTexture1D,deviceProp.maxTexture2D[0],deviceProp.maxTexture2D[1]
            ,deviceProp.maxTexture3D[0],deviceProp.maxTexture3D[1],deviceProp.maxTexture3D[2]);
    printf("  Max Layered Texture Size (dim) x layers       1D=(%d) x %d,2D=(%d,%d) x %d\n",
            deviceProp.maxTexture1DLayered[0],deviceProp.maxTexture1DLayered[1],
            deviceProp.maxTexture2DLayered[0],deviceProp.maxTexture2DLayered[1],
            deviceProp.maxTexture2DLayered[2]);
    printf("  Total amount of constant memory               %lu bytes\n",
            deviceProp.totalConstMem);
    printf("  Total amount of shared memory per block:      %lu bytes\n",
            deviceProp.sharedMemPerBlock);
    printf("  Total number of registers available per block:%d\n",
            deviceProp.regsPerBlock);
    printf("  Wrap size:                                    %d\n",deviceProp.warpSize);
    printf("  Maximun number of thread per multiprocesser:  %d\n",
            deviceProp.maxThreadsPerMultiProcessor);
    printf("  Maximun number of thread per block:           %d\n",
            deviceProp.maxThreadsPerBlock);
    printf("  Maximun size of each dimension of a block:    %d x %d x %d\n",
            deviceProp.maxThreadsDim[0],deviceProp.maxThreadsDim[1],deviceProp.maxThreadsDim[2]);
    printf("  Maximun size of each dimension of a grid:     %d x %d x %d\n",
            deviceProp.maxGridSize[0],
	    deviceProp.maxGridSize[1],
	    deviceProp.maxGridSize[2]);
    printf("  Maximu memory pitch                           %lu bytes\n",deviceProp.memPitch);
    printf("----------------------------------------------------------\n");
    printf("Number of multiprocessors:                      %d\n", deviceProp.multiProcessorCount);
    printf("Total amount of constant memory:                %4.2f KB\n",
	deviceProp.totalConstMem/1024.0);
    printf("Total amount of shared memory per block:        %4.2f KB\n",
     deviceProp.sharedMemPerBlock/1024.0);
    printf("Total number of registers available per block:  %d\n",
    deviceProp.regsPerBlock);
    printf("Warp size                                       %d\n", deviceProp.warpSize);
    printf("Maximum number of threads per block:            %d\n", deviceProp.maxThreadsPerBlock);
    printf("Maximum number of threads per multiprocessor:  %d\n",
	deviceProp.maxThreadsPerMultiProcessor);
    printf("Maximum number of warps per multiprocessor:     %d\n",
	deviceProp.maxThreadsPerMultiProcessor/32);
    return EXIT_SUCCESS;
   
}

既然是并发编程，那么另一个需要考虑的就是同步的问题了，不可能代码一直异步到最后的。当需要进行同步的时候执行：cudaDeviceSynchronize()即可。从前面的代码可以看到。

尽管同步就是一个简单的等待，但是可能会因此产生错误。线程块中的线程束以未定义的顺序运行，CUDA因此设置了一个函数来标记同步点：

__device__ void __syncthreads(void);

当上述函数被调用时，一个线程块内的所有线程都得等待线程块内所有线程到这个同步点，在此之前的内存访问将会在这个同步点之后对其他线程可见。但是这种强制线程空闲的操作可能会导致性能的下降。

此外，如果代码涉及到竞争的时候，很有可能由于互相需要资源而死锁。总而言之，不同块、不同块之间的线程不会同步，只有一个线程束中的会同步。为此，保证安全的方法最好就是在核函数结束执行全局同步点并进行下一个核函数。

避免分支分化

在完成上述内容之后就需要一个简单的例子来说明前面所讲述的内容，这里采用的例子是归约求和，相较于串行的循环求和，可以使用并发编程实现求和，尽管是简单的求和，依旧要考虑如何组织使得整体的效率最高，显然有且不仅有下面的两种组织方式：

当然，考虑到空间连续性的问题，我们自然会觉得3-19的性能更高一些，接下来就是详细的核函数编写过程了。首先是一个简单的串行的递归求和，我们基于此函数进行修改和并行的实现：

int recursiveReduce(int *data, int const size)
{
	// terminate check
	if (size == 1) return data[0];
	// renew the stride
	int const stride = size / 2;
	if (size % 2 == 1)
	{
		for (int i = 0; i < stride; i++)
		{
			data[i] += data[i + stride];
		}
		data[0] += data[size - 1];
	}
	else
	{
		for (int i = 0; i < stride; i++)
		{
			data[i] += data[i + stride];
		}
	}
	// call
	return recursiveReduce(data, stride);
}

虽然可以简单用循环来实现加法，但是对于并行计算问题中的一类问题：归约问题，就是仿照上面的形式，而且是并行算法中的一个关键运算。

接下来就是对相邻元素进行归约的算法：

__global__ void reduceNeighbored(int * g_idata,int * g_odata,unsigned int n) 
{
	//set thread ID
	unsigned int tid = threadIdx.x;
	//boundary check
	if (tid >= n) return;
	//convert global data pointer to the local pointer of this block
	int *idata = g_idata + blockIdx.x*blockDim.x;
	//in-place reduction in global memory
	for (int stride = 1; stride < blockDim.x; stride *= 2)
	{
		if ((tid % (2 * stride)) == 0)
		{
			idata[tid] += idata[tid + stride];
		}
		//synchronize within block
		__syncthreads();
	}
	//write result for this block to global mem
	if (tid == 0)
		g_odata[blockIdx.x] = idata[0];

}

其具体的实现倒也不复杂，首先是获取线程id，并获取全局数据指针。要注意的是这里用的都是一维的块和一维的线程，所以相对来说简单了好多。相较于前面的递归归约，在并行算法中就不需要进行递归了，只需要使用循环来不断减少步进即可。不过一个有趣的问题在于这段代码只用了偶数线程做操作，因此感觉性能会稍微弱一点。每次都有接近$1-\frac{1}{stride}$比例的线程在等待。所以这种设计是简单但是显然是拉胯的。接下来就是对这部分代码进行修正使其性能提升：


__global__ void reduceNeighboredLess(int * g_idata,int *g_odata,unsigned int n)
{
	unsigned int tid = threadIdx.x;
	unsigned idx = blockIdx.x*blockDim.x + threadIdx.x;
	// convert global data pointer to the local point of this block
	int *idata = g_idata + blockIdx.x*blockDim.x;
	if (idx > n)
		return;
	//in-place reduction in global memory
	for (int stride = 1; stride < blockDim.x; stride *= 2)
	{
		//convert tid into local array index
		int index = 2 * stride *tid;
		if (index < blockDim.x)
		{
			idata[index] += idata[index + stride];
		}
		__syncthreads();
	}
	//write result for this block to global men
	if (tid == 0)
		g_odata[blockIdx.x] = idata[0];
}

前面的操作和上一个的类似，但是接下来就有区别了，相较于线程1对多映射数据，这里把tid转化为对应的索引进行处理，也因此不存在分化问题，从而消除了分化的情况（但是它的代码实现还是挺优雅的可恶）。

尽管消除了偶数线程不工作的情况，但是上述代码仍然不是最优的，考虑这样一种情况，数据量很大而同时只能运行512个线程的块是，前8个线程束进行归约的时候剩下的8个线程束无事可做，而第二轮前4个线程束工作的时候还有12个线程束不工作。这样是消除了块间的分化，是没什么问题的。但是在最后5轮中，每一轮的线程总数小于线程束大小时，分化就又出现了（16个线程干活而另外16+480个线程束在休息）从而提高了延迟。

在前面，我们提到了程序的空间一致性，因此尽可能调度相邻的数据进行处理，然而由于一个SM上的数据是共享的，在满足能放下的情况下似乎没有必要考虑相邻配对所带来的性能提升，因此是否存在一种更加高效的并行方法来减少分化呢？还真有：

__global__ void reduceInterleaved(int * g_idata, int *g_odata, unsigned int n)
{
	unsigned int tid = threadIdx.x;
	unsigned idx = blockIdx.x*blockDim.x + threadIdx.x;
	// convert global data pointer to the local point of this block
	int *idata = g_idata + blockIdx.x*blockDim.x;
	if (idx >= n)
		return;
	//in-place reduction in global memory
	for (int stride = blockDim.x/2; stride >0; stride >>=1)
	{
		
		if (tid <stride)
		{
			idata[tid] += idata[tid + stride];
		}
		__syncthreads();
	}
	//write result for this block to global men
	if (tid == 0)
		g_odata[blockIdx.x] = idata[0];
}

这次是一个交错进行数据交互的过程，其他的与之前的类似，也就是依然会有上面的线程分化问题，但是这一次结果如下：

（不得不说现代编译器优化还是太夸张了，不把o3优化关掉gpu都打不过CPU的）

展开循环

展开循环是一个在CPU串行就十分有效的方法，通过如下操作展开循环并减少循环轮次：

for(int i=0;i < 100; i++){
    r += data[i]
}
// 展开后
for(int i=0;i < 100; i+=2){
    r += data[i];
    r += data[i+1];
}

但这种循环展开为什么会带来性能提升，这一点可能和编译器优化有关。由于在每个循环中，语句的读写是独立的，因此CPU可以同时发出内存操作。但是在每次判断继续循环条件时次数减半（减少分支预测了）https://lazybing.github.io/blog/2019/04/17/loop-unroll/ （详细看这篇文章）

对于CUDA并行编程，由于判断的减少可以带来分支分化的减少，这时候相比CPU编程时的编译器优化，就更加需要手动操作循环来降低分支的次数了。

对于归约求和的例子，对于每个核函数，其实只处理了一个数据块，那么是否能拆开，让一个线程处理两个数据块，答案是可以的，但是这里给出的例子倒是有点稍微不同，仅仅是在进行归约前先进行一次：

__global__ void reduceUnroll2(int * g_idata,int * g_odata,unsigned int n)
{
	//set thread ID
	unsigned int tid = threadIdx.x;
	unsigned int idx = blockDim.x*blockIdx.x*2+threadIdx.x;
	//boundary check
	if (tid >= n) return;
	//convert global data pointer to the
	int *idata = g_idata + blockIdx.x*blockDim.x*2;
	if(idx+blockDim.x<n)
	{
		g_idata[idx]+=g_idata[idx+blockDim.x];

	}
	__syncthreads();
	//in-place reduction in global memory
	for (int stride = blockDim.x/2; stride>0 ; stride >>=1)
	{
		if (tid <stride)
		{
			idata[tid] += idata[tid + stride];
		}
		//synchronize within block
		__syncthreads();
	}
	//write result for this block to global mem
	if (tid == 0)
		g_odata[blockIdx.x] = idata[0];

}

通过上述简单的修改，可以进一步的降低运行的时间。当然，既然可以两个两个展开，也就可以4、8个展开（但是对应的grid数量得依次减少2/4/8倍，虽然前面没有详细提到，不过现在看来block*grid应该和每个核函数要处理的线程块*数据量相同），最后的结果如下：

但是看得出来，后面还有个怪东西的运行效率更高，这就是前面提到的如果小于32个线程的时候，每次都要等待所带来的影响，既然有了循环展开手段，因此就可以利用这个技术进行处理了处理方式也很简单，对于最后的6个迭代直接一个大展开即可：

if(tid<32)
{
	volatile int *vmem = idata;
	vmem[tid]+=vmem[tid+32];
	vmem[tid]+=vmem[tid+16];
	vmem[tid]+=vmem[tid+8];
	vmem[tid]+=vmem[tid+4];
	vmem[tid]+=vmem[tid+2];
	vmem[tid]+=vmem[tid+1];
}

要注意的是这里出现了一个新的关键字：volatile，该关键字说明接下来的赋值，每次都会将vmem[tid]的值写回全局内存（如果不佳则可能使用由编译器优化而保存在寄存器上的缓存）

接下来，最后一步，由于在部分情况下（比如说我们的例子中），循环条件是可知的，对于Fermi、Kepler架构，由于线程块最大为1024为已知，因此根据已知的的线程块，而代码中归约核函数的循环迭代次数实际上和这个1024是有关的，因此实际上可以完全展开而不适用for循环，最终结果：

template <unsigned int iBlockSize>
__global__ void reduceCompleteUnroll(int * g_idata,int * g_odata,unsigned int n)
{
	//set thread ID
	unsigned int tid = threadIdx.x;
	unsigned int idx = blockDim.x*blockIdx.x*8+threadIdx.x;
	//boundary check
	if (tid >= n) return;
	//convert global data pointer to the
	int *idata = g_idata + blockIdx.x*blockDim.x*8;
    // 相当于每个block分配了8*idx个数量的数据，先对8个数据进行归约
	if(idx+7 * blockDim.x<n)
	{
		int a1=g_idata[idx];
		int a2=g_idata[idx+blockDim.x];
		int a3=g_idata[idx+2*blockDim.x];
		int a4=g_idata[idx+3*blockDim.x];
		int a5=g_idata[idx+4*blockDim.x];
		int a6=g_idata[idx+5*blockDim.x];
		int a7=g_idata[idx+6*blockDim.x];
		int a8=g_idata[idx+7*blockDim.x];
		g_idata[idx]=a1+a2+a3+a4+a5+a6+a7+a8;

	}
	__syncthreads();
	//in-place reduction in global memory
    /*
	if(blockDim.x>=1024 && tid <512)
		idata[tid]+=idata[tid+512];
	__syncthreads();
	if(blockDim.x>=512 && tid <256)
		idata[tid]+=idata[tid+256];
	__syncthreads();
	if(blockDim.x>=256 && tid <128)
		idata[tid]+=idata[tid+128];
	__syncthreads();
	if(blockDim.x>=128 && tid <64)
		idata[tid]+=idata[tid+64];
	__syncthreads();
	*/
    // 接下来只剩下1024个数据了，因此可以简单的进行划分，线程束数量逐步下降。
	if(iBlockSize>=1024 && tid <512)
		idata[tid]+=idata[tid+512];
	__syncthreads();
	if(iBlockSize>=512 && tid <256)
		idata[tid]+=idata[tid+256];
	__syncthreads();
	if(iBlockSize>=256 && tid <128)
		idata[tid]+=idata[tid+128];
	__syncthreads();
	if(iBlockSize>=128 && tid <64)
		idata[tid]+=idata[tid+64];
	__syncthreads();
	//write result for this block to global mem
	if(tid<32)
	{
		volatile int *vmem = idata;
		vmem[tid]+=vmem[tid+32];
		vmem[tid]+=vmem[tid+16];
		vmem[tid]+=vmem[tid+8];
		vmem[tid]+=vmem[tid+4];
		vmem[tid]+=vmem[tid+2];
		vmem[tid]+=vmem[tid+1];

	}

	if (tid == 0)
		g_odata[blockIdx.x] = idata[0];

}

这里的代码与上面的代码实际上还有一点不一样，加了一个模板函数（该特性在CUDA编程中同样支持），相较于使用blockDim.x做大小比较使用模板函数可以在编译的时候就决定是否进行该循环从而进一步实现性能的优化，至此，这个归并求和的优化就算到头了。

总而言之，CUDA编程可以通过避免分支分化和展开循环来初步增加性能。考虑到现在的内容还是整本书1/4的部分，因此后面的优化估计也会更加复杂了。

动态执行

前面的代码都是静态执行，很简单，就是用CPU启动核函数执行，这是在各种条件都获取的比较充分的情况下才能做到的，实际上很多地方为了扩展性会不得不采用更加复杂的实现。因此也就出现了动态执行有用的场合。

动态执行的理解很简单，用核函数执行核函数。

注意，前面没有用到的grid到这里估计要用到咯。

关于资源，父网格和子网格共享全局和常量的内存存储，但是局部和共享内存不同。

首先拿一段代码简要说明一下：

#include <cuda_runtime.h>
#include <stdio.h>
__global__ void nesthelloworld(int iSize,int iDepth)
{
    unsigned int tid=threadIdx.x;
    printf("depth : %d blockIdx: %d,threadIdx: %d\n",iDepth,blockIdx.x,threadIdx.x);
    if (iSize==1)
        return;
    int nthread=(iSize>>1);
    if (tid==0 && nthread>0)
    {
        nesthelloworld<<<1,nthread>>>(nthread,++iDepth);
        printf("-----------> nested execution depth: %d\n",iDepth);
    }

}

int main(int argc,char* argv[])
{
    int size=8;
    int block_x=4;
    dim3 block(block_x,1);
    dim3 grid((size-1)/block.x+1,1);
    nesthelloworld<<<grid,block>>>(size,0);
    cudaGetLastError();
    cudaDeviceReset();
    return 0;
}

相当简单的代码，前面的部分就不用详细说了，在代码内，每次进行线程数量为原来一半的一次调用并增加深度。（但是要注意一点，由于这个动态运行是需要一个动态运行库支持的，因此在编译时需要加上参数-lcudadevrt，此外，为了在输出时效果为同步进行的，还需要使用参数--relocatable-device-code true，当然，在正常运行时是不需要使用这个参数的），运行结果如下：

到这里又重新迷惑了。接下来依次解释前面就该迷惑的点了（毕竟抽象的层次太复杂了）。

关于block到底是什么。block中定义了线程的数量，也就是一个块中有多少个线程工作，上一节中就是一个block采用了1024个线程工作，同时一个block也是资源共享的一个单位
关于grid是什么。前面由于对grid的使用不是很清楚，大多数情况下用来扩充block的倍数从而实现更大的运算，实际上在这一节，一个grid可以看成是一个block的执行单元，一个block中的所有线程放在一个SM上运行。举个例子：grid(3,1)，block(1024,1)就是开3个SM，每个SM里面跑的block由1024个线程运行，实际上1024还要除以32才是真正运行的单元
那么接下来，为什么会出现上面的结果（这里着重讲的blockIdx都是0），当然也很简单，因为每个父网格里面执行的子网格（或者说嵌套执行的子网格）的grid都是1，所以只开了一个block。仅此而已。

既然提到了这个嵌套的操作，在进行CPU的归约求和的时候我们使用了嵌套的实现，那么是否可以将上面的归约求和的嵌套形式也应用到CUDA中，当然可以！而且也简单了好多：

__global__ void  gpuRecursiveReduce(int *g_idata, int *g_odata, unsigned int isize) {
    unsigned int tid = threadIdx.x;
    
    int *idata = g_idata + blockIdx.x*blockDim.x;
    int *odata = &g_odata[blockIdx.x];
    
    if (isize==2 && tid==0) {
        g_odata[blockIdx.x] = idata[0] + idata[1];
        return ;
    }
    
    int istride = isize >> 1;
    if(istride > 1 && tid < istride) {
        idata[tid] += idata[tid+istride];
    }
    __syncthreads();
    if(tid == 0){
        gpuRecursiveReduce<<<1, istride>>>(idata, odata, istride);
    	cudaDeviceSynchronize();
    }
    __syncthreads()
}
// updated:
 if(istride > 1 && tid < istride) {
        idata[tid] += idata[tid+istride];
    if(tid == 0){
        gpuRecursiveReduce<<<1, istride>>>(idata, odata, istride);
    }
}

但是结果其实相当不好，即使相比于cpu也慢了不少。

由于有2048个线程块，每个线程块执行8次递归，所以总共计算一下，创造了16384个子线程块，同时中间还涉及到了同步的相关操作，从而导致性能极低。但是考虑到子网格和父网格的内存是完全相同的，因此子网格启动前线程块内部的同步是不需要的。

但是性能还是很差，毕竟子网格启动还是需要一定时间的。所以接下来要考虑的是减少子网格的调用：

__global__ void  gpuRecursiveReduce(int *g_idata, int *g_odata, int istride, int const iDim) {
    unsigned int tid = threadIdx.x;
    int *idata = g_idata + blockIdx.x*blockDim.x;
    
    if (iStride == 1 && threadIdx.x == 0){
        g_odata[blockIdx.x] = idata[0] + idata[1];
        return ;
    }
    
	idata[idx] += idata[idx+iStride]
    if(tid == 0 && blockIdx.x==0){
        gpuRecursiveReduce<<<gridDim.x, iStride/2>>>(g_idata, g_odata, iStride/2, iDim);
    }
}

相比于上面的操作，这里的代码就更加简单了，而且也更加清楚了，这里的核心代码在于在创建子网格时的条件和调用的过程，可以看出只有第一个block的第一个线程才能创造子网格，而且块大小变为原来的一半。

这里不管是博客还是书中都解释的比较含糊，考虑到这里可能并不是重点，咱且不多详细说明，暂且认为相较于网格创建的开销，块创建的开销极大同时由于一些内存的问题导致了性能的提升，当然，尽管如此，这种实现的性能依然不是最佳的。

简而言之，如果需要使用动态并行方式实现并发编程，需要考虑因素就更多了：同步正确性、同步次数、内核嵌套层数、调用次数。整体来说是更加复杂的。

全局内存

其实对于什么语言，内存管理都是重要的一栏（虽然对于现代高度抽象的语言另说）

CPU内存部分的管理已经学的非常差不多了，接下来就是对CUDA内存的