CUDA编程——常用存储器的分类与介绍

CUDA的存储器从物理上可分为两类：

• 板载显存（On-board memory）
• 片上内存（On-chip memory）

其中板载显存主要包括全局内存（global memory）、本地内存（local memory）、常量内存（constant memory）和纹理内存（texture memory），而片上内存主要包括寄存器（register）和共享内存（shared memory）。它们的主要特点如下表所列：

片上内存的读写效率通常比板载显存更快，而寄存器又是所有显存类型中最快的存储器。本文我们将分别介绍这些存储器的应用场景及用法。

01

寄存器与本地内存

从代码实现上看，寄存器变量与本地内存变量的定义方式是一样的，它们都是定义于CUDA核函数中的变量，很像C/C++函数中的局部变量定义，比如以下核函数代码中的变量A、index、tmp都属于寄存器变量或本地变量：

__global__ void cuda_kernel(float *a, float b, float *c, int row, int col)
{
    int x = threadIdx.x + blockDim.x * blockIdx.x;  //col
    int y = threadIdx.y + blockDim.y * blockIdx.y;  //row
    
    if(x < col && y < row)
    {
        float A[200];
        int index = y*col + x;
        float tmp = a[index]*b[index] + b[index]*b[index];
        c[index] = tmp*tmp;
    }
}

那么定义于CUDA核函数中的变量，什么时候是寄存器变量，什么时候是本地变量呢？

通常以下三种情况下定义的变量为本地变量，其余情况则是寄存器变量：

• 在编译阶段编译器无法确定数组的值，这种情况下该数组是本地变量，也即其数据内容存储在本地内存。
• 如果数组或结构体占用内存空间很大，则系统将其分配到本地内存，也即本地变量。
• 寄存器空间是很小的，如果核函数中定义了很多变量，那些超过寄存器空间限制的变量则被分配到本地内存，也即本地变量。

由以上可知，一般不会在核函数中定义太多或者太大的变量，不然系统自动将超出寄存器限制的变量分配到本地内存，影响程序运行效率。

02

共享内存

共享内存的主要特点在于“共享”，也即同一个线程块中的所有线程都可以对这一块存储进行读写操作，所以“共享”是针对同一个线程块中所有线程而言的。一旦共享内存被定义并指定大小，系统将给所有线程块都分配相同大小的共享内存，比如定义一个大小为8 bytes的unsigned char型共享内存，那么所有线程块都会被分配一个8 bytes的unsigned char型共享内存。

前文我们已经详细介绍共享内存的特点与应用，此处不再重复：

CUDA加速——共享内存介绍及其应用

03

常量内存

常量内存在device端（GPU端）只读，在host端（CPU端）可读可写，通常情况下使用__constant__修饰的变量，其数据存储于常量内存，并且该变量为全局变量，对同一个.cu文件中且定义于其后面的所有核函数都可见。比如以下代码中，变量A的值存储于常量内存，kernel1和kernel2都可以使用A，但kernel0因定义在A前面而不能使用A。

__global__ void kernel0()
{


}


__constant__ float A[128];


__global__ void kernel1()
{


}


__global__ void kernel2()
{


}

常量内存在device端只读，因此只能在host端对其初始化和修改，通过调用cudaMemcpyToSymbol函数实现，比如以下代码：

__constant__ float A[10];


void init_constant(void)
{
     float B[10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};
     //数组A定义在device端的常量内存，数组B定义在host端，使用数组B的内容初始化数组A
     cudaMemcpyToSymbol(A, B, 10*sizeof(float));  
}

什么情况下使用常量内存呢？

首先我们来讲一下warp和half-warp的概念。

从软件的角度来看，一个block包含的最大线程数通常为512或1024，多个block包含的所有线程都是并行执行的。然而从硬件的角度来看并不是这样，硬件上把每个block中的所有线程分成每32个一组的线程束，一个线程束称为一个warp。同一个warp中的线程才是真正意义上的并行执行，并且它们使用各自的数据执行相同的处理指令。不同warp的执行由系统调度，所以不同的warp不一定并行执行。

