并行计算与计算机架构

并行计算

并行计算涉及两个不同的技术领域：

• 计算机架构（硬件）
• 并行程序设计（软件）

一个生产工具，一个用工具产生各种不同应用。
硬件主要的目标就是为软件提供更快的计算速度，更低的性能功耗比，硬件结构上支持更快的并行。
软件的主要目的是使用当前的硬件压榨出最高的性能，给应用提供更稳定快速的计算结果。
我们传统的计算机结构一般是哈佛体系结构（后来演变出冯·诺依曼结构）主要分成三部分：

• 内存（指令内存，数据内存）
• 中央处理单元（控制单元和算术逻辑单元）
• 输入、输出接口

写并行和串行的最大区别就是，写串行程序可能不需要学习不同的硬件平台，但是写并行程序就需要对硬件有一定的了解。

并行性

写并行程序主要是分解任务，我们一般把一个程序看成是指令和数据的组合，当然并行也可以分为这两种：

• 指令并行
• 数据并行

我们的任务更加关注数据并行，所以我们的主要任务是分析数据的相关性，哪些可以并行，哪些不能不行。我们研究的是大规模数据计算，计算过程比较单一（不同的数据基本用相同的计算过程）但是数据非常多，所以我们主要是数据并行，分析好数据的相关性，决定了我们的程序设计。CUDA 非常适合数据并行程序设计，第一步就是把数据依据线程进行划分

1. 块划分，把一整块数据切成小块，每个小块随机的划分给一个线程，每个块的执行顺序随机

2. 周期划分，线程按照顺序处理相邻的数据块，每个线程处理多个数据块，比如我们有五个线程，线程 1 执行块 1，线程 2 执行块 2……线程 5 执行块 5，线程 1 执行块 6

下面是示意图，注意颜色相同的块使用的同一个线程，从执行顺序上看如下：

下面是数据集的划分上看：

不同的数据划分严重影响程序性能，所以针对不同的问题和不同计算机结构，我们要通过和理论和试验共同来决定最终最优的数据划分。

计算机架构

Flynn’s Taxonomy

划分不同计算机结构的方法有很多，广泛使用的一种被称为佛林分类法 Flynn’s Taxonomy，它根据指令和数据进入 CPU 的方式分类，分为以下四类：

分别以数据和指令进行分析：

• 单指令单数据 SISD（传统串行计算机，386）
• 单指令多数据 SIMD（并行架构，比如向量机，所有核心指令唯一，但是数据不同，现在 CPU 基本都有这类的向量指令）
• 多指令单数据 MISD（少见，多个指令围殴一个数据）
• 多指令多数据 MIMD（并行架构，多核心，多指令，异步处理多个数据流，从而实现空间上的并行，MIMD 多数情况下包含 SIMD，就是 MIMD 有很多计算核，计算核支持 SIMD）

为了提高并行的计算能力，我们要从架构上实现下面这些性能提升：

• 降低延迟
• 提高带宽
• 提高吞吐量

延迟是指操作从开始到结束所需要的时间，一般用微秒计算，延迟越低越好。
带宽是单位时间内处理的数据量，一般用 MB/s 或者 GB/s 表示。
吞吐量是单位时间内成功处理的运算数量，一般用 gflops 来表示（十亿次浮点计算），吞吐量和延迟有一定关系，都是反应计算速度的，一个是时间除以运算次数，得到的是单位次数用的时间–延迟，一个是运算次数除以时间，得到的是单位时间执行次数–吞吐量。

根据内存划分

计算机架构也可以根据内存进行划分：

1. 分布式内存的多节点系统
2. 共享内存的多处理器系统

第一个更大，通常叫做集群，就是一个机房好多机箱，每个机箱都有内存处理器电源等一系列硬件，通过网络互动，这样组成的就是分布式。

第二个是单个主板有多个处理器，他们共享相同的主板上的内存，内存寻址空间相同，通过 PCIe 和内存互动。

多个处理器可以分多片处理器，和单片多核（众核 many-core)，也就是有些主板上挂了好多片处理器，也有的是一个主板上就一个处理器，但是这个处理器里面有几百个核。
GPU 就属于众核系统。当然现在 CPU 也都是多核的了，但是他们还是有很大区别的：

• CPU 适合执行复杂的逻辑，比如多分支，其核心比较重（复杂）
• GPU 适合执行简单的逻辑，大量的数据计算，其吞吐量更高，但是核心比较轻（结构简单）

异构计算与 CUDA

异构计算

异构计算，首先必须了解什么是异构，不同的计算机架构就是异构。
GPU 本来的任务是做图形图像的，也就是把数据处理成图形图像，图像有个特点就是并行度很高，基本上一定距离以外的像素点之间的计算是独立的，所以属于并行任务。
GPU 之前是不可编程的，或者说不对用户开放的，人家本来是做图形计算控制显示器的，之后 CUDA 诞生使得 GPU 可以实现编程。

x86 CPU+GPU 的这种异构应该是最常见的，也有 CPU+FPGA，CPU+DSP 等各种各样的组合，CPU+GPU 在每个笔记本或者台式机上都能找到。当然超级计算机大部分也采用异构计算的方式来提高吞吐量。异构架构虽然比传统的同构架构运算量更大，但是其应用复杂度更高，因为要在两个设备上进行计算，控制，传输，这些都需要人为干预，而同构的架构下，硬件部分自己完成控制，不需要人为设计。

异构架构

运行程序的时候，CPU 像是一个控制者，指挥 GPU 完成工作后进行汇总，和下一步工作安排，所以 CPU 我们可以把它看做一个指挥者，主机端，host，而完成大量计算的 GPU 是我们的计算设备，device。

上面这张图能大致反应 CPU 和 GPU 的架构不同。

• 左图：一个四核 CPU 一般有四个 ALU，ALU 是完成逻辑计算的核心，也是我们平时说四核八核的核，控制单元，缓存也在片上，DRAM 是内存，一般不在片上，CPU 通过总线访问内存。

• 右图：GPU，绿色小方块是 ALU，我们注意红色框内的部分 SM，这一组 ALU 共用一个 Control 单元和 Cache，这个部分相当于一个完整的多核 CPU，但是不同的是 ALU 多了，control 部分变小，可见计算能力提升了，控制能力减弱了

  所以对于控制（逻辑）复杂的程序，一个 GPU 的 SM 是没办法和 CPU 比较的，但是对于逻辑简单，数据量大的任务，GPU 更高效，并且一个 GPU 有好多个 SM，而且越来越多。

CPU 和 GPU 之间通过 PCIe 总线连接，用于传递指令和数据，这部分也是后面要讨论的性能瓶颈之一。
一个异构应用包含两种以上架构，所以代码也包括不止一部分：

• 主机代码
• 设备代码

主机代码在主机端运行，被编译成主机架构的机器码，设备端的在设备上执行，被编译成设备架构的机器码，所以主机端的机器码和设备端的机器码是隔离的，自己执行自己的，没办法交换执行。
主机端代码主要是控制设备，完成数据传输等控制类工作，设备端主要的任务就是计算。
因为当没有 GPU 的时候 CPU 也能完成这些计算，只是速度会慢很多，所以可以把 GPU 看成 CPU 的一个加速设备。
NVIDIA 目前的计算平台（不是架构）有：

• Tegra
• Geforce
• Quadro
• Tesla

每个平台针对不同的应用场景，比如 Tegra 用于嵌入式，Geforce 是我们平时打游戏用到。

上面是根据应用场景分类的几种平台。

衡量 GPU 计算能力的主要靠下面两种容量特征：

• CUDA 核心数量（越多越好）
• 内存大小（越大越好）

相应的也有计算能力的性能指标:

• 峰值计算能力
• 内存带宽

nvidia 自己有一套描述 GPU 计算能力的代码，其名字就是“计算能力”，主要区分不同的架构，早期架构的计算能力不一定比新架构的计算能力强

这里的 Tesla 架构，与上面的 Tesla 平台不同，不要混淆，一个是平台名字，一个是架构名字

范例

CPU 和 GPU 相互配合，各有所长，各有所短

• 低并行逻辑复杂的程序适合用 CPU
• 高并行逻辑简单的大数据计算适合 GPU

一个程序可以进行如下分解，串行部分和并行部分：

CPU 和 GPU 线程的区别：

1. CPU 线程是重量级实体，操作系统交替执行线程，线程上下文切换花销很大
2. GPU 线程是轻量级的，GPU 应用一般包含成千上万的线程，多数在排队状态，线程之间切换基本没有开销。
3. CPU 的核被设计用来尽可能减少一个或两个线程运行时间的延迟，而 GPU 核则是大量线程，最大幅度提高吞吐量

CUDA：一种异构计算平台

CUDA 平台不是单单指软件或者硬件，而是建立在 Nvidia GPU 上的一整套平台，并扩展出多语言支持

CUDA C 是标准 ANSI C 语言的扩展，扩展出一些语法和关键字来编写设备端代码，而且 CUDA 库本身提供了大量 API 来操作设备完成计算。

对于 API 也有两种不同的层次，一种相对较高层，一种相对底层。

• CUDA 驱动 API
• CUDA 运行时 API

驱动 API 是低级的 API，使用相对困难，运行时 API 是高级 API，使用简单，其实现基于驱动 API。
这两种 API 是互斥的，也就是你只能用一个，两者之间的函数不可以混合调用，只能用其中的一个库。

一个 CUDA 应用通常可以分解为两部分，

• CPU 主机端代码
• GPU 设备端代码

CUDA nvcc 编译器会自动分离你代码里面的不同部分，如图中主机代码用 C 写成，使用本地的 C 语言编译器编译，设备端代码，也就是核函数，用 CUDA C 编写，通过 nvcc 编译，链接阶段，在内核程序调用或者明显的 GPU 设备操作时，添加运行时库。

注意：核函数是我们后面主要接触的一段代码，就是设备上执行的程序段

nvcc 是从 LLVM 开源编译系统为基础开发的。

CUDA 工具箱提供编译器，数学库，调试优化等工具，当然 CUDA 的文档是相当完善的，可以去查阅，当然在我们基本了解基础结构的情况下，直接上来看文档会变得机械化。

Hello World

1
2
3
4
5
6
7
8
9
10

#include <stdio.h>
__global__ void hello_world(void) {
    printf("GPU: Hello world!\n");
}
int main(int argc, char **argv) {
    printf("CPU: Hello world!\n");
    hello_world<<<1, 10>>>();
    cudaDeviceReset(); // if no this line ,it can not output hello world from gpu
    return 0;
}

简单介绍其中几个关键字

1

__global__

是告诉编译器这个是个可以在设备上执行的核函数

1

hello_world<<<1,10>>>();

这句话 C 语言中没有<<<>>>是对设备进行配置的参数，也是 CUDA 扩展出来的部分。

1

cudaDeviceReset();

这句话如果没有，则不能正常的运行，因为这句话包含了隐式同步，GPU 和 CPU 执行程序是异步的，核函数调用后成立刻会到主机线程继续，而不管 GPU 端核函数是否执行完毕，所以上面的程序就是 GPU 刚开始执行，CPU 已经退出程序了，所以我们要等 GPU 执行完了，再退出主机线程。
一般 CUDA 程序分成下面这些步骤：

1. 分配 GPU 内存
2. 拷贝内存到设备
3. 调用 CUDA 内核函数来执行计算
4. 把计算完成数据拷贝回主机端
5. 内存销毁

CPU 与 GPU 的编程主要区别在于对 GPU 架构的熟悉程度，理解机器的结构是对编程效率影响非常大的一部分，了解你的机器，才能写出更优美的代码，而目前计算设备的架构决定了局部性将会严重影响效率。

数据局部性分两种

• 空间局部性
• 时间局部性

这个两个性质告诉我们，当一个数据被使用，其附近的数据将会很快被使用，当一个数据刚被使用，则随着时间继续其被再次使用的可能性降低，数据可能被重复使用。

CUDA 中有两个模型是决定性能的：

• 内存层次结构
• 线程层次结构

CUDA C 写核函数的时候我们只写一小段串行代码，但是这段代码被成千上万的线程执行，所有线程执行的代码都是相同的，CUDA 编程模型提供了一个层次化的组织线程，直接影响 GPU 上的执行顺序。

CUDA 抽象了硬件实现：

1. 线程组的层次结构
2. 内存的层次结构
3. 障碍同步

这些都是我们后面要研究的，线程，内存是主要研究的对象，我们能用到的工具相当丰富，NVIDIA 为我们提供了：

• Nvidia Nsight 集成开发环境
• CUDA-GDB 命令行调试器
• 性能分析可视化工具
• CUDA-MEMCHECK 工具
• GPU 设备管理工具

CUDA 编程模型概述

CUDA 编程模型概述

CUDA 编程模型为应用和硬件设备之间的桥梁，所以 CUDA C 是编译型语言，不是解释型语言，OpenCL 就有点类似于解释型语言，通过编译器和链接，给操作系统执行（操作系统包括 GPU 在内的系统），下面的结构图片能形象的表现他们之间的关系：

其中 Communication Abstraction 是编程模型和编译器，库函数之间的分界线。
GPU 中大致可以分为：

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

等几个关键部分。
以上这些理论同时也适用于其他非 CPU+GPU 异构的组合。
下面我们会说两个我们 GPU 架构下特有几个功能：

• 通过组织层次结构在 GPU 上组织线程的方法
• 通过组织层次结构在 GPU 上组织内存的方法

从宏观上我们可以从以下几个环节完成 CUDA 应用开发：

1. 领域层
2. 逻辑层
3. 硬件层

第一步就是在领域层（也就是你所要解决问题的条件）分析数据和函数，以便在并行运行环境中能正确，高效地解决问题。
当分析设计完程序就进入了编程阶段，我们关注点应转向如何组织并发进程，这个阶段要从逻辑层面思考。
CUDA 模型主要的一个功能就是线程层结构抽象的概念，以允许控制线程行为。这个抽象为并行编程提供了良好的可扩展性（这个扩展性后面有提到，就是一个 CUDA 程序可以在不同的 GPU 机器上运行，即使计算能力不同）。
在硬件层上，通过理解线程如何映射到机器上，能充分帮助我们提高性能。

CUDA 编程结构

一个异构环境，通常有多个 CPU 多个 GPU，他们都通过 PCIe 总线相互通信，也是通过 PCIe 总线分隔开的。所以我们要区分一下两种设备的内存：

• 主机：CPU 及其内存
• 设备：GPU 及其内存

注意这两个内存从硬件到软件都是隔离的（CUDA6.0 以后支持统一寻址），我们目前先不研究统一寻址，我们现在还是用内存来回拷贝的方法来编写调试程序，以巩固大家对两个内存隔离这个事实的理解。

一个完整的 CUDA 应用可能的执行顺序如下图：

