GPU与cuda

GPU

1. GPU的历史

  NVIDIA公司在1999年8月31日发布GeForce 256图形处理芯片时首先提出GPU的概念。
  GPU之所以被称为图形处理器，最主要的原因是因为它可以进行几乎全部与计算机图形有关的数据运算，而这些在过去是CPU的专利。
  目前，计算机图形学正处于前所未有的发展时期。近年来，GPU技术以令人惊异的速度在发展。渲染速率每6个月就翻一番。性能自99年，多年来翻番了十倍百倍，也就是（2的10次方比2）提高了上千倍！与此同时，不仅性能得到了提高，计算质量和图形编程的灵活性也逐渐得以改善。

1.1 NVidia GPU架构发展史

2008 - Tesla
  Tesla最初是给计算处理单元使用的，应用于早期的CUDA系列显卡芯片中，并不是真正意义上的普通图形处理芯片。

2010 - Fermi
  Fermi是第一个完整的GPU计算架构。首款可支持与共享存储结合纯cache层次的GPU架构，支持ECC的GPU架构。

2012 - Kepler
  Kepler相较于Fermi更快，效率更高，性能更好。

2014 - Maxwell
  其全新的立体像素全局光照 (VXGI) 技术首次让游戏 GPU 能够提供实时的动态全局光照效果。基于 Maxwell 架构的 GTX 980 和 970 GPU 采用了包括多帧采样抗锯齿 (MFAA)、动态超级分辨率 (DSR)、VR Direct 以及超节能设计在内的一系列新技术。

2016 - Pascal
  Pascal 架构将处理器和数据集成在同一个程序包内，以实现更高的计算效率。1080系列、1060系列基于Pascal架构

2017 - Volta
  Volta 配备640 个Tensor 核心，每秒可提供超过100 兆次浮点运算(TFLOPS) 的深度学习效能，比前一代的Pascal 架构快5 倍以上。

2018 - Turing
  Turing 架构配备了名为 RT Core 的专用光线追踪处理器，能够以高达每秒 10 Giga Rays 的速度对光线和声音在 3D 环境中的传播进行加速计算。Turing 架构将实时光线追踪运算加速至上一代 NVIDIA Pascal™ 架构的 25 倍，并能以高出 CPU 30 多倍的速度进行电影效果的最终帧渲染。2060系列、2080系列显卡也是跳过了Volta直接选择了Turing架构。

2. GPU的结构

  GPU全称是Graphic Processing Unit－－图形处理器，其最大的作用就是进行各种绘制计算机图形所需的运算，包括顶点设置、光影、像素操作等。GPU实际上是一组图形函数的集合，而这些函数由硬件实现。以前，这些工作都是有CPU配合特定软件进行的，GPU从某种意义上讲就是为了在图形处理过程中充当主角而出现的。
  GPU由多个streaming-multiprocessors (SMs)组成，它们通过crossbar内部互联网络共享L2 Cache和DRAM控制器。一个SM包含多个scalar processor cores (SPs) 和两种其他类型的功能单元（the Double-Precision Units (DPUs) for double-precision (DP) floating-point calculations and the Special-Function Units (SFUs) for processing transcendental functions and texture-fetching interpolations），还包含register files (RFs), load- store units (LSUs), scratchpad memory (i.e., shared memory), and various caches (i.e., instruction cache, constant cache, texture/read-only cache, L1 cache) for on-chip data caching。

2.1 功能单元

计算

• Scalar-Processor (SP): SM的主要基础处理器，实现基本的整数、浮点数计算，比较和类型转换等操作。每个SP包含一个单精度浮点数处理单元FPU和一个整数算数/逻辑单元ALU。都是流水线的。
  
• Special-Function-Unit (SFU): SFU实现了快速透明的函数计算（sin，cos等）和平面属性插值。
• Double-Precision-Unit (DPU): 专门的双精度计算单元。
• Load-Store-Unit (LSU): 从内存读写数据的单元。包含独立的计算单元来快速计算内存请求的来源和目的地址。

内存

