记录自己阅读《Professional CUDA C Programming》这本书学习CUDA编程的一些知识，同时供大家参考。

主要参考文献：

①谭升大佬的博客应该查询过CUDA编程的同学都应该有所了解，该博客将《Professional CUDA C Programming》这本书中的知识点进行了浓缩。

《Professional CUDA C Programming》相关资料整理汇总：

本教材的所有学习笔记都是按照书本中的章节进行划分的，所以文章中的所有源码在我整理的仓库对应章节都能够轻松找到，欢迎大家收藏，fork：

彩色电子版《Professional CUDA C Programming》：

文章和教材中各章节的代码样例：

教材中各章节的习题答案：

正文：

二、CUDA编程模型

1、CUDA编程模型概述

CUDA编程模型提供了一个计算机架构抽象来作为应用程序和可用硬件之间的桥梁。

同时CUDA编程模型利用GPU架构的计算能力提供了两种层次结构：

• 访问内存的层次结构：Local、Share、Global、Const、Texture
• 组织线程的层次结构：Grid、Block、Thread

（1）CUDA编程结构

①host与device

一个异构环境由多个CPU和GPU组成，每个CPU和GPU的内存由一条PCI-Express总线分隔开。正如之前提到过的在CUDA中将CPU视为host端，GPU视为device端。在后续的学习中也将这两清晰的区分开：

• host：CPU及其host内存，相关变量用h_前缀作为标识
• device：GPU及其device内存，相关变量用d_前缀作为标识

从CUDA 6.0开始，NVIDIA提供了名为“统一寻址”的改进，可以使用单一指针访问CPU和GPU的内存，但是为了更加清晰的了解和区分开host和device端分割的异构思想，目前还是先学习下如何分别为host和device端分配内存以及如何在CPU和GPU之间拷贝数据，通过手动管理内存可以加深对于底层硬件调度的理解。

②kernel核函数

又一个新的关键词，kernel核函数是真正在GPU上运行的代码，目前就简单的理解为在GPU上会有非常多的线程会同时一个kernel函数，同时使用之前学过的数据划分的方式来让每个kernel函数计算不同的数据。

抽象一次kernel函数运行过程：

从图中需要关注3个点：

• 相同计算：所有的thread执行的是同一个计算逻辑，也就是kernel函数
• 不同数据：每个thread作用的数据块是不相同的，这也就是之前学过的数据划分
• 并行计算：所有的thread是并行计算的，不需要排队，可以并行获取数据和计算

而CUDA编程的核心其实也就是如何合理的划分数据并且针对数据结构编写高效的kernel函数

③CUDA程序流程

一个典型的CUDA程序实现流程遵循以下模式：

• 计算内存申请：在CPU和GPU上申请用于计算的数据内存
• CPU->GPU搬移：把数据从CPU内存拷贝到GPU内存
• GPUkernel计算：调用kernel函数对存储在GPU内存中的数据进行计算
• GPU->CPU回迁：将计算结果从GPU内存拷贝回CPU内存
• 计算内存释放：将CPU和GPU上申请的内存进行释放

需要重点关注的是kernel一旦启动，管理权立刻返回给host，CUDA编程模型是异步的。

下面就依据CUDA程序流程来逐一讲解相应的内存管理机制和函数。

（2）内存管理

①内存管理函数

host端的代码实际上就是标准的C/C++语言，直接使用标准C函数即可，相应的CUDA也提供了在GPU上申请内存的函数，并且他们为了使用和记忆的方便，将函数名统一为了cuda + c中标准函数名：

重点关注cudaMalloc和cudaMemcpy两个函数：

• cudaMalloc：执行GPU内存分配
  此函数负责向GPU分配一段连续的内存，并以devPtr的形式返回指向所分配内存的指针。
  cudaError_t cudaMalloc ( void** devPtr, size_t size );
