CUDA编程指南阅读笔记

转自http://blog.csdn.net/csgxy123/article/category/1500471

随着多核CPU和众核GPU的到来，并行编程已经获得了业界愈来愈多的重视，CPU-GPU异构程序可以极大提升现有计算机系统的运算性能，对于科学计算等运算密集型程序有着很是重要的意义。这一系列文章是根据《CUDA C语言编程指南》来整理的，该指南是NVIDIA公司提供的CUDA学习资料，介绍了CUDA编程最基本最核心的概念，是学习CUDA必不可少的阅读材料。

初学CUDA，笔记错误之处在所不免，还请发现问题的诸位读者不吝赐教。

1. 什么是CUDA？

      CUDA全称是Compute Unified Device Architecture，中文名称即统一计算设备架构，它是NVIDIA公司提出了一种通用的并行计算平台和编程模型。使用CUDA，咱们能够开发出同时在CPU和GPU上运行的通用计算程序，更加高效地利用现有硬件进行计算。为了简化并行计算学习，CUDA为程序员提供了一个类C语言的开发环境以及一些其它的如FORTRAN、DirectCOmpute、OpenACC的高级语言/编程接口来开发CUDA程序。

2. CUDA编程模型如何扩展？

     咱们知道，不一样的GPU拥有不一样的核心数目，在核心较多的系统上CUDA程序运行的时间较短，而在核心较少的系统上CUDA程序的执行时间较多。那么，CUDA是如何作到的呢？

      并行编程的中心思想是分而治之：将大问题划分为一些小问题，再把这些小问题交给相应的处理单元并行地进行处理。在CUDA中，这一思想便体如今它的具备两个层次的问题划分模型。一个问题能够首先被粗粒度地划分为若干较小的子问题，CUDA使用被称为块（Block）的单元来处理它们，每一个块都由一些CUDA线程组成，线程是CUDA中最小的处理单元，将这些较小的子问题进一步划分为若干更小的细粒度的问题，咱们即可以使用线程来解决这些问题了。对于一个普通的NVIDIA GPU，其CUDA线程数目一般能达到数千个甚至更多，所以，这样的问题划分模型即可以成倍地提高计算机的运算性能。

      GPU是由多个流水多处理器构成的，流水处理器以块（Block）为基本调度单元，所以，对于流水处理器较多的GPU，它一次能够处理的块（Block）更多，从而运算速度更快，时间更短。而反之对于流水处理器较少的GPU，其运算速度便会较慢。这一原理能够经过下图形象地看出来：

3. CUDA基本概念（上）

      本节将介绍CUDA的一些基本的编程概念，该节用到的例子来自于CUDA Sample中的VectorAdd项目。

3.1 内核（Kernels）

      CUDA C是C语言的一个扩展，它容许程序员定义一种被称为内核函数（Kernel Functions）的C函数，内核函数运行在GPU上，一旦启动，CUDA中的每个线程都将会同时并行地执行内核函数中的代码。

      内核函数使用关键字__global__来声明，运行该函数的CUDA线程数则经过<<<...>>>执行配置语法来设置。（参见章节"C语言扩展"），每个执行内核函数的线程都由一个惟一的线程ID，这一ID能够经过在内核函数中访问threadIdx变量来获得。

      下面经过一些示例代码来展现刚刚提到的这些概念该如何应用在编程中：

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

 int main() {
     ...
     // Kernel invocation with N threads
     VecAdd<<<1, N>>>(A, B, C);
    ...
 }

在上面的代码中，N个线程将会并行地同时执行加法运算。

