CUDA编程-part2

全局内存

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

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

CPU内存层次中L1和L2缓存都是不可编程的，具体的实现都是由硬件和系统决定的。但是对于CUDA编程，为了性能上的优化，CUDA内存模型直接提出了一堆可以用的编程内存类型，比较清楚的一幅图：

接下来就需要依次描述每个的作用了（挺折腾的）

寄存器：一旦涉及到寄存器，那就是最快的。寄存器被每个线程所私有，和核函数的生命周期相同。但是并非一个线程对应一个寄存器，在Fermi架构下一个线程最多可以用63个寄存器。当然也因此如果代码所需要的数据超过了寄存器数量，就会使用内存来代替，寄存器溢出从而降低性能。
共享内存：使用__shared__修饰符修饰的变量放在共享内存中。从图中可以看出共享内存在块内是共享的，而且性能较强，有点类似L1缓存但是可编程。也因此可以用于进行线程间的通信等操作。SM中的共享内存和L1缓存都使用的64KB的片上内存，是静态划分的，但是同样可以修改。
本地内存：符合存储在寄存器条件但是被溢出到的地方，通常会有如下变量：未知索引本地数组、较大的结构体或者数组、不满足核函数寄存器限定条件的变量。虽然叫作本地内存，但是从硬件上来看和全局内存是一个东西，而且具备较高的性能，可以理解成L1、L2缓存。
常量内存：利用__constant__修饰的变量，而且必须要全局空间和核函数之外进行声明，且只能占用64KB的常量内存，且线程只能读取数据。听起来比较死板，但是实际上有有用的用途，比如存储数学公式中的悉数
纹理内存：是一个特别优化的内存部分，在此不详细说明（在图形计算上有用）
全局内存：全局内存大但是慢，同样可以被所有的SM访问。可以是静态或者动态的，静态的使用__device__修饰（不同于CPU的static）。鉴于其全局而且可以被读写，因此很容易导致读写的不一致问题。

鉴于全局内存是使用比较频繁的部分，对其也有着较大的优化空间，对于全局内存的请求是通过内存事务进行的，内存事务的大小通常为32、64、128字节（有点像CPU的寻址，必须对齐），在线程束进行内存加载或者保存时，就会通过内存事务来调度数据，因此其需要满足的传输数量取决于：跨线程的内存地址分布和事务内存的对齐方式。

书中这里描述的比较硬核，实际上感觉这里想要说明的是内存的数据空间连续性使得调度频率高的数据尽可能比较接近这样的？

CUDA自然也有缓存，同样也不可编程，主要就是为了保证访存快而设计的。

接下来是一个静态声明全局变量的例子：

#include <cuda_runtime.h>
#include <stdio.h>
__device__ float devData;
__global__ void checkGlobalVariable()
{
    printf("Device: The value of the global variable is %f\n",devData);
    devData+=2.0;
}
int main()
{
    float value=3.14f;
    cudaMemcpyToSymbol(devData,&value,sizeof(float));
    printf("Host: copy %f to the global variable\n",value);
    checkGlobalVariable<<<1,1>>>();
    cudaMemcpyFromSymbol(&value,devData,sizeof(float));
    printf("Host: the value changed by the kernel to %f \n",value);
    cudaDeviceReset();
    return EXIT_SUCCESS;
}
/*
Host: copy 3.140000 to the global variable
Device: The value of the global variable is 3.140000
Host: the value changed by the kernel to 5.140000 
*/

非常简单的例子，核函数试图修改全局变量的值但是木大了。上述代码简单但是有几个比较核心的地方1、cudaMemcpyToSymbol是CUDA运行时的API所以可以直接用2、devData是一个标识符，不是全局内存的变量地址，不能使用&获取地址，获取地址需要通过cudaGetSymbolAddress函数操作3、核函数中的devData是一个全局内存的变量。

内存管理

在最CUDA编程最开始就提到了一些CUDA设备内存的调度函数，要注意的一点是由于内存的分配和释放代价较大，因此还需要尽可能地重利用设备内存。

尽管GPU的带宽非常夸张，但是实际上由于是通过PCIE通道进行的，因此带宽远远达不这个程度PCIe3 只有8G而PCIe4也只有16G

除了之前提到的诸如cudaMemcpy这样的函数，还有一些可以直接从主机端获取内存的函数可以调用。要注意的是，由于主机内存是可分页、虚拟的，因此可能存在奇怪的引用操作。因此CUDA的函数cudaMallocHost从可分页的内存中固定一部分来提供设备端使用，并且需要使用对应的函数进行释放cudaFreeHost：