顾名思义，half-warp就是半个warp，也即同一个warp中的16个线程。

接下来我们介绍使用常量内存的好处，主要有两个：

• 针对一个half-warp包含的16个线程，GPU只需要对常量内存执行一次读操作，这16个线程就都能获取到数据，而不需要执行16次读操作。
• 常量内存具有缓存（cache）机制，如果GPU访问过常量内存的一个地址，则将该地址的值缓存，再次访问该地址时直接从cache中获取其值，不需要再进行一次读操作，因此速度会快很多。

由以上可知，在每个half-warp中16个线程需要访问相同内存地址的情况下，使用常量内存能够大大提升效率，但是如果每个half-warp中16个线程需要访问不同的内存地址，这种情况下则不适合使用常量内存。

04

全局内存

全局内存是GPU上容量最大的存储器，可达到10 GB，所以CUDA编程时通常把较大的数据存储在全局内存，因此全局内存也是所有GPU存储器中最常用的存储器。下面我们分别介绍全局内存的申请、拷贝、使用和释放。

• 全局内存的申请

通过调用cudaMalloc函数可以方便地申请全局内存，不过需要注意该函数第三个参数单位是byte，因此如果数据类型不是char/unsigned char类型，那么第三个参数需要将数据长度再乘以数据类型所占字节数（sizeof）：

//定义指针
unsigned char *A;
int *B;
float *C;
//定义数据长度
const int data_len = 128;
//申请全局内存
cudaMalloc((void**)&A, data_len * sizeof(unsigned char));
cudaMalloc((void**)&B, data_len * sizeof(int));
cudaMalloc((void**)&C, data_len * sizeof(float));

• 全局内存的拷贝

通常调用cudaMemcpy函数将数据从host端内存拷贝到device端全局内存，或者从device端全局内存拷贝到host端内存。同样，需要注意该函数第三个参数单位也是byte，因此第三个参数需要将数据长度再乘以数据类型所占字节数：

float A[8] = {0, 1, 2, 3, 4, 5, 6, 7};
float *C;


cudaMalloc((void**)&C, 8 * sizeof(float));
//将数据从host端内存拷贝到device端全局内存
cudaMemcpy(C, A, 8 * sizeof(float), cudaMemcpyHostToDevice);
//将数据从device端全局内存拷贝到host端内存
cudaMemcpy(A, C, 8 * sizeof(float), cudaMemcpyDeviceToHost);

注意以上代码中决定拷贝方向的是第四个参数，如果是从host拷贝到device，那么该参数是cudaMemcpyHostToDevice，如果是device到host，那么该参数是cudaMemcpyDeviceToHost。

• 全局内存的使用和释放

一个典型的CUDA并行任务流程是这样的：

(1) 把数据从host端拷贝到device端（通常是拷贝到device端的全局内存）。

(2) 在device端开启多线程并行处理数据。

(3) 待开启的所有线程处理数据完毕，将最后处理结果从device端再拷贝回host端。

所以全局内存通常在CUDA核函数中使用，通常使用线程id号来索引全局内存中的数据，使每个线程与存储数据的全局内存地址一一对应。下面我们举一个简单例子来说明全局内存的使用。

假如有两个相同尺寸的矩阵A和矩阵B，且A、B都是float型数据矩阵，现在要使用CUDA并行计算A、B相同坐标点数据的平方差。

首先，是CPU实现代码，循环遍历所有点计算平方差即可：

void CPU_cal(Mat A, Mat B)
{
  Mat C(A.size(), CV_32FC1);
  
  for(int i = 0; i < A.rows; i++)  //行遍历
  {
    float *pA = A.ptr<float>(i);
    float *pB = B.ptr<float>(i);
    float *pC = C.ptr<float>(i);
    
    for(int j = 0; j < A.cols; j++)  //列遍历
    {
      pC[j] = (pA[j] - pB[j])*(pA[j] - pB[j]);
    }
  }


}