3.2 线程层次（Thread Hierarchy）

      CUDA的每个线程都有其线程ID，线程的ID信息由变量threadIdx给出。threadIdx是CUDA C语言的内建变量，一般它用一个三维数组来表示。使用三维数组的方便之处在于能够很方便地表示一维、二维和三维线程索引，进而方便地表示一维、二维和三维线程块（thread block）。这样，不管是数组、矩阵仍是体积的计算，均可以很容易地使用CUDA进行运算。

      线程的索引与线程ID之间存在着直接的换算关系，对于一个索引为(x, y, z)的线程来讲：

      一、若是线程块（block）是一维的，则线程ID = x

      二、若是线程块是二维的，假设块尺寸为（Dx，Dy），那么线程ID = x + y * Dx

      三、若是线程块是三维的，设其尺寸为（Dx，Dy，Dz），那么线程ID = x + y * Dx + z * Dx * Dy

      下面的例子展现了两个NxN矩阵相加的CUDA实现：

 // Kernel definition
 __global__ void MatAdd(float A[N][N], float B[N][N], float C[N][N]) {
     int i = threadIdx.x;
     int j = threadIdx.y;
     C[i][j] = A[i][j] + B[i][j];
 }

 int main() {
     ...
     // Kernel invocation with one block of N * N * 1 threads
     int numBlocks = 1;
     dim3 threadsPerBlock(N, N);
     MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
     ...
 }

     每一个线程块（block）中的线程数量是有限制的，由于依据前面所说，同一线程块（block）中的全部线程都会被分配到同一个处理器核上运行，共享有限的存储资源，所以对于当前的GPU，线程块所能包含的最大线程数目为1024。

      上面的例子中numBlocks表明线程块的数量，这里的值为1。在通常的CUDA程序中，这个值一般大于1，也就是说将会有多个线程块被分配到多个处理器核中同时进行处理，这样就大大提升了程序的并行性。

      在CUDA中，线程块包含在线程格（grid）当中，线程格能够是一维、二维或者三维的，线程格的尺寸通常根据待处理数据的规模或者处理器的数量来指定。线程格中所包含的线程块数目一般远远大于GPU处理器核心的数目。下图展现了线程格（grid）、线程块（block）以及线程（thread）之间的关系：

      内核函数的调用能够简化为kernel<<<A，B>>>(parameters)，在尖括号中，A表明线程格（grid）的尺寸，它能够是三维的，用类型dim3表示，也能够是一维的，用int类型表示。B表明线程块（block）的尺寸，它与A相似，也可分别用dim3或int类型表示。

      在内核函数内部，CUDA为咱们内建了一些变量用于访问线程格、线程块的尺寸和索引等信息，它们是：

      1. gridDim：表明线程格（grid）的尺寸，gridDim.x为x轴尺寸，gridDim.y、gridDim.z相似。拿上图来讲，它的gridDim.x = 3，gridDim.y = 2，gridDim.z = 1。

      2. blockIdx：表明线程块（block）在线程格（grid）中的索引值，拿上图来讲，Block（1,1）的索引值为：blockIdx.x = 1，blockIdx.y = 1。

      3. blockDim：表明线程块（block）的尺寸，blockDIm.x为x轴尺寸，其它依此类推。拿上图来讲，注意到Block（1,1）包含了4 * 3个线程，所以blockDim.x = 4, blockDim.y = 3。

      4. threadIdx：线程索引，前面章节已经详细探讨过了，这里再也不赘述。

      明白了这些变量的含义，那么下面的矩阵加法程序便不难理解了：

 // Kernel definition
 __global__ void MatAdd(float A[N][N], float B[N][N], float C[N][N]) {
     int i = blockIdx.x * blockDim.x + threadIdx.x;
     int j = blockIdx.y * blockDim.y + threadIdx.y;
     if (i < N && j < N)
         C[i][j] = A[i][j] + B[i][j];
 }

 int main() {
     ...
     // Kernel invocation
     dim3 threadsPerBlock(16, 16);
     dim3 numBlocks(N / threadsPerBlock.x, N / threadsPerBlock.y);
     MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
     ...
 }

      在上面的程序中，线程块（block）的尺寸是16x16，这是CUDA编程中一个很是广泛的选择。线程格（grid）包含了足够多的线程块（block）来进行计算。

      线程块（block）是独立执行的，在执行的过程当中线程块之间互不干扰，所以它们的执行顺序是随机的。

      同一线程块中的线程能够经过访问共享内存（shared memory）或者经过同步函数__syncthreads()来协调合做。这些概念将在之后的章节中详细解释。

3. CUDA基本概念（下）

3.3 内存层次（Memory Hierarchy）

      在GPU上CUDA线程能够访问到的存储资源有不少，每一个CUDA线程拥有独立的本地内存（local Memory）；每个线程块（block）都有其独立的共享内存（shared memory），共享内存对于线程块中的每一个线程都是可见的，它与线程块具备相同的生存时间；同时，还有一片称为全局内存（global memory）的区域对全部的CUDA线程都是可访问的。

      除了上述三种存储资源之外，CUDA还提供了两种只读内存空间：常量内存（constant memory）和纹理内存（texture memory），同全局内存相似，全部的CUDA线程均可以访问它们。对于一些特殊格式的数据，纹理内存提供多种寻址模式以及数据过滤方法来操做内存。这两类存储资源主要用于一些特殊的内存使用场合。

      一个程序启动内核函数之后，全局内存、常量内存以及纹理内存将会一直存在直到该程序结束。下面是CUDA的内存层次图：

3.4 异构编程（Heterogeneous Programming）

      CUDA的异构编程模型假定CUDA线程都运行在一个可被看作CPU协处理器的芯片上，这就使得CUDA内核函数能够和CPU端C程序的运行并行运行，从而加快程序的运行效率。为了达到这个效果，CUDA程序须要管理两大块由DRAM构成的内存区域：CPU端能够访问到的主机内存（host memory）以及GPU端供CUDA内核访问到的设备内存（device memory），设备内存主要由全局内存、常量内存以及纹理内存构成。如今，CUDA程序的运行机制便很明了了：CPU端代码生成原始数据，经过CUDA运行时函数库将这些原始数据传输到GPU上，在CPU端启动CUDA内核函数进行运算，而后将运算结果从设备端传输到主机端，计算任务便完成了。

4. CUDA C语言编程接口

    异构程序设计跟传统的串行程序设计差异是很大的，学习起来也是很是不容易的。NVIDIA很是够意思，为了简化CUDA的学习曲线，它采用了绝大多数程序员都熟悉的C语言做为其根基，CUDA C是NVIDIA为程序员提供的一类编程接口，它其实是一个C语言的扩展，在C的基础上增长了一些新的语法和变量，而且提供了功能丰富的库函数，方便程序员使用GPU进行异构计算。

        除了前面章节提到的CUDA最基本、最核心的概念之外，CUDA C呈现给程序员的接口主要由两大类API构成，它们分别是CUDA Runtime API和CUDA Driver API，Runtime API其实是对于Driver API的封装，其目的天然是方便程序员的代码编写工做。Driver API为用户提供了更细一层的控制手段，经过它能够控制诸如CUDA Contexts（一种相似主机进程的概念）以及CUDA Modules（相似主机动态加载库的概念）等更加底层的CUDA模块。