• Register Files (RF): GPU的寄存器很大，出于吞吐量的考虑，划分成banks。因此相比于CPU寄存器，GPU寄存器的latency更大，而且有更多潜在的bank conflicts。
• Local Memory (LM): local memory不是物理空间，而是global memory的一部分。它是线程私有的，主要用来临时的spilling，比如寄存器溢出，或者数组在Kernel里声明了，但编译器无法获得他准确的indexing。Local memory在Fremi和Kepler中可以被L1和L2 cached，但是在Maxwell和Pascal中只能被L2 cached。寄存器溢出到local memory会造成巨大的性能下降（引起了更多的指令和内存拥挤）尤其是cache miss时。
• Shared Memory (HM，SMEM): 又称为scratchpad memory，是片上的存储，SM里的所有单元共享。它可以用来作为同一个thread blocks里不同线程之间快速数据交换的通讯接口。由于是片上存储，带宽很大，访问延迟很小。因此，把global/local memory的访问迁移到SMEM上的优化是被编程手册推荐的。为了获得更高的带宽，SMEM被划分成banks，这样就可以并行的访问（寄存器文件和L2 cache也类似）。但是如果在一次内存请求中，访问的两个地址落在了同一个bank，就造成了bank conflict，此时请求需要串行化，从而降低了SMEM的性能。
• Global Memory (GM): 又称device memory，GPU片下内存，GPU主存。它是GPU性能的主要瓶颈。GM可达到的吞吐量取决于：(1)GM物理带宽限制，即pin数、线长度和DRAM的物理特性。因此在kelper以前增长缓慢，但是Pascal使用了3D栈内存技术，性能获得了巨大提升。(2) 访存请求合并。LSU一开始会计算每个warp的目标地址。在memory fetch之前，有一个特殊的地址合并硬件来检查同一个warp里地址是否是连续分布的。
• Constant Memory (CM) / Constant Cache: 常量内存用来存储在Kernel执行期间没有改变的数据。在所有GPU上都是64KB的off-chip。和local memory类似，都是Global memory的一部分。不被L1或L2 cached，有专门的constant cache。每个SM上8/10KB的常量cache被专门设计，以便 the data of a single memory address can be broadcast to all
  threads across the warp at a time。
• Texture Memory ™ / Texture Cache: 又称surface memory也在全局内存上。被texture cache cached。Texture cache针对2D的空间局部性进行了优化。

缓存

• L1 Data-Cache：在Fermi上首次被提出来。SM的私有L1 Cache和SMEM共享片上存储，他们的大小是可以配置的，Fermi上16/48 or 48/16，Kepler上32/32 or 48/16。L1 cache line是128B。L1 data cache缓存global memory读和local memory的读写，并且是non-corherent。local memory主要用来寄存器溢出、函数调用和自动变量。当L1 cache 被用来缓存global memory时是只读的，当被用来缓存local memory时也是可写的。从Maxwell开始，传统的L1 cache被统一到了Texture cache里。
• L2 Cache： 它被用来缓存各种类型的memory access，且和宿主机的CPU memory保持一致性（？）。采取写回策略。

NC and ROP

• Interconnection Network (NC): crossbar network.它允许多个SM和L2 banks之间同时通信。crossbar NC包含一条地址线和两条数据线。地址线是单向的，从SM到L2 banks，而两条数据线则是双向的，存储SM和L2 banks。因此，通讯是点到点的。一个内存请求队列（memory-request queue, MRQ）和一个bank 加载队列（ bank load queue, BLQ）分别对应了一个SM和L2 bank。当load请求从SM中的LSU产生时，它会首先缓存在MRQ里，然后通过NC被分发到目的BLQ中。在BLQ中等待一会，然后这个请求将被L2 banks处理。crossbar network在同时多个连接时有很高的切换代价，特别是当访问请求随机且数据量大。
• Raster Operation Processor (ROP)。

2.2 几种架构特性简介