• cudaMemcpy：执行CPU和GPU以及他们之间的内存拷贝
  此函数主要负责host与device之间的数据拷贝，就函数参数来说是将src地址开始的count长度的内存拷贝到dst地址上。
  cudaError_t cudaMemcpy ( void* dst, const void* src, size_t count, cudaMemcpyKind kind );
  cudaMemcpyKind有四种，并且从名字就能看出它们的作用：
  
• cudaMemcpyHostToHost：从host拷贝到host
• cudaMemcpyHostToDevice：从host拷贝到device
• cudaMemcpyDeviceToHost：从device拷贝到host
• cudaMemcpyDeviceToDevice：从device拷贝到device

关于这个函数有一个非常需要注意的点：此函数以同步方式执行，因为在cudaMemcpy函数返回以及传输操作完成之前host端是阻塞的，必须等到拷贝完毕才会继续执行后续计算。

同时注意除了内核启动之外的CUDA调用都会返回一个错误枚举类型cuda_Error_t，如果GPU内存分配成功则返回cudaSuccess否则返回cudaErrorMemoryAllocation。可以使用以下CUDA运行时函数将错误代码转化为可读的错误消息：

char* cudaGetErrorString(cudaError_t error)

②内存层次结构

CUDA编程模型从GPU架构中抽象出了一个内存层次结构，第4、5章会详细介绍，这里先讲解一个由全局内存（global memory）和共享内存（shared memory）组成的简化版本。

实际上最常接触到的也就是global memory和shared memory，并且CUDA中是有global和shared这两个关键字的，所以之后为了方便理解代码，在文字中也都以英文来进行记录。

③内存管理样例

下面以一个经典的数组相加的例子来学习如何在host和device端进行数据传输，以及如何使用CUDA C编写kernel函数。数组a的第i个元素与数组b的第i个元素相加得到数组c的第i个元素，线性和并行的区别如图：

完整可以运行代码请移步代码仓库：

• 内存申请：

int nElem = 1024;
size_t nBytes = nElem * sizeof(float); // 设置数组的长度以及需要的内存大小
​
float *h_A, *h_B, *h_C;
h_A = (float *)malloc(nByte);          // 在host上分配A、B、C三个数组的内存并获取各自的地址
h_B = (float *)malloc(nByte);
host_res = (float *)malloc(nByte);     // 用于存放host求和的结果
device_res = (float *)malloc(nByte);   // 用于存放device求和的结果
​
initialData(h_A, nElem);               // 初始化h_A和h_B的数组的数值
initialData(h_B, nElem);
​
float *d_A, *d_B, *d_C;
cudaMalloc((float**)&d_A, nBytes);     // 在device上分配A、B、C三个数组的内存并获取各自的地址
cudaMalloc((float**)&d_B, nBytes);
cudaMalloc((float**)&d_C, nBytes);

• 内存拷贝：

将数据从host内存拷贝到device的global memory中，cudaMemcpyHostToDevice指明了拷贝方向

cudaMemcpy(d_A, h_A, nBytes, cudaMemcpyHostToDevice); // 将host端数组A和B的值拷贝到device端
cudaMemcpy(d_B, h_B, nBytes, cudaMemcpyHostToDevice);

• kernel计算：

这里先给出kernel函数的定义，后续介绍如何启动kernel函数

__global__ void sumArraysOnGPU(float *d_A, float *d_B, float *d_C) {
  int i = threadIdx.x;      // 这里假设每个线程处理一个数据
  d_C[i] = d_A[i] + d_B[i]; // A和B的第i个数据进行求和并存放到C的第i个位置
}

• 内存回拷贝：

将device的计算结果拷贝回host端的d_res结果数组中

cudaMemcpy(device_res, d_C, nBytes, cudaMemcpyDeviceToHost);

这其中包含了两个重要的点：

强制同步：之前提到host启动了kernel之后控制前就会立刻回到host这里，host和device是异步的，而这里的cudaMemcpy函数是一个同步函数，因此host端需要等到device中的kernel计算完毕并且将结果从d_C成功拷贝回host端的device_res之后才会进行后续的计算。