4.1 NVCC编译器

        任何一种程序设计语言都须要相应的编译器将其编译为二进制代码，进而在目标机器上获得执行。对于异构计算而言，这一过程与传统程序设计语言是有一些区别的。为何？由于CUDA它本质上不是一种语言，而是一种异构计算的编程模型，使用CUDA C写出的代码须要在两种体系结构彻底不一样的设备上执行：一、CPU；二、GPU。所以，CUDA C的编译器所作的工做就有点略多了。一方面，它须要将源代码中运行在GPU端的代码编译获得能在CUDA设备上运行的二进制程序。另外一方面，它也须要将源代码中运行在CPU端的程序编译获得能在主机CPU上运行的二进制程序。最后，它须要把这两部分有机地结合起来，使得两部分代码可以协调运行。

        CUDA C为咱们提供了这样的编译器，它即是NVCC。严格意义上来说，NVCC并不能称做编译器，NVIDIA称其为编译器驱动（Compiler Driver），本节咱们暂且使用编译器来描述NVCC。使用nvcc命令行工具咱们能够简化CUDA程序的编译过程，NVCC编译器的工做过程主要能够划分为两个阶段：离线编译（Offline Compilation）和即时编译（Just-in-Time Compilation）。

离线编译（Offline Compilation）

    下面这幅图简单说明了离线编译的过程：

        在CUDA源代码中，既包含在GPU设备上执行的代码，也包括在主机CPU上执行的代码。所以，NVCC的第一步工做即是将两者分离开来，这一过程结束以后：

        1. 运行于设备端的代码将被NVCC工具编译为PTX代码（GPU的汇编代码）或者cubin对象（二进制GPU代码）；

        2. 运行于主机端的代码将被NVCC工具改写，将其中的内核启动语法（如<<<...>>>）改写为一系列的CUDA Runtime函数，并利用外部编译工具（gcc for linux，或者vc compiler for windows）来编译这部分代码，以获得运行于CPU上的可执行程序。

        完事以后，NVCC将自动把输出的两个二进制文件连接起来，获得异构程序的二进制代码。

即时编译（Just-in-time Compile）

        任何在运行时被CUDA程序加载的PTX代码都会被显卡的驱动程序进一步编译成设备相关的二进制可执行代码。这一过程被称做即时编译（just-in-time compilation）。即时编译增长了程序的装载时间，可是也使得编译好的程序能够重新的显卡驱动中得到性能提高。同时到目前为止，这一方法是保证编译好的程序在还未问世的GPU上运行的惟一解决方案。

        在即时编译的过程当中，显卡驱动将会自动缓存PTX代码的编译结果，以免屡次调用同一程序带来的重复编译开销。NVIDIA把这部分缓存称做计算缓存（compute cache），当显卡驱动升级时，这部分缓存将会自动清空，以使得程序可以自动得到新驱动为即时编译过程带来的性能提高。

        有一些环境变量能够用来控制即时编译过程：

        1. 设置CUDA_CACHE_DISABLE为1将会关闭缓存功能

        2. CUDA_CACHE_MAXSIZE变量用于指定计算缓存的字节大小，默认状况下它的值是32MB，它最大能够被设置为4GB。任何大于缓存最大值得二进制代码将不会被缓存。在须要的状况下，一些旧的二进制代码可能被丢弃以腾出空间缓存新的二进制代码。

        3. CUDA_CACHE_PATH变量用于指定计算缓存的存储目录地址，它的缺省值以下：

        4. 设置CUDA_FORCE_PTX_JIT为1会强制显卡驱动忽略应用程序中的二进制代码而且即时编译程序中的嵌入PTX代码。若是一个内核函数没有嵌入的PTX代码，那么它将会装载失败。该变量能够用来确认程序中存在嵌入的PTX代码。同时，使用即时编译（just-in-time Compilation）技术也可确保程序的向前兼容性。

4.2 兼容性

      一、二进制兼容性

      二进制代码是设备相关的，使用NVCC编译器编译时，若指定-code选项，则会编译产生目标设备的二进制cubin对象。例如，编译时使用-code=sm_13会产生适用于计算能力1.3的二进制代码。二进制代码在CUDA计算设备上具备小版本的向前兼容性，可是在大版本上不具有兼容性。也就是说，对于计算能力X.y的硬件，使用-code=sm_Xy编译后，程序可以运行于计算能力X.z（其中z>=y）的硬件上，但不能运行在计算能力M.n（M!=X）的硬件上。

      二、PTX代码兼容性

      不一样计算能力的设备所支持的PTX指令条数是不一样的，一些PTX指令只在拥有较高计算能力的设备上被支持。例如，全局内存（global Memory）的原子操做指令只能用于计算能力不小于1.1的设备；双精度浮点运算指令只能用于计算能力不小于1.3的设备。在将C语言编译为PTX代码时，NVCC使用-arch编译选项指定PTX代码目标设备的计算能力。所以，要想使用双精度运算，编译时必须使用选项-arch=sm_13（或使用更高的计算能力），不然NVCC会自动将双精度操做降级为单精度操做。

      为某一特定设备产生的PTX代码，在运行时老是可以被具备更高计算能力的设备JIT编译为可执行的二进制代码。

      三、应用程序兼容性

      执行CUDA程序有两种方式，一种是直接加载编译好的CUDA二进制代码运行，另外一种是首先加载程序中的PTX代码，再执行JIT编译获得二进制的设备可执行文件，而后运行。特别须要注意的是，为了让程序运行具备更高计算能力的将来设备上，必须让程序加载PTX代码。

      事实上，在一个CUDA C程序中能够嵌入不止一个版本的PTX/二进制代码。那么，具体执行时哪个版本的PTX或者二进制代码会获得执行呢？答案是：最兼容的那个版本。例如编译一个名为x.cu的CUDA源代码：

