随着多核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)的线程来说：
1、如果线程块（block）是一维的，则线程ID = x
2、如果线程块是二维的，假设块尺寸为（Dx，Dy），那么线程ID = x + y * Dx
3、如果线程块是三维的，设其尺寸为（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.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写出的代码需要在两种体系结构完全不同的设备上执行：1、CPU；2、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 兼容性
1、二进制兼容性
二进制代码是设备相关的，使用NVCC编译器编译时，若指定-code选项，则会编译产生目标设备的二进制cubin对象。例如，编译时使用-code=sm_13会产生适用于计算能力1.3的二进制代码。二进制代码在CUDA计算设备上具有小版本的向前兼容性，但是在大版本上不具备兼容性。也就是说，对于计算能力X.y的硬件，使用-code=sm_Xy编译后，程序能够运行于计算能力X.z（其中z>=y）的硬件上，但不能运行在计算能力M.n（M!=X）的硬件上。
2、PTX代码兼容性
不同计算能力的设备所支持的PTX指令条数是不同的，一些PTX指令只在拥有较高计算能力的设备上被支持。例如，全局内存（global Memory）的原子操作指令只能用于计算能力不小于1.1的设备；双精度浮点运算指令只能用于计算能力不小于1.3的设备。在将C语言编译为PTX代码时，NVCC使用-arch编译选项指定PTX代码目标设备的计算能力。因此，要想使用双精度运算，编译时必须使用选项-arch=sm_13（或使用更高的计算能力），否则NVCC会自动将双精度操作降级为单精度操作。
为某一特定设备产生的PTX代码，在运行时总是能够被具有更高计算能力的设备JIT编译为可执行的二进制代码。
3、应用程序兼容性
执行CUDA程序有两种方式，一种是直接加载编译好的CUDA二进制代码运行，另一种是首先加载程序中的PTX代码，再执行JIT编译得到二进制的设备可执行文件，然后运行。特别需要注意的是，为了让程序运行具有更高计算能力的未来设备上，必须让程序加载PTX代码。
事实上，在一个CUDA C程序中可以嵌入不止一个版本的PTX/二进制代码。那么，具体执行时哪一个版本的PTX或者二进制代码会得到执行呢？答案是：最兼容的那个版本。例如编译一个名为x.cu的CUDA源代码：
将会产生兼容计算能力1.1硬件的二进制代码（第一排的-gencode选项）以及兼容计算能力1.1设备的PTX和二进制代码，这些代码都将会嵌入到编译后的目标文件中。
主机端将会产生一些额外的代码，在程序运行时，这些代码会自动决定装载哪一个版本的代码来执行。对于上面的例子：
- 计算能力1.0的设备运行该程序将会装载1.0版本的二进制代码
- 计算能力1.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的简化形式。更多详尽的内容请参阅该手册。
4、C/C++兼容性
NVCC编译器前端使用C++语法啊规则来处理CUDA源文件。在主机端，CUDA支持完整的C++语法；而在设备端，只有部分C++语法是被支持的。这方面更为详尽的讨论请参见《CUDA C程序设计指南》的C/C++语言支持章节。
5、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 CUDA C Runtime
4.3.3 共享内存（Shared Memory）
共享内存是CUDA设备中非常重要的一个存储区域，有效地使用共享内存可以充分利用CUDA设备的潜能，极大提升程序性能。那么，共享内存有哪些特点呢？
1、共享内存（shared Memory）是集成在GPU处理器芯片上的（on-chip），因此相比于存在于显存颗粒中的全局内存（global Memory）和本地内存（local Memory），它具有更高的传输带宽，一般情况下，共享内存的带宽大约是全局内存带宽的7-10倍。
2、共享内存的容量很小。根据NVIDIA官方文档的说法，在计算能力1.x的设备中，每一个流多处理器（Streaming Multiprocessor）上的共享内存容量为16KB。对于计算能力2.x、3.0及3.5的设备该参数为48KB。因此共享内存是稀有资源。
3、共享内存在物理上被划分为很多块，每一块被称为一个存储体（bank）。在同一时刻，CUDA设备可以同时访问多个存储体。因此，如果一次针对共享内存的访存操作需要读取n个地址，而这n个地址恰好分布在n个不同的存储体（bank）中，那么只需要一个存取周期就可以完成n个地址的访存任务了。对于计算能力1.x的设备，共享内存被平均划分为16个存储体。而对于计算能力2.x、3.0及3.5的设备此参数为32。在共享内存中，相邻两块32bit的数据分别属于相邻的两个存储体。存储体每两个时钟周期可以传输32位数据。
4、共享内存既可以静态分配，也可以动态分配。
从共享内存的这些特点中我们可以看出，它实际上相当于一个程序员可以操控的缓存（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可以被划分为若干个较小的子方阵C
sub，我们使用一个线程块（thread block）来负责某一子方阵的计算，线程块中的每一个线程（thread）正好负责子方阵C
sub中一个元素的计算。这样划分后，任何一个结果子方阵C
sub'（尺寸为block_size * block_size）都是与该方阵具有相同行索引的尺寸为A.width * block_size的A的子矩阵A
sub和与该方阵具有相同列索引的尺寸为block_size * B.height的B的子矩阵B
sub相乘所得到。
为了匹配设备的计算资源，两个子矩阵Asub和Bsub被划分为尽可能多的分离的维度为block_size的子方阵，Csub的值便是这些子矩阵相乘后相加所得到的结果。子矩阵乘法的执行顺序都是首先将它们从全局内存（global memory）拷贝到共享内存（shared memory）（线程块中的每一个线程正好负责方阵一个元素的拷贝），然后由线程自己完成相应元素的计算任务，利用寄存器存储局部结果，最后将寄存器的内容与新得到的计算结果依此累加起来得到最终运算结果并将其传输到全局内存（global memory）中。
通过使用这种分治的计算策略，共享内存得到了很好的利用，采用这种方案计算完成时全局内存中矩阵A被访问的次数为B.width / block_size，矩阵B被访问的次数为A.height / block_size，很明显，这为我们节省了非常多的全局内存带宽。优化后的矩阵计算示意图如下所示：
为了提升计算效率，我们为类型Matrix增加了一个成员变量stride。__device__函数用来获得和设置子矩阵的元素。下面是优化后的代码：
http://blog.csdn.net/csgxy123/article/details/9704461
// 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提供了一些异步函数。异步是指设备在完成功能执行之前就将控制权交还给主机线程，以便主机线程继续执行。这些函数有：
1、内涵启动（Kernel Launches）；
2、同一设备内存中两个地址块之间的数据传输；
3、从主机到设备的小于或等于64KB的一次数据块传输；
4、使用Async前缀函数进行的数据传输；
5、内存置位函数调用（Memory set functions calls）。
程序员可以通过设置环境变量CUDA_LAUNCH_BLOCKING来开启或者关闭内核启动（Kernel Launch）的异步功能。但是这一方法仅限于调试，在任何产品代码中不应当关闭异步内核启动。
内核启动在下面这些情况下则是同步的：
1、应用程序通过调试器或者内存检查器运行在计算能力为1.x的设备上。
2、硬件计数器信息正被性能分析器收集。
将内核启动与数据传输重叠起来
对于一些计算能力等于或高于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）
应用程序通过流来管理并行。一个流是一个顺次执行的命令序列。不同的流之间并行执行，没有固定的执行顺序。
1、流的创建与销毁
定义一个流的过程通常包括：创建一个流对象，然后指定它为内核启动或者主机设备间数据传输的流参数。下面的一段代码创建了两个流并且在锁页内存中分配了一块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]);2、默认流（Default stream）
在内核启动或者数据拷贝过程中如果不指定流，或者设置流参数为0，则相应的指令将会运行在默认流上，它们也因此而顺次执行。
3、明同步（Explicit Synchronization）
在CUDA中有很多种方式可以用来同步流的执行：
cudaDeviceSynchronize()函数使得主机端线程阻塞直到所有流中的指令执行完成。
cudaStreamSynchronize()函数将一个流对象作为输入参数，用以等待指定流中的所有指令执行完成。
cudaStreamWaitEvent()函数将一个流对象和一个事件作为输入参数，它将延迟该函数调用后在指定流中所有新加入的命令的执行直到指定的事件完成为止。流参数可以为0，在该情形下所有流中的任何新加入的指令都必须等待指定事件的发生，然后才可以执行。
cudaStreamQuery()函数为应用程序提供了一个检测指定流中之前指令是否执行完成的方法。
为了避免同步带来的性能下降，所有上述同步函数最好用于计时目的或者分离错误的内核执行或数据拷贝。
4、暗同步（Implicit Synchronization）
如果任何一个流中正在执行以下操作，那么其它流是不能与其并行运行的：
a. 分配锁页内存空间
b. 设备内存分配
c. 设备内存置位
d. 同一设备两个不同地址间正在进行数据拷贝
e. 默认流中有指令正在执行
f. L1/shared内存配置的转换
对于支持并行内核执行并且计算能力3.0或以下的设备来说，任何一个需要检查依赖性以确定流内核启动是否完成的操作：
a. 只有当前CUDA上下文中所有流中所有之前的内核启动之后才能够启动执行。
b. 将会阻塞所有当前CUDA上下文中的任意流中新加入的内核调用直到内核检查完成。
需要进行依赖性检查的操作包括执行检查的内核启动所在流中的其它指令以及任何在该流上对cudaStreamQuery()函数的调用。因此，应用程序可以遵照以下指导原则来提升潜在并行性：
（1）所有非依赖操作应当比依赖性操作提前进行
（2）任何类型的同步越迟越好
5、重叠行为（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]);6、回调函数
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的调用，直接的或者间接的都是不可以的。
http://blog.csdn.net/csgxy123/article/details/9704461
扫一扫
6122
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