从 host 的串行到调用核函数（核函数被调用后控制马上归还主机线程，也就是在第一个并行代码执行时，很有可能第二段 host 代码已经开始同步执行了)。

我们接下来的研究层次是：

• 内存

• 线程

• 核函数

  • 启动核函数
  • 编写核函数
  • 验证核函数

• 错误处理

内存管理

内存管理在传统串行程序是非常常见的，寄存器空间，栈空间内的内存由机器自己管理，堆空间由用户控制分配和释放，CUDA 程序同样，只是 CUDA 提供的 API 可以分配管理设备上的内存，当然也可以用 CDUA 管理主机上的内存，主机上的传统标准库也能完成主机内存管理。

下面表格有一些主机 API 和 CUDA C 的 API 的对比：

我们先研究最关键的一步，这一步要走总线的

1
2

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

这个函数是内存拷贝过程，可以完成以下几种过程（cudaMemcpyKind kind）

• cudaMemcpyHostToHost
• cudaMemcpyHostToDevice
• cudaMemcpyDeviceToHost
• cudaMemcpyDeviceToDevice

如果函数执行成功，则会返回 cudaSuccess 否则返回 cudaErrorMemoryAllocation

使用下面这个指令可以吧上面的错误代码翻译成详细信息：

1

char* cudaGetErrorString(cudaError_t error)

内存是分层次的，下图可以简单地描述，但是不够准确，后面我们会详细介绍每一个具体的环节：

解释下内存管理部分的代码：

1

cudaMalloc((float**)&a_d,nByte);

分配设备端的内存空间，为了区分设备和主机端内存，我们可以给变量加后缀或者前缀 h表示 host，d表示 device

线程管理

当内核函数开始执行，如何组织 GPU 的线程就变成了最主要的问题了，我们必须明确，一个核函数只能有一个 grid，一个 grid 可以有很多个块，每个块可以有很多的线程，这种分层的组织结构使得我们的并行过程更加自如灵活：

一个线程块 block 中的线程可以完成下述协作：

• 同步
• 共享内存

不同块内线程不能相互影响！他们是物理隔离的！

接下来就是给每个线程一个编号了，我们知道每个线程都执行同样的一段串行代码，那么怎么让这段相同的代码对应不同的数据呢？首先第一步就是让这些线程彼此区分开，才能对应到相应从线程，使得这些线程也能区分自己的数据。如果线程本身没有任何标记，那么没办法确认其行为。
依靠下面两个内置结构体确定线程标号：

• blockIdx（线程块在线程网格内的位置索引）
• threadIdx（线程在线程块内的位置索引）

注意这里的 Idx 是 index 的缩写，这两个内置结构体基于 uint3 定义，包含三个无符号整数的结构，通过三个字段来指定：

• blockIdx.x
• blockIdx.y
• blockIdx.z
• threadIdx.x
• threadIdx.y
• threadIdx.z

上面这两个是坐标，当然我们要有同样对应的两个结构体来保存其范围，也就是 blockIdx 中三个字段的范围 threadIdx 中三个字段的范围：

• blockDim
• gridDim

他们是 dim3 类型(基于 uint3 定义的数据结构)的变量，也包含三个字段 x,y,z.

• blockDim.x
• blockDim.y
• blockDim.z

网格和块的维度一般是二维和三维的，也就是说一个网格通常被分成二维的块，而每个块常被分成三维的线程。

注意：dim3 是手工定义的，主机端可见。uint3 是设备端在执行的时候可见的，不可以在核函数运行时修改，初始化完成后 uint3 值就不变了。他们是有区别的！这一点必须要注意。

网格和块的维度存在几个限制因素，块大小主要与可利用的计算资源有关，如寄存器共享内存。
分成网格和块的方式可以使得我们的 CUDA 程序可以在任意的设备上执行。

核函数概述

核函数就是在 CUDA 模型上诸多线程中运行的那段串行代码，这段代码在设备上运行，用 NVCC 编译，产生的机器码是 GPU 的机器码。所以我们写 CUDA 程序就是写核函数，第一步我们要确保核函数能正确的运行产生正切的结果，第二优化 CUDA 程序的部分，无论是优化算法，还是调整内存结构，线程结构都是要调整核函数内的代码，来完成这些优化的。

启动核函数

启动核函数，通过的以下的 ANSI C 扩展出的 CUDA C 指令：

1

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

其标准 C 的原型就是 C 语言函数调用

1

function_name(argument list);

这三个尖括号<<<grid,block>>>内是对设备代码执行的线程结构的配置（或者简称为对内核进行配置），也就是我们上一篇中提到的线程结构中的网格，块。回忆一下上文，我们通过 CUDA C 内置的数据类型 dim3 类型的变量来配置 grid 和 block（上文提到过：在设备端访问 grid 和 block 属性的数据类型是 uint3 不能修改的常类型结构，这里反复强调一下）。
通过指定 grid 和 block 的维度，我们可以配置：

• 内核中线程的数目
• 内核中使用的线程布局

我们可以使用 dim3 类型的 grid 维度和 block 维度配置内核，也可以使用 int 类型的变量，或者常量直接初始化：

1

kernel_name<<<4,8>>>(argument list);

上面这条指令的线程布局是：

我们的核函数是同时复制到多个线程执行的，上文我们说过一个对应问题，多个计算执行在一个数据，肯定是浪费时间，所以为了让多线程按照我们的意愿对应到不同的数据，就要给线程一个唯一的标识.

由于设备内存是线性的（基本市面上的内存硬件都是线性形式存储数据的）我们观察上图，可以用 threadIdx.x 和 blockIdx.x 来组合获得对应的线程的唯一标识.

接下来我们就是修改代码的时间了，改变核函数的配置，产生运行出结果一样，但效率不同的代码：

1. 一个块：

1

kernel_name<<<1,32>>>(argument list);

1. 32 个块

1

kernel_name<<<32,1>>>(argument list);

上述代码如果没有特殊结构在核函数中，执行结果应该一致，但是有些效率会一直比较低。

上面这些是启动部分，当主机启动了核函数，控制权马上回到主机，而不是主机等待设备完成核函数的运行,想要主机等待设备端执行可以用下面这个指令：

1

cudaError_t cudaDeviceSynchronize(void);

这是一个显式方法，对应的也有隐式方法，隐式方法就是不明确说明主机要等待设备端，而是设备端不执行完，主机没办法进行，比如内存拷贝函数：

1
2

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

所有 CUDA 核函数的启动都是异步的，这点与 C 语言是完全不同的

编写核函数

声明核函数有一个比较模板化的方法：

1

__global__ void kernel_name(argument list);

注意：声明和定义是不同的，这点 CUDA 与 C 语言是一致的

在 C 语言函数前没有的限定符global ,CUDA C 中还有一些其他我们在 C 中没有的限定符，如下：

这里有个特殊的情况就是有些函数可以同时定义为 __device__ 和 __host__ ，这种函数可以同时被设备和主机端的代码调用，主机端代码调用函数很正常，设备端调用函数与 C 语言一致，但是要声明成设备端代码，告诉 nvcc 编译成设备机器码，同时声明主机端设备端函数，那么就要告诉编译器，生成两份不同设备的机器码

Kernel 核函数编写有以下限制

1. 只能访问设备内存
2. 必须有 void 返回类型
3. 不支持可变数量的参数
4. 不支持静态变量
5. 显示异步行为

并行程序中经常的一种现象：把串行代码并行化时对串行代码块 for 的操作，也就是把 for 并行化。

1
2
3
4
5
6
7
8
9
10

// 串行：
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];
}

验证核函数

在开发阶段，每一步都进行验证是绝对高效的，比把所有功能都写好，然后进行测试这种过程效率高很多，同样写 CUDA 也是这样的每个代码小块都进行测试，看起来慢，实际会提高很多效率。

CUDA 小技巧，当我们进行调试的时候可以把核函数配置成单线程的：

kernel_name<<<1,1>>>(argument list)

错误处理

所有编程都需要对错误进行处理，早期的编码错误，编译器会帮我们搞定，内存错误也能观察出来，但是有些逻辑错误很难发现，而且 CUDA 基本都是异步执行的，当错误出现的时候，不一定是哪一条指令触发的，这一点非常头疼；这时候我们就需要对错误进行防御性处理了，例如如下这个宏：

1
2
3
4
5
6
7
8
9

#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);                                                         \
        }                                                                    \
    }

就是获得每个函数执行后的返回结果，然后对不成功的信息加以处理，CUDA C 的 API 每个调用都会返回一个错误代码，这个代码我们就可以好好利用。

给核函数计时

用 CPU 计时

使用 cpu 计时的方法是测试时间的一个常用办法，我们在写 C 程序的时候最多使用的计时方法是：

1
2
3
4
5

clock_t start, finish;
start = clock();
// 要测试的部分
finish = clock();
duration = (double)(finish - start) / CLOCKS_PER_SEC;

其中 clock()是个关键的函数，“clock 函数测出来的时间为进程运行时间，单位为滴答数(ticks)”；字面上理解 CLOCKS_PER_SEC 这个宏，就是每秒钟多少 clocks，在不同的系统中值可能不同。

必须注意的是，并行程序这种计时方式有严重问题！如果想知道具体原因，可以查询 clock 的源代码（c 语言标准函数）

1
2
3
4
5
6

#include <sys/time.h>
double cpuSecond(){
  struct timeval tp;
  gettimeofday(&tp,NULL);
  return((double)tp.tv_sec+(double)tp.tv_usec*1e-6);
}

gettimeofday 是 linux 下的一个库函数，创建一个 cpu 计时器，从 1970 年 1 月 1 日 0 点以来到现在的秒数，需要头文件 sys/time.h

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

#include "freshman.h"
#include <cuda_runtime.h>
#include <stdio.h>
__global__ void sumArraysGPU(float *a, float *b, float *res, int N) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N)
        res[i] = a[i] + b[i];
}
int main(int argc, char **argv) {
    // timer
    double iStart, iElaps;
    iStart = cpuSecond();
    sumArraysGPU<<<grid, block>>>(a_d, b_d, res_d, nElem);
    cudaDeviceSynchronize();
    iElaps = cpuSecond() - iStart;
}

主要分析计时这段，首先 iStart 是 cpuSecond 返回一个秒数，接着执行核函数，核函数开始执行后马上返回主机线程，所以我们必须要加一个同步函数等待核函数执行完毕，如果不加这个同步函数，那么测试的时间是从调用核函数，到核函数返回给主机线程的时间段，而不是核函数的执行时间，加上了

1

cudaDeviceSynchronize();

函数后，计时是从调用核函数开始，到核函数执行完并返回给主机的时间段，下面图大致描述了执行过程的不同时间节点：

我们可以大概分析下核函数启动到结束的过程：

1. 主机线程启动核函数
2. 核函数启动成功
3. 控制返回主机线程
4. 核函数执行完成
5. 主机同步函数侦测到核函数执行完

我们要测试的是 2~4 的时间，但是用 CPU 计时方法，只能测试 1~5 的时间，所以测试得到的时间偏长。
接着我们调整下我们的参数，来看看不同线程维度对速度的影响，看看计时能不能反映出来点问题，这里我们考虑一维线程模型

• 2 的幂次数据量 1<<24，16 兆数据：
  • 每个块 256 个线程
    
  • 每个块 512 个线程
    
  • 每个块 1024 个线程
• 2 的非幂次数据量 (1<<24)+1，16 兆加一个数据：
  • 每个块 256 个线程
  • 每个块 512 个线程
  • 每个块 1024 个线程

这个可以使用一点小技巧，比如只传输可完整切割数据块，然后剩下的 1，2 个使用 cpu 计算

用 nvprof 计时

CUDA 5.0 后有一个工具叫做 nvprof 的命令行分析工具。
nvprof 的用法如下：

1

nvprof [nvprof_args] <application>[application_args]

于是我们执行命令得到
工具不仅给出了 kernel 执行的时间，比例，还有其他 cuda 函数的执行时间。
nvprof 这个强大的工具给了我们优化的目标，分析数据可以得出我们重点工作要集中在哪部分。

理论界限最大化

得到了实际操作值，我们需要知道的是我们能优化的极限值是多少，也就是机器的理论计算极限，这个极限我们永远也达不到，但是我们必须明确的知道，比如理论极限是 2 秒，我们已经从 10 秒优化到 2.01 秒了，基本就没有必要再继续花大量时间优化速度了，而应该考虑买更多的机器或者更新的设备。
各个设备的理论极限可以通过其芯片说明计算得到，比如说：

• Tesla K10 单精度峰值浮点数计算次数：745MHz 核心频率 x 2GPU/芯片 x（8 个多处理器 x 192 个浮点计算单元 x 32 核心/多处理器） x 2 OPS/周期 =4.58 TFLOPS
• Tesla K10 内存带宽峰值： 2GPU/芯片 x 256 位 x 2500 MHz 内存时钟 x 2 DDR/8 位/字节 = 320 GB/s
• 指令比：字节 4.58 TFLOPS/320 GB/s =13.6 个指令： 1 个字节

组织并行线程

这一篇我们就详细介绍每一个线程是怎么确定唯一的索引，然后建立并行计算，并且不同的线程组织形式是怎样影响性能的：

• 二维网格二维线程块
• 一维网格一维线程块
• 二维网格一维线程块

使用块和线程建立矩阵索引

多线程的优点就是每个线程处理不同的数据计算，那么怎么分配好每个线程处理不同的数据，而不至于多个不同的线程处理同一个数据，或者避免不同的线程没有组织的乱访问内存。

如果多线程不能按照组织合理的干活，那么就相当于一群没训练过的哈士奇拉雪橇，往不同的方向跑，那么是没办法前进的，必须有组织，有规则的计算才有意义。

我们的线程模型前面已经有个大概的介绍，但是下图可以非常形象的反映线程模型，不过注意硬件实际的执行和存储不是按照图中的模型来的，大家注意区分：

这里(ix,iy)就是整个线程模型中任意一个线程的索引，或者叫做全局地址，局部地址当然就是(threadIdx.x,threadIdx.y)了，当然这个局部地址目前还没有什么用处，他只能索引线程块内的线程，不同线程块中有相同的局部索引值。
图中的横坐标就是：

ix=threadIdx.x+blockIdx.x×blockDim.x

纵坐标是：

iy=threadIdx.y+blockIdx.y×blockDim.y

这样我们就得到了每个线程的唯一标号，并且在运行时 kernel 是可以访问这个标号的。前面讲过 CUDA 每一个线程执行相同的代码，也就是异构计算中说的多线程单指令，如果每个不同的线程执行同样的代码，又处理同一组数据，将会得到多个相同的结果，显然这是没意义的，为了让不同线程处理不同的数据，CUDA 常用的做法是让不同的线程对应不同的数据，也就是用线程的全局标号对应不同组的数据。
设备内存或者主机内存都是线性存在的，比如一个二维矩阵 (8×6)(8×6)，存储在内存中是这样的：