将会产生兼容计算能力1.1硬件的二进制代码（第一排的-gencode选项）以及兼容计算能力1.1设备的PTX和二进制代码，这些代码都将会嵌入到编译后的目标文件中。

      主机端将会产生一些额外的代码，在程序运行时，这些代码会自动决定装载哪个版本的代码来执行。对于上面的例子：

• 计算能力1.0的设备运行该程序将会装载1.0版本的二进制代码
• 计算能力1.一、1.2或者1.3的设备运行该程序将会装载1.1版本的二进制代码
• 计算能力2.0或者更高的设备运行该程序将会装载1.1版本的PTX代码进而对其进行JIT编译获得相应设备的二进制代码

      同时，x.cu还能够在程序中使用一些特殊的宏来改变不一样设备的代码执行路径。例如，对于计算能力1.1的设备而言，宏__CUDA_ARCH__等于110，在程序中能够对该宏的值进行判断，而后分支执行程序。

      NVCC用户手册列出了不少-arch，-code和-gencode等编译选项的简化书写形式。例如，-arch=sm_13就是-arch=compute_13 -code=compute13, sm_13的简化形式。更多详尽的内容请参阅该手册。

      四、C/C++兼容性

      NVCC编译器前端使用C++语法啊规则来处理CUDA源文件。在主机端，CUDA支持完整的C++语法；而在设备端，只有部分C++语法是被支持的。这方面更为详尽的讨论请参见《CUDA C程序设计指南》的C/C++语言支持章节。

      五、64位兼容性

      64位版本的nvcc编译器将设备代码编译为64位模式，即指针是64位的。运行64位设备代码的先决条件是主机端代码必须也使用64位模式进行编译。一样，32位版本的nvcc将设备代码编译为32位模式，这些代码也必须与相应的32位主机端代码相配合方能运行。

      32位nvcc编译器可使用-m64编译选项将设备代码编译为64位模式。同时64位nvcc编译器也可以使用-m32编译选项将设备代码编译为32位模式。

4.3 CUDA C Runtime

        CUDA C Runtime使用cudart动态连接库实现（cudart.dll或者cudart.so），运行时中全部的入口函数都以cuda为前缀。

4.3.1 初始化

        CUDA C Runtime函数库没有明确的初始化函数，在程序第一次调用Runtime库函数时它会自动初始化。所以，在记录Runtime函数调用时间和理解程序中第一个Runtime调用返回的错误代码时，须要将初始化考虑在内。

        在初始化期间，Runtime将会为系统中每个设备建立一个CUDA上下文（相似CPU中进程的数据结构），这个上下文是设备的基本上下文，它被程序中全部的主机线程所共享。建立过程在后台运行，而且，Runtime将隐藏基本上下文使之对Runtime API这一层次的程序员不可见。

        当一个主机线程调用cudaDeviceReset()函数时，它将会销毁线程当前控制设备的基本上下文。也就是说，当线程下一次调用runtime函数时将会重启初始化，一个新的CUDA基本上下文将被建立出来。

4.3.2 设备内存

        正如前面异构计算章节所讲，CUDA编程模型假定系统是由主机和设备构成的，它们分别具备本身独立的内存空间。Runtime负责设备内存的分配，回收，拷贝以及在主机和设备间传输数据的工做。

        设备内存能够有两种分配方式：线性内存或者CUDA数组

        CUDA数组是一块不透明的内存空间，它主要被优化用于纹理存取。

        线性内存空间与平时咱们访问的内存相似，对于计算能力1.x的设备来讲，它存在于一个32位的地址空间。对于更高计算能力的设备而言，它存在于一个40位的地址空间中。所以，单独分配的实体可使用指针来相互应用。

        咱们一般使用cudaMalloc()函数分配线性内存空间，使用cudaFree()函数释放线性内存空间，使用cudaMemcpy()函数在主机和设备之间传输数据。下面是CUDA Vector Add代码示例的一些片断：

 // Device code
 __global__ void VecAdd(float *A, float *B, float *C, int N) {
     int i = blockDim.x * blockIdx.x + threadIdx.x;
     if (i < N)
         C[i] = A[i] + B[i];
 }

 // Host code
 int main() {
     int N = ...;
     size_t size = N * sizeof(float);

     // Allocate input vectors h_A and h_B in host memory
     float *h_A = (float*)malloc(size);
     float *h_B = (float*)malloc(size);

     // Initialize input vectors
     ...

     // Allocate vectors in device memory
     float *d_A, *d_B, *d_C;
     cudaMalloc(&d_A, size);
     cudaMalloc(&d_B, size);
     cudaMalloc(&d_C, size);

     // Copy vectors from host memory to device memory
     cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
     cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);

     // Invoke kernel
     int threadsPerBlock = 256;
     int blocksPerGrid = (N +threadsPerBlock - 1) / threadsPerBlock;
     VecAdd<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);

     // Copy result from device memory to host Memory
     cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);

     // Free device memory
     cudaFree(d_A);
     cudaFree(d_B);
     cudaFree(d_C);

     // Free host memory
     ...
 }