上面的概念没有代码可能会有些混乱：

  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;
  // pine memory malloc
  CHECK(cudaMallocHost((float**)&a_d,nByte));
  CHECK(cudaMallocHost((float**)&b_d,nByte));
  CHECK(cudaMallocHost((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));

想要强调的就是数据通过CUDA的API调用而不是直接使用malloc之后转移（固定分页的内存同样需要进行数据转移）。

理论上设备和主机之间的内存是不能直接互相访问的，但是为了更加方便，这时候就出现了一个给两者通用的内存：零拷贝内存。其具有以下特点：

扩充设备内存
避免数据显式传输
提高PCIe传输率

本来以为零拷贝内存是一个简单小块的内容但是实际上还是比较复杂的。函数如下：

cudaError_t cudaHostAlloc()(void **pHost, size_t count, unsigned int flags);其中flags有四种：

cudaHostAllocDefalt该flag用途和cudaMallocHost一致
cudaHostAllocPortable该函数返回可以被CUDA上下文使用的固定内存（全局了）
cudaHostAllocWriteCombine可以在某些系统配置上实现更快地传输，但是不容易在大多数设备上读取，可以用于缓冲区
cudaHostAllocMapped零拷贝内存最明显的表示，实现主机和设备间映射的内存

由于零拷贝内存走的是PCIe通道，因此尽管有优势，但是如果频繁读写的话是会带来较大的通讯代价的。

这里理论上会有一个例子说明，但是占空间暂且不说了。零拷贝内存还有一个有用的地方在于一种异构的计算架构：集成架构。在该架构下CPU和GPU是在一个芯片的因此共享内存，这种情况下零拷贝内存性能更高（不走PCIe了）。但是问题在于现在的高性能计算基本上都是离散架构的，这里说一下图一乐。

但是无论如何优化，不可避免地一个问题就是代码复杂太多了，一个数据要从主机端根据数据特性、调用形式来进行优化，相当复杂，也因此在CUDA4.0引入了UVA统一虚拟寻址。相较于前面复杂的内存拷贝过程，UVA可以直接将cudaHostAlloc返回的指针直接传递给核函数进行使用（注意不是cudaMalloc函数）。

然而优化还没有停止，由于UVA只能应用于cudaHostAlloc这样的零拷贝内存，但是更多场合下我们似乎更经常使用cudaMalloc来实现内存的分配，因此，在CUDA6.0版本出现了统一内存寻址。统一内存寻址依赖于UVA但是又不同于UVA、UVA为系统所有处理器提供了一个单一的虚拟内存地址空间，但是无法自动将数据从物理位置进行转移，但是统一内存寻址做到了。其使用方式有两种：

__device__ __managed__ int y;
cudaError_t cudaMallocManaged(void **devPtr, size_t size, unsigned int flags=0)

类似cudaMalloc但是该函数可以自动实现数据传输和重复指针消除的功能，当然剩下的特性稍后在说明。

注意，这一部分详细说明了几种内存模式和对应的调用，同时涉及到了一些软件上的优化，后续会更加清楚地理解如何应用（虽然重要的好像主要是cudaMalloc和cudaHostAlloc）

内存访问模式

CUDA核函数的运行数据的访问最开始都是从全局内存的，且极容易收到内存带宽的限制，因此也是需要对内存访问进行优化的。对于CUDA而言，其特性前面提到特殊一点的就是指令必须以线程束为单位进行发布和执行，这也可以应用到存储操作上，在读取时，每个线程都提出了一个包含请求地址的单一内存访问请求，因此不同的内存分布可能会影响到整体的性能，具体的内容接下来开始看。

在CPU的内存层次结构我们知道不管中间怎么操作，最后都是会需要将数据从内存传输到缓存和寄存器中的：

如图所示，对全局内存的访问都是需要通过二级缓存的（取决于GPU的架构，也有一些会通过一级缓存），如果两级都通过了，那么得到的就是128B的内存事务，如果只使用了二级缓存，则得到32B的内存事务。这也不难理解，对于一级缓存，如果每个线程发起了4B的请求，那么合起来就是128B，正好与缓存行和设备内存段大小相契合（注意，在架构那一章我们看图可以知道二级缓存是所有SM共享的）。

与CPU访存类似，对于不是128B的数据访存是会产生浪费的，而且如果没有和线程对齐，则一次访问可能会带来更多的内存事务，下图所示的访存事务就导致了多次的访存。

所以最快的访问方式就是通过一级缓存之后进行读写了（注意前面那张图还有只读缓存和常量缓存，但是那两个需要显式指明访问地址，暂且不多说明，常见的还是一级缓存）

接下来就是对加载情况以及实际代码导致存在的问题进行详细举例了，但是这部分和CPU访存还是高度一致的，因此不用复杂的说明了。

但是有一些需要考虑的细节：在不经过一级缓存的时候得到的是32B的一个内存事务，则一次128B的随机访问如果正好是4个块的话性能是不会变化的。

接下来就是只读缓存了，这个功能本来是留给纹理内存加载使用的，但是对于3.5以上的GPU也能用来代替一级缓存了。由于只读缓存的粒度为32B，因此在分散读取的时候性能会略高于一级缓存。这里没有详细的说明了，说明不太好用。

关于全局内存的写入和访问类似，同样需要尽可能对齐使得为32或128的整数倍从而使得每次每次都由四段长事务或者一段长事务组成。

说到访存，我们在之前的例子里面用到了数组来存储和使用，如果使用了合理的映射关系那么可以实现非常好用而且高性能的加载，然而单单使用数组只是简单的情况，在进行实际应用的时候，我们还需要使用结构体数组或者反过来数组结构体。

woc这里书里面的命名太拉胯了实在不好区分，这里采用：数组的结构体（对应结构体数组SoA）和结构体的数组（对应数组结构体AoS）来区分。两者的内存排列如下：

至于这种访存任务下哪个性能高就不难理解了。

在了解了CUDA的内存访问机制之后就是如何应用了。在上一章我们使用了诸多手段优化指令吞吐量，这里要优化的就是内存访问了。前一章使用了展开技术来实现指令吞吐的优化，这里同样可以使用展开技术进行优化：

__global__ void sumArraysGPU(float*a,float*b,float*res,int offset,int n)
{
  //int i=threadIdx.x;
  int i=blockIdx.x*blockDim.x*4+threadIdx.x;
  int k=i+offset;
  if(k+3*blockDim.x<n)
  {
      res[i]=a[k]+b[k];
      res[i+blockDim.x]=a[k+blockDim.x]+b[k+blockDim.x];
      res[i+blockDim.x*2]=a[k+blockDim.x*2]+b[k+blockDim.x*2];
      res[i+blockDim.x*3]=a[k+blockDim.x*3]+b[k+blockDim.x*3];
  }

}

仅仅是较小的改变就能使得性能提高3倍，这主要是有由于该核函数是I/O密集的，且访存优先级较高，理论上来说，展开其实不影响访存的次数，只影响并发执行的数量，也因此带来了巨大的性能提升：

是否展开并没有改变加载和存储效率。从某个角度来看实际上还是对并行指令的优化，因为内存优化的核心其实还是对齐和使用正确的指令，也就是接下来这章节使用的例子：矩阵转置。

矩阵转置问题

对于转置问题，由于矩阵的存储依然是线性的，因此实际上的转置变换如下：

不难看出，转置前后必然有一个过程是无法实现对齐的内存事务的。接下来首先是一个直觉上的感知，以下那种实现的速度更快：

不显然，如果禁用了一级缓存，那么两者的性能相同，但是如果没有禁用，那么就是第二个快一点了，这里取一张图说明一下：

画成这样大概率就理解了，写的时候没有什么额外操作，但是在读的时候是有缓存的，所以如果一次读了128B的事务，在后面的使用过程中是仍然会用到后面的数据的，从而减少访存次数。

在完成矩阵转置之前首先测试一下根据上述判断仅仅通过简单的复制操作得到的结果：

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

通过上述代码可以大致得到我们计算得到的带宽的上下限。接下来就是矩阵转置算法的实现了，倒也并不复杂：

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

完成大致的搭建之后就是使用展开的方法继续进行优化了，代码如下，实现起来和前面的是一致的，但是同样要注意在展开之后修改块的数量：

__global__ void transformNaiveColUnroll(float * MatA,float * MatB,int nx,int ny)
{
    int ix=threadIdx.x+blockDim.x*blockIdx.x*4;
    int iy=threadIdx.y+blockDim.y*blockIdx.y;
    int idx_row=ix+iy*nx;
    int idx_col=ix*ny+iy;
    if (ix<nx && iy<ny)
    {
        MatB[idx_row]=MatA[idx_col];
        MatB[idx_row+1*blockDim.x]=MatA[idx_col+ny*1*blockDim.x];
        MatB[idx_row+2*blockDim.x]=MatA[idx_col+ny*2*blockDim.x];
        MatB[idx_row+3*blockDim.x]=MatA[idx_col+ny*3*blockDim.x];
    }
}

从结构的角度来看其实上面的优化已经到头了，但是实际上还有一些接近算法层面的优化──对角转置。

在启用核函数的时候块被分配的SM是由块ID确定的，但是由于不同SM执行的速度不同，有些还没运行的块就被分配到其他的SM上，从而使得本来连续的块ID变不连续了。虽然控制不了整体的调度顺序，但是可以用块坐标来进行数据的读取。

__global__ void transformNaiveRowDiagonal(float * MatA,float * MatB,int nx,int ny)
{
    int block_y=blockIdx.x;
    int block_x=(blockIdx.x+blockIdx.y)%gridDim.x;
    int ix=threadIdx.x+blockDim.x*block_x;
    int iy=threadIdx.y+blockDim.y*block_y;
    int idx_row=ix+iy*nx;
    int idx_col=ix*ny+iy;
    if (ix<nx && iy<ny)
    {
      MatB[idx_col]=MatA[idx_row];
    }
}

但是这种方法为什么有效？说起来就更加复杂了，这里的优化涉及到了DRAM的并行访问。

发送给全局内存的请求是有DRAM分区实现的，这里拿博客的解释来说明一下吧，对于DRAM的每一个分区，尽管读写看上去似乎没有什么区别，但是过多的对一个分区进行读写会导致一种排队的现象从而使得该区的负载较大，从而降低性能（一开始启动了4个块1，2，3，4但是访问的都是第一个DRAM分区大致这样的）。如果使用了对角坐标系，则对内存的访问会稍微均匀一点从而减少排队现象。

最后一种方案就是修改块的尺寸来增加并行性了，这就属于一种简单的微调了，整体也并不复杂：

图中的第二个就是瘦块，增加存储在线程块中的连续元素的数量来提高性能。

到这里整体的矩阵转置就全部优化结束了。（真复杂啊）

当然还有一个使用统一内存的方式实现的矩阵加法，但是并不复杂，是对API的学习（毕竟这个统一内存不需要考虑拷贝来拷贝去）但是性能会稍微下降。令一点要注意的是在设备到主机间传输数据的时候，由于统一内存方式是一个主机内存地址和设备内存地址的一个映射这样的，但是主机访问设备内存的时候是需要通过一个页面故障来进行的（怎么感觉是翻译的问题，应该是中断）。

共享内存和常量内存

接下来是最后两个内存的讲解，其实比较奇怪，前面已经详细说明了一个全局内存和缓存的一些相关关系，这两个似乎只有很少的出场机会。实际上，对于跨全局内存以及非合并的内存访问在有些情况下还是会遇到的，算法层面可能无法完全抹除非合并内存访问问题，但是跨全局内存的问题还是可以使用共享内存以及常量内存解决的。

概述

共享内存其实是在SM缓存内部的，在图中称为SMEM的部分，用于一个线程块内部的线程共享，也因此可以减少对全局变量的访问。而且这部分缓存是可以程序显式管理的，所以也可以被描述为可编程管理的缓存。但是这个东西的使用并不是简单的必然带来高效，如果线程块的32个线程顺序访问共享内存的，则同样会带来较大的性能丢失，这时候就没有使用共享内存的意义了，也因此需要对其进行算法和程序上的优化。而且，由于共享内存是块内的线程共享的，一个核函数对共享内存的需求越多，整体处于活跃状态的线程块也是会下降的。

共享内存的变量定义也并不复杂：__shared__。但是不同的域内使用这样的定义效果也不同。核函数内的声明，则其作用域只在该内核内，而在核函数外进行声明，则该变量的作用域对于所有核函数都是全局的。

其次，共享内存是一个可以在编译阶段不指定具体大小的内存空间，可以使用extern关键字进行声明，但是由于编译位置，在调用的时候，还是需要让核函数知道到底需要多少空间的，这时候就需要用上<<<>>>中的第三个参数了：

kernel«<grid, block, isize * sizeof(int)»>(…)

要注意，这种动态声明只能用于一维数组。

对于内存访问的优化我们是需要考虑延迟和带宽的，但是对于共享内存，由于其是放在SM内因此可以用来隐藏内存延迟。

共享内存的设计非常有趣，为了获得高内存带宽，共享内存被分为32个大小相同的内存模型，每个称为存储体，可以被同时访问，在这种情况下，显然数据在不同的存储体的时候性能会最大化，在同一个存储体的数据会被多次访问从而得到更多的内存事务。

既然知道这个道理，那应该如何进行内存的访问呢，对于不同版本的计算设备，其存储体的宽度不太相同（2.x是4B而3.x是8B），因此共享内存地址到存储体索引的映射可以按照如下公式进行计算：

$$存储体索引=(字节地址/4字节)%32存储体$$

对于4B的访问可能会导致冲突，但是对于8B的一个字访问，由于取出的是4B的内容，因此访问一个地址不一定会导致冲突：

如图所示5-9是一个双向冲突而5-10是一个三向冲突。

既然为了尽可能减少冲突，那么一个有效的手段就是把数据进行填充，把原本属于一个存储块的数据放到不同的存储块中（但是实际上由于不同的架构的存储块位宽不同，可能会导致冲突）。在很前面的地方就说明了缓存的设计：

cudaError_t cudaDeviceSetCacheConfig(cudaFuncCache cacheConfig);

这里就可以用到了：

cudaFuncCachePreferNone: no preference(default)
cudaFuncCachePreferShared: prefer 48KB shared memory and 16 KB L1 cache
cudaFuncCachePreferL1: prefer 48KB L1 cache and 16 KB shared memory
cudaFuncCachePreferEqual: prefer 32KB L1 cache and 32 KB shared memory

可以根据使用情况来决定使用那种配置了（共享内存用的多就2，缓存用得多就3）。L1缓存是GPU控制存储和删除的，因此共享内存的使用算是给了用户更大的控制空间。（除了上面使用固定的设置以外也可以指定不同的核函数使用不同的方案cudaFuncSetCacheConfig）。对于CPU的内存管理往往可以使用一些FIFO等算法实现，但是由于GPU的访存相当复杂，其不可预知导致我们尽可能使用自己的方式显示管理。

最后是同步的相关概念，虽然不知道为什么在这里说明。同步的两个方法就是障碍和内存栅栏。在灵活使用同步之前的一道门槛是理解CUDA的弱排序内存模型。这个内存模型也很容易理解：即写入数据的顺序和这些数据在源代码的访问数据不一致，换句话说，核函数内连续两个内存访问指令，如果独立，其不一定哪个先被执行。所以必须引入同步手段防止不一致的结果。也因此有了__syncthreads()用于线程块的阻塞（如果想要实现线程块之间的同步，建议从核函数启动的地方出发）和内存栅栏。

但是内存栅栏的层次分为块、网格和内存栅栏三个层次：

__threadfence_block用于线程块内的栅栏，使得调用线程对共享内存和全局内存的读写结果对同一块的所有线程可见（这个没必要，本身就是可见的而且不执行线程同步）
__threadfence挂起调用的线程直到所有写操作对相同网格内的所有线程可见
__threadfence_system这个更是离谱，对整个系统可见

最后是volatile关键字修饰变量防止编译器优化把数据缓存到寄存器或内存中，防止写回时不一致的场景。（说到底都是为了保持数据的一致性）

共享内存的数据布局

共享内存从前面的截图可以看出来居然是二维的方形内存，因此其灵活性也是大大增加了。

那么既然是这样定义的，接下来就可以实践一下了，首先思考一下下面哪种方式性能高：

tile[threadIdx.y][threadIdx.x]
tile[threadIdx.x][threadIdx.y]

由于前面知道了一个bank的访问是会排队的，所以显然一行一行的访问更快，在这种情况下，合理的方式就是第一种存取方案。直觉上的是合理的，但是往往还是得来点实验说明一下，也就是接下来要讲的了。

我们使用这样一个例子进行说明：在核函数内仅进行两个操作：将全局内存按照一定的索引（至于哪种索引就是接下来要实践的内容了）写到二维共享内存中、从共享内存中按照某种索引再写到全局内存里。

__global__ void setRowReadRow(int * out)
{
    __shared__ int tile[BDIMY][BDIMX];
    unsigned int idx=threadIdx.y*blockDim.x+threadIdx.x;

    tile[threadIdx.y][threadIdx.x]=idx;
    __syncthreads();
    out[idx]=tile[threadIdx.y][threadIdx.x];
}
__global__ void setColReadCol(int * out)
{
    __shared__ int tile[BDIMY][BDIMX];
    unsigned int idx=threadIdx.y*blockDim.x+threadIdx.x;

    tile[threadIdx.x][threadIdx.y]=idx;
    __syncthreads();
    out[idx]=tile[threadIdx.x][threadIdx.y];
}

__global__ void setRowReadCol(int * out)
{
    __shared__ int tile[BDIMY][BDIMX];
    unsigned int idx=threadIdx.y*blockDim.x+threadIdx.x;

    tile[threadIdx.y][threadIdx.x]=idx;
    __syncthreads();
    out[idx]=tile[threadIdx.x][threadIdx.y];
}

一个一个拉出来批判，首先第一个是没问题的，因为连续的线程有连续的threadIdx.x值所以是连续访问的。而对于第二个，会按列主序访问，因此在Fermi架构下会有32路冲突而在Kepler有16路冲突（64位一次两个数据）

接下来一个示例类似，换一下行列的问题，在读的时候会导致16路冲突。

由于前面提到了使用填充数组来处理对存储体的冲突访问，所以这里实践一下。

__global__ void setRowReadColIpad(int * out)
{
    __shared__ int tile[BDIMY][BDIMX+IPAD];
    unsigned int idx=threadIdx.y*blockDim.x+threadIdx.x;

    tile[threadIdx.y][threadIdx.x]=idx;
    __syncthreads();
    out[idx]=tile[threadIdx.x][threadIdx.y];
}

上述代码的每行增加了一个元素，使得列元素分布在不同的存储体从而读写无冲突了，非常微妙倒是。对于静态的处理已经很舒服，但是在动态的时候就会稍微复杂一点了：

__global__ void setRowReadColRectDynPad(int * out)
{
    extern __shared__ int tile[];
    unsigned int idx=threadIdx.y*blockDim.x+threadIdx.x;
    unsigned int icol=idx%blockDim.y;
    unsigned int irow=idx/blockDim.y;
    unsigned int row_idx=threadIdx.y*(IPAD+blockDim.x)+threadIdx.x;
    unsigned int col_idx=icol*(IPAD+blockDim.x)+irow;
    tile[row_idx]=idx;
    __syncthreads();
    out[idx]=tile[col_idx];
}
setRowReadColDynIpad<<<grid,block,(BDIMX+IPAD)*BDIMY*sizeof(int)>>>(out);

从二维到一维转换时对于每一行必须跳过填充的内存空间。对应的需要在启动核函数的时候指定填充的共享内存的大小。与前面静态填充的效果一样。整体的效果如下：

减少全局内存访问

对全局内存的访问减少最简单的就是使用缓存（但是一般是不可编程的所以合理的手段就是使用共享内存）这里用一个前面归约求和的例子说明吧。

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

    smem[tid]=idata[tid];
	__syncthreads();
	//in-place reduction in global memory
	if(blockDim.x>=1024 && tid <512)
		smem[tid]+=smem[tid+512];
	__syncthreads();
	if(blockDim.x>=512 && tid <256)
		smem[tid]+=smem[tid+256];
	__syncthreads();
	if(blockDim.x>=256 && tid <128)
		smem[tid]+=smem[tid+128];
	__syncthreads();
	if(blockDim.x>=128 && tid <64)
		smem[tid]+=smem[tid+64];
	__syncthreads();
	//write result for this block to global mem
	if(tid<32)
	{
		volatile int *vsmem = smem;
		vsmem[tid]+=vsmem[tid+32];
		vsmem[tid]+=vsmem[tid+16];
		vsmem[tid]+=vsmem[tid+8];
		vsmem[tid]+=vsmem[tid+4];
		vsmem[tid]+=vsmem[tid+2];
		vsmem[tid]+=vsmem[tid+1];

	}

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

}