• Fermi第一次出现了两级cache的设计。
• Kelper有巨大的计算能力，因为它们每个SM有最多的CUDA cores。
• 从Maxwell开始注重功耗。
• Pascal则因为其3D的memory架构拥有了快速的half-percision(16 bits)的计算。
• GTX是桌面发型版，Telsa面向HPC，Jetson则是嵌入式设备。

2.3 具体的几种微观架构

NVIDIA Tesla

Tesla微观架构总览图如上。下面将阐述它的特性和概念：

• 拥有7组TPC（Texture/Processor Cluster，纹理处理簇）
• 每个TPC有两组SM（Stream Multiprocessor，流多处理器）
• 每个SM包含：
    6个SP（Streaming Processor，流处理器）
    2个SFU（Special Function Unit，特殊函数单元）
    L1缓存、MT Issue（多线程指令获取）、C-Cache（常量缓存）、共享内存
    除了TPC核心单元，还有与显存、CPU、系统内存交互的各种部件。

NVidia Fermi架构

• 拥有16个SM
• 每个SM：
    2个Warp（线程束）
    两组共32个Core
    16组加载存储单元（LD/ST）
    4个特殊函数单元（SFU）
• 每个Warp：
    16个Core
    Warp编排器（Warp Scheduler）
    分发单元（Dispatch Unit）
• 每个Core：
    1个FPU（浮点数单元）
    1个ALU（逻辑运算单元）

NVidia Maxwell
  采用了Maxwell的GM204，拥有4个GPC，每个GPC有4个SM，对比Tesla架构来说，在处理单元上有了很大的提升。

NVidia Kepler
  Kepler除了在硬件有了提升，有了更多处理单元之外，还将SM升级到了SMX。SMX是改进的架构，支持动态创建渲染线程（下图），以降低延迟。

3. GPU执行模型

  GPU的运算速度如此之快，主要得益于GPU是对图形实时渲染量身定制的，具有两点主要特征：超长流水线与并行计算。GPU中使用SIMT模型。

3.1 SIMT

• single-instruction-multiple-threads，SIMT，单指令多线程。

• thread blocks or Cooperative-Thread-Arrays (CTAs)

  一个Kernel函数包含很多同时运行的轻量级GPU线程，而这些线程被划分成多个thread blocks组。当Kernel开始运行时，CTAs被分配到不同的SM上去，也可能其中几个CTA被分配到相同的SM上去。

  CTA里面的threads则进一步分组管理，这些分组叫warp，warp内遵循lockstep规则，即所有线程同步执行。在SM里，warp是基本的调度、执行和读写cache/memory的单元。

  如果一个warp里的线程在某个点（if-else等）出现了分歧(warp divergence)，则所有的分支顺序执行。执行if的时候，else的线程被挂起。在分歧（if-else）结束以后，再继续lockstep的执行。warp divergence造成了巨大的overhead。

  GPU支持multi-issuing和multi-dispatching。在执行期间，dual- or quad- warp scheduler选择两个或四个ready状态的warps分调度到不同的功能单元上（SPs，SFUs，当然大多数指令都是在SP上完成的）。

3.2 线程级别的映射

线程以warp为单位被映射到SP、SFU或者DPU上。

• CTA被映射到SM上。
• grid映射到GPU device。

3.3 编译

  PTX stands for Parallel-Thread-Execution，是一个中间级的汇编代码。而cubin二进制则已经指定了架构。

  Shader-Assembly (SASS)，真正的机器汇编，由cubin文件经过cuobjdump工具转换而来。目前没有官方的sass to cubin的工具。

cuda

1. CUDA简介

  CUDA：Compute Unified Device Architecture，统一计算设备架构，CUDA是一种由NVIDIA推出的通用并行计算架构，该架构使GPU (Graphics processing unit) 能够解决复杂的计算问题。它包含了CUDA指令集架构（ISA）以及GPU内部的并行计算引擎。

  2006年，NVIDIA公司发布了CUDA，CUDA是建立在NVIDIA的CPUs上的一个通用并行计算平台和编程模型，基于CUDA编程可以利用GPUs的并行计算引擎来更加高效地解决比较复杂的计算难题。近年来，GPU最成功的一个应用就是深度学习领域，基于GPU的并行计算已经成为训练深度学习模型的标配。

  CUDA是NVIDIA公司所开发的GPU编程模型，它提供了GPU编程的简易接口，基于CUDA编程可以构建基于GPU计算的应用程序。CUDA提供了对其它编程语言的支持，如C/C++，Python，Fortran等语言，