片断展现了设备内存的分配，传输以及回收过程。

        除了上面展现的方法，咱们还可使用cudaMallocPitch()和cudaMalloc3D()函数来分配线性内存。这些函数可以确保分配的内存知足设备内存访问的对齐要求，对于行地址的访问以及多维数组间的数据传输提供高性能保证，所以很是适合对于二维和三维数组内存空间的分配。下面的代码片断展现了分配和使用尺寸为width x height的二维数组的技术：

 // Host code
 int width = 64, height = 64;
 float *devPtr;
 size_t pitch;
 cudaMallocPitch(&devPtr, &pitch, width * sizeof(float), height);
 MyKernel<<<100, 512>>>(devPtr, pitch, width, height);

 // Device code
 __global__ void MyKernel(float* devPtr, size_t pitch, int width, int height) {
     for (int r = 0; r < height; ++r) {
         float* row = (float*)((char*)devPtr + r * pitch);
         for (int c = 0; c < width; ++c) {
             float element = row[c];
         }
     }
 }

        下面的代码片断展现了一个尺寸为width x height x depth的三维数组的分配和使用方法：

 // Host code
 int width = 64, height = 64, depth = 64;
 cudaExtent extent = make_cudaExtent(width * sizeof(float), height, depth);
 cudaPitchedPtr devPitchedPtr;
 cudaMalloc3D(&devPitchedPtr, extent);
 MyKernel<<<100, 512>>>(devPitchedPtr, width, height, depth);

 // Device code
 __global__ void MyKernel(cudaPitchedPtr devPitchedPtr, int width, int height, int depth) {
     char* devPtr = devPitchedPtr.ptr;
     size_t pitch = devPitchedPtr.pitch;
     size_t slicePitch = pitch * height;
     for (int z = 0; z < depth; ++z) {
         char* slice = devPtr + z * slicePitch;
         for (int y = 0; y < height; ++y) {
             float* row = (float*)(slice + y * pitch);
             for (int x = 0; x < width; ++x)
                 float element = row[x];
         }
     }
 }

更多详细的内容请查阅参考手册。

        下面的代码示例展现了多种使用Runtime API访问全局变量的技术：

 __constant__ float constData[256];
 float data[256];
 cudaMemcpyToSymbol(constData, data, sizeof(data));
 cudaMemcpyFromSymbol(data, constData, sizeof(data));

 __device__ float devData;
 float value = 3.14f;
 cudaMemcpyToSymbol(devData, &value, sizeof(float));

 __device__ float* devPointer;
 float* ptr;
 cudaMalloc(&ptr, 256 * sizeof(float));
 cudaMemcpyToSymbol(devPointer, &ptr, sizeof(ptr));

使用cudaGetSymbolAddress()函数能够得到被声明存储在全局内存中的变量地址。为了得到分配内存的大小，可使用cudaGetSymbolSize()函数。

4.3.3 共享内存（Shared Memory）

        共享内存是CUDA设备中很是重要的一个存储区域，有效地使用共享内存能够充分利用CUDA设备的潜能，极大提高程序性能。那么，共享内存有哪些特色呢？

        一、共享内存（shared Memory）是集成在GPU处理器芯片上的（on-chip），所以相比于存在于显存颗粒中的全局内存（global Memory）和本地内存（local Memory），它具备更高的传输带宽，通常状况下，共享内存的带宽大约是全局内存带宽的7-10倍。

        二、共享内存的容量很小。根据NVIDIA官方文档的说法，在计算能力1.x的设备中，每个流多处理器（Streaming Multiprocessor）上的共享内存容量为16KB。对于计算能力2.x、3.0及3.5的设备该参数为48KB。所以共享内存是稀有资源。

        三、共享内存在物理上被划分为不少块，每一块被称为一个存储体（bank）。在同一时刻，CUDA设备能够同时访问多个存储体。所以，若是一次针对共享内存的访存操做须要读取n个地址，而这n个地址刚好分布在n个不一样的存储体（bank）中，那么只须要一个存取周期就能够完成n个地址的访存任务了。对于计算能力1.x的设备，共享内存被平均划分为16个存储体。而对于计算能力2.x、3.0及3.5的设备此参数为32。在共享内存中，相邻两块32bit的数据分别属于相邻的两个存储体。存储体每两个时钟周期能够传输32位数据。

        四、共享内存既能够静态分配，也能够动态分配。

        从共享内存的这些特色中咱们能够看出，它实际上至关于一个程序员能够操控的缓存（cache），下面，咱们使用矩阵乘法的例子来讲明如何有效使用共享内存。

        首先，咱们使用最直观的方法来完成矩阵乘法C = A x B：读取A的每一行和B的每一列，顺次完成计算任务。矩阵乘法的示意图以下所示：