禁止混用：切忌不要不要不要混用host和device端的内存，因为host和device内存之间是物理隔离的，它两之间隔了一坐PCIe大桥。在host端依据device端变量的地址去访问数据是错误并且危险的，反之亦然，因此当出现混用的时候运行是会报错的。当然CUDA 6.0有了统一寻址机制，但我个人觉得将host和device显著的分隔开更能够体现异构架构。

device_res = d_C; // 通过这种方式直接将一个device端的地址赋给host端的变量是完全错误的

• 内存释放

计算结束，我们将动态申请的内存全部释放掉，防止内存泄漏：

free(h_A), free(h_B), free(h_C);
cudaFree(d_A), cudaFree(d_B), cudaFree(d_C);

（3）线程管理

kernel函数在host端启动后，它的执行会移动到device端，此时device中会产生大量的线程thread并且每个thread都执行由kernel函数指定的语句。

①两层线程层次结构

CUDA同样抽象出了线程层次结构来便于组织线程，这是一个两层的线程层次结构，由线程网络grid和线程块block构成：

从图中我们可以看出，一个kernel函数对应了device内存上的一个grid，而一个grid中有很多个block，每个block中又有很多的thread，最终每个thread来执行kernel函数中的语句。下面给出书中的定义：

• grid：一个kernel所启动的所有thread统称为一个grid，同一个grid中的所有线程共享相同的global memory
• block：一个grid由多个block构成，一个block包含一组thread，同一个block内的thread可以进行如下协作：
  
• 同步
• 共享内存：同一个block中的所有thread，共享相同的shared memory

不同block中的线程不能协作。

②关键变量：索引Idx

线程具有唯一的一组标识符变量：

• blockIdx：thread所在block在grid中的索引block Index，包括3个无符号整数
  
• blockIdx.x：block在grid中的x坐标
• blockIdx.y：block在grid中的y坐标
• blockIdx.z：block在grid中的z坐标

• threadIdx：thread在blcok中的索引thread Index，包括3个无符号整数
  
• threadIdx.x：thread在block中的x坐标
• threadIdx.y：thread在block中的y坐标
• threadIdx.z：thread在block中的z坐标

这些变量是kernel函数需要预初始化的内置变量，是基于unit3定义的CUDA内置向量类型。当执行一个kernel函数时，CUDA会为每个thread分配坐标变量blockIdx和threadIdx，基于这两个变量我们就可以确定唯一的thread。

③关键变量：维度Dim

能够自由获取thread的blockIdx和threadIdx坐标，自然也需要对grid和block的布局维度进行规划和设置，以①中布局为例：

• gridDim：grid的维度，设置一个grid中block的布局，①中设置的布局为：
  dim3 grid(3, 2); // 设置x方向3，y方向2的block结构，注意x代表的是横向向右的计数，y代表纵向向下的计数
  
• gridDim.x：grid中x方向block的个数，①中为3
• gridDim.y：grid中y方向block的个数，①中为2
• gridDim.z：grid中z方向block的个数，①中为1

• blockDim：block的维度，设置一个block中thread的布局，①中设置的布局为：
  dim3 block(5, 3); // 设置x方向5，y方向3的thread结构
  
• blockDim.x：block中x方向thread的个数，①中为5
• blockDim.y：block中y方向thread的个数，①中为3
• blockDim.z：block中z方向thread的个数，①中为1

它们都是dim3类型的变量，是基于unit3定义的整形向量，用于表示xyz三个维度的大小，未指定的维度初始化为1。

④打印索引和维度信息

注意自己通过dim3定义的block和grid实际上是一个可以随意更名的变量，想要叫什么名字都可以，而通过一下语句设置了block和grid的布局后CUDA为每个thread分配的坐标变量blockIdx和threadIdx是不可以改名的：

dim3 grid_def(3, 2);   // 自定义的dim3变量，想叫啥叫啥，可以把这里的grid_def改成任意的xxx，这里的block_def同理
dim3 block_def(5, 3); 
<<<grid_def, block_def>>>; // 启动kernel函数，设定了kernel函数中grid和block的结构