相较于上一个版本主要的变化就是smem的使用，但是尽管只是变换了一下内存，时间就嘎嘎又提升了接近两倍。从结果来看很明显的减少了对全局内存的访问。

然而，在前面的核函数中仅仅每个线程块处理一个数据块，但是现在可以再进一步，一次展开四个线程块从而：增加并行I/O提高吞吐量、减少存储事务。（实际上这部分代码在第三章说明过了，我就说为什么这么怪）

合并的全局内存访问

共享内存的另一个重大作用是合并全局内存的访问。前面的矩阵转置问题读是可以合并的但是写是交叉进行的，因此会大幅度影响带宽，因此这里就可以尝试先用共享内存进行转置之后对全局内存进行合并和写操作。

但是问题有好多，首先共享内存是有存储片的，这样按列获取性能应该很差，而且需要数据全部读到共享内存中才能进行？接下来从代码开始理解：

__global__ void transformSmem(float * in,float* out,int nx,int ny)
{
	__shared__ float tile[BDIMY][BDIMX];
	unsigned int ix,iy,transform_in_idx,transform_out_idx;
	ix=threadIdx.x+blockDim.x*blockIdx.x;
    iy=threadIdx.y+blockDim.y*blockIdx.y;
	transform_in_idx=iy*nx+ix;

	unsigned int bidx,irow,icol;
	bidx=threadIdx.y*blockDim.x+threadIdx.x;
	irow=bidx/blockDim.y;
	icol=bidx%blockDim.y;


	ix=blockIdx.y*blockDim.y+icol;
	iy=blockIdx.x*blockDim.x+irow;


	transform_out_idx=iy*ny+ix;

	if(ix<nx&& iy<ny)
	{
		tile[threadIdx.y][threadIdx.x]=in[transform_in_idx];
		__syncthreads();
		out[transform_out_idx]=tile[icol][irow];

	}

}