接着是GPU实现代码，开启多线程并行计算每个点的平方差：

/*
CUDA核函数
*/
__global__ void GPU_cal_kernel(float *A_cuda, float *B_cuda, float *C_cuda, int row, int col)
{
  //线程的x方向id
  int x = threadIdx.x + blockDim.x * blockIdx.x;  //col
  //线程的y方向id
  int y = threadIdx.y + blockDim.y * blockIdx.y;  //row
  
  if(x < col && y < row)
  {
    //将线程的二维id转换为全局内存的一维地址索引，并保存到寄存器变量index
    int index = y * col + x;
    
    //根据一一对应关系，线程(x, y)对应全局内存地址index = y * col + x
    //也即线程(x, y)负责处理全局内存地址index保存的数据
    //使用index来索引A_cuda、B_cuda，就相当于从全局内存的index地址读取数据
    //将从A_cuda、B_cuda读取的index地址数据相减，并把差值保存到寄存器变量diff
    float diff = A_cuda[index] - B_cuda[index];
    
    //寄存器变量diff保存了差值，因此计算平方的时候可直接使用diff中保存的值
    //得到平方值之后，再将结果保存到全局内存C_cuda的index位置
    //使用index索引C_cuda，并对其赋值，相当于对全局内存C_cuda的index地址进行写操作
    C_cuda[index] = diff * diff;
  }
}


/*
调用以上核函数
*/
void GPU_cal(Mat A, Mat B)
{
  float *A_cuda, *B_cuda, *C_cuda;
  //计算数据长度，注意如果是byte长度还需再乘以sizeof(float)
  const int data_len = A.rows * A.cols;
  //申请全局内存
  cudaMalloc((void**)&A_cuda, data_len * sizeof(float));
  cudaMalloc((void**)&B_cuda, data_len * sizeof(float));
  cudaMalloc((void**)&C_cuda, data_len * sizeof(float));
  //将数据从host内存拷贝到device全局内存
  cudaMemcpy(A_cuda, (float *)A.data, data_len * sizeof(float), cudaMemcpyHostToDevice);
  cudaMemcpy(B_cuda, (float *)B.data, data_len * sizeof(float), cudaMemcpyHostToDevice);
  //定义线程块、线程块中的线程都为二维索引
  dim3 cuda_Block(16, 16);   //每个线程块有16*16个线程
  int M = (A.cols + cuda_Block.x - 1) / cuda_Block.x;
  int N = (A.rows + cuda_Block.y - 1) / cuda_Block.y;
  dim3 cuda_Grid(M, N);  //线程网格总共有M*N个线程块
  //调用核函数并行处理
  GPU_cal_kernel<<<cuda_Grid, cuda_Block>>>(A_cuda, B_cuda, C_cuda, A.rows, A.cols);
  
  Mat C(A.size(), CV_32FC1);
  //将并行计算结果从device全局内存拷贝到host内存
  cudaMemcpy((float *)C.data, C_cuda, data_len * sizeof(float), cudaMemcpyDeviceToHost);
  //释放申请的全局内存
  cudaFree(A_cuda);
  cudaFree(B_cuda);
  cudaFree(C_cuda);


}

需注意，对于申请的全局内存，如果不再使用，必须调用cudaFree函数将其释放。而且以上代码中线程块、线程的索引都是二维的，其实它们的索引还可以是一维或三维，以后我们再详细讲怎么通过线程块、线程的一维或二维或三维索引来确定其对应的全局内存地址。

此外，针对全局内存的访问有合并、对齐的说法，只有在合并、对齐的情况下才能高效地访问全局内存，在下篇文章我们再详细探讨这个话题。

05

纹理内存