1.1 GPU和CPU

  GPU并不是一个独立运行的计算平台，而需要与CPU协同工作，可以看成是CPU的协处理器，因此当我们在说GPU并行计算时，其实是指的基于CPU+GPU的异构计算架构。在异构计算架构中，GPU与CPU通过PCIe总线连接在一起来协同工作，CPU所在位置称为为主机端（host），而GPU所在位置称为设备端（device）。

  GPU包括更多的运算核心，其特别适合数据并行的计算密集型任务，如大型矩阵运算，而CPU的运算核心较少，但是其可以实现复杂的逻辑运算，因此其适合控制密集型任务。另外，CPU上的线程是重量级的，上下文切换开销大，但是GPU由于存在很多核心，其线程是轻量级的。因此，基于CPU+GPU的异构计算平台可以优势互补，CPU负责处理逻辑复杂的串行程序，而GPU重点处理数据密集型的并行计算程序，从而发挥最大功效。

与CPU相比，GPU的优点缺点如下
优点

• 更大的内存带宽
• 更多的执行单元，虽然频率比CPU低
• 价格更低

缺点

• 运算单元多，只适合高度并行化的工作
• 对于具有高度分支的程序，效率会比较差
  

  GPU特别适合大量并行的数据运算（高运算密度）。由于对每个数据进行相同的操作，所以对复杂的流控制需求较低，并且因为处理许多数据单元并且具有高运算密度，内存读取延时可以通过运算来隐藏，CPU采用的是高速缓存cache来缩减延时 （latency ）

1.2 可伸缩Scaleable的编程模式

核心是三个关键的抽象：

• 线程组的层次结构
• 共享内存
• 障碍同步
  

2. CUDA编程模型

2.1 CUDA编程模型中的概念

  CUDA编程模型是一个异构模型，需要CPU和GPU协同工作。在CUDA中，host和device是两个重要的概念，我们用host指代CPU及其内存，而用device指代GPU及其内存。
  CUDA程序中既包含host程序，又包含device程序，它们分别在CPU和GPU上运行。同时，host与device之间可以进行通信，这样它们之间可以进行数据拷贝。典型的CUDA程序的执行流程如下：

1. 分配host内存，并进行数据初始化；
2. 分配device内存，并从host将数据拷贝到device上；
3. 调用CUDA的核函数在device上完成指定的运算；
4. 将device上的运算结果拷贝到host上；
5. 释放device和host上分配的内存。

2.1.1 Kernel

  kernel是CUDA中一个重要的概念，kernel是在device上线程中并行执行的函数，核函数用__global__符号声明，在调用时需要用<<<grid, block>>>来指定kernel要执行的线程数量。

  在CUDA中，每一个线程都要执行核函数，并且每个线程会分配一个唯一的线程号thread ID，这个ID值可以通过核函数的内置变量threadIdx来获得。

  由于GPU实际上是异构模型，所以需要区分host和device上的代码，在CUDA中是通过函数类型限定词开区别host和device上的函数，主要的三个函数类型限定词如下：

• global：在device上执行，从host中调用（一些特定的GPU也可以从device上调用），返回类型必须是void，不支持可变参数参数，不能成为类成员函数。注意用__global__定义的kernel是异步的，这意味着host不会等待kernel执行完就执行下一步。

• device：在device上执行，单仅可以从device中调用，不可以和__global__同时用。

• host：在host上执行，仅可以从host上调用，一般省略不写，不可以和__global__同时用，但可和__device__，此时函数会在device和host都编译。

// 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); 
    ... 
}

2.1.2 Thread Hierarchy

  要深刻理解kernel，必须要对kernel的线程层次结构有一个清晰的认识。首先GPU上很多并行化的轻量级线程。kernel在device上执行时实际上是启动很多线程，一个kernel所启动的所有线程称为一个网格（grid），同一个网格上的线程共享相同的全局内存空间，grid是线程结构的第一层次。