下面是矩阵乘法的CUDA C主要实现代码：

 // Matrices are stored in row-major order:
 // M(row, col) = *(M.elements + row * M.width + col)
 typedef struct {
     int width;
     int height;
     float *elements;
 } Matrix;

 // Thread block size
 #define BLOCK_SIZE 16

 // Forward declaration of the matrix multiplication kernel
 __global__ void MatMulKernel(const Matrix, const Matrix, Matrix);

 // Matrix multiplication - Host code
 // Matrix dimensions are assumed to be multiples of BLOCK_SIZE
 void MatMul(const Matrix A, const Matrix B, Matrix C) {
     // Load A and B to device memory
     Matrix d_A;
     d_A.width = A.width; d_A.height = A.height;
     size_t size = A.width * A.height * sizeof(float);
     cudaMalloc(&d_A.elements, size);
     cudaMemcpy(d_A.elements, A.elements, size, cudaMemcpyHostToDevice);
     Matrix d_B;
     d_B.width = B.width; d_B.height = B.height;
     size = B.width * B.height * sizeof(float);
     cudaMalloc(&d_B.elements, size);
     cudaMemcpy(d_B.elements, B.elements, size, cudaMemcpyHostToDevice);

     // Allocate C in device memory
     Matrix d_C;
     d_C.width = C.width; d_C.height = C.height;
     size = C.width * C.height * sizeof(float);
     cudaMalloc(&d_C.elements, size);

     // Invoke kernel
     dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
     dim3 dimGrid(B.width / dimBlock.x, A.height / dimBlock.y);
     MatMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C);

     // Read C from device memory
     cudaMemcpy(C.elements, d_c.elements, size, cudaMemcpyDeviceToHost);

     // Free device memory
     cudaFree(d_A.elements);
     cudaFree(d_B.elements);
     cudaFree(d_C.elements);
 }

 // Matrix multiplication kernel called by MatMul()
 __global__ void MatMulKernel(Matrix A, Matrix B, Matrix C) {
     // Each thread computes one element of C
     // by accumulating results into Cvalue
     float Cvalue = 0;
     int row  = blockIdx.y * blockDim.y + threadIdx.y;
     int col = blockIdx.x * blockDim.x + threadIdx.xl
     for (int e = 0; e < A.width; ++e)
         Cvalue += A.elements[row * A.width + e] * B.elements[e * B.width + col];
     C.elements[row * C.width + col] = Cvalue;
 }

能够看出，为了计算矩阵C的任何一个元素，程序都须要从全局内存（global memory）中得到矩阵A的一行和矩阵B的一列。所以，完成这一计算矩阵A被读取了B.width次，矩阵B被读取了A.height次。
        如今咱们来使用共享内存（shared memory）实现矩阵乘法。假设矩阵C能够被划分为若干个较小的子方阵Csub，咱们使用一个线程块（thread block）来负责某一子方阵的计算，线程块中的每个线程（thread）正好负责子方阵Csub中一个元素的计算。这样划分后，任何一个结果子方阵Csub'（尺寸为block_size * block_size）都是与该方阵具备相同行索引的尺寸为A.width * block_size的A的子矩阵Asub和与该方阵具备相同列索引的尺寸为block_size * B.height的B的子矩阵Bsub相乘所获得。

        为了匹配设备的计算资源，两个子矩阵Asub和Bsub被划分为尽量多的分离的维度为block_size的子方阵，Csub的值即是这些子矩阵相乘后相加所获得的结果。子矩阵乘法的执行顺序都是首先将它们从全局内存（global memory）拷贝到共享内存（shared memory）（线程块中的每个线程正好负责方阵一个元素的拷贝），而后由线程本身完成相应元素的计算任务，利用寄存器存储局部结果，最后将寄存器的内容与新获得的计算结果依此累加起来获得最终运算结果并将其传输到全局内存（global memory）中。

        经过使用这种分治的计算策略，共享内存获得了很好的利用，采用这种方案计算完成时全局内存中矩阵A被访问的次数为B.width / block_size，矩阵B被访问的次数为A.height / block_size，很明显，这为咱们节省了很是多的全局内存带宽。优化后的矩阵计算示意图以下所示：

        为了提高计算效率，咱们为类型Matrix增长了一个成员变量stride。__device__函数用来得到和设置子矩阵的元素。下面是优化后的代码：

 // Matrices are stored in row-major order;
 // M(row, col) = *(M.elements + row * M.stride + col)
 typedef struct {
     int width;
     int height;
     int stride;
     float* elements;
 } Matrix;

 // Get a matrix element
 __device__ float GetElement(const Matrix A, int row, int col) {
     return A.elements[row * A.stride + col];
 }

 // Set a matrix element
 __device__ void SetElement(Matrix A, int row, int col, float value) {
     A.elements[row * A.stride + col] = value;
 }

 // Get the BLOCK_SIZExBLOCK_SIZE sub-matrix Asub of A that is
 // located col sub-matrices to the right and row sub-matrices down
 // from the upper-left corner of A
 __device__ Matrix GetSubMatrix(Matrix A, int row, int col) {
     Matrix Asub;
     Asub.width = BLOCK_SIZE;
     Asub.height = BLOCK_SIZE;
     Asub.stride = A.stride;
     Asub.elements = &A.elements[A.stride * BLOCK_SIZE * row + BLOCK_SIZE * col];
     return Asub;
 }

 // Thread block size
 #define BLOCK_SIZE 16

 // Forward declaration of the matrix multiplication kernel
 __global__ void MatMulKernel(const Matrix, const Matrix, Matrix);

 // Matrix multiplication - Host code
 // Matrix dimensions are assumed to be multiples of BLOCK_SIZE
 void MatMul(const Matrix A, const Matrix B, Matrix C) {
     // Load A and B to device memory
     Matrix d_A;
     d_A.width = d_A.stride = A.width;
     d_A.height = A.height;
     size_t size = A.width * A.height * sizeof(float);
     cudaMalloc(&d_A.elements, size);
     cudaMemcpy(d_A.elements, A.elements, size, cudaMemcpyHostToDevice);
     Matrix d_B;
     d_B.width = d_B.stride = B.width;
     d_B.height = B.height;
     size = B.width * B.height * sizeof(float);
     cudaMalloc(&d_B.elements, size);
     cudaMemcpy(d_B.elements, B.elements, size, cudaMemcpyHostToDevice);

     // Allocate C in device memory
     Matrix d_C;
     d_C.width = d_C.stride = C.width;
     d_C.height = C.height;
     size = C.width * C.height * sizeof(float);
     cudaMalloc(&d_C.elements, size);

     // Invoke kernel
     dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
     dim3 dimGrid(B.width / dimBlock.x, A.height / dimBlock.y);
     MatMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C);

     // Read C from device memory
     cudaMemcpy(C.elements, d_C.elements, size, cudaMemcpyDeviceToHost);

     // Free device memory
     cudaFree(d_A.elements);
     cudaFree(d_B.elements);
     cudaFree(d_C.elements);
 }

 // Matrix multiplication kernel called by MatMul()
 __global__ void MatMulKernel(Matrix A, Matrix B, Matrix C) {
     // Block row and column
     int blockRow = blockIdx.y;
     int blockCol = blockIdx.x;

     // Each thread block computes one sub-matrix Csub of C
     Matrix Csub = GetSubMatrix(C, blockRow, blockCol);

     // Each thread computes one element of Csub
     // by accumulating results into Cvalue
     float Cvalue = 0;

     // Thread row and column within Csub
     int row = threadIdx.y;
     int col = threadIdx.x;

     // Look over all the sub-matrices of A and B that are required to compute Csub
     // Multiply each pair of sub-matrices together and accumulate the results
     for (int m = 0; m < (A.width / BLOCK_SIZE); ++m) {
         // Get sub-matrix Asub of A
         Matrix Asub = GetSubMatrix(A, blockRow, m);

         // Get sub-matrix Bsub of B
         Matrix Bsub = GetSubMatrix(B, m, blockCol);

         // Shared memory used to store Asub and Bsub respectively
         __shared__ float As[BLOCK_SIZE][BLOCK_SIZE];
         __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];

         // Load Asub and Bsub from device memory to shared memory
         // Each thread loads one element of each sub-matrix
         As[row][col] = GetElement(Asub, row, col);
         Bs[row][col] = GetElement(Bsub, row, col);

         // Synchronize to make sure the sub-matrices are loaded
         // before starting the computation
         __syncthreads();

         // Multiply Asub and Bsub together
         for (int e = 0; e < BLOCK_SIZE; ++e)
             Cvalue += As[row][e] * Bs[e][col];

         // Synchronize to make sure that the preceding computation is done before
         // loading two new sub-matrices of A and B in the next iteration
         __syncthreads();
     }

     // Write Csub to device memory
     // Each thread writes one element
     SetElement(Csub, row, col, Cvalue);
 }