纹理内存是GPU中的一种只读存储器，其使用方式为将某一段全局内存绑定到纹理内存，这段全局内存通常的表现形式为一维CUDA数组/全局内存、二维或三维CUDA数组，然后通过读取纹理内存（也称为纹理拾取）来获取全局内存的数据。相比全局内存的访问要求对齐、合并，纹理内存对非对齐访问和随机访问具有良好的加速效果。

• 一维纹理

一维纹理可以绑定到CUDA数组，也可以直接绑定到全局内存。下面举一个简单的例子来介绍一维纹理内存的使用。

首先是一维纹理的定义，纹理内存通常定义为全局变量：

//float表示数据类型
//cudaTextureType1D、cudaTextureType2D、cudaTextureType3D分别表示一维、二维、三维
//cudaReadModeElementType表示只读模式
//tex_1D为定义的纹理内存变量
texture<float, cudaTextureType1D, cudaReadModeElementType> tex_1D;

其次，是纹理内存的绑定，这里我们直接把全局内存绑定到纹理：

const int data_len = 128;
const int data_size = data_len * sizeof(float);
//初始化host端数组
float *data_host = (float *)malloc(data_size);
for(int i = 0; i < data_len; i++)
{
  data_host[i] = i; 
}
//申请device端全局内存
float *data_device;
cudaMalloc((void**)&data_device, data_size);
//将数据从host端拷贝到device端全局内存
cudaMemcpy(data_device, data_host, data_size, cudaMemcpyHostToDevice);
//参数一表示以bytes为单位的偏移量，也即绑定到纹理的全局内存的起始偏移地址
//参数二为纹理内存变量
//参数三为全局内存地址变量
cudaBindTexture(0, tex_1D, data_device);  

接着是在核函数中纹理拾取，通过调用tex1Dfetch函数实现：

//功能：将纹理内存的数据拷贝到全局内存A_cuda
__global__ void cuda_kernel(float *A_cuda, int data_len)
{
  //线程id
  int x = threadIdx.x + blockDim.x * blockIdx.x;
  
  if(x < data_len)
  {
    //使用线程id来索引全局内存A_cuda和纹理内存tex_1D
    A_cuda[x] = tex1Dfetch(tex_1D, x);
  }
}

最后是纹理内存的解绑，函数执行完毕之后需要对纹理内存进行解绑：

cudaUnbindTexture(tex_1D);

完整代码：

texture<float, cudaTextureType1D, cudaReadModeElementType> tex_1D;


//功能：将纹理内存的数据拷贝到全局内存A_cuda
__global__ void cuda_kernel(float *A_cuda, int data_len)
{
  //线程id
  int x = threadIdx.x + blockDim.x * blockIdx.x;
  
  if(x < data_len)
  {
    //使用线程id来索引全局内存A_cuda和纹理内存tex_1D
    A_cuda[x] = tex1Dfetch(tex_1D, x);
  }
}


void cuda_copy_data(float *data_host_dst)
{
  const int data_len = 128;
  const int data_size = data_len * sizeof(float);
  
  //初始化host端数组
  float *data_host = (float *)malloc(data_size);
  for(int i = 0; i < data_len; i++)
  {
    data_host[i] = i; 
  }
  
  //申请device端全局内存
  float *data_device, *data_dst;
  cudaMalloc((void**)&data_device, data_size);
  cudaMalloc((void**)&data_dst, data_size);
  //将数据从host端拷贝到device端全局内存
  cudaMemcpy(data_device, data_host, data_size, cudaMemcpyHostToDevice);
  
  //参数一表示以bytes为单位的偏移量，也即绑定到纹理的全局内存的起始偏移地址
  //参数二为纹理内存变量
  //参数三为全局内存地址变量
  cudaBindTexture(0, tex_1D, data_device);  
  
  dim3 tex_Block(16);  //每个block有16个线程
  //总共有((data_len + 15) / 16)个block
  dim3 tex_Grid((data_len + tex_Block.x - 1) / tex_Block.x);
  
  //调用核函数
  cuda_kernel<<<tex_Grid, tex_Block>>>(data_dst, data_len);
  
  //将数据从device端拷贝到host端
  cudaMemcpy(data_host_dst, data_dst, data_size, cudaMemcpyDeviceToHost);
  
  //纹理解绑
  cudaUnbindTexture(tex_1D);
  //释放全局内存
  cudaFree(data_device);
  cudaFree(data_dst);
  //释放host内存
  free(data_host);


}