代码中最多的实际上是索引的转换部分，核心的地方就是if语句里面的三行，首先把全局内存的数据放到了共享内存中，之后按列进行数据读出，但是确实，读取一列的时候确实会带来存储体冲突（会带来16路冲突，尽管如此性能还是提升了）

但是对于共享存储有个非常有效的解决方案──使用填充共享内存（至于添加1还是2得取决于存储块的位宽），稍微修改一下就好了：

__shared__ float tile[BDIMY][BDIMX+IPAD];

通过错开访问和存储，直接一波消除了所有的存储体冲突。

最后就是使用一些展开技术了，相较于一个线程使用一个数据块，不如一下处理两个数据元素：


__global__ void transformSmemUnrollPad(float * in,float* out,int nx,int ny)
{
	__shared__ float tile[BDIMY*(BDIMX*2+IPAD)];


	unsigned int ix,iy,transform_in_idx,transform_out_idx;
	ix=threadIdx.x+blockDim.x*blockIdx.x*2;
    iy=threadIdx.y+blockDim.y*blockIdx.y;
	transform_in_idx=iy*nx+ix;

	unsigned int bidx,irow,icol;
	bidx=threadIdx.y*blockDim.x+threadIdx.x;
	irow=bidx/blockDim.y;
	icol=bidx%blockDim.y;


	unsigned int ix2=blockIdx.y*blockDim.y+icol;
	unsigned int iy2=blockIdx.x*blockDim.x*2+irow;


	transform_out_idx=iy2*ny+ix2;

	if(ix+blockDim.x<nx&& iy<ny)
	{
		unsigned int row_idx=threadIdx.y*(blockDim.x*2+IPAD)+threadIdx.x;
		tile[row_idx]=in[transform_in_idx];
		tile[row_idx+BDIMX]=in[transform_in_idx+BDIMX];
		__syncthreads();
		unsigned int col_idx=icol*(blockDim.x*2+IPAD)+irow;
        out[transform_out_idx]=tile[col_idx];
		out[transform_out_idx+ny*BDIMX]=tile[col_idx+BDIMX];

	}

}