这里只是为了体会用于设置grid和block而自定义的dim3变量grid_def和block_def与真实grid和block的差异，防止理解出现混淆，但实际上理解了之后大多数情况都是直接用grid和block来命名这两个设置维度的变量。

可以用下面的代码分别打印自定义grid_def和block_def以及启动kernel函数预初始化的grid和block变量：

printf("grid_def.x %d grid_def.y %d grid_def.z %d\n", grid_def.x, grid_def.y, grid_def.z);
printf("block_def.x %d block_def.y %d block_def.z %d\n", block_def.x, block_def.y, block_def.z);
__global__ void printIndex(void)
{
  printf("threadIdx:(%d,%d,%d) blockIdx:(%d,%d,%d) blockDim:(%d,%d,%d) gridDim(%d,%d,%d)\n",
         threadIdx.x, threadIdx.y, threadIdx.z, blockIdx.x, blockIdx.y, blockIdx.z, 
         blockDim.x, blockDim.y, blockDim.z, gridDim.x, gridDim.y, gridDim.z);
}

打印结果为：

可以看到每个thread有自己的唯一坐标，所有的thread都有相同的gridDim和blockDim。

⑤确定grid和block布局

对于一个给定的数据大小，确定grid和block布局的一般步骤为：

• 确定block布局
• 根据数据大小和block布局计算出grid布局

而想要确定block布局，通常需要考虑：

• 内核的性能特性
• GPU资源的限制

为了展示数据大小、block布局和grid布局之间的关系，打印如下数据：

 // 定义数据大小
 int nElem = 1024;
​
 // 定义grid和block布局
 dim3 block (1024);
 dim3 grid ((nElem+block.x-1)/block.x);
 printf("grid.x %d block.x %d \n",grid.x, block.x);
​
 // 重新设置block布局
 block.x = 512;
 grid.x = (nElem+block.x-1)/block.x;
 printf("grid.x %d block.x %d \n",grid.x, block.x);
​
 // 重新设置block布局
 block.x = 256;
 grid.x = (nElem+block.x-1)/block.x;
 printf("grid.x %d block.x %d \n",grid.x, block.x);
​
 // 重新设置block布局
 block.x = 128;
 grid.x = (nElem+block.x-1)/block.x;
 printf("grid.x %d block.x %d \n",grid.x, block.x);

打印结果为：

由于数据大小是固定的，所以block布局变化时grid布局也相应发生改变，实际上设置合适的block布局对于计算性能影响很大，但是这里先不进行介绍，知道有这个概念即可。

（4）启动一个CUDA核函数

①调用语句

首先大家肯定对于C/C++中的函数调用语句不陌生：

function_name (argument list);

那么其实CUDA C也是采用了一个非常类型的调用语句，只不过需要通过<<< >>>指定配置的grid和block布局：

kernel_name <<<grid, block>>>(argument list); // 依据自定义的dim3类型grid和block设置启动的grid和block布局

• kernel_name其实就是kernel函数的名字，和function_name没有任何差别
• <<<grid, block>>>就是用于指定启动的grid和block布局
• argument list是参数列表，也没有任何区别

所以可以发现还是非常容易记忆的，只需要加一个<<<grid, block>>>就可以了

②不同的grid和block划分

通过设置<<<grid, block>>>我们可以设置：

• kernel中thread的数目
• kernel中使用的thread布局

同一个block中的线程之间可以相互协作，不同block中的线程不能协作，对于同一个问题可以使用不同的grid和block布局来组织thread，并且不同的布局会显著影响计算的性能。

这里先假设有32个数据用于计算，将每8个元素分为一个block，需要启动4个block：

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

如果把32个元素都放到一个block，那么只会得到一个block：

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

如果每个block只有一个元素，那么会有32个block：

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