我们要做管理的就是：

• 线程和块索引（来计算线程的全局索引）
• 矩阵中给定点的坐标（ix,iy）
• (ix,iy)对应的线性内存的位置

线性位置的计算方法是：

idx=ix+iy∗nx

我们上面已经计算出了线程的全局坐标，用线程的全局坐标对应矩阵的坐标，也就是说，线程的坐标(ix,iy)对应矩阵中(ix,iy)的元素，这样就形成了一一对应，不同的线程处理矩阵中不同的数据。
我们接下来的代码来输出每个线程的标号信息：

1
2
3
4
5
6
7
8
9

__global__ void printThreadIndex(float *A, const int nx, const int ny) {
    int ix = threadIdx.x + blockIdx.x * blockDim.x;
    int iy = threadIdx.y + blockIdx.y * blockDim.y;
    unsigned int idx = iy * nx + ix;
    printf("thread_id(%d,%d) block_id(%d,%d) coordinate(%d,%d)"
           "global index %2d ival %2d\n",
           threadIdx.x, threadIdx.y,
           blockIdx.x, blockIdx.y, ix, iy, idx, A[idx]);
}

二维矩阵加法

我们利用上面的线程与数据的对应完成了下面的核函数：

1
2
3
4
5
6
7
8

__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];
    }
}

下面我们调整不同的线程组织形式，测试一下不同的效率并保证得到正确的结果，但是什么时候得到最好的效率是后面要考虑的，我们要做的就是用各种不同的相乘组织形式得到正确结果

不同网格组合运行结果

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108

#include "freshman.h"
#include <cuda_runtime.h>
#include <stdio.h>
void sumMatrix2D_CPU(float *MatA, float *MatB, float *MatC, int nx, int ny) {
    float *a = MatA;
    float *b = MatB;
    float *c = MatC;
    for (int j = 0; j < ny; j++) {
        for (int i = 0; i < nx; i++) {
            c[i] = a[i] + b[i];
        }
        c += nx;
        b += nx;
        a += nx;
    }
}
__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("starting...\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;
}

运行结果：

总结

用不同的线程组织形式会得到正确结果，但是效率有所区别，观察结果没有多大差距，但是明显比 CPU 快了很多，而且最主要的是我们本文用不同的线程组织模式都得到了正确结果，并且

• 改变执行配置（线程组织）能得到不同的性能
• 传统的核函数可能不能得到最好的效果
• 一个给定的核函数，通过调整网格和线程块大小可以得到更好的效果

GPU 设备信息

API 查询 GPU 信息

在软件内查询信息，用到如下代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

#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("  Maximum number of thread per multiprocesser:  %d\n",
           deviceProp.maxThreadsPerMultiProcessor);
    printf("  Maximum number of thread per block:           %d\n",
           deviceProp.maxThreadsPerBlock);
    printf("  Maximum size of each dimension of a block:    %d x %d x %d\n",
           deviceProp.maxThreadsDim[0], deviceProp.maxThreadsDim[1], deviceProp.maxThreadsDim[2]);
    printf("  Maximum size of each dimension of a grid:     %d x %d x %d\n",
           deviceProp.maxGridSize[0],
           deviceProp.maxGridSize[1],
           deviceProp.maxGridSize[2]);
    printf("  Maximum 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;
}

运行的效果如下：

这里面很多参数是我们后面要介绍的，而且每一个都对性能有影响：

1. CUDA 驱动版本
2. 设备计算能力编号
3. 全局内存大小（1.95G）
4. GPU 主频
5. GPU 带宽
6. L2 缓存大小
7. 不同维度下的纹理维度最大值
8. 层叠纹理维度最大值
9. 常量内存大小
10. 块内共享内存大小
11. 块内寄存器大小
12. 线程束大小
13. 每个处理器硬件处理的最大线程数
14. 每个块处理的最大线程数
15. 块的最大尺寸
16. 网格的最大尺寸
17. 最大连续线性内存

上面这些都是后面要用到的关键参数，这些会严重影响我们的效率。后面会一一说到，不同的设备参数要按照不同的参数来使得程序效率最大化，所以我们必须在程序运行前得到所有我们关心的参数。

NVIDIA-SMI

nvidia-smi 是 nvidia 驱动程序内带的一个工具，可以返回当前环境的设备信息：

nvitop更好用

这个命令可以加各种参数，当然参数你要查文档查文档查文档：
利用下面这个参数可以精简上面那么一大堆的信息，而这些可以在脚本中帮我们得到设备信息，比如我们可以写通用程序时在编译前执行脚本来获取设备信息，然后在编译时固化最优参数，这样程序运行时就不会被查询设备信息的过程浪费资源。
也就是我们可以用以下两种方式编写通用程序：

• 运行时获取设备信息：
  • 编译程序
  • 启动程序
  • 查询信息，将信息保存到全局变量
  • 功能函数通过全局变量判断当前设备信息，优化参数
  • 程序运行完毕
• 编译时获取设备信息
  • 脚本获取设备信息
  • 编译程序，根据设备信息调整固化参数到二进制机器码
  • 运行程序
  • 程序运行完毕

详细信息使用

1

nvidia-smi -q -i

CUDA 执行模型概述

这一篇开始我们开始接近 CUDA 最核心的部分，就是有关硬件和程序的执行模型，用 CUDA 的目的其实说白了就是为计算速度快，所以压榨性能，提高效率其实就是 CUDA 学习的最终目的。

概述

CUDA 执行模型揭示了 GPU 并行架构的抽象视图，再设计硬件的时候，其功能和特性都已经被设计好了，然后去开发硬件，如果这个过程模型特性或功能与硬件设计有冲突，双方就会进行商讨妥协，知道最后产品定型量产，功能和特性算是全部定型，而这些功能和特性就是变成模型的设计基础，而编程模型又直接反应了硬件设计，从而反映了设备的硬件特性。
比如最直观的一个就是内存，线程的层次结构帮助我们控制大规模并行，这个特性就是硬件设计最初设计好，然后集成电路工程师拿去设计，定型后程序员开发驱动，然后在上层可以直接使用这种执行模型来控制硬件。
所以了解 CUDA 的执行模型，可以帮助我们优化指令吞吐量，和内存使用来获得极限速度。

GPU 架构概述

GPU 架构是围绕一个流式多处理器（SM）的扩展阵列搭建的。通过复制这种结构来实现 GPU 的硬件并行。

上图包括关键组件：

• CUDA 核心

• 共享内存/一级缓存

• 寄存器文件

• 加载/存储单元

• 特殊功能单元

• 线程束调度器

SM

GPU 中每个 SM 都能支持数百个线程并发执行，每个 GPU 通常有多个 SM，当一个核函数的网格被启动的时候，多个 block 会被同时分配给可用的 SM 上执行。

注意: 当一个 blcok 被分配给一个 SM 后，他就只能在这个 SM 上执行了，不可能重新分配到其他 SM 上了，多个线程块可以被分配到同一个 SM 上。

在 SM 上同一个块内的多个线程进行线程级别并行，而同一线程内，指令利用指令级并行将单个线程处理成流水线。

线程束

CUDA 采用单指令多线程 SIMT 架构管理执行线程，不同设备有不同的线程束大小，但是到目前为止基本所有设备都是维持在 32，也就是说每个 SM 上有多个 block，一个 block 有多个线程（可以是几百个，但不会超过某个最大值），但是从机器的角度，在某时刻 T，SM 上只执行一个线程束，也就是 32 个线程在同时同步执行，线程束中的每个线程执行同一条指令，包括有分支的部分，这个我们后面会讲到，

SIMD vs SIMT

单指令多数据的执行属于向量机，比如我们有四个数字要加上四个数字，那么我们可以用这种单指令多数据的指令来一次完成本来要做四次的运算。这种机制的问题就是过于死板，不允许每个分支有不同的操作，所有分支必须同时执行相同的指令，必须执行没有例外。
相比之下单指令多线程 SIMT 就更加灵活了，虽然两者都是将相同指令广播给多个执行单元，但是 SIMT 的某些线程可以选择不执行，也就是说同一时刻所有线程被分配给相同的指令，SIMD 规定所有人必须执行，而 SIMT 则规定有些人可以根据需要不执行，这样 SIMT 就保证了线程级别的并行，而 SIMD 更像是指令级别的并行。
SIMT 包括以下 SIMD 不具有的关键特性：

1. 每个线程都有自己的指令地址计数器
2. 每个县城都有自己的寄存器状态
3. 每个线程可以有一个独立的执行路径

而上面这三个特性在编程模型可用的方式就是给每个线程一个唯一的标号（blckIdx,threadIdx），并且这三个特性保证了各线程之间的独立

从概念上讲，32 是 SM 以 SIMD 方式同时处理的工作粒度，这句话这么理解，可能学过后面的会更深刻的明白，一个 SM 上在某一个时刻，有 32 个线程在执行同一条指令，这 32 个线程可以选择性执行，虽然有些可以不执行，但是他也不能执行别的指令.

CUDA 编程的组件与逻辑

下图从逻辑角度和硬件角度描述了 CUDA 编程模型对应的组件。

SM 中共享内存，和寄存器是关键的资源，线程块中线程通过共享内存和寄存器相互通信协调。

寄存器和共享内存的分配可以严重影响性能！

因为 SM 有限，虽然我们的编程模型层面看所有线程都是并行执行的，但是在微观上看，所有线程块也是分批次的在物理层面的机器上执行，线程块里不同的线程可能进度都不一样，但是同一个线程束内的线程拥有相同的进度。
并行就会引起竞争，多线程以未定义的顺序访问同一个数据，就导致了不可预测的行为，CUDA 只提供了一种块内同步的方式，块之间没办法同步！
同一个 SM 上可以有不止一个常驻的线程束，有些在执行，有些在等待，他们之间状态的转换是不需要开销的。

Fermi 架构

Fermi 架构是第一个完整的 GPU 架构，所以了解这个架构是非常有必要的。

Fermi 架构逻辑图如上，具体数据如下:

1. 512 个加速核心，CUDA 核
2. 每个 CUDA 核心都有一个全流水线的整数算术逻辑单元 ALU，和一个浮点数运算单元 FPU
3. CUDA 核被组织到 16 个 SM 上
4. 6 个 384-bits 的 GDDR5 的内存接口
5. 支持 6G 的全局机载内存
6. GigaThread 疫情，分配线程块到 SM 线程束调度器上
7. 768KB 的二级缓存，被所有 SM 共享

而 SM 则包括下面这些资源：

• 执行单元（CUDA 核）
• 调度线程束的调度器和调度单元
• 共享内存，寄存器文件和一级缓存

每个多处理器 SM 有 16 个加载/存储单元所以每个时钟周期内有 16 个线程（半个线程束）计算源地址和目的地址

特殊功能单元 SFU 执行固有指令，如正弦，余弦，平方根和插值，SFU 在每个时钟周期内的每个线程上执行一个固有指令。
每个 SM 有两个线程束调度器，和两个指令调度单元，当一个线程块被指定给一个 SM 时，线程块内的所有线程被分成线程束，两个线程束选择其中两个线程束，在用指令调度器存储两个线程束要执行的指令
像第一张图上的显示一样，每 16 个 CUDA 核心为一个组，还有 16 个加载/存储单元或 4 个特殊功能单元。当某个线程块被分配到一个 SM 上的时候，会被分成多个线程束，线程束在 SM 上交替执行：

上面曾经说过，每个线程束在同一时间执行同一指令，同一个块内的线程束互相切换是没有时间消耗的。
Fermi 上支持同时并发执行内核。并发执行内核允许执行一些小的内核程序来充分利用 GPU，如图：

Kepler 架构

Kepler 架构作为 Fermi 架构的后代，有以下技术突破：

• 强化的 SM
• 动态并行
• Hyper-Q 技术

技术参数也提高了不少，比如单个 SM 上 CUDA 核的数量，SFU 的数量，LD/ST 的数量等：

kepler 架构的最突出的一个特点就是内核可以启动内核了，这使得我们可以使用 GPU 完成简单的递归操作，流程如下。

Hyper-Q 技术主要是 CPU 和 GPU 之间的同步硬件连接，以确保 CPU 在 GPU 执行的同时做更多的工作。Fermi 架构下 CPU 控制 GPU 只有一个队列，Kepler 架构下可以通过 Hyper-Q 技术实现多个队列如下图。

计算能力概览：

使用 Profile 进行优化（Profile-Driven Optimization）

中文翻译的这个标题是配置文件驱动优化，驱动这个词在这里应该是个动词，或者翻译的人直接按照字面意思翻译的，其实看完内容以后的意思是根据 profile 这个文件内的信息对程序进行优化。
性能分析通过以下方法来进行：

1. 应用程序代码的空间(内存)或时间复杂度
2. 特殊指令的使用
3. 函数调用的频率和持续时间

程序优化是建立在对硬件和算法过程理解的基础上的，如果都不了解，靠试验，那么这个结果可想而知。理解平台的执行模型也就是硬件特点，是优化性能的基础。
开发高性能计算程序的两步：

1. 保证结果正确，和程序健壮性
2. 优化速度

Profile 可以帮助我们观察程序内部。

• 一个原生的内核应用一般不会产生最佳效果，也就是我们基本不能一下子就写出最好最快的内核，需要通过性能分析工具分析性能。找出性能瓶颈
• CUDA 将 SM 中的计算资源在该 SM 中的多个常驻线程块之间进行分配，这种分配方式可能导致一些资源成为性能限制因素，性能分析工具可以帮我们找出来这些资源是如何被使用的
• CUDA 提供了一个硬件架构的抽象。它能够让用户控制线程并发。性能分析工具可以检测和优化，并肩优化可视化

总结起来一句话，想优化速度，先学好怎么用性能分析工具。

• nvvp
• nvprof

限制内核性能的主要包括但不限于以下因素

• 存储带宽
• 计算资源
• 指令和内存延迟

理解线程束执行的本质

我们前面已经大概的介绍了 CUDA 执行模型的大概过程，包括线程网格，线程束，线程间的关系，以及硬件的大概结构，例如 SM 的大概结构，而对于硬件来说，CUDA 执行的实质是线程束的执行，因为硬件根本不知道每个块谁是谁，也不知道先后顺序，硬件(SM)只知道按照机器码跑，而给他什么，先后顺序，这个就是硬件功能设计的直接体现了。
从外表来看，CUDA 执行所有的线程，并行的，没有先后次序的，但实际上硬件资源是有限的，不可能同时执行百万个线程，所以从硬件角度来看，物理层面上执行的也只是线程的一部分，而每次执行的这一部分，就是我们前面提到的线程束。