也从中看出来了，其实难点就两点：

利用书中的技术和论断进行优化以及
索引的映射

常量内存

和全局内存一样位于DRAM，但是有一个专用的片上缓存，和一级缓存和共享内存一样。但是大小有64KB。

相较于前面的各种内存又是横着读又是竖着读或者对齐操作，常量内存的最优访问模式就是所有线程访问相同的地址。至此可以简单列一个表格说明每种内存的最佳访问模式：

寄存器	共享内存	全局内存	常量内存	本地内存
多访问就好	对齐	经过一级缓存后合成一个访存事务	访问一个地方	无法被管理咯

代码实现的例子是一个九点模板：

__constant__ float coef[TEMP_RADIO_SIZE];//if in midle of the program will be error
....
__global__ void stencil_1d(float * in,float * out)
{
    __shared__ float smem[BDIM+2*TEMP_RADIO_SIZE];
    int idx=threadIdx.x+blockDim.x*blockIdx.x;
    int sidx=threadIdx.x+TEMP_RADIO_SIZE;
    smem[sidx]=in[idx];

    if (threadIdx.x<TEMP_RADIO_SIZE)

    {
        if(idx>TEMP_RADIO_SIZE)
            smem[sidx-TEMP_RADIO_SIZE]=in[idx-TEMP_RADIO_SIZE];
        if(idx<gridDim.x*blockDim.x-BDIM)
            smem[sidx+BDIM]=in[idx+BDIM];

    }

    __syncthreads();
    if (idx<TEMP_RADIO_SIZE||idx>=gridDim.x*blockDim.x-TEMP_RADIO_SIZE)
        return;
    float temp=.0f;
    #pragma unroll
    for(int i=1;i<=TEMP_RADIO_SIZE;i++)
    {
        temp+=coef[i-1]*(smem[sidx+i]-smem[sidx-i]);
    }
    out[idx]=temp;
    //printf("%d:GPU :%lf,\n",idx,temp);
}
int main(int argc,char** argv)
{
...
float templ_[]={-1.0,-2.0,2.0,1.0};
CHECK(cudaMemcpyToSymbol(coef,templ_,TEMP_RADIO_SIZE*sizeof(float)));
...
}