这里再次强调kernel函数的调用与host线程是异步的，kernel函数调用后，控制权立刻返回给host端，但是可以使用以下函数强制host端程序等待所有kernel函数执行结束：

cudaError_t cudaDeviceSynchronize(void);

（5）编写核函数

学会了如何调用和启动kernel函数，我们可以正式的学一下如何编写kernel函数：

①核函数-kernel function

kernel函数是在device端执行的函数，在其中需要为一个thread规定要进行的计算以及要处理的数据。当kernel被调用时，许多不同的CUDAthread并行执行同一个计算任务，但是处理的数据通常是不同的。以下使用__global__声明定义kernel函数：

__global__ void kernel_name(argument list);

kernel函数的返回类型必须为void！！！

②函数类型限定符

CUDA C中提供了集中函数类型限定符，能够指定一个函数在host执行还是在device执行，以及可以被host还是device调用

__device__和__host__可以一起使用，表示函数可以同时在device和host端进行编译。。

③kernel函数固定限制

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

④kernel函数举例

同样是两个长度为N的向量A和B相加得到C的例子：

void sumArraysOnHost(float *A, float *B, float *C, const int N) { // host端函数
    for (int i = 0; i < N; i++) {
        C[i] = A[i] + B[i];
    }
}
​
__global__ void sumArraysOnGPU(float *A, float *B, float *C) {    // device端函数
    int i = threadIdx.x;
    C[i] = A[i] + B[i];
}

可以看到通过kernel核函数中的并行计算，我们去除了host端函数中的for循环，不严谨的假设计算第i个位置的和需要的时间为t，那么在host端函数中需要用n*t时间，而在device端函数中只需要用到t时间，因为kernel函数是并行执行的，再用第一章绘制的图解释一下区别：

（6）验证核函数

kernel函数编写好了之后我们需要将host端的计算结果和device端的计算结果进行比较，因此编写如下checkResult函数：

void checkResult(float *hostRef, float *gpuRef, const int N) {
 double epsilon = 1.0E-8;
 int match = 1;
 for (int i = 0; i < N; i++) {
    if (abs(hostRef[i] - gpuRef[i]) > epsilon) {
        match = 0;
        printf("Arrays do not match!\n");
        printf("host %5.2f gpu %5.2f at current %d\n",
        hostRef[i], gpuRef[i], i);
        break;
    }
 }
 if (match) printf("Arrays match.\n\n");
 return;
}

除了这样进行逐一比较之外，其实还有一个非常好用的方法来验证kernel函数中的计算逻辑是否正确：

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

通过设置kernel函数启动的布局为1个block和1个thread，那么相当于我们的kernel函数也是在模拟串行计算

（7）处理错误

由于许多CUDA调用是异步的，所以有时可能很难确定某个错误是由哪一步程序引起的，因此我们可以定义一个错误处理宏来封装所有的CUDA API调用，简化错误检查的过程：

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

比如在内存拷贝的时候封装，检查内存拷贝是否正常：

CHECK(cudaMemcpy(d_C, gpuRef, nBytes, cudaMemcpyHostToDevice));

也可以在kernel函数调用后检查kernel函数是否出现错误：

kernel_function<<<grid, block>>>(argument list);
CHECK(cudaDeviceSynchronize());

CHECK(cudaDeviceSynchronize())会阻塞host线程的运行直到device端所有的请求任务都结束并且确保最后的kernel函数启动部分不会出错，当然这只会用于调试过程，因为其会阻塞host端线程，成为全局屏障。

（8）编译和执行

将所有之前的代码进行整合然后编译运行，这一部分就不贴源码了，大家请到gitee仓库进行查看：

2、给核函数记时

由于CUDA是用来进行高性能并行运算的，因此了解kernel函数执行的时间是优化计算的基础，本节介绍两种简单的计算kernel函数执行时间的方法，第六章会学习如何使用CUDA特定的记时程序。

（1）用CPU计时器记时

①记时方法