线程束和线程块

线程束是 SM 中基本的执行单元，当一个网格被启动（网格被启动，等价于一个内核被启动，每个内核对应于自己的网格），网格中包含线程块，线程块被分配到某一个 SM 上以后，将分为多个线程束，每个线程束一般是 32 个线程，在一个线程束中，所有线程按照单指令多线程 SIMT 的方式执行，每一步执行相同的指令，但是处理的数据为私有的数据，下图反应的就是逻辑，实际，和硬件的图形化

线程块是个逻辑产物，因为在计算机里，内存总是一维线性存在的，所以执行起来也是一维的访问线程块中的线程，但是我们在写程序的时候却可以以二维三维的方式进行，原因是方便我们写程序，比如处理图像或者三维的数据，三维块就会变得很直接，很方便。
在块中，每个线程有唯一的编号（可能是个三维的编号），threadIdx。
网格中，每个线程块也有唯一的编号(可能是个三维的编号)，blockIdx
那么每个线程就有在网格中的唯一编号。
当一个线程块中有 128 个线程的时候，其分配到 SM 上执行时，会分成 4 个块：

1
2
3
4

warp0: thread  0,........thread31
warp1: thread 32,........thread63
warp2: thread 64,........thread95
warp3: thread 96,........thread127

当编号使用三维编号时，x 位于最内层，y 位于中层，z 位于最外层，想象下 c 语言的数组，如果把上面这句话写成 c 语言，假设三维数组 t 保存了所有的线程，那么(threadIdx.x,threadIdx.y,threadIdx.z)表示为

1

t[z][y][x];

计算出三维对应的线性地址是：

tid=threadIdx.x+threadIdx.y×blockDim.x+threadIdx.z×blockDim.x×blockDim.y

上面的公式可以借助 c 语言的三维数组计算相对地址的方法理解

一个线程块包含多少个线程束呢？

WarpsPerBlock=ceil(ThreadsPerBlockwarpSize)

线程束和线程块，一个是硬件层面的线程集合，一个是逻辑层面的线程集合，我们编程时为了程序正确，必须从逻辑层面计算清楚，但是为了得到更快的程序，硬件层面是我们应该注意的。

线程束分化

线程束被执行的时候会被分配给相同的指令，处理各自私有的数据，在 CUDA 中支持 C 语言的控制流，比如 if…else, for ,while 等，CUDA 中同样支持，但是如果一个线程束中的不同线程包含不同的控制条件，那么当我们执行到这个控制条件是就会面临不同的选择。
这里要讲一下 CPU 了，当我们的程序包含大量的分支判断时，从程序角度来说，程序的逻辑是很复杂的，因为一个分支就会有两条路可以走，如果有 10 个分支，那么一共有 1024 条路走，CPU 采用流水线话作业，如果每次等到分支执行完再执行下面的指令会造成很大的延迟，所以现在处理器都采用分支预测技术，而 CPU 的这项技术相对于 gpu 来说高级了不止一点点，而这也是 GPU 与 CPU 的不同，设计初衷就是为了解决不同的问题。CPU 适合逻辑复杂计算量不大的程序，比如操作系统，控制系统，GPU 适合大量计算简单逻辑的任务，所以被用来算数。
如下一段代码：

1
2
3
4
5
6

if (con){
    //do something
}
else{
    //do something
}

假设这段代码是核函数的一部分，那么当一个线程束的 32 个线程执行这段代码的时候，如果其中 16 个执行 if 中的代码段，而另外 16 个执行 else 中的代码块，同一个线程束中的线程，执行不同的指令，这叫做线程束的分化。
我们知道在每个指令周期，线程束中的所有线程执行相同的指令，但是线程束又是分化的，所以这似乎是相悖的，但是事实上这两个可以不矛盾。
解决矛盾的办法就是每个线程都执行所有的 if 和 else 部分，当一部分 con 成立的时候，执行 if 块内的代码，这些 con 不成立的线程等待，线程束分化会产生严重的性能下降。条件分支越多，并行性削弱越严重。
注意线程束分化研究的是一个线程束中的线程，不同线程束中的分支互不影响。
执行过程如下：

因为线程束分化导致的性能下降就应该用线程束的方法解决，根本思路是避免同一个线程束内的线程分化。

我们能控制线程束内线程行为的原因是线程块中线程分配到线程束是有规律的而不是随机的。这就使得我们根据线程编号来设计分支是可以的，补充说明下，当一个线程束中所有的线程都执行 if 或者，都执行 else 时，不存在性能下降；只有当线程束内有分歧产生分支的时候，性能才会急剧下降。

线程束内的线程是可以被我们控制的，那么我们就把都执行 if 的线程塞到一个线程束中，或者让一个线程束中的线程都执行 if，另外线程都执行 else 的这种方式可以将效率提高很多。

下面这个 kernel 可以产生一个比较低效的分支：

1
2
3
4
5
6
7
8
9
10
11

__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;
}

这种情况下我们假设只配置一个 x=64 的一维线程块，那么只有两个个线程束，线程束内奇数线程（threadIdx.x 为奇数）会执行 else，偶数线程执行 if，分化很严重。
但是如果我们换一种方法，得到相同但是错乱的结果 C，这个顺序其实是无所谓的，因为我们可以后期调整。那么下面代码就会很高效

1
2
3
4
5
6
7
8
9
10
11

__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;
}

第一个线程束内的线程编号 tid 从 0 到 31，tid/warpSize 都等于 0，那么就都执行 if 语句。
第二个线程束内的线程编号 tid 从 32 到 63，tid/warpSize 都等于 1，执行 else
线程束内没有分支，效率较高。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119

#include "freshman.h"
#include <cuda_runtime.h>
#include <stdio.h>
#include <stdlib.h>

__global__ void warmup(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;
    }
    // printf("%d %d %f \n",tid,warpSize,a+b);
    c[tid] = a + b;
}
__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;
}
__global__ void mathKernel3(float *c) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float a = 0.0;
    float b = 0.0;
    bool ipred = (tid % 2 == 0);
    if (ipred) {
        a = 100.0f;
    } else {
        b = 200.0f;
    }
    c[tid] = a + b;
}

int main(int argc, char **argv) {
    int dev = 0;
    cudaDeviceProp deviceProp;
    cudaGetDeviceProperties(&deviceProp, dev);
    printf("%s using Device %d: %s\n", argv[0], dev, deviceProp.name);

    // set up data size
    int size = 64;
    int blocksize = 64;
    if (argc > 1)
        blocksize = atoi(argv[1]);
    if (argc > 2)
        size = atoi(argv[2]);
    printf("Data size %d ", size);

    // set up execution configuration
    dim3 block(blocksize, 1);
    dim3 grid((size + block.x - 1) / block.x, 1);
    printf("Execution Configure (block %d grid %d)\n", block.x, grid.x);

    // allocate gpu memory
    float *C_dev;
    size_t nBytes = size * sizeof(float);
    float *C_host = (float *)malloc(nBytes);
    cudaMalloc((float **)&C_dev, nBytes);

    // run a warmup kernel to remove overhead
    double iStart, iElaps;
    cudaDeviceSynchronize();
    iStart = cpuSecond();
    warmup<<<grid, block>>>(C_dev);
    cudaDeviceSynchronize();
    iElaps = cpuSecond() - iStart;

    printf("warmup	   <<<%4d,%4d>>>elapsed %lf sec \n", grid.x, block.x, iElaps);

    // run kernel 1
    iStart = cpuSecond();
    mathKernel1<<<grid, block>>>(C_dev);
    cudaDeviceSynchronize();
    iElaps = cpuSecond() - iStart;
    printf("mathKernel1<<<%4d,%4d>>>elapsed %lf sec \n", grid.x, block.x, iElaps);
    cudaMemcpy(C_host, C_dev, nBytes, cudaMemcpyDeviceToHost);
    // run kernel 2
    iStart = cpuSecond();
    mathKernel2<<<grid, block>>>(C_dev);
    cudaDeviceSynchronize();
    iElaps = cpuSecond() - iStart;
    printf("mathKernel2<<<%4d,%4d>>>elapsed %lf sec \n", grid.x, block.x, iElaps);

    // run kernel 3
    iStart = cpuSecond();
    mathKernel3<<<grid, block>>>(C_dev);
    cudaDeviceSynchronize();
    iElaps = cpuSecond() - iStart;
    printf("mathKernel3<<<%4d,%4d>>>elapsed %lf sec \n", grid.x, block.x, iElaps);

    cudaFree(C_dev);
    free(C_host);
    cudaDeviceReset();
    return EXIT_SUCCESS;
}

事件和指标

上面我们提到了事件 event，事件是可计算的活动，比如这个分支就是一个可以计算的活动，对应一个在内核执行期间被搜集的硬件计数器。
指标是内核的特征，有一个或多个事件计算得到。

资源分配

我们前面提到过，每个 SM 上执行的基本单位是线程束，也就是说，单指令通过指令调度器广播给某线程束的全部线程，这些线程同一时刻执行同一命令，当然也有分支情况，线程束中一类是已经激活的，也就是说这类线程束其实已经在 SM 上准备就绪了，只是没轮到他执行，这时候他的状态叫做阻塞，还有一类可能分配到 SM 了，但是还没上到片上，这类我称之为未激活线程束。
而每个 SM 上有多少个线程束处于激活状态，取决于以下资源：

• 程序计数器
• 寄存器
• 共享内存

线程束一旦被激活来到片上，那么他就不会再离开 SM 直到执行结束。
每个 SM 都有 32 位的寄存器组，每个架构寄存器的数量不一样，其存储于寄存器文件中，为每个线程进行分配，同时，固定数量的共享内存，在线程块之间分配。
一个 SM 上被分配多少个线程块和线程束取决于 SM 中可用的寄存器和共享内存，以及内核需要的寄存器和共享内存大小。

这是一个平衡问题，当 kernel 占用的资源较少，那么更多的线程（这时线程越多线程束也就越多）处于活跃状态，相反则线程越少。
关于寄存器资源的分配：

关于共享内存的分配：

上面讲的主要是线程束，如果从逻辑上来看线程块的话，可用资源的分配也会影响常驻线程块的数量。
特别是当 SM 内的资源没办法处理一个完整块，那么程序将无法启动。

以下是资源列表：

当寄存器和共享内存分配给了线程块，这个线程块处于活跃状态，所包含的线程束称为活跃线程束。
活跃的线程束又分为三类：

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

当 SM 要执行某个线程束的时候，执行的这个线程束叫做选定的线程束，准备要执行的叫符合条件的线程束，如果线程束不符合条件还没准备好就是阻塞的线程束。
满足下面的要求，线程束才算是符合条件的：

• 32 个 CUDA 核心可以用于执行
• 执行所需要的资源全部就位

Kepler 活跃的线程束数量从开始到结束不得大于 64，可以等于。
任何周期选定的线程束小于等于 4。
由于计算资源是在线程束之间分配的，且线程束的整个生命周期都在片上，所以线程束的上下文切换是非常快速的。
下面我们介绍如何通过大量的活跃的线程束切换来隐藏延迟

延迟隐藏

延迟隐藏，延迟是什么，就是当你让计算机帮你算一个东西的时候计算需要用的时间
所以最大化是要最大化硬件，尤其是计算部分的硬件满跑，都不闲着的情况下利用率是最高的，总有人闲着，利用率就会低很多，即最大化功能单元的利用率。利用率与常驻线程束直接相关。
硬件中线程调度器负责调度线程束调度，当每时每刻都有可用的线程束供其调度，这时候可以达到计算资源的完全利用，以此来保证通过其他常驻线程束中发布其他指令的，可以隐藏每个指令的延迟。

与其他类型的编程相比，GPU 的延迟隐藏及其重要。对于指令的延迟，通常分为两种：

• 算术指令
• 内存指令

算数指令延迟是一个算术操作从开始，到产生结果之间的时间，这个时间段内只有某些计算单元处于工作状态，而其他逻辑计算单元处于空闲。
内存指令延迟很好理解，当产生内存访问的时候，计算单元要等数据从内存拿到寄存器，这个周期是非常长的。
延迟：

• 算术延迟 10~20 个时钟周期
• 内存延迟 400~800 个时钟周期

下图就是阻塞线程束到可选线程束的过程逻辑图：

其中线程束 0 在阻塞两段时间后恢复可选模式，但是在这段等待时间中，SM 没有闲置。
那么至少需要多少线程，线程束来保证最小化延迟呢？
little 法则给出了下面的计算公式

所需线程束=延迟 × 吞吐量

注意带宽和吞吐量的区别，带宽一般指的是理论峰值，最大每个时钟周期能执行多少个指令，吞吐量是指实际操作过程中每分钟处理多少个指令。

这个可以想象成一个瀑布，像这样，绿箭头是线程束，只要线程束足够多，吞吐量是不会降低的：

下面表格给出了 Fermi 和 Kepler 执行某个简单计算时需要的并行操作数：

另外有两种方法可以提高并行：

• 指令级并行(ILP): 一个线程中有很多独立的指令
• 线程级并行(TLP): 很多并发地符合条件的线程

同样，与指令周期隐藏延迟类似，内存隐藏延迟是靠内存读取的并发操作来完成的。

需要注意的是，指令隐藏的关键目的是使用全部的计算资源，而内存读取的延迟隐藏是为了使用全部的内存带宽，内存延迟的时候，计算资源正在被别的线程束使用，所以我们不考虑内存读取延迟的时候计算资源在做了什么，这两种延迟我们看做两个不同的部门但是遵循相同的道理。

我们的根本目的是把计算资源，内存读取的带宽资源全部使用满，这样就能达到理论的最大效率。
同样下表根据 Little 法则给出了需要多少线程束来最小化内存读取延迟，不过这里有个单位换算过程，机器的性能指标内存读取速度给出的是 GB/s 的单位，而我们需要的是每个时钟周期读取字节数，所以要用这个速度除以频率，例如 C 2070 的内存带宽是 144 GB/s 化成时钟周期： 144GB/s1.566GHz=92B/t,这样就能得到单位时间周期的内存带宽了。
得出下表的数据

需要说明的是这个速度不是单个 SM 的而是整个 GPU 设备的，内存带宽是 GPU 设备的而不是针对一个 SM 的。
Fermi 需要并行的读取 74 的数据才能让 GPU 带宽满载，如果每个线程读取 4 个字节，我们大约需要 18500 个线程，大约 579 个线程束才能达到这个峰值。
所以，延迟的隐藏取决于活动的线程束的数量，数量越多，隐藏的越好，但是线程束的数量又受到上面的说的资源影响。所以这里就需要寻找最优的执行配置来达到最优的延迟隐藏。