异步并行执行

主机和设备间并行执行

        为了支持主机和设备的并行执行，CUDA提供了一些异步函数。异步是指设备在完成功能执行以前就将控制权交还给主机线程，以便主机线程继续执行。这些函数有：

        一、内涵启动（Kernel Launches）；

        二、同一设备内存中两个地址块之间的数据传输；

        三、从主机到设备的小于或等于64KB的一次数据块传输；

        四、使用Async前缀函数进行的数据传输；

        五、内存置位函数调用（Memory set functions calls）。

程序员能够经过设置环境变量CUDA_LAUNCH_BLOCKING来开启或者关闭内核启动（Kernel Launch）的异步功能。可是这一方法仅限于调试，在任何产品代码中不该当关闭异步内核启动。

        内核启动在下面这些状况下则是同步的：

        一、应用程序经过调试器或者内存检查器运行在计算能力为1.x的设备上。

        二、硬件计数器信息正被性能分析器收集。

将内核启动与数据传输重叠起来

        对于一些计算能力等于或高于1.1的设备，它们能够将内核启动任务和锁页内存到设备内存的数据传输任务并行执行。应用程序能够检查设备属性中的asyncEngineCount项来肯定设备是否支持这一功能。当该项值大于0时表明设备支持这一层次的并行。对于计算能力1.x的设备，该功能不支持经过cudaMallocPitch()函数分配的CUDA数组或2D数组。

并行内核执行

        一些计算能力2.x或更高的设备能够同时并行执行多个内核函数。应用程序能够检查设备属性中的concurrentKernels项来肯定设备是否支持这一功能，值为1表明支持。运算能力3.5的设备在同一时刻可以并行执行的最大内核函数数量为32，运算能力小于3.5的硬件则最多支持同时启动16个内核函数的执行。同时须要注意的是，在一个CUDA上下文中的内核函数不能与另外一个CUDA上下文中的内核函数同时执行。使用不少纹理内存或者大量本地内存的内核函数也极可能没法与其它内核函数并行执行。

并行数据传输

        一些计算能力为2.x或更高的设备能够将锁页内存到设备内存的数据传输和设备内存到锁页内存的数据传输并行执行。应用程序可检查设备属性中的asyncEngineCount项来肯定这一功能的支持程度，等于2时表示支持。

流（Streams）

        应用程序经过流来管理并行。一个流是一个顺次执行的命令序列。不一样的流之间并行执行，没有固定的执行顺序。

一、流的建立与销毁

        定义一个流的过程一般包括：建立一个流对象，而后指定它为内核启动或者主机设备间数据传输的流参数。下面的一段代码建立了两个流而且在锁页内存中分配了一块float类型的数组hostPtr：

 cudaStream_t stream[2];
 for (int i = 0; i < 2; ++i)
     cudaStreamCreate(&stream[i]);
 float *hostPtr;
 cudaMallocHost(&hostPtr, 2 * size);

下面的代码定义了每个流的行为：从主机端拷贝数据到设备端，内核启动，从设备端拷贝数据到主机端：

 for (int i = 0; i < 2; ++i) {
     cudaMemcpyAsync(inputDevPtr + i * size, hostPtr + i * size, size, cudaMemcpyHostToDevice, stream[i]);
     MyKernel<<<100, 512, 0, stream[i]>>>(outputDevPtr + i * size, inputDevPtr + i * size, size);
     cudaMemcpyAsync(hostPtr + i * size, outputDevPtr + i * size, size, cudaMemcpyDeviceToHost, stream[i]);
 }