使用gettimeofday系统调用来封装一个CPU计时器函数cpuSecond()，以获取系统的时钟时间，其返回的是自1970年1月1日零点以来到现在的秒数，其包含在sys/time.h头文件中：

double cpuSecond() {
    struct timeval tp;
    gettimeofday(&tp,NULL);
    return ((double)tp.tv_sec + (double)tp.tv_usec*1.e-6);
}

之后可以直接使用cpuSecond()函数来测试kernel函数的耗时：

double iStart = cpuSecond();                 // kernel函数开始前的时间
kernel_name<<<grid, block>>>(argument list);
cudaDeviceSynchronize();                     // 等待所有GPU thread运行结束才算kernel函数计算结束
double iElaps = cpuSecond() - iStart;        // kernel函数结束的时间-kernel开始前的时间≈kernel函数运行的时间

由于kernel和host程序是异步的，所以需要使用cudaDeviceSynchronize()强制同步等到所有GPU thread运行结束再记时。

②测试用例

设置一个有16M个元素的大向量进行测试：

int nElem = 1<<24;

由于我们设置先设置的block再设置的grid，所以经常会出现数据大小无法整除blockDim的情况，肯定不能丢掉最后的余数，那么就最后一个block自然就无法装满：

此时我们遇到的问题就是创建的thread线程总数大于向量元素总数，对于同样一段kernel最后的一个block中的尾部thread是没有可处理数据的，因此我们需要加以判断，防止他们访问了错误的地址导致计算错误：

__global__ void sumArraysOnGPU(float *A, float *B, float *C, const int N) {
 int i = blockIdx.x * blockDim.x + threadIdx.x; // 计算当前thread对应的数据索引，参加下方图解进行理解
 if (i < N) C[i] = A[i] + B[i]; // 判断当前thread对应的数据有没有超过向量元素总数
}

运行代码样例：

• 第0次实验：block(1024)+grid(16384)，16384个block的一维grid，每个block中有1024个thread

• 第1次实验：block(512)+grid(32768)，32768个block的一维grid，每个block中有512个thread

• 第2次实验：block(256)+grid(65536)，65536个block的一维grid，每个block中有256个thread

报错信息：block的总数超过了一维grid的限制

从上述实验可以看出两点：

• 布局影响性能：不断缩小block中thread的大小可以发现计算的时间逐步下降，这印证了之前学过的不同布局对性能有影响
• 布局存在限制： 调整布局最关键的就是了解设备对于grid和block维度的限制，其最大尺寸取决于设备，对于Fermi设备，每个block最大的thread数为1024，且grid的xyz三个方向上的维度最大值为65535

（2）用nvprof工具记时

①记时方法

自CUDA 5.0起，NVIDIA提供了一个名为nvprof的命令行分析工具，可以帮助从应用程序的CPU和GPU活动详情中获取时间线信息，其中包括kernel执行、内存传输以及CUDA API调用，其用法如下：

$ nvprof [nvprof_args] <application> [application_args]

可以使用以下语句打印nvprof的帮助信息：

$ nvprof --help

通过nvprof测试：

$ nvprof ./sumArraysOnGPU-timer

打印结果为：

其中前半部分为程序的输出，后半部分为nvprof的输出，可以看到其显示出了整个程序流程中各个计算部分的耗时，从结果中可以发现以下几点：

• CPU计时器与nvprof统计不一致：CPU计时器显示kernel函数耗时为3.26ms，而nvprof显示kernel函数耗时为2.90ms，nvprof的更加准确，因为CPU计时器还包含了来自nvprof运行的额外时间。
• 数据传输比kernel执行更耗时：可以到在host端和device端之间数据传输的时间之和是kernel函数运行的时间的几十倍

绘图可以更显著的看出各个阶段的耗时，对于这种计算时间少于数据传输时间的程序，应尽量减少host和device之间的传输。

②理论界限最大化

在进行程序优化时，学会计算理论界限非常重要，通过对计算流程进行分析可以得出当前程序的最高性能，以此对当前程序的优化程度进行判断。比如理论极限是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个字节