那么我们怎么样确定一个线程束的下界呢，使得当高于这个数字时 SM 的延迟能充分的隐藏，其实这个公式很简单，也很好理解，就是 SM 的计算核心数乘以单条指令的延迟，
比如 32 个单精度浮点计算器，每次计算延迟 20 个时钟周期，那么我需要最少 32x20 =640 个线程使设备处于忙碌状态。

占用率

占用率是一个 SM 种活跃的线程束的数量，占 SM 最大支持线程束数量的比，

最大 64 个线程束每个 SM。
CUDA 工具包中提供一个叫做 UCDA 占用率计算器的电子表格，填上相关数据可以帮你自动计算网格参数：

上面我们已经明确内核使用寄存器的数量会影响 SM 内线程束的数量，nvcc 的编译选项也有手动控制寄存器的使用。
也可以通过调整线程块内线程的多少来提高占用率，当然要合理不能太极端：

• 小的线程块：每个线程块中线程太少，会在所有资源没用完就达到了线程束的最大要求
• 大的线程块：每个线程块中太多线程，会导致每个 SM 中每个线程可用的硬件资源较少

同步

并发程序对同步非常有用，比如 pthread 中的锁，openmp 中的同步机制，这么做的主要目的是避免内存竞争
CUDA 同步这里只讲两种：

• 线程块内同步
• 系统级别

块级别的就是同一个块内的线程会同时停止在某个设定的位置，用__syncthread();这个函数完成

这个函数只能同步同一个块内的线程，不能同步不同块内的线程，想要同步不同块内的线程，就只能让核函数执行完成，控制程序交换主机，这种方式来同步所有线程。

内存竞争是非常危险的，一定要非常小心，这里经常出错。

可扩展性

可扩展性其实是相对于不同硬件的，当某个程序在设备 1 上执行的时候时间消耗是 T 当我们使用设备 2 时，其资源是设备 1 的两倍，我们希望得到 T/2 的运行速度，这种性质是 CUDA 驱动部分提供的特性，目前来说 Nvidia 正在致力于这方面的优化，如下图：

并行性表现

本文的主要内容就是进一步理解线程束在硬件上执行的本质过程，结合上几篇关于执行模型的学习，本文相对简单，通过修改核函数的配置，来观察核函数的执行速度，以及分析硬件利用数据，分析性能，调整核函数配置是 CUDA 开发人员必须掌握的技能，本篇只研究对核函数的配置是如何影响效率的（也就是通过网格，块的配置来获得不同的执行效率。）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92

#include "freshman.h"
#include <cuda_runtime.h>
#include <stdio.h>
void sumMatrix2D_CPU(float *MatA, float *MatB, float *MatC, int nx, int ny) {
    float *a = MatA;
    float *b = MatB;
    float *c = MatC;
    for (int j = 0; j < ny; j++) {
        for (int i = 0; i < nx; i++) {
            c[i] = a[i] + b[i];
        }
        c += nx;
        b += nx;
        a += nx;
    }
}
__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) {
    initDevice(0);
    int nx = 1 << 13;
    int ny = 1 << 13;
    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 = argc > 2 ? atoi(argv[1]) : 32;
    int dimy = argc > 2 ? atoi(argv[2]) : 32;

    double iStart, iElaps;
    // cpu compute
    iStart = cpuSecond();
    sumMatrix2D_CPU(A_host, B_host, C_host, nx, ny);
    iElaps = cpuSecond() - iStart;
    printf("CPU Execution Time elapsed %f sec\n", iElaps);
    // warm up
    //  2d block and 2d grid
    dim3 block_0(32, 32);
    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());
    printf("Warm Up \n");

    // 2d block and 2d grid
    dim3 block(dimx, dimy);
    dim3 grid((nx - 1) / block.x + 1, (ny - 1) / block.y + 1);
    iStart = cpuSecond();
    sumMatrix<<<grid, block>>>(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.x, grid.y, block.x, block.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;
}

可见我们用两个 8192×8192 的矩阵相加来测试我们效率。
注意一下这里的 GPU 内存，一个矩阵是 213×213×22=228 字节 也就是 256M，三个矩阵就是 768M 因为我们的 GPU 内存就是 2G 的，所以我们没办法进行更大的矩阵计算了。

用 nvprof 检测活跃的线程束

对比性能要控制变量，上面的代码只用两个变量，也就是块的 x 和 y 的大小，所以，调整 x 和 y 的大小来产生不同的效率，我们先来看看结果,数据结果如下

当块大小超过硬件的极限，并没有报错，而是返回了错误值，所以需要使用CHECK来辅助检查，另外，每个机器执行此代码效果可能不一样，所以大家要根据自己的硬件分析数据。

活跃线程束比例高的未必执行速度快，但实际上从原理出发，应该是利用率越高效率越高，但是还受到其他因素制约。

活跃线程束比例的定义是：每个周期活跃的线程束的平均值与一个 sm 支持的线程束最大值的比。

用 nvprof 检测内存操作

下面我们继续用 nvprof 来看看内存利用率如何。
首先使用:

1

nvprof --metrics gld_throughput ./simple_sum_matrix

来看一下内核的内存读取效率：

可以看出虽然第一种配置的线程束活跃比例不高，但是吞吐量最大所以可见吞吐量和线程束活跃比例一起都对最终的效率有影响。
接着我们看看全局加载效率，全局效率的定义是：被请求的全局加载吞吐量占所需的全局加载吞吐量的比值（全局加载吞吐量），也就是说应用程序的加载操作利用了设备内存带宽的程度；注意区别吞吐量和全局加载效率的区别，这个在前面我们已经解释过吞吐量了，忘了的同学回去看看。

1

nvprof --metrics gld_efficiency ./simple_sum_matrix

结果中在当前机器上进行测试所有的利用率都是 100% ，可见 CUDA 对核函数进行了优化，在 M2070 上 使用以前的 CUDA 版本，并没有如此高的加载效率，有效加载效率是指在全部的内存请求中（当前在总线上传递的数据）有多少是我们要用于计算的。

书上说如果线程块中内层的维度（blockDim.x）过小，小于线程束会影响加载效率，但是目前来看，不存在这个问题了。
随着硬件的升级，以前的一些问题，可能就不是问题了，当然对付老的设备，这些技巧还是很有用的。

增大并行性

上面说 “线程块中内层的维度（blockDim.x）过小” 是否对现在的设备还有影响，我们来看一下下面的试验：

通过这个表我们发现，最快的还是第一个，块最小的反而获得最高的效率，这里与书上的结果又不同了，我再想书上的数据量大可能会影响结果，当数据量大的时候有可能决定时间的因素会发生变化，但是一些结果是可以观察到

• 尽管（64，4） 和 （128，2） 有同样大小的块，但是执行效率不同，说明内层线程块尺寸影响效率。
• 最后的块参数无效
• 第一种方案速度最快

我们调整块的尺寸，还是为了增加并行性，或者说增加活跃的线程束，我们来看看线程束的活跃比例，得到如下数据:

可见最高的利用率没有最高的效率。没有任何一个因素可以直接左右最后的效率，一定是大家一起作用得到最终的结果，多因一效的典型例子，于是在优化的时候，我们应该首先保证测试时间的准确性，客观性，以及稳定性，说实话，我们上面的时间测试方法并不那么稳定，更稳定方法应该是测几次的平均时间，来降低人为误差。

总结

指标与性能

• 大部分情况，单一指标不能优化出最优性能
• 总体性能直接相关的是内核的代码本质（内核才是关键）
• 指标与性能之间选择平衡点
• 从不同的角度寻求指标平衡，最大化效率
• 网格和块的尺寸为调节性能提供了一个不错的起点

避免分支分化

并行规约问题

在串行编程中，我们最常见的一个问题就是一组特别多数字通过计算变成一个数字，比如加法，也就是求这一组数据的和，或者乘法，这种计算当有如下特点的时候，我们可以用并行归约的方法处理他们：

1. 结合性
2. 交换性

对应的加法或者乘法就是交换律和结合律，在我们的数学分析系列已经详细的介绍了加法和乘法的结合律和交换律的证明。所以对于所有有这两个性质的计算，都可以使用归约式计算。
为什么叫归约，归约是一种常见的计算方式（串并行都可以），归约的归有递归的意思，约就是减少，这样就很明显了，每次迭代计算方式都是相同的（归），从一组多个数据最后得到一个数（约）。
归约的方式基本包括如下几个步骤：

1. 将输入向量划分到更小的数据块中
2. 用一个线程计算一个数据块的部分和
3. 对每个数据块的部分和再求和得到最终的结果。

数据分块保证我们可以用一个线程块来处理一个数据块。
一个线程处理更小的块，所以一个线程块可以处理一个较大的块，然后多个块完成整个数据集的处理。
最后将所有线程块得到的结果相加，就是结果，这一步一般在 cpu 上完成。

归约问题最常见的加法计算是把向量的数据分成对，然后用不同线程计算每一对元素，得到的结果作为输入继续分成对，迭代的进行，直到最后一个元素。
成对的划分常见的方法有以下两种：

1. 相邻配对：元素与他们相邻的元素配对
   
2. 交错配对：元素与一定距离的元素配对
   

图中将两种方式表现的很清楚了，我们可以用代码实现以下。
首先是 cpu 版本实现交错配对归约计算的代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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);
}

这个加法运算可以改成任何满足结合律和交换律的计算。比如乘法，求最大值等。下面我们就来通过不同的配对方式，不同的数据组织来看 CUDA 的执行效率。

并行规约中的分化

有判断条件的地方就会产生分支，比如 if 和 for 这类关键词。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

__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
    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];
}

这里面唯一要注意的地方就是同步指令

1

__syncthreads();

改善并行规约的分化

1

if ((tid % (2 * stride)) == 0)

这个条件判断给内核造成了极大的分支，如图所示：

可见线程利用率是非常低的，因为这些线程在一个线程束，所以，只能等待，不能执行别的指令。
对于上面的低利用率，我们想到了下面这个方案来解决：

注意橙色圆形内的标号是线程符号，这样的计算线程的利用率是高于原始版本的，核函数如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

__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

int index = 2 * stride * tid;

这一步保证 index 能够向后移动到有数据要处理的内存位置，而不是简单的用 tid 对应内存地址，导致大量线程空闲。
首先我们保证在一个块中前几个执行的线程束是在接近满跑的，而后半部分线程束基本是不需要执行的，当一个线程束内存在分支，而分支都不需要执行的时候，硬件会停止他们调用别人，这样就节省了资源，所以高效体现在这，如果还是所有分支不满足的也要执行，即便整个线程束都不需要执行的时候，那这种方案就无效了，还好现在的硬件比较只能，于是，我们执行后得到如下结果

交错配对的规约

上面的套路是修改线程处理的数据，使部分线程束最大程度利用数据，接下来采用同样的思想，但是方法不同，接下来我们使用的方法是调整跨度，也就是我们每个线程还是处理对应的内存的位置，但内存对不是相邻的了，而是隔了一定距离的：
示意图是最好的描述方法：

我们依然把上图当做一个完整的线程块，那么前半部分的线程束依然是最大负载在跑，后半部分的线程束也是啥活不干

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

__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];
}

reduceInterleaved 内存效率是最低的，但是线程束内分化却是最小的，而书中说 reduceInterleaved 的优势在内存读取，而非线程束分化，此处需要查看机器码，确定两个内核的实际不同

展开循环

GPU 喜欢确定的东西，像前面讲解执行模型和线程束的时候，明确的指出，GPU 没有分支预测能力，所有每一个分支他都是执行的，所以在内核里尽量别写分支，分支包括啥，包括 if 当然还有 for 之类的循环语句。
如果你不知道到为啥 for 算分支语句我给你写个简单到不能运行的例子：

1

for (int i=0;i<tid;i++){}

如果上面这段代码出现在内核中，就会有分支，因为一个线程束第一个线程和最后一个线程 tid 相差 32（如果线程束大小是 32 的话） 那么每个线程执行的时候，for 终止时完成的计算量都不同，这就有人要等待，这也就产生了分支。

循环展开是一个尝试通过减少分支出现的频率和循环维护指令来优化循环的技术。不止并行算法可以展开，传统串行代码展开后效率也能一定程度的提高，因为省去了判断和分支预测失败所带来的迟滞。

1
2
3

for (int i=0;i<100;i++){
    a[i]=b[i]+c[i];
}

这个是最传统的写法，这个写法在各种 c++教材上都能看到，不做解释，如果我们进行循环展开呢？

1
2
3
4
5
6

for (int i=0;i<100;i+=4){
    a[i+0]=b[i+0]+c[i+0];
    a[i+1]=b[i+1]+c[i+1];
    a[i+2]=b[i+2]+c[i+2];
    a[i+3]=b[i+3]+c[i+3];
}

这样做的好处，从串行较多来看是减少了条件判断的次数。但是如果你把这段代码拿到机器上跑，其实看不出来啥效果，因为现代编译器把上述两个不同的写法，编译成了类似的机器语言，也就是，我们这不循环展开，编译器也会帮我们做。

不过值得注意的是：目前 CUDA 的编译器还不能帮我们做这种优化，人为的展开核函数内的循环，能够非常大的提升内核性能

在 CUDA 中展开循环的目的还是那两个：

1. 减少指令消耗

2. 增加更多的独立调度指令来提高性能

被添加到 CUDA 流水线上，是非常受欢迎的，因为其能最大限度的提高指令和内存带宽。
下面我们就在前面归约的例子上继续挖掘性能，看看是否能得到更高的效率。

展开的归约

我们的内核函数 reduceInterleaved 核函数中每个线程块只处理对应那部分的数据，我们现在的一个想法是能不能用一个线程块处理多块数据，其实这是可以实现的，如果在对这块数据进行求和前（因为总是一个线程对应一个数据）使用每个线程进行一次加法，从别的块取数据，相当于先做一个向量加法，然后再归约，这样将会用一句指令，完成之前一般的计算量。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

__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];
}

这里面的第二句，第四句，在确定线程块对应的数据的位置的时候有个乘 2 的偏移量，

这就是第二句，第四句指令的意思，我们只处理红色的线程块，而旁边白色线程块我们用

1
2
3

if(idx+blockDim.x<n){
    g_idata[idx]+=g_idata[idx+blockDim.x];
}

处理掉了，注意我们这里用的是一维线程，也就是说，我们用原来的一半的块的数量，而每一句只添加一小句指令的情况下，完成了原来全部的计算量，这个效果应该是客观的，所以我们来看一下效果之前先看一下调用核函数的部分：

1
2
3
4
5
6
7
8
9
10
11
12
13