而网格又可以分为很多线程块（block），一个线程块里面包含很多线程，这是第二个层次。

  线程两层组织结构如下图所示，这是一个gird和block均为2-dim的线程组织。grid和block都是定义为dim3类型的变量，dim3可以看成是包含三个无符号整数（x，y，z）成员的结构体变量，在定义时，缺省值初始化为1。因此grid和block可以灵活地定义为1-dim，2-dim以及3-dim结构，对于图中结构（主要水平方向为x轴），定义的grid和block如下所示，kernel在调用时也必须通过执行配置<<<grid, block>>>来指定kernel所使用的线程数及结构。

dim3 grid(3, 2);
dim3 block(5, 3);
kernel_fun<<< grid, block >>>(prams...);

2.1.3 索引和ID

  所以，一个线程需要两个内置的坐标变量（blockIdx，threadIdx）来唯一标识，它们都是dim3类型变量，其中blockIdx指明线程所在grid中的位置，而threaIdx指明线程所在block中的位置。

  一个线程块上的线程是放在同一个流式多处理器（SM)上的，但是单个SM的资源有限，这导致线程块中的线程数是有限制的，现代GPUs的线程块可支持的线程数可达1024个。有时候，我们要知道一个线程在blcok中的全局ID，此时就必须还要知道block的组织结构，这是通过线程的内置变量blockDim来获得。它获取线程块各个维度的大小。

• 对于一个1-dim的block，一维索引直接对应线程ID
• 对于一个2-dim的block(Dx, Dy)的块， 线程索引为 (x, y)，对应ID为(x + y Dx);
• 三维(Dx, Dy, Dz)块, 线程索引为(x, y, z) ，对应ID为 (x + y Dx+ z Dx Dy).
  另外线程还有内置变量gridDim，用于获得网格块各个维度的大小。

2.1.4 同步和原子

同步 synchronization
  线程通信和资源共享：
  通过同步来相互协调，调用固有函数__syncthreads()
  只在block层次产生作用，类似于C/C++中的 barrier()函数

原子操作 Atomic Operation
  执行读-修改-写的原子操作，在全局或共享内存空间中串行操作

2.2 CUDA的内存模型

  每个线程有自己的私有本地内存（Local Memory），而每个线程块有包含共享内存（Shared Memory）,可以被线程块中所有线程共享，其生命周期与线程块一致。此外，所有的线程都可以访问全局内存（Global Memory）。还可以访问一些只读内存块：常量内存（Constant Memory）和纹理内存（Texture Memory）。内存结构涉及到程序优化。

2.3 GPU硬件的核心组件SM

Streaming Multiprocessors
  A GPU is built around a scalable array of multithreaded Streaming Multiprocessors (SMs).

• Streaming Multiprocessors (SMs)
    SM的核心组件包括CUDA核心，共享内存，寄存器等，SM可以并发地执行数百个线程，并发能力就取决于SM所拥有的资源数。当一个kernel被执行时，它的gird中的线程块被分配到SM上，一个线程块只能在一个SM上被调度。SM一般可以调度多个线程块，这要看SM本身的能力。那么有可能一个kernel的各个线程块被分配多个SM，所以grid只是逻辑层，而SM才是执行的物理层。

  SM采用的是SIMT (Single-Instruction, Multiple-Thread，单指令多线程)架构，，一个线程内部通过指令流水进行指令级别的并行，通过硬件多线程进行线程级别的并行。

  基本的执行单元是线程束（wraps)，线程束包含32个线程。一个GPU中有多个SM,每个SM有多个core（processor），但是只有一个指令单元，同时只能够执行完全相同的指令集。

  这些线程同时执行相同的指令，但是每个线程都包含自己的指令地址计数器和寄存器状态，也有自己独立的执行路径。所以尽管线程束中的线程同时从同一程序地址执行，但是可能具有不同的行为，比如遇到了分支结构，一些线程可能进入这个分支，但是另外一些有可能不执行，它们只能死等，因为GPU规定线程束中所有线程在同一周期执行相同的指令，线程束分化会导致性能下降。当线程块被划分到某个SM上时，它将进一步划分为多个线程束，因为这才是SM的基本执行单元，但是一个SM同时并发的线程束数是有限的。这是因为资源限制，SM要为每个线程块分配共享内存，而也要为每个线程束中的线程分配独立的寄存器。所以SM的配置会影响其所支持的线程块和线程束并发数量。

  总之，就是网格和线程块只是逻辑划分，一个kernel的所有线程其实在物理层是不一定同时并发的。所以kernel的grid和block的配置不同，性能会出现差异，这点是要特别注意的。还有，由于SM的基本执行单元是包含32个线程的线程束，所以block大小一般要设置为32的倍数。