看上去还是相当方便使用的。在介绍内存的时候我们把只读缓存一笔带过了，当时说是为了纹理流水线但是可以用来使用的内容，对于只读缓存，其与L1缓存的好处在于更适合使用分散的读取（不方便用在相同的地址读取）

只读缓存相较于常量缓存可以存储较大的内容（虽然实际上更小）。一般很少自己手动用到这个东西，只有需要显式控制以及代码足够复杂乃至无法编译器判断的时候才会显式的需要手动管理只读缓存。这里我的选择是开摆，但是要知道两个关键字是干什么的：

__ldg(&input[idx])用来解引用并设置从只读缓存中获取。const __restrict__ *xxx用来说明通过只读缓存访问。不过这里既然说到了就说明其和常量内存必然是有共同性的：

两者都是只读的，但是只读缓存适合分散读取而常量缓存适合统一读取。

线程束洗牌操作

这个名字听起来很酷诶，在计算能力大于3.0的设备上可以实现束内线程的数据交流。对应设计的API也相对多一点而且很酷诶！

int __shfl(int var, int srcLane, int width=warpSize)

上式指令的返回值是var，通过由srcLane确定的同一线程束中的线程传递给__shfl，而srcLane的含意变化取决于width。举个简单的例子：

int y = shfl(x,3,16);

表示线程0~15从线程3获得x的值，16~31从线程19获得x的值。接下来简单说明一下剩下几个常用的API和一个实际的例子：

int __shfl_up(int var, unsigned int delta, int width=warpSize)

int __shfl_down(int var, unsigned int delta, int width=warpSize)

int __shfl_xor(int var, int laneMask, int width=warpSize)

这里就不对每一种API都做详细测试了，这里把前面的并行归约再拉出来优化一遍。在最后的优化中我们把最后的32个线程的归约进行了完全的展开，现在对于每个线程束可以直接进行合并（虽然听起来更复杂了，但是由于不需要读取内存了，性能会更好一点）

__inline__ __device__ int warpReduce(int localSum)
{
    localSum += __shfl_xor(localSum, 16);
    localSum += __shfl_xor(localSum, 8);
    localSum += __shfl_xor(localSum, 4);
    localSum += __shfl_xor(localSum, 2);
    localSum += __shfl_xor(localSum, 1);

    return localSum;
}
__global__ void reduceShfl(int * g_idata,int * g_odata,unsigned int n)
{
	//set thread ID
    __shared__ int smem[DIM];
	unsigned int idx = blockDim.x*blockIdx.x+threadIdx.x;
	//convert global data pointer to the

	int mySum=g_idata[idx];
	int laneIdx=threadIdx.x%warpSize;
	int warpIdx=threadIdx.x/warpSize;

	mySum=warpReduce(mySum);

	if(laneIdx==0)
		smem[warpIdx]=mySum;
	__syncthreads();
	mySum=(threadIdx.x<DIM)?smem[laneIdx]:0;
	if(warpIdx==0)
		mySum=warpReduce(mySum);
	if(threadIdx.x==0)
		g_odata[blockIdx.x]=mySum;

}

原文还在这里扯了一小段解释，感觉倒是没那么复杂。但是对索引的理解还是得深刻一点啧啧啧。

流和并发