// kernel 1:reduceUnrolling2
CHECK(cudaMemcpy(idata_dev, idata_host, bytes, cudaMemcpyHostToDevice));
CHECK(cudaDeviceSynchronize());
iStart = cpuSecond();
reduceUnroll2<<<grid.x / 2, block>>>(idata_dev, odata_dev, size);
cudaDeviceSynchronize();
iElaps = cpuSecond() - iStart;
cudaMemcpy(odata_host, odata_dev, grid.x * sizeof(int), cudaMemcpyDeviceToHost);
gpu_sum = 0;
for (int i = 0; i < grid.x / 2; i++)
    gpu_sum += odata_host[i];
printf("reduceUnrolling2            elapsed %lf ms gpu_sum: %d<<<grid %d block %d>>>\n",
       iElaps, gpu_sum, grid.x / 2, block.x);

这里需要注意由于合并了一半的线程块，这里的网格个数都要对应的减少一半，来看效率
我们上面框里有 2，4，8 三种尺度的展开，分别是一个块计算 2 个块，4 个块和 8 个块的数据。执行效率是和内存吞吐量是呈正相关的

完全展开的归约

接着我们的目标是最后那 32 个线程，因为归约运算是个倒金字塔，最后的结果是一个数，所以每个线程最后 64 个计算得到一个数字结果的过程，没执行一步，线程的利用率就降低一倍，因为从 64 到 32，然后 16。。这样到 1 的，我们现在像个办法，展开最后的 6 步迭代（64，32，16，8，4，2，1）使用下面的核函数来展开最后 6 步分支计算：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

__global__ void reduceUnrollWarp8(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;
    // unrolling 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
    for (int stride = blockDim.x / 2; stride > 32; stride >>= 1) {
        if (tid < stride) {
            idata[tid] += idata[tid + stride];
        }
        // synchronize within block
        __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];
}

在 unrolling8 的基础上，我们对于 tid 在 [0,32] 之间的线程用这个代码展开

1
2
3
4
5
6
7

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 int 类型变量我们先不说，我们先把最后这个展开捋顺一下，当只剩下最后下面三角部分，从 64 个数合并到一个数，首先将前 32 个数，按照步长为 32，进行并行加法，前 32 个 tid 得到 64 个数字的两两和，存在前 32 个数字中
然后这 32 个数加上步长为 16 的变量，理论上，这样能得到 16 个数，这 16 个数的和就是最后这个块的归约结果，但是根据上面 tid<32 的判断条件线程 tid 16 到 31 的线程还在运行，但是结果已经没意义了，这一步很重要（这一步可能产生疑惑的另一个原因是既然是同步执行，会不会比如线程 17 加上了线程 33 后写入 17 号的内存了，这时候 1 号才来加 17 号的结果，这样结果就不对了，因为我们的 CUDA 内核从内存中读数据到寄存器，然后进行加法都是同步进行的，也就是 17 号线程和 1 号线程同时读 33 号和 17 号的内存，这样 17 号即便在下一步修改，也不影响 1 号线程寄存器里面的值了），虽然 32 以内的 tid 的线程都在跑，但是没进行一步，后面一半的线程结果将没有用途了，
这样继续计算，得到最后的一个有效的结果就是 tid[0]。
上面这个过程有点复杂，但是我们自己好好想一想，从硬件取数据，到计算，每一步都分析一下，就能得到实际的结果。

volatile int 类型变量是控制变量结果写回到内存，而不是存在共享内存，或者缓存中，因为下一步的计算马上要用到它，如果写入缓存，可能造成下一步的读取会读到错误的数据

1
2

vmem[tid]+=vmem[tid+32];
vmem[tid]+=vmem[tid+16];

模板函数的归约

根据上面展开最后 64 个数据，我们可以直接就展开最后 128 个，256 个，512 个，1024 个。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

__global__ void reduceCompleteUnrollWarp8(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;
    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();
    // 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];
}

内核代码如上，这里用到了 tid 的大小，和最后 32 个没用到 tid 不同的是，这些如果计算完整会有一半是浪费的，而最后 32 个已经是线程束最小的大小了，所以无论后面的数据有没有意义，那些进程都不会停。

模板函数的归约

我们看上面这个完全展开的函数，

1
2
3
4
5
6
7
8
9
10
11
12

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();

这一步比较应该是多于的，因为 blockDim.x 自内核启动时一旦确定，就不能更改了，所以模板函数帮我解决了这个问题，当编译时编译器会去检查，blockDim.x 是否固定，如果固定，直接可以删除掉内核中不可能的部分也就是上半部分，下半部分是要执行的。

小结

我们总结下本文四步优化的指标对比：

动态并行

到目前为止，我们所有的内核都是在主机线程中调用的，那么我们肯定会想，是否我们可以在内核中调用内核，这个内核可以是别的内核，也可以是自己，那么我们就需要动态并行了，这个功能在早期的设备上是不支持的。
动态并行的好处之一就是能让复杂的内核变得有层次，坏处就是写出来的程序更复杂，因为并行行为本来就不好控制。
动态并行的另一个好处是等到执行的时候再配置创建多少个网格，多少个块，这样就可以动态的利用 GPU 硬件调度器和加载平衡器了，通过动态调整，来适应负载。并且在内核中启动内核可以减少一部分数据传输消耗。

嵌套执行

前面我们大费周章的其实也就只学了网格，块和启动配置，以及一些线程束的知识，现在我们要做的是从内核中启动内核。
内核中启动内核，和 cpu 并行中有一个相似的概念，就是父线程和子线程。子线程由父线程启动，但是到了 GPU，这类名词相对多了些，比如父网格，父线程块，父线程，对应的子网格，子线程块，子线程。子网格被父线程启动，且必须在对应的父线程，父线程块，父网格结束之前结束。所有的子网格结束后，父线程，父线程块，父网格才会结束。

上图清晰地表明了父网格和子网格的使用情况，一种典型的执行方式：

主机启动一个网格（也就是一个内核）-> 此网格（父网格）在执行的过程中启动新的网格（子网格们）->所有子网格们都运行结束后-> 父网格才能结束，否则要等待

如果调用的线程没有显示同步启动子网格，那么运行时保证，父网格和子网格隐式同步。
图中显式的同步了父网格和子网格，通过设置栅栏的方法。
父网格中的不同线程启动的不同子网格，这些子网格拥有相同的父线程块，他们之间是可以同步的。线程块中所有的线程创建的所有子网格完成之后，线程块执行才会完成。如果块中的所有线程在子网格完成前退出，那么子网格隐式同步会被触发。隐式同步就是虽然没用同步指令，但是父线程块中虽然所有线程都执行完毕，但是依旧要等待对应的所有子网格执行完毕，然后才能退出。
前面我们讲过隐式同步，比如 cudaMemcpy 就能起到隐式同步的作用，但是主机内启动的网格，如果没有显式同步，也没有隐式同步指令，那么 cpu 线程很有可能就真的退出了，而你的 gpu 程序可能还在运行，这样就非常尴尬了。父线程块启动子网格需要显示的同步，也就是说不通的线程束需要都执行到子网格调用那一句，这个线程块内的所有子网格才能依据所在线程束的执行，一次执行。
接着是最头疼的内存，内存竞争对于普通并行就很麻烦了，现在对于动态并行，更麻烦，主要的有下面几点：

1. 父网格和子网格共享相同的全局和常量内存。
2. 父网格子网格有不同的局部内存
3. 有了子网格和父网格间的弱一致性作为保证，父网格和子网格可以对全局内存并发存取。
4. 有两个时刻父网格和子网格所见内存一致：子网格启动的时候，子网格结束的时候
5. 共享内存和局部内存分别对于线程块和线程来说是私有的
6. 局部内存对线程私有，对外不可见。

在 GPU 上嵌套 Hello World

为了研究初步动态并行，我们先来写个 Hello World 进行操作，代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

#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 = 64;
    int block_x = 2;
    dim3 block(block_x, 1);
    dim3 grid((size - 1) / block.x + 1, 1);
    nesthelloworld<<<grid, block>>>(size, 0);
    cudaGetLastError();
    cudaDeviceReset();
    return 0;
}

这就是完成可执行代码，编译的命令与之前有些不同，工程中使用 cmake 管理，但是本程序没有纳入其中，而是使用了一个单独的 makefile

1

nvcc -arch=sm_35  nested_Hello_World.cu -o nested_Hello_World -lcudadevrt --relocatable-device-code true

-lcudadevrt –relocatable-device-code true 是前面没有的，这两个指令是动态并行需要的一个库，relocatable-device-code 表示生成可重新定位的代码。

可见，当多层调用子网格的时候，同一家的（就是用相同祖宗线程的子网）是隐式同步的，而不同宗的则是各跑各的。

内存模型概述

CUDA 采用的内存模型，结合了主机和设备内存系统，展现了完整的内存层次模型，其中大部分内存我们可以通过编程控制，来使我们的程序性能得到优化。

内存层次结构的优点

程序具有局部性特点，包括：

1. 时间局部性
2. 空间局部性

时间局部性，就是一个内存位置的数据某时刻被引用，那么在此时刻附近也很有可能被引用，随时间流逝，该数据被引用的可能性逐渐降低。
空间局部性，如果某一内存位置的数据被使用，那么附近的数据也有可能被使用。
现代计算机的内存结构主要如下：

这个内存模型在程序局部性原则成立的时候有效。学习过串行编程的人也应该知道内存模型，速度最快的是寄存器，他能和 cpu 同步的配合，接着是缓存，在 CPU 片上，然后是主存储器，现在常见的就是内存条，显卡上也有内存芯片，然后是硬盘，这些内存设备的速度和容量相反，越快的越小，越慢的越大。
局部性是个非常有趣的事情，首先局部性的产生并不是因为设备的原因，而是程序从一开始被编写就有这个特征，与生俱来，所以当我们发现此特征后，就开始设计满足此特征硬件结构，也就是内存模型，当内存模型设计成如上结构的时候，如果你想写快速高效的程序，就要让自己的程序局部性足够好，所以这就进入了一个死循环，最后为了追求高效率，设备将越来越优化局部性，而程序也会越来越局部化。
总结下最后一层（硬盘磁带之类的）的特点：

• 每个比特位的价格要更低
• 容量要更高
• 延迟较高
• 处理器访问频率低

CPU 和 GPU 的主存都是采用 DRAM——动态随机存取存储器，而低延迟的内存，比如一级缓存，则采用 SRAM——静态随机存取存储器。虽然底层的存储器延迟高，容量大，但是其中有数据被频繁使用的时候，就会向更高一级的层次传输，比如我们运行程序处理数据的时候，程序第一步就是把硬盘里的数据传输到主存里面。
GPU 和 CPU 的内存设计有相似的准则和模型。但他们的区别是：CUDA 编程模型将内存层次结构更好的呈献给开发者，让我们显示的控制其行为。

CUDA 内存模型

对于程序员来说，分类内存的方法有很多中，但是对于我们来说最一般的分法是：

• 可编程内存
• 不可编程内存

对于可编程内存，如字面意思，你可以用你的代码来控制这组内存的行为；相反的，不可编程内存是不对用户开放的，也就是说其行为在出厂后就已经固化了，对于不可编程内存，我们能做的就是了解其原理，尽可能的利用规则来加速程序，但对于通过调整代码提升速度来说，效果很一般。
CPU 内存结构中，一级二级缓存都是不可编程（完全不可控制）的存储设备。
另一方面，CUDA 内存模型相对于 CPU 来说那是相当丰富了，GPU 上的内存设备有：

• 寄存器
• 共享内存
• 本地内存
• 常量内存
• 纹理内存
• 全局内存

上述各种都有自己的作用域，生命周期和缓存行为。CUDA 中每个线程都有自己的私有的本地内存；线程块有自己的共享内存，对线程块内所有线程可见；所有线程都能访问读取常量内存和纹理内存，但是不能写，因为他们是只读的；全局内存，常量内存和纹理内存空间有不同的用途。对于一个应用来说，全局内存，常量内存和纹理内存有相同的生命周期。下图总结了上面这段话，后面的大篇幅文章就是挨个介绍这些内存的性质和使用的。

寄存器

寄存器无论是在 CPU 还是在 GPU 都是速度最快的内存空间，但是和 CPU 不同的是 GPU 的寄存器储量要多一些，而且当我们在核函数内不加修饰的声明一个变量，此变量就存储在寄存器中，但是 CPU 运行的程序有些不同，只有当前在计算的变量存储在寄存器中，其余在主存中，使用时传输至寄存器。在核函数中定义的有常数长度的数组也是在寄存器中分配地址的。
寄存器对于每个线程是私有的，寄存器通常保存被频繁使用的私有变量，注意这里的变量也一定不能使共有的，不然的话彼此之间不可见，就会导致大家同时改变一个变量而互相不知道，寄存器变量的生命周期和核函数一致，从开始运行到运行结束，执行完毕后，寄存器就不能访问了。
寄存器是 SM 中的稀缺资源，Fermi 架构中每个线程最多 63 个寄存器。Kepler 结构扩展到 255 个寄存器，一个线程如果使用更少的寄存器，那么就会有更多的常驻线程块，SM 上并发的线程块越多，效率越高，性能和使用率也就越高。
那么问题就来了，如果一个线程里面的变量太多，以至于寄存器完全不够呢？这时候寄存器发生溢出，本地内存就会过来帮忙存储多出来的变量，这种情况会对效率产生非常负面的影响，所以，不到万不得已，一定要避免此种情况发生。
为了避免寄存器溢出，可以在核函数的代码中配置额外的信息来辅助编译器优化，比如：

1
2
3

__global__ void
__lauch_bounds__(maxThreadaPerBlock,minBlocksPerMultiprocessor)
kernel(...) {}

这里面在核函数定义前加了一个 关键字 lauch_bounds，然后他后面对应了两个变量：

1. maxThreadaPerBlock：线程块内包含的最大线程数，线程块由核函数来启动

2. minBlocksPerMultiprocessor：可选参数，每个 SM 中预期的最小的常驻内存块参数。

   注意，对于一定的核函数，优化的启动边界会因为不同的结构而不同

   也可以在编译选项中加入

1

-maxrregcount=32

来控制一个编译单元里所有核函数使用的最大数量。

本地内存

核函数中符合存储在寄存器中但不能进入被核函数分配的寄存器空间中的变量将存储在本地内存中，编译器可能存放在本地内存中的变量有以下几种：

• 使用未知索引引用的本地数组
• 可能会占用大量寄存器空间的较大本地数组或者结构体
• 任何不满足核函数寄存器限定条件的变量

本地内存实质上是和全局内存一样在同一块存储区域当中的，其访问特点——高延迟，低带宽。对于 2.0 以上的设备，本地内存存储在每个 SM 的一级缓存，或者设备的二级缓存上。

共享内存

在核函数中使用如下修饰符的内存，称为共享内存：

1

__share__

每个 SM 都有一定数量的由线程块分配的共享内存，共享内存是片上内存，跟主存相比，速度要快很多，也即是延迟低，带宽高。其类似于一级缓存，但是可以被编程。