当然这非常的高深，说实在的我也看不懂，需要深入学习之后再来理解。

3、组织并行线程

从之前一维的例子可以看出，如果使用了合适的grid和block布局来组织thread能够对kernel函数性能产生很大的影响，本节使用一个二维矩阵加法的例子来进一步说明这一点。对于矩阵运算，传统做法就是使用一个二维的grid和二维block布局来组织线程，但是这种方法无法获得最佳性能，因此本节尝试使用以下三种布局来学习更多关于grid和block的用法：

• 由二维block构成的二维grid
• 由一维block构成的一维grid
• 由一维block构成的二维grid

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

①构建thread与矩阵数据的关系

一个矩阵通常是以行优先的方式存储在global memory中的一段连续内存中，以一个3行5列的矩阵举例：

之前学习过一维grid和一维block中线程与数据idx的对应关系，而对于二维grid和二维block以及矩阵来说需要管理3种索引：

• 线程和块索引：每个线程的threadIdx和blockIdx
• 矩阵点的坐标：每个线程对应在矩阵中的坐标ix、iy、iz
• 全局内存偏移：每个线程对应数据在global memory上距离起始地址的偏移量，其实也就是索引

第一步：将每个线程的threadIdx和blockIdx隐射到矩阵坐标上：

ix = threadIdx.x + blockIdx.x * blockDim.x;
iy = threadIdx.y + blockIdx.y * blockDim.y;

第二步：将矩阵坐标隐射到global memory的索引上：

idx = iy * nx + ix;

以矩阵（9， 12）举例，启动布局为<<<grid = (4, 3), block = (3, 3)>>>：：

②打印thread信息

通过printThreadInfo函数可以输出关于每个线程的threadIdx、blockIdx、ix、iy、idx、对应元素值，源码参考：

这些thread信息之间的关系也可以通过书中这张图进行理解，当然我觉得我上面画的图也挺好理解的哈哈哈：

（2）使用二维网格和二维块对矩阵求和

现在开始正式进入编码kernel函数阶段，分别求和函数的host和device版本：

void sumMatrixOnHost (float *A, float *B, float *C, const int nx, const int ny) {
    float *ia = A;
    float *ib = B;
    float *ic = C;
    for (int iy=0; iy<ny; iy++) {
        for (int ix=0; ix<nx; ix++) {
            ic[ix] = ia[ix] + ib[ix];
        }
        ia += nx; ib += nx; ic += nx;
    }
}
__global__ void sumMatrixOnGPU2D(float *MatA, float *MatB, float *MatC, int nx, int ny) {
    unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
    unsigned int iy = threadIdx.y + blockIdx.y * blockDim.y;
    unsigned int idx = iy * nx + ix;
    if (ix < nx && iy < ny) MatC[idx] = MatA[idx] + MatB[idx];
}

源码请参考：

这边直接给出不同grid和block布局下kernel的运行时间：

可以看到实验1几乎比实验0性能翻了一倍，但是进一步减小block布局后的实验2并没有比实验1性能更高，在第3章中将会学习为什么不同的配置会影响kernel函数的性能。

（3）使用一维网络和一维块对矩阵求和

依然以矩阵（9， 12）举例，启动布局为<<<grid = (4), block = (3)>>>：

需要注意的是，由于是一维grid和一维block所以实际启动的thread只有图中绿色的格子，下面的灰色部分是矩阵的相应坐标，需要通过在kernel函数中添加for循环来实现y方向的访问，device端kernel函数如下：

__global__ void sumMatrixOnGPU1D(float *MatA, float *MatB, float *MatC, int nx, int ny) {
    unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x; // 计算当前thread对应矩阵的ix
    if (ix < nx ) {
        for (int iy = 0; iy < ny; iy++) {                    // 循环访问矩阵y方向的坐标iy
            int idx = iy * nx + ix;                          // 计算global memory上数据的idx
            MatC[idx] = MatA[idx] + MatB[idx];
        }
    }
}