• 二维纹理

通常将存储一张二维图像的全局内存绑定到二维纹理，在核函数中可高效地随机访问二维纹理数据。

以全局变量的方式定义二维纹理：

texture<float, cudaTextureType2D, cudaReadModeElementType>  tex_src;

定义CUDA数组，并申请CUDA数组内存：

//声明数据类型为float
cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();
//定义CUDA数组
cudaArray *cuArray_src;
//分配大小为col_c*row_c的CUDA数组
int col_c = 512;
int row_c = 512;
cudaMallocArray(&cuArray_src, &channelDesc, col_c, row_c);  

设置纹理内存参数，并将CUDA数组绑定到纹理内存：

//寻址方式
//cudaAddressModeWrap--循环寻址，如果超出最大地址则转成从最小地址开始
//cudaAddressModeClamp--钳位寻址，如果超出最大地址则访问最大地址
tex_src1.addressMode[0] = cudaAddressModeWrap;
tex_src1.addressMode[1] = cudaAddressModeWrap;
//是否对纹理坐标归一化
tex_src1.normalized = false;
//纹理的滤波模式：
//cudaFilterModePoint--最邻近插值
//cudaFilterModeLinear--双线性插值
tex_src1.filterMode = cudaFilterModePoint;
//纹理绑定，将CUDA数组绑定到纹理tex_src
cudaBindTextureToArray(&tex_src, cuArray_src, &channelDesc);  

将数据从host内存拷贝到CUDA数组：

Mat M_float = Mat::zeros(row_c, col_c, CV_32FC1);
cudaMemcpyToArray(cuArray_src, 0, 0, (float *)M_float.data, row_c*col_c*sizeof(float), cudaMemcpyHostToDevice);

在CUDA核函数中调用tex2D函数进行纹理拾取：

//tex_src--要拾取的纹理内存
//x--纹理内存的x坐标
//y--纹理内存的y坐标
//功能：将纹理内存中(x,y)坐标处的数据加载到寄存器变量d
float d = tex2D(tex_src, x, y);

纹理解绑并释放CUDA数组：

cudaUnbindTexture(tex_src);
cudaFreeArray(cuArray_src);

• 三维纹理

有时候需要使用CUDA处理多帧的图像，把多帧时间序列的图像传入到GPU中，此时就可以把保存图像的全局内存绑定到二维纹理内存（将每帧图像展开拼接为一行），核函数通过纹理拾取来访问输入的图像数据。不过二维纹理内存的宽是有限制的：

cudaMallocArray函数的第三个参数img_size为宽，也即每帧图像的总数据个数。然而二维纹理对宽是有限制的，如果宽超过64K就会出错。所以如果每帧图像的大小超过了64K，则不能使用二维纹理内存，这时候可以使用三维纹理内存。

  cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();
  cudaMallocArray((cudaArray**)&arr_mat_x, &channelDesc, img_size, IIR_N_X); 

三维纹理的使用与一维、二维纹理区别较大，下面将详细说明三维纹理内存的使用。

以全局变量的方式定义三维纹理：

texture<float, cudaTextureType3D, cudaReadModeElementType>  tex_mat;    //定义为3D类型的纹理内存

定义三维的CUDA数组：