前面的关注点可能只有内核级别的开发，分别从编程模型、执行模型和内存模型三个角度进行了优化，接下来就是对网格级别的并发，多个内核在统一设备上同时执行使得设备的利用率更高。

CUDA流

cuda流是一系列异步的CUDA操作（实际上由于资源的问题也可能会导致同步）

所有的CUDA操作都在一个流中显式或者隐式的运行，其主要有两种类型：

隐式声明的流（空流）
显式声明的流（非空流）

如果我们没有特别声明一个流，那么我们的所有操作是在默认的空流中完成的，我们前面的所有例子都是在默认的空流中进行的。

空流是没办法管理的，因为他连个名字都没有，似乎也没有默认名，所以当我们想控制流，非空流是非常必要的。

基于流的异步内核启动和数据传输支持以下类型的粗粒度并发

重叠主机和设备计算
重叠主机计算和主机设备数据传输
重叠主机设备数据传输和设备计算
并发设备计算（多个设备）

CUDA编程和普通的C++不同的就是，我们有两个“可运算的设备”也就是CPU和GPU这两个东西，这种情况下，他们之间的同步并不是每一步指令都互相通信执行进度的，设备不知道主机在干啥，主机也不是完全知道设备在干啥。但是数据传输是同步的，也就是主机要等设备接收完数据才干别的。

前面用的cudaMemcpy就是个同步操作，我们还提到过隐式同步——从设备复制结果数据回主机，要等设备执行完。当然数据传输有异步版本：

cudaError_t cudaMemcpyAsync(void* dst, const void* src, size_t count,cudaMemcpyKind kind, cudaStream_t stream = 0);

值得注意的就是最后一个参数，stream表示流，一般情况设置为默认流，这个函数和主机是异步的，执行后控制权立刻归还主机，当然我们需要声明一个非空流：

cudaError_t cudaStreamCreate(cudaStream_t* pStream);

主机虚拟内存中分配的数据在物理内存中是随时可能被移动的，我们必须确保其在整个生存周期中位置不变，这样在异步操作中才能准确的转移数据，否则如果操作系统移动了数据的物理地址，那么我们的设备可能还是回到之前的物理地址取数据，这就会出现未定义的错误。

在非空流中执行内核需要在启动核函数的时候加入一个附加的启动配置：

kernel_name<<<grid, block, sharedMemSize, stream>>>(argument list);

pStream参数就是附加的参数，使用目标流的名字作为参数，比如想把核函数加入到a流中，那么这个stream就变成a。

前面我们为一个流分配资源，当然后面就要回收资源，回收方式：

cudaError_t cudaStreamDestroy(cudaStream_t stream);

这个回收函数很有意思，由于流和主机端是异步的，你在使用上面指令回收流的资源的时候，很有可能流还在执行，这时候，这条指令会正常执行，但是不会立刻停止流，而是等待流执行完成后，立刻回收该流中的资源。这样做是合理的也是安全的。当然，我们可以查询流执行的怎么样了，下面两个函数就是帮我们查查我们的流到哪了：

cudaError_t cudaStreamSynchronize(cudaStream_t stream);
cudaError_t cudaStreamQuery(cudaStream_t stream);

这两条执行的行为非常不同，cudaStreamSynchronize会阻塞主机，直到流完成。cudaStreamQuery则是立即返回，如果查询的流执行完了，那么返回cudaSuccess否则返回cudaErrorNotReady。

下面这段示例代码就是典型多个流中调度CUDA操作的常见模式：

for (int i = 0; i < nStreams; i++) {
    int offset = i * bytesPerStream;
    cudaMemcpyAsync(&d_a[offset], &a[offset], bytePerStream, streams[i]);
    kernel<<grid, block, 0, streams[i]>>(&d_a[offset]);
    cudaMemcpyAsync(&a[offset], &d_a[offset], bytesPerStream, streams[i]);
}
for (int i = 0; i < nStreams; i++) {
    cudaStreamSynchronize(streams[i]);
}

相较于单核的进行，显然使用多核进行可以提升运行的速度。虽然但是由于是用到了PCIE进行数据传输所以还是需要等待数据传输结束之后才能进行新的操作。

由于架构的不同，不同架构的最大并行数量也有所不同，其中Fermi支持16路并发而Kepler支持32路并发。

然而尽管从逻辑上所有的流可以同时执行，但是由于硬件的限制，绝大多数情况下并不能实现十分充分的并行，因此还需要考虑如何调度流。首先是Fermi架构上的虚假依赖的例子：

对于Fermi架构而言，由于只有一个硬件工作队列（所以实际上是虚假的并发），而且会检测依赖关系，因此尽管A、P、X三个任务序列不会互相依赖但是仍然无法充分并发进行。

当然上述的解决方案也很简单：多开几个工作队列就好了，这也就是后续架构的Hyper-Q技术，Kelper架构的32个硬件工作队列可以同时执行多个流从而实现流的并发，减少虚假依赖。

在计算能力3.5以上的设备就支持优先级了：

优先级只影响核函数，不影响数据传输，高优先级的流可以占用低优先级的工作。下面函数创建一个有指定优先级的流

cudaError_t cudaStreamCreateWithPriority(cudaStream_t* pStream, unsigned int flags,int priority);

不同的设备有不同的优先级等级，下面函数可以查询当前设备的优先级分布情况：

cudaError_t cudaDeviceGetStreamPriorityRange(int *leastPriority, int *greatestPriority);

leastPriority表示最低优先级（整数，远离0）

greatestPriority表示最高优先级（整数，数字较接近0）

如果设备不支持优先级返回0

CUDA事件不同于我们前面介绍的内存事务，不要搞混，事件也是软件层面上的概念。事件的本质就是一个标记，它与其所在的流内的特定点相关联。可以使用时间来执行以下两个基本任务：

同步流执行
监控设备的进展

流中的任意点都可以通过API插入事件以及查询事件完成的函数，只有事件所在流中其之前的操作都完成后才能触发事件完成。默认流中设置事件，那么其前面的所有操作都完成时，事件才出发完成。