完整源码参考：

（4）使用二维网络和一维块对矩阵求和

老规矩还是以9 * 12的矩阵为例：

device中kernel函数为：

__global__ void sumMatrixOnGPUMix(float *MatA, float *MatB, float *MatC, int nx, int ny) {
    unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x; // thread对应矩阵中数据的ix
    unsigned int iy = blockIdx.y;                            // thread对应矩阵中数据的iy
    unsigned int idx = iy*nx + ix;                           // thread对应矩阵中数据在global memory中的idx
    if (ix < nx && iy < ny) MatC[idx] = MatA[idx] + MatB[idx];
}

可以看出这个函数和sumMatrixOnGPU2D非常的类似，但是通过配置grid和block省去了一次整数乘法和一次整数加法运算：

unsigned int iy = threadIdx.y + blockIdx.y * blockDim.y; // sumMatrixOnGPU2D计算iy
unsigned int iy = blockIdx.y;                            // sumMatrixOnGPUMix计算iy

完整源码参考：

（5）总结

通过以上三种布局策略：

• 由二维block构成的二维grid
• 由一维block构成的一维grid
• 由一维block构成的二维grid

得到不同的kernel运行时间如下：

通过矩阵加法的例子可以看出：

• 改变执行配置（grid和block布局）对kernel有影响
• 传统的kernel函数实现一般不能获得最佳性能
• 对于一个给定的kernel函数，尝试使用不同的grid和block布局可以获得更好的性能

4、设备管理

这一节主要介绍了CUDA中提供的几个用于查询设备信息的API，看看就好，需要用的时候再去查就行了

（1）使用运行时API查询GPU信息

参考源码：

• 查询GPU设备的所有信息：
  cudaError_t cudaGetDeviceProperties(cudaDeviceProp* prop, int device);
• 查看GPU属性
  cudaDeviceProp结构体返回GPU设备的属性，可以通过以下网址查看其内容：
  

（2）确定最优GPU

在多GPU并且每个GPU不同的情况下， 使用以下函数选择性能最好的GPU来运行kernel函数：

int numDevices = 0;
cudaGetDeviceCount(&numDevices);
if (numDevices > 1) {
    int maxMultiprocessors = 0, maxDevice = 0;
    for (int device = 0; device < numDevices; device++) {
        cudaDeviceProp props;
        cudaGetDeviceProperties(&props, device);
        if (maxMultiprocessors < props.multiProcessorCount) {
            maxMultiprocessors = props.multiProcessorCount;
            maxDevice = device;
        }
    }
    cudaSetDevice(maxDevice);
}

（3）使用nvidia-smi查询GPU信息

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

支持以下一些参数：

• MEMORY：用于获取GPU显存的使用情况。
• UTILIZATION：用于获取GPU的利用率，包括GPU内核利用率、内存利用率和总线利用率等。
• ECC：用于获取GPU ECC（Error-Correcting Code，纠错码）的信息，包括错误数和纠错数等。
• TEMPERATURE：用于获取GPU的温度信息。
• POWER：用于获取GPU的功耗信息。
• CLOCK：用于获取GPU的时钟信息，包括GPU核心时钟、内存时钟和shader时钟等。
• COMPUTE：用于获取GPU的计算能力（Compute Capability）。
• PIDS：用于获取进程ID（PID）和相关GPU的显存使用情况。
• PERFORMANCE：用于获取GPU的性能信息。
• SUPPORTED_CLOCKS：用于获取GPU支持的时钟速率信息。
• PAGE_RETIREMENT：用于获取GPU页面退役（Page Retirement）的信息，包括页面退役率等。
• ACCOUNTING：用于获取GPU资源使用的信息，包括GPU使用时间和GPU资源使用率等。

（4）在运行时设置设备

当有多设备时可以设置CUDA_VISIBLE_DEVICES来指定用于运行的GPU设备，GPU设备编号从0开始：

export CUDA_VISIBLE_DEVICES=2, 3