使用共享内存的时候一定要注意，不要因为过度使用共享内存，而导致 SM 上活跃的线程束减少，也就是说，一个线程块使用的共享内存过多，导致更过的线程块没办法被 SM 启动，这样影响活跃的线程束数量。

共享内存在核函数内声明，生命周期和线程块一致，线程块运行开始，此块的共享内存被分配，当此块结束，则共享内存被释放。
因为共享内存是块内线程可见的，所以就有竞争问题的存在，也可以通过共享内存进行通信，当然，为了避免内存竞争，可以使用同步语句：

1

void __syncthreads();

此语句相当于在线程块执行时各个线程的一个障碍点，当块内所有线程都执行到本障碍点的时候才能进行下一步的计算，这样可以设计出避免内存竞争的共享内存使用程序、

注意，__syncthreads();频繁使用会影响内核执行效率。

SM 中的一级缓存，和共享内存共享一个 64k 的片上内存（不知道现在的设备有没有提高），他们通过静态划分，划分彼此的容量，运行时可以通过下面语句进行设置：

1

cudaError_t cudaFuncSetCacheConfig(const void * func,enum cudaFuncCache);

这个函数可以设置内核的共享内存和一级缓存之间的比例。cudaFuncCache 参数可选如下配置：

1
2
3
4

cudaFuncCachePreferNone//无参考值，默认设置
cudaFuncCachePreferShared//48k共享内存，16k一级缓存
cudaFuncCachePreferL1// 48k一级缓存，16k共享内存
cudaFuncCachePreferEqual// 32k一级缓存，32k共享内存

Fermi 架构支持前三种，后面的设备都支持。

常量内存

常量内存驻留在设备内存中，每个 SM 都有专用的常量内存缓存，常量内存使用：__constant__修饰，常量内存在核函数外，全局范围内声明，对于所有设备，只可以声明 64k 的常量内存，常量内存静态声明，并对同一编译单元中的所有核函数可见。
叫常量内存，显然是不能被修改的，这里不能被修改指的是被核函数修改，主机端代码是可以初始化常量内存的，不然这个内存谁都不能改就没有什么使用意义了，常量内存，被主机端初始化后不能被核函数修改，初始化函数如下：

1

cudaError_t cudaMemcpyToSymbol(const void* dst,const void *src,size_t count);

同 cudaMemcpy 的参数列表相似，从 src 复制 count 个字节的内存到 dst 里面，也就是设备端的常量内存。多数情况下此函数是同步的，也就是会马上被执行。
当线程束中所有线程都从相同的地址取数据时，常量内存表现较好，比如执行某一个多项式计算，系数都存在常量内存里效率会非常高，但是如果不同的线程取不同地址的数据，常量内存就不那么好了，因为常量内存的读取机制是——一次读取会广播给所有线程束内的线程。

纹理内存

纹理内存驻留在设备内存中，在每个 SM 的只读缓存中缓存，纹理内存是通过指定的缓存访问的全局内存，只读缓存包括硬件滤波的支持，它可以将浮点插入作为读取过程中的一部分来执行，纹理内存是对二维空间局部性的优化。
总的来说纹理内存设计目的应该是为了 GPU 本职工作显示设计的，但是对于某些特定的程序可能效果更好，比如需要滤波的程序，可以直接通过硬件完成。

全局内存

GPU 上最大的内存空间，延迟最高，使用最常见的内存，global 指的是作用域和生命周期，一般在主机端代码里定义，也可以在设备端定义，不过需要加修饰符，只要不销毁，是和应用程序同生命周期的。全局内存对应于设备内存，一个是逻辑表示，一个是硬件表示
全局内存可以动态声明，或者静态声明，可以用下面的修饰符在设备代码中静态的声明一个变量：

1

__device__

我们前面声明的所有的在 GPU 上访问的内存都是全局内存，或者说到目前为止我们还没对内存进行任何优化。
因为全局内存的性质，当有多个核函数同时执行的时候，如果使用到了同一全局变量，应注意内存竞争。
全局内存访问是对齐，也就是一次要读取指定大小（32，64，128）整数倍字节的内存，所以当线程束执行内存加载/存储时，需要满足的传输数量通常取决与以下两个因素：

1. 跨线程的内存地址分布
2. 内存事务的对齐方式。

一般情况下满足内存请求的事务越多，未使用的字节被传输的可能性越大，数据吞吐量就会降低，换句话说，对齐的读写模式使得不需要的数据也被传输，所以，利用率低到时吞吐量下降。1.1 以下的设备对内存访问要求非常严格（为了达到高效，访问受到限制）因为当时还没有缓存，现在的设备都有缓存了，所以宽松了一些。

GPU 缓存

与 CPU 缓存类似，GPU 缓存不可编程，其行为出厂时已经设定好了。GPU 上有 4 种缓存：

1. 一级缓存
2. 二级缓存
3. 只读常量缓存
4. 只读纹理缓存

每个 SM 都有一个一级缓存，所有 SM 公用一个二级缓存。一级二级缓存的作用都是被用来存储本地内存和全局内存中的数据，也包括寄存器溢出的部分。Fermi，Kepler 以及以后的设备，CUDA 允许我们配置读操作的数据是使用一级缓存和二级缓存，还是只使用二级缓存。
与 CPU 不同的是，CPU 读写过程都有可能被缓存，但是 GPU 写的过程不被缓存，只有加载会被缓存！
每个 SM 有一个只读常量缓存，只读纹理缓存，它们用于设备内存中提高来自于各自内存空间内的读取性能。

CUDA 变量声明总结

用表格进行总结：

设备存储器的重要特征：

静态全局内存

CPU 内存有动态分配和静态分配两种类型，从内存位置来说，动态分配在堆上进行，静态分配在站上进行，在代码上的表现是一个需要 new，malloc 等类似的函数动态分配空间，并用 delete 和 free 来释放。在 CUDA 中也有类似的动态静态之分，我们前面用的都是要 cudaMalloc 的，所以对比来说就是动态分配，我们今天来个静态分配的，不过与动态分配相同是，也需要显式的将内存 copy 到设备端，我们用下面代码来看一下程序的运行结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

#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;
}

运行结果

这个唯一要注意的就是，这一句

1

cudaMemcpyToSymbol(devData,&value,sizeof(float));

函数原型说的是第一个应该是个 void*，但是这里写了一个__device__ float devData;变量，这个说到底还是设备上的变量定义和主机变量定义的不同，设备变量在代码中定义的时候其实就是一个指针，这个指针指向何处，主机端是不知道的，指向的内容也不知道，想知道指向的内容，唯一的办法还是通过显式的办法传输过来：

1

cudaMemcpyFromSymbol(&value,devData,sizeof(float));

这里需要注意的只有这点：

1. 在主机端，devData 只是一个标识符，不是设备全局内存的变量地址

2. 在核函数中，devData 就是一个全局内存中的变量。

   主机代码不能直接访问设备变量，设备也不能访问主机变量，这就是 CUDA 编程与 CPU 多核最大的不同之处

1

cudaMemcpy(&value,devData,sizeof(float));

是不可以的！这个函数是无效的！就是你不能用动态 copy 的方法给静态变量赋值！
如果你死活都要用 cudaMemcpy，只能用下面的方式：

1
2
3

float *dptr=NULL;
cudaGetSymbolAddress((void**)&dptr,devData);
cudaMemcpy(dptr,&value,sizeof(float),cudaMemcpyHostToDevice);

主机端不可以对设备变量进行取地址操作！这是非法的！
想要得到 devData 的地址可以用下面方法：

1
2

float *dptr=NULL;
cudaGetSymbolAddress((void**)&dptr,devData);

当然也有一个例外，可以直接从主机引用 GPU 内存——CUDA 固定内存。后面我们会研究这部分。
CUDA 运行时 API 能访问主机和设备变量，但这取决于你给正确的函数是否提供了正确的参数，使用运行时 API，如果参数填错，尤其是主机和设备上的指针，结果是无法预测的。

内存管理

CUDA 编程的目的是给我们的程序加速，我们在控制硬件，控制硬件的语言属于底层语言，比如 C 语言，最头疼的就是管理内存，python，php 这些语言有自己的内存管理机制，c 语言的内存管理机制——程序员管理。CUDA 是 C 语言的扩展，内存方面基本集成了 C 语言的方式，由程序员控制 CUDA 内存，当然，这些内存的物理设备是在 GPU 上的，而且与 CPU 内存分配不同，CPU 内存分配完就完事了，GPU 还涉及到数据传输，主机和设备之间的传输。
接下来我们要了解的是：

• 分配释放设备内存
• 在主机和设备间传输内存

为达到最优性能，CUDA 提供了在主机端准备设备内存的函数，并且显式地向设备传递数据，显式的从设备取回数据。

内存分配和释放

1

cudaError_t cudaMalloc(void ** devPtr,size_t count)

这个函数用过很多次了，唯一要注意的是第一个参数，是指针的指针，一般的用法是首先我们生命一个指针变量，然后调用这个函数：

1
2

float * devMem=NULL;
cudaError_t cudaMalloc((float**) devMem, count)

devMem 是一个指针，定义时初始化指向 NULL，这样做是安全的，避免出现野指针

cudaMalloc 函数要修改 devMem 的值，所以必须把他的指针传递给函数，如果把 devMem 当做参数传递，经过函数后，指针的内容还是 NULL。
不知道这个解释有没有听明白，通俗的讲，如果一个参数想要在函数中被修改，那么一定要传递他的地址给函数，如果只传递本身，函数是值传递的，不会改变参数的值。
当分配完地址后，可以使用下面函数进行初始化：

1

cudaError_t cudaMemset(void * devPtr,int value,size_t count)

用法和 Memset 类似，但是注意，这些被我们操作的内存对应的物理内存都在 GPU 上。
当分配的内存不被使用时，使用下面语句释放程序。

1

cudaError_t cudaFree(void * devPtr)

注意这个参数一定是前面 cudaMalloc 类的函数（还有其他分配函数）分配到空间，如果输入非法指针参数，会返回 cudaErrorInvalidDevicePointer 错误，如果重复释放一个空间，也会报错。目前为止，套路基本和 C 语言一致。

注:设备内存的分配和释放非常影响性能。

内存传输

下面介绍点 C 语言没有的，C 语言的内存分配完成后就可以直接读写了，但是对于异构计算，这样是不行的，因为主机线程不能访问设备内存，设备线程也不能访问主机内存，这时候我们要传送数据了：

1

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

这个函数我们前面也反复用到，注意这里的参数是指针，而不是指针的指针，第一个参数 dst 是目标地址，第二个参数 src 是原始地址，然后是拷贝的内存大小，最后是传输类型，传输类型包括以下几种：

• cudaMemcpyHostToHost
• cudaMemcpyHostToDevice
• cudaMemcpyDeviceToHost
• cudaMemcpyDeviceToDevice

四种方式，都写在字面上来，唯一有点问题的就是有个 host 到 host，不知道为啥存在，估计很多人跟我想法一样，可能后面有什么高级的用法。
这个例子也不用说了，前面随便找个有数据传输的都有这两步：从主机到设备，然后计算，最后从设备到主机。
代码省略，来张图：

GPU 的内存采用的 DDR5 制式，2011 三星才做出来 DDR4 的主机内存，但是 GPU 却一直在使用 DDR5，GPU 的内存理论峰值带宽非常高，对于 Fermi C2050 有 144GB/s，这个值估计现在的 GPU 应该都超过了，CPU 和 GPU 之间通信要经过 PCIe 总线，总线的理论峰值要低很多——8GB/s 左右，也就是说所，管理不当，算到半路需要从主机读数据，那效率瞬间全挂在 PCIe 上了。

CUDA 编程需要减少主机和设备之间的内存传输。

固定内存

主机内存采用分页式管理，通俗的说法就是操作系统把物理内存分成一些“页”，然后给一个应用程序一大块内存，但是这一大块内存可能在一些不连续的页上，应用只能看到虚拟的内存地址，而操作系统可能随时更换物理地址的页（从原始地址复制到另一个地址）但是应用是不会差觉得，但是从主机传输到设备上的时候，如果此时发生了页面移动，对于传输操作来说是致命的，所以在数据传输之前，CUDA 驱动会锁定页面，或者直接分配固定的主机内存，将主机源数据复制到固定内存上，然后从固定内存传输数据到设备上：

上图左边是正常分配内存，传输过程是：锁页-复制到固定内存-复制到设备
右边时分配时就是固定内存，直接传输到设备上。
下面函数用来分配固定内存：

1

cudaError_t cudaMallocHost(void ** devPtr,size_t count)

分配 count 字节的固定内存，这些内存是页面锁定的，可以直接传输到设备的,这样就是的传输带宽变得高很多。
固定的主机内存释放使用：

1

cudaError_t cudaFreeHost(void *ptr)

我们可以测试一下固定内存和分页内存的传输效率，代码如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

#include "freshman.h"
#include <cuda_runtime.h>
#include <stdio.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 = blockIdx.x * blockDim.x + threadIdx.x;
    res[i] = a[i] + b[i];
}
int main(int argc, char **argv) {
    int dev = 0;
    cudaSetDevice(dev);

    int nElem = 1 << 14;
    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;
    // 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));

    dim3 block(1024);
    dim3 grid(nElem / block.x);
    sumArraysGPU<<<grid, block>>>(a_d, b_d, res_d);
    printf("Execution configuration<<<%d,%d>>>\n", grid.x, block.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);
    cudaFreeHost(a_d);
    cudaFreeHost(b_d);
    cudaFreeHost(res_d);

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

    return 0;
}

结论：

• 固定内存的释放和分配成本比可分页内存要高很多，但是传输速度更快，所以对于大规模数据，固定内存效率更高。
• 尽量使用流来使内存传输和计算之间同时进行。

零拷贝内存

截止到目前，我们所接触到的内存知识的基础都是：主机直接不能访问设备内存，设备不能直接访问主机内存。对于早期设备，这是肯定的，但是后来，一个例外出现了——零拷贝内存。
GPU 线程可以直接访问零拷贝内存，这部分内存在主机内存里面，CUDA 核函数使用零拷贝内存有以下几种情况：

• 当设备内存不足的时候可以利用主机内存
• 避免主机和设备之间的显式内存传输
• 提高 PCIe 传输率

前面我们讲，注意线程之间的内存竞争，因为他们可以同时访问同一个内存地址，现在设备和主机可以同时访问同一个设备地址了，所以，我们要注意主机和设备的内存竞争——当使用零拷贝内存的时候。
零拷贝内存是固定内存，不可分页。可以通过以下函数创建零拷贝内存：

1

cudaError_t cudaHostAlloc(void ** pHost,size_t count,unsigned int flags)

最后一个标志参数，可以选择以下值：

• cudaHostAllocDefalt