cudaArray *arr_mat;   //定义CUDA数组
cudaExtent extent;   //定义图像的尺寸和帧数结构体
cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();  //定义数据类型为float
//创建extent时，在旧版本中宽度w以字节为单位，即必须乘上sizeof(DTYPE)，新版本已经不用乘以sizeof(DTYPE)，否则反而会出错!文档和函数说明都没有改过来，这里是坑！
extent.width = col;     //每帧图像的列数，这里不需要再乘以sizeof(float)
extent.height = row;   //每帧图像的行数
extent.depth = picnum;   //图像的总帧数
//创建picnum帧row*col的存储空间
cudaMalloc3DArray((cudaArray**)&arr_mat, &channelDesc, extent);   

将数据从host内存拷贝到CUDA数组中，其中pic为Mat类型的vector数组，其包含了picnum帧图像：

vector<Mat> pic;

cudaMemcpy3DParms HostToDev = {0};   //定义数据传输的结构体
HostToDev.dstArray = arr_mat;    //指定数据传输的目标地址为cuda数组
HostToDev.extent = make_cudaExtent(col, row, 1);      //创建extent时，在旧版本中宽度w以字节为单位，即必须乘上sizeof(DTYPE)，新版本已经不用乘以sizeof(DTYPE)，否则反而会出错!文档和函数说明都没有改过来，这里是坑！  
HostToDev.kind = cudaMemcpyHostToDevice;  //定义传输方向为CPU到GPU显存
HostToDev.srcPos = make_cudaPos(0, 0, 0);   //定义数据传输的源地址的偏移量(w, h, img_index)
for(int i = 0; i < picnum; i++)   //拷贝多帧图像到cuda数组
{
  //指定数据传输的源地址，注意这里的第二个参数需要乘以数据类型所占的字节数
  HostToDev.srcPtr = make_cudaPitchedPtr((void *)pic[i].data, col*sizeof(float), col, row);        
  HostToDev.dstPos = make_cudaPos(0, 0, i);   //指定目标地址的偏移量，分别为x，y，z地址
  cudaMemcpy3D(&HostToDev);     //根据以上设置的参数实行拷贝
}

设置三维纹理参数，并将CUDA数组绑定到三维纹理：

tex_mat.normalized = 0;   //索引地址不归一化
//filterMode：滤波模式。仅对绑定 CUDA 数组的纹理有效。当使用浮点型的坐标寻址纹理时，将根据设定返回不同类型的值。设定可以有：cudaFilterModePoint 和 cudaFilterModeLinear。分别表示最近邻插值和线性插值
tex_mat.filterMode = cudaFilterModePoint;    
tex_mat.addressMode[0] = cudaAddressModeClamp;  //寻址模式，即如何处理越界的纹理坐标。可设置：cudaAddressModeClamp 和 cudaAddressModeWrap。Clamp 即钳位模式，Wrap 为循环模式。循环模式只支持归一化的纹理坐标
tex_mat.addressMode[1] = cudaAddressModeClamp;
tex_mat.addressMode[2] = cudaAddressModeClamp;
tex_mat.channelDesc = channelDesc;   //描述纹理返回值类型，同cuda数组部分的内容
cudaBindTextureToArray(tex_mat, (cudaArray *)arr_mat, channelDesc);   //绑定纹理内存

在核函数中调用tex3D函数执行纹理拾取：

//后面三个参数分别是x，y，z坐标
tex3D(tex_mat, x, y, z); 

最后是纹理解绑和释放CUDA数组（这个与二维纹理一样）：

cudaUnbindTexture(tex_mat);
cudaFreeArray(arr_mat);

• 纹理内存的硬件插值功能

纹理内存具有硬件插值功能，包括最邻近插值和双线性插值这两种插值方式。如果纹理拾取时输入的访问坐标地址是浮点数，纹理内存将自动根据设置插值方式对浮点坐标进行插值，然后返回插值结果。这个插值过程不需要开发者来实现，是硬件自动完成的，开发者只需要设置好插值方式为最邻近插值或者双线性插值即可，因此可以节省很多计算时间。

这个我们前文已经介绍，详细请参考：

CUDA纹理内存--硬件插值功能的应用