事件就像一个个路标，其本身不执行什么功能，就像我们最原始测试c语言程序的时候插入的无数多个printf一样。

事件的声明如下：

cudaEvent_t event;

同样声明完后要分配资源：

cudaError_t cudaEventCreate(cudaEvent_t* event);

回收事件的资源

cudaError_t cudaEventDestroy(cudaEvent_t event);

如果回收指令执行的时候事件还没有完成，那么回收指令立即完成，当事件完成后，资源马上被回收。

事件的一个主要用途就是记录事件之间的时间间隔。

事件通过下面指令添加到CUDA流：

cudaError_t cudaEventRecord(cudaEvent_t event, cudaStream_t stream = 0);

在流中的事件主要左右就是等待前面的操作完成，或者测试指定流中操作完成情况，下面和流类似的事件测试指令（是否出发完成）会阻塞主机线程知道事件被完成。

cudaError_t cudaEventSynchronize(cudaEvent_t event);

同样，也有异步版本：

cudaError_t cudaEventQuery(cudaEvent_t event);

这个不会阻塞主机线程，而是直接返回结果和stream版本的类似。

另一个函数用在事件上的是记录两个事件之间的时间间隔：

cudaError_t cudaEventElapsedTime(float* ms, cudaEvent_t start, cudaEvent_t stop);

这个函数记录两个事件start和stop之间的时间间隔，单位毫秒，两个事件不一定是同一个流中。这个时间间隔可能会比实际大一些，因为cudaEventRecord这个函数是异步的，所以加入时间完全不可控，不能保证两个事件之间的间隔刚好是两个事件之间的。一段简单的记录事件时间间隔的代码

// create two events
cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
// record start event on the default stream
cudaEventRecord(start);
// execute kernel
kernel<<<grid, block>>>(arguments);
// record stop event on the default stream
cudaEventRecord(stop);
// wait until the stop event completes
cudaEventSynchronize(stop);
// calculate the elapsed time between two events
float time;
cudaEventElapsedTime(&time, start, stop);
// clean up the two events
cudaEventDestroy(start);
cudaEventDestroy(stop);

这段代码显示，我们的事件被插入到空流中，设置两个事件作为标记，然后记录他们之间的时间间隔。

cudaEventRecord是异步的，所以间隔不准，这是特别要注意的。

然而，并发归并发，我们依然需要进行同步从而减少各种问题和损耗，从主机的角度而言，CUDA操作可以分为两类：

内存相关操作
内核启动

内核启动总是异步的，虽然某些内存是同步的，但是他们也有异步版本。

前面我们提到了流的两种类型：

异步流（非空流）
同步流（空流/默认流）

没有显式声明的流式默认同步流，程序员声明的流都是异步流，异步流通常不会阻塞主机，同步流中部分操作会造成阻塞，主机等待，什么都不做，直到某操作完成。

非空流并不都是非阻塞的，其也可以分为两种类型：

阻塞流
非阻塞流

虽然正常来讲，非空流都是异步操作，不存在阻塞主机的情况，但是有时候可能被空流中的操作阻塞。如果一个非空流被声明为非阻塞的，那么没人能阻塞他，如果声明为阻塞流，则会被空流阻塞。

有点晕，就是非空流有时候可能需要在运行到一半和主机通信，这时候我们更希望他能被阻塞，而不是不受控制，这样我们就可以自己设定这个流到底受不受控制，也就是是否能被阻塞，下面我们研究如何使用这两种流。

cudaStreamCreate创建的是阻塞流，意味着里面有些操作会被阻塞，直到空流中默写操作完成。

空流不需要显式声明，而是隐式的，他是阻塞的，跟所有阻塞流同步。下面这个过程很重要：

当操作A发布到空流中，A执行之前，CUDA会等待A之前的全部操作都发布到阻塞流中，所有发布到阻塞流中的操作都会挂起，等待，直到在此操作指令之前的操作都完成，才开始执行。

有点复杂，因为这涉及到代码编写的过程和执行的过程，两个过程混在一起说，肯定有点乱，我们来个例子压压惊就好了：

kernel_1<<<1, 1, 0, stream_1>>>();
kernel_2<<<1, 1>>>();
kernel_3<<<1, 1, 0, stream_2>>>();

上面这段代码，有三个流，两个有名字的，一个空流，我们认为stream_1和stream_2是阻塞流，空流是阻塞的，这三个核函数都在阻塞流上执行，具体过程是，kernel_1被启动，控制权返回主机，然后启动kernel_2，但是此时kernel_2 不会并不会马山执行，他会等到kernel_1执行完毕，同理启动完kernel_2 控制权立刻返回给主机，主机继续启动kernel_3,这时候kernel_3 也要等待，直到kernel_2执行完，但是从主机的角度，这三个核都是异步的，启动后控制权马上还给主机。

然后我们就想创建一个非阻塞流，因为我们默认创建的是阻塞版本：

cudaError_t cudaStreamCreateWithFlags(cudaStream_t* pStream, unsigned int flags);

第二个参数就是选择阻塞还是非阻塞版本：

cudaStreamDefault;// 默认阻塞流
cudaStreamNonBlocking: //非阻塞流，对空流的阻塞行为失效。

如果前面的stream_1和stream_2声明为非阻塞的，那么上面的调用方法的结果是三个核函数同时执行。

既然是并行的操作，那么仅仅进行阻塞自然不是很灵活，可以考虑使用同步来实现复杂的并行操作。

前面提到了好多同步的函数，但是这些都是显式的同步手段，还有一些隐式的同步手段，大多数都是在内存操作上的：

锁页主机内存分布
设备内存分配
设备内存初始化
同一设备两地址之间的内存复制
一级缓存，共享内存配置修改

这个太细了，估计只有在真正编程的时候才会遇到这些问题了。