• cudaHostAllocPortable

• cudaHostAllocWriteCombined

• cudaHostAllocMapped

  cudaHostAllocDefalt 和 cudaMallocHost 函数一致，cudaHostAllocPortable 函数返回能被所有 CUDA 上下文使用的固定内存，cudaHostAllocWriteCombined 返回写结合内存，在某些设备上这种内存传输效率更高。cudaHostAllocMapped 产生零拷贝内存。

  注意，零拷贝内存虽然不需要显式的传递到设备上，但是设备还不能通过 pHost 直接访问对应的内存地址，设备需要访问主机上的零拷贝内存，需要先获得另一个地址，这个地址帮助设备访问到主机对应的内存，方法是：

  1	cudaError_t cudaHostGetDevicePointer(void ** pDevice,void * pHost,unsigned flags);

pDevice 就是设备上访问主机零拷贝内存的指针了！此处 flag 必须设置为 0。
零拷贝内存可以当做比设备主存储器更慢的一个设备。频繁的读写，零拷贝内存效率极低，这个非常容易理解，因为每次都要经过 PCIe，速度肯定慢。
我们把实验对比结果写在一个表里面：

这是通过观察运行时间得到的，当然也可以通过我们上面的 nvprof 得到内核执行时间：

直接上数据，但是这种比较方法有点问题，因为零拷贝内存在执行内核的时候相当于还执行了内存传输工作，所以我觉得应该把内存传输也加上，那样看速度就基本差不多了，但是如果常规内存完成传输后可以重复利用，那又是另一回事了。

但是零拷贝内存也有例外的时候，比如当 CPU 和 GPU 继承在一起的时候，他们的物理内存公用的，这时候零拷贝内存，效果相当不错。但是如果离散架构，主机和设备之间通过 PCIe 连接，那么零拷贝内存将会非常耗时。

统一虚拟寻址

设备架构 2.0 以后，Nvida 又有新创意，他们搞了一套称为同一寻址方式（UVA）的内存机制，这样，设备内存和主机内存被映射到同一虚拟内存地址中。如图

UVA 之前，我们要管理所有的设备和主机内存，尤其是他们的指针，零拷贝内存尤其麻烦。有了 UVA 再也不用怕，一个人一个名，走到哪里都能用，通过 UVA，cudaHostAlloc 函数分配的固定主机内存具有相同的主机和设备地址，可以直接将返回的地址传递给核函数。
前面的零拷贝内存，可以知道以下几个方面：

• 分配映射的固定主机内存
• 使用 CUDA 运行时函数获取映射到固定内存的设备指针
• 将设备指针传递给核函数

有了 UVA，可以不用上面的那个获得设备上访问零拷贝内存的函数了：

1

cudaError_t cudaHostGetDevicePointer(void ** pDevice,void * pHost,unsigned flags);

UVA 来了以后，此函数基本失业了。
试验，代码：

1
2
3
4
5
6
7
8
9
10
11
12

float *a_host, *b_host, *res_d;
CHECK(cudaHostAlloc((float **)&a_host, nByte, cudaHostAllocMapped));
CHECK(cudaHostAlloc((float **)&b_host, nByte, cudaHostAllocMapped));
CHECK(cudaMalloc((float **)&res_d, nByte));
res_from_gpu_h = (float *)malloc(nByte);

initialData(a_host, nElem);
initialData(b_host, nElem);

dim3 block(1024);
dim3 grid(nElem / block.x);
sumArraysGPU<<<grid, block>>>(a_host, b_host, res_d);

UVA 代码主要就是差个获取指针，UVA 可以直接使用主机端的地址。

统一内存寻址

Nvidia CUDA6.0 的时候又来了个统一内存寻址，提出的目的也是为了简化内存管理统一内存中创建一个托管内存池（CPU 上有，GPU 上也有），内存池中已分配的空间可以通过相同的指针直接被 CPU 和 GPU 访问，底层系统在统一的内存空间中自动的进行设备和主机间的传输。数据传输对应用是透明的，大大简化了代码。
搞个内存池，这部分内存用一个指针同时表示主机和设备内存地址，依赖于 UVA 但是是完全不同的技术。
统一内存寻址提供了一个“指针到数据”的编程模型，概念上类似于零拷贝，但是零拷贝内存的分配是在主机上完成的，而且需要互相传输，但是统一寻址不同。
托管内存是指底层系统自动分配的统一内存，未托管内存就是我们自己分配的内存，这时候对于核函数，可以传递给他两种类型的内存，已托管和未托管内存，可以同时传递。

托管内存可以是静态的，也可以是动态的，添加 managed 关键字修饰托管内存变量。静态声明的托管内存作用域是文件

托管内存分配方式：

1

cudaError_t cudaMallocManaged(void ** devPtr,size_t size,unsigned int flags=0)

CUDA6.0 中设备代码不能调用 cudaMallocManaged，只能主机调用，所有托管内存必须在主机代码上动态声明，或者全局静态声明

内存访问模式

多数 GPU 程序容易受到内存带宽的限制，所以最大程度的利用全局内存带宽，提高全局加载效率，是调控内核函数性能的基本条件。如果不能正确调控全局内存使用，那么优化方案可能收效甚微。
CUDA 执行模型告诉我们，CUDA 执行的基本单位是线程束，所以，内存访问也是以线程束为基本单位发布和执行的，存储也一致。我们本文研究的就是这一个线程束的内存访问，不同线程的内存请求，其目标位置的不同，可以产生非常多种情况。所以本篇就是研究这些不同情况的，以及如何实现最佳的全局内存访问。

注意：访问可以是加载，也可以是存储。

对齐与合并访问

全局内存通过缓存实现加载和存储的过程如下图

全局内存是一个逻辑层面的模型，我们编程的时候有两种模型考虑：一种是逻辑层面的，也就是我们在写程序的时候（包括串行程序和并行程序），写的一维（多维）数组，结构体，定义的变量，这些都是在逻辑层面的；一种是硬件角度，就是一块 DRAM 上的电信号，以及最底层内存驱动代码所完成数字信号的处理。
L1 表示一级缓存，每个 SM 都有自己 L1，但是 L2 是所有 SM 公用的，除了 L1 缓存外，还有只读缓存和常量缓存，这个我们后面会详细介绍。
核函数运行时需要从全局内存（DRAM)中读取数据，只有两种粒度，这个是关键的：

• 128 字节
• 32 字节

解释下“粒度”，可以理解为最小单位，也就是核函数运行时每次读内存，哪怕是读一个字节的变量，也要读 128 字节，或者 32 字节，而具体是到底是 32 还是 128 还是要看访问方式：

• 使用一级缓存
• 不使用一级缓存

对于 CPU 来说，一级缓存或者二级缓存是不能被编程的，但是 CUDA 是支持通过编译指令停用一级缓存的。如果启用一级缓存，那么每次从 DRAM 上加载数据的粒度是 128 字节，如果不适用一级缓存，只是用二级缓存，那么粒度是 32 字节。
还要强调一下 CUDA 内存模型的内存读写，我们现在讨论的都是单个 SM 上的情况，多个 SM 只是下面我们描述的情形的复制：SM 执行的基础是线程束，也就是说，当一个 SM 中正在被执行的某个线程需要访问内存，那么，和它同线程束的其他 31 个线程也要访问内存，这个基础就表示，即使每个线程只访问一个字节，那么在执行的时候，只要有内存请求，至少是 32 个字节，所以不使用一级缓存的内存加载，一次粒度是 32 字节而不是更小。
在优化内存的时候，我们要最关注的是以下两个特性

• 对齐内存访问
• 合并内存访问

我们把一次内存请求——也就是从内核函数发起请求，到硬件响应返回数据这个过程称为一个内存事务（加载和存储都行）。
当一个内存事务的首个访问地址是缓存粒度（32 或 128 字节）的偶数倍的时候：比如二级缓存 32 字节的偶数倍 64，128 字节的偶数倍 256 的时候，这个时候被称为对齐内存访问，非对齐访问就是除上述的其他情况，非对齐的内存访问会造成带宽浪费。
当一个线程束内的线程访问的内存都在一个内存块里的时候，就会出现合并访问。
对齐合并访问的状态是理想化的，也是最高速的访问方式，当线程束内的所有线程访问的数据在一个内存块，并且数据是从内存块的首地址开始被需要的，那么对齐合并访问出现了。为了最大化全局内存访问的理想状态，尽量将线程束访问内存组织成对齐合并的方式，这样的效率是最高的。下面看一个例子。

• 一个线程束加载数据，使用一级缓存，并且这个事务所请求的所有数据在一个 128 字节的对齐的地址段上（对齐的地址段是我自己发明的名字，就是首地址是粒度的偶数倍，那么上面这句话的意思是，所有请求的数据在某个首地址是粒度偶数倍的后 128 个字节里），具体形式如下图，这里请求的数据是连续的，其实可以不连续，但是不要越界就好。
  
  上面蓝色表示全局内存，下面橙色是线程束要的数据，绿色就是我称为对齐的地址段。
• 如果一个事务加载的数据分布在不一个对齐的地址段上，就会有以下两种情况：
  1. 连续的，但是不在一个对齐的段上，比如，请求访问的数据分布在内存地址 1~128，那么 0~127 和 128~255 这两段数据要传递两次到 SM
  2. 不连续的，也不在一个对齐的段上，比如，请求访问的数据分布在内存地址 0~63 和 128~191 上，明显这也需要两次加载。
     
     上图就是典型的一个线程束，数据分散开了，thread0 的请求在 128 之前，后面还有请求在 256 之后，所以需要三个内存事务，而利用率，也就是从主存取回来的数据被使用到的比例，只有 128/128×3/128 的比例。这个比例低会造成带宽的浪费，最极端的表现，就是如果每个线程的请求都在不同的段，也就是一个 128 字节的事务只有 1 个字节是有用的，那么利用率只有 1/128

这里总结一下内存事务的优化关键：用最少的事务次数满足最多的内存请求。事务数量和吞吐量的需求随设备的计算能力变化。

全局内存读取

SM 加载数据，根据不同的设备和类型分为三种路径：

1. 一级和二级缓存
2. 常量缓存
3. 只读缓存

常规的路径是一级和二级缓存，需要使用常量和只读缓存的需要在代码中显式声明。但是提高性能，主要还是要取决于访问模式。
控制全局加载操作是否通过一级缓存可以通过编译选项来控制，当然比较老的设备可能就没有一级缓存。
编译器禁用一级缓存的选项是：

1

-Xptxas -dlcm=cg

编译器启用一级缓存的选项是：

1

-Xptxas -dlcm=ca

当一级缓存被禁用的时候，对全局内存的加载请求直接进入二级缓存，如果二级缓存缺失，则由 DRAM 完成请求。
每次内存事务可由一个两个或者四个部分执行，每个部分有 32 个字节，也就是 32，64 或者 128 字节一次（注意前面我们讲到是否使用一级缓存决定了读取粒度是 128 还是 32 字节，这里增加的 64 并不在此情况，所以需要注意）。

启用一级缓存后，当 SM 有全局加载请求会首先通过尝试一级缓存，如果一级缓存缺失，则尝试二级缓存，如果二级缓存也没有，那么直接 DRAM。在有些设备上一级缓存不用来缓存全局内存访问，而是只用来存储寄存器溢出的本地数据，比如 Kepler 的 K10,K20。

内存加载可以分为两类：

• 缓存加载
• 没有缓存的加载

内存访问有以下特点：

• 是否使用缓存：一级缓存是否介入加载过程
• 对齐与非对齐的：如果访问的第一个地址是 32 的倍数（前面说是 32 或者 128 的偶数倍，这里似乎产生了矛盾，为什么我现在也很迷惑）
• 合并与非合并，访问连续数据块则是合并的

缓存加载

下面是使用一级缓存的加载过程，图片表达很清楚，我们只用少量文字进行说明：

1. 对齐合并的访问，利用率 100%
   
2. 对齐的，但是不是连续的，每个线程访问的数据都在一个块内，但是位置是交叉的，利用率 100%
   
3. 连续非对齐的，线程束请求一个连续的非对齐的，32 个 4 字节数据，那么会出现，数据横跨两个块，但是没有对齐，当启用一级缓存的时候，就要两个 128 字节的事务来完成
   
4. 线程束所有线程请求同一个地址，那么肯定落在一个缓存行范围（缓存行的概念没提到过，就是主存上一个可以被一次读到缓存中的一段数据。），那么如果按照请求的是 4 字节数据来说，使用一级缓存的利用率是 4/128=3.125%
   
5. 比较坏的情况，前面提到过最坏的，就是每个线程束内的线程请求的都是不同的缓存行内，这里比较坏的情况就是，所有数据分布在 N 个缓存行上，其中 1≤N≤32，那么请求 32 个 4 字节的数据，就需要 N 个事务来完成，利用率也是 1N
   

CPU 和 GPU 的一级缓存有显著的差异，GPU 的一级缓存可以通过编译选项等控制，CPU 不可以，而且 CPU 的一级缓存是的替换算法是有使用频率和时间局部性的，GPU 则没有。

没有缓存的加载

没有缓存的加载是指的没有通过一级缓存，二级缓存则是不得不经过的。
当不使用一级缓存的时候，内存事务的粒度变为 32 字节，更细粒度的好处是提高利用律，这个很好理解，比如你每次喝水只能选择一瓶大瓶 500ml 的或则一个小瓶的 250ml，当你非常渴的时候需要 400ml 水分，喝大瓶的，比较方便，因为如果喝小瓶的一瓶不够，还需要再喝一瓶，此时大瓶的方便.但如果你需要 200ml 的水分的时候，小瓶的利用率就高很多。细粒度的访问就是用小瓶喝水，虽然体积小，但是每次的利用率都高了不少，针对上面使用缓存的情况 5，可能效果会更好。
继续我们的图解：

1. 对齐合并访问 128 字节，不用说，还是最理想的情况，使用 4 个段，利用率 100%100%
   
2. 对齐不连续访问 128 字节，都在四个段内，且互不相同，这样的利用率也是 100%100%
   
3. 连续不对齐，一个段 32 字节，所以，一个连续的 128 字节的请求，即使不对齐，最多也不会超过五个段，所以利用率是 45=80%45=80% ,如果不明白为啥不能超过 5 个段，请注意前提是连续的，这个时候不可能超过五段
   
4. 所有线程访问一个 4 字节的数据，那么此时的利用率是 432=12.5%432=12.5%
   
5. 最欢的情况，所有目标数据分散在内存的各个角落，那么需要 N 个内存段， 此时与使用一级缓存的作比较也是有优势的因为 N×128 还是要比 N×32 大不少，这里假设 N 不会因为 128128 还是 3232 而变的，而实际情况，当使用大粒度的缓存行的时候， N 有可能会减小