这部分代码中有一点须要注意：为了并行化数据拷贝和内核执行，主机端内存必须分配为锁页（page-locked）内存。

要销毁一个流须要调用函数cudaStreamDestroy()

 for (int i = 0; i < 2; ++i)
     cudaStreamDestroy(stream[i]);

cudaStreamDestroy()函数等待以前流中的指令序列运行完成，而后销毁指定流，将控制权返还给主机端。

二、默认流（Default stream）

        在内核启动或者数据拷贝过程当中若是不指定流，或者设置流参数为0，则相应的指令将会运行在默认流上，它们也所以而顺次执行。

三、明同步（Explicit Synchronization）

        在CUDA中有不少种方式能够用来同步流的执行：

        cudaDeviceSynchronize()函数使得主机端线程阻塞直到全部流中的指令执行完成。

       cudaStreamSynchronize()函数将一个流对象做为输入参数，用以等待指定流中的全部指令执行完成。

       cudaStreamWaitEvent()函数将一个流对象和一个事件做为输入参数，它将延迟该函数调用后在指定流中全部新加入的命令的执行直到指定的事件完成为止。流参数能够为0，在该情形下全部流中的任何新加入的指令都必须等待指定事件的发生，而后才能够执行。

       cudaStreamQuery()函数为应用程序提供了一个检测指定流中以前指令是否执行完成的方法。

        为了不同步带来的性能降低，全部上述同步函数最好用于计时目的或者分离错误的内核执行或数据拷贝。

四、暗同步（Implicit Synchronization）

        若是任何一个流中正在执行如下操做，那么其它流是不能与其并行运行的：

        a. 分配锁页内存空间

        b. 设备内存分配

        c. 设备内存置位

        d. 同一设备两个不一样地址间正在进行数据拷贝

        e. 默认流中有指令正在执行

        f. L1/shared内存配置的转换

        对于支持并行内核执行而且计算能力3.0或如下的设备来讲，任何一个须要检查依赖性以肯定流内核启动是否完成的操做：

        a. 只有当前CUDA上下文中全部流中全部以前的内核启动以后才可以启动执行。

        b. 将会阻塞全部当前CUDA上下文中的任意流中新加入的内核调用直到内核检查完成。

        须要进行依赖性检查的操做包括执行检查的内核启动所在流中的其它指令以及任何在该流上对cudaStreamQuery()函数的调用。所以，应用程序能够遵守如下指导原则来提高潜在并行性：

        （1）全部非依赖操做应当比依赖性操做提早进行

        （2）任何类型的同步越迟越好

五、重叠行为（Overlapping Behavior）

        两个流间重叠行为的数量取决于如下几个因素：

        （1）每一个流中命令发出的次序

        （2）设备是否支持内核启动与数据传输并行

        （3）设备是否支持多内核并行启动

        （4）设备是否支持多数据传输并行

        例如，在不支持并行数据传输的设备上，“流的建立与销毁”章节中代码样例中的操做就不能并行，由于在stream[0]中发出设备端到主机端的数据拷贝后，stream[1]又发出主机端到设备端的数据拷贝命令，这两个命令式不能重叠执行的。假设设备支持数据传输与内核启动并行，那么以下代码：

 for (int i = 0; i < 2; ++i)
     cudaMemcpyAsync(inputDevPtr + i * size, hostPtr + i * size, size, cudaMemcpyHostToDevice, stream[i]);
 for (int i = 0; i < 2; ++i)
     MyKernel<<<100, 512, 0, stream[i]>>>(outputDevPtr + i * size, inputDevPtr + i * size, size);
 for (int i = 0; i < 2; ++i)
     cudaMemcpyAsync(hostPtr + i * size, outputDevPtr + i * size, size, cudaMemcpyDeviceToHost, stream[i]);

可将stream[0]的内核启动和stream[1]从主机端到设备端的数据拷贝重叠起来并行执行。

六、回调函数

        CUDA运行时提供了cudaStreamAddCallback()函数以在流中的任意位置插入一个回调函数点。回调函数运行于主机端，若是在默认流中插入回调函数，那么它将等待全部其它流中的命令执行完成以后才会开始执行。

        下面的代码展现了回调函数技术的应用：

 void CUDART_CB MyCallback(cudaStream_t stream, cudaError_t status, void **data) {
     printf("Inside callback %d\n", (int)data);
 }
 ...
 for (int i = 0; i < 2; ++i) {
     cudaMemcpyAsync(devPtrIn[i], hostPtr[i], size, cudaMemcpyHostToDevice, stream[i]);
     MyKernel<<<100, 512, 0, stream[i]>>>(devPtrOut[i], devPtrIn[i], size);
     cudaMemcpyAsync(hostPtr[i], devPtrOut[i], size, cudaMemcpyDeviceToHost, stream[i]);
     cudaStreamAddCallback(stream[i], MyCallback, (void**)i, 0);
 }

上面的代码定义了两个流的操做，每一个流都完成一次主机端到设备端的数据拷贝，一次内核启动，一次设备端到主机端的数据拷贝，最后增长了一个加入回调函数的操做。当设备端代码运行到回调函数点的时候，设备将控制权交还给主机端，主机端运行完成之后再将控制权返还给设备端，而后设备端继续运行。

        值得注意的是，在一个回调函数中，必定不能进行任何CUDA API的调用，直接的或者间接的都是不能够的。