• Warp
    Multiprocessor执行线程块时，把他们以32个并行线程为一组，称为warps,每个warp由warp scheduler调度执行，warp中的每个线程在相同的程序地址处开始执行，但他们有自己的指令地址计数器和寄存器状态。

  当他们执行相同的指令时可以达到最大的效率，但当因为由于依靠数据决定的条件分支产生了分歧（warp divergence），warp连续执行每个分支路径，禁止不相关的线程，执行完成后，又回归到相同的执行路径。不同warps执行互相独立。

  多线程的CUDA程序具有自动适应性，它被分解为相互独立的线程块，可以以任意的顺序执行，不管是并行还是串行，所以一个编译好的CUDA程序可以在任意数目的多处理器上运行，只有运行的系统需要知道实际的处理器数量。

3. CUDA编程实例

3.1 向量加法

// 两个向量加法kernel，grid和block均为一维
__global__ void add(float* x, float * y, float* z, int n)
{
    // 获取全局索引
    int index = threadIdx.x + blockIdx.x * blockDim.x;
    // 步长
    int stride = blockDim.x * gridDim.x;
    for (int i = index; i < n; i += stride)
    {
        z[i] = x[i] + y[i];
    }
}

int main()
{
    int N = 1 << 20;
    int nBytes = N * sizeof(float);
    // 申请host内存
    float *x, *y, *z;
    x = (float*)malloc(nBytes);
    y = (float*)malloc(nBytes);
    z = (float*)malloc(nBytes);

    // 初始化数据
    for (int i = 0; i < N; ++i)
    {
        x[i] = 10.0;
        y[i] = 20.0;
    }

    // 申请device内存
    float *d_x, *d_y, *d_z;
    cudaMalloc((void**)&d_x, nBytes);
    cudaMalloc((void**)&d_y, nBytes);
    cudaMalloc((void**)&d_z, nBytes);

    // 将host数据拷贝到device
    cudaMemcpy((void*)d_x, (void*)x, nBytes, cudaMemcpyHostToDevice);
    cudaMemcpy((void*)d_y, (void*)y, nBytes, cudaMemcpyHostToDevice);
    // 定义kernel的执行配置
    dim3 blockSize(256);
    dim3 gridSize((N + blockSize.x - 1) / blockSize.x);
    // 执行kernel
    add << < gridSize, blockSize >> >(d_x, d_y, d_z, N);

    // 将device得到的结果拷贝到host
    cudaMemcpy((void*)z, (void*)d_z, nBytes, cudaMemcpyHostToDevice);

    // 检查执行结果
    float maxError = 0.0;
    for (int i = 0; i < N; i++)
        maxError = fmax(maxError, fabs(z[i] - 30.0));
    std::cout << "最大误差: " << maxError << std::endl;

    // 释放device内存
    cudaFree(d_x);
    cudaFree(d_y);
    cudaFree(d_z);
    // 释放host内存
    free(x);
    free(y);
    free(z);

    return 0;
}

在上面的实现中，我们需要单独在host和device上进行内存分配，并且要进行数据拷贝。CUDA 6.0引入统一内存（Unified Memory）来避免这种麻烦，简单来说就是统一内存使用一个托管内存来共同管理host和device中的内存，并且自动在host和device中进行数据传输。CUDA中使用cudaMallocManaged函数分配托管内存：

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

// 使用统一内存管理
int main()
{
    int N = 1 << 20;
    int nBytes = N * sizeof(float);

    // 申请托管内存
    float *x, *y, *z;
    cudaMallocManaged((void**)&x, nBytes);
    cudaMallocManaged((void**)&y, nBytes);
    cudaMallocManaged((void**)&z, nBytes);

    // 初始化数据
    for (int i = 0; i < N; ++i)
    {
        x[i] = 10.0;
        y[i] = 20.0;
    }

    // 定义kernel的执行配置
    dim3 blockSize(256);
    dim3 gridSize((N + blockSize.x - 1) / blockSize.x);
    // 执行kernel
    add << < gridSize, blockSize >> >(x, y, z, N);

    // 同步device 保证结果能正确访问
    cudaDeviceSynchronize();
    // 检查执行结果
    float maxError = 0.0;
    for (int i = 0; i < N; i++)
        maxError = fmax(maxError, fabs(z[i] - 30.0));
    std::cout << "最大误差: " << maxError << std::endl;

    // 释放内存
    cudaFree(x);
    cudaFree(y);
    cudaFree(z);

    return 0;
}

3.2 矩阵乘法

定义矩阵的结构体：

// 矩阵类型，行优先，M(row, col) = *(M.elements + row * M.width + col)
struct Matrix
{
    int width;
    int height;
    float *elements;
};

定义了两个辅助的__device__函数分别用于获取矩阵的元素值和为矩阵元素赋值

// 获取矩阵A的(row, col)元素
__device__ float getElement(Matrix *A, int row, int col)
{
	return A->elements[row * A->width + col];
}

// 为矩阵A的(row, col)元素赋值
__device__ void setElement(Matrix *A, int row, int col, float value)
{
	A->elements[row * A->width + col] = value;
}

// 矩阵相乘kernel，2-D，每个线程计算一个元素
__global__ void matMulKernel(Matrix *A, Matrix *B, Matrix *C)
{
	float Cvalue = 0.0;
	int row = threadIdx.y + blockIdx.y * blockDim.y;
	int col = threadIdx.x + blockIdx.x * blockDim.x;
	for (int i = 0; i < A->width; ++i)
	{
		Cvalue += getElement(A, row, i) * getElement(B, i, col);
	}
	setElement(C, row, col, Cvalue);
}

采用统一内存编写矩阵相乘的测试实例

int main()
{
    int width = 1 << 10;
    int height = 1 << 10;
    Matrix *A, *B, *C;
    // 申请托管内存
    cudaMallocManaged((void**)&A, sizeof(Matrix));
    cudaMallocManaged((void**)&B, sizeof(Matrix));
    cudaMallocManaged((void**)&C, sizeof(Matrix));
    int nBytes = width * height * sizeof(float);
    cudaMallocManaged((void**)&A->elements, nBytes);
    cudaMallocManaged((void**)&B->elements, nBytes);
    cudaMallocManaged((void**)&C->elements, nBytes);

    // 初始化数据
    A->height = height;
    A->width = width;
    B->height = height;
    B->width = width;
    C->height = height;
    C->width = width;
    for (int i = 0; i < width * height; ++i)
    {
        A->elements[i] = 1.0;
        B->elements[i] = 2.0;
    }

    // 定义kernel的执行配置
    dim3 blockSize(32, 32);
    dim3 gridSize((width + blockSize.x - 1) / blockSize.x, 
        (height + blockSize.y - 1) / blockSize.y);
    // 执行kernel
    matMulKernel << < gridSize, blockSize >> >(A, B, C);


    // 同步device 保证结果能正确访问
    cudaDeviceSynchronize();
    // 检查执行结果
    float maxError = 0.0;
    for (int i = 0; i < width * height; ++i)
        maxError = fmax(maxError, fabs(C->elements[i] - 2 * width));
    std::cout << "最大误差: " << maxError << std::endl;

    return 0;
}