CUDA - 如何在CUDA C/C++内核中高效地访问全局内存

原文链接：How to Access Global Memory Efficiently in CUDA C/C++ Kernels

---

在前两篇文章（第一篇，第二篇）中，我们研究了如何在主机和设备之间高效地移动数据。在这个CUDA C/C++系列的第六篇文章中，我们将讨论如何从内核中高效地访问设备内存，特别是全局内存（global memory）。

CUDA设备上有几种内存，每种内存都有不同的作用域、生存期和缓存行为。到目前为止，在本系列文章中，我们已经使用了全局内存，它位于设备DRAM中，用于主机和设备之间的传输以及内核的数据输入源和输出位置。这里的名称全局(global) 指的是作用域，因为从主机和设备都可以访问和修改它。全局内存可以在全局（变量）作用域中使用__device__声明说明符进行声明，如以下代码段的第一行所示；也可以使用cudaMalloc()动态开辟全局内存并分配给常规C指针变量，如第7行所示。开辟的全局内存在应用程序的生存期内持续存在。根据设备的计算能力，全局内存有可能被高速缓存在芯片上。

__device__ int globalArray[256];
 
void foo()
{
    ...
    int *myDeviceMemory = 0;
    cudaError_t result = cudaMalloc(&myDeviceMemory, 256 * sizeof(int));
    ...
}

在我们讨论全局内存访问性能之前，我们需要完善对CUDA执行模型（CUDA execution model）的理解。我们已经讨论了如何将线程(threads)分组为线程块(thread blocks)，这些线程块被分配给设备上的multiprocessor。在执行过程中，线程被更精细地分组为束（warps）。GPU上的multiprocessor以SIMD（Single Instruction Multiple Data, 单指令多数据）方式为每个warp执行指令。所有当前具有CUDA功能的GPU的warp size（实际上是SIMD width）是32个线程。

全局内存合并（coalesce）

将线程分组为warp不仅与计算有关，还与全局内存访问有关。设备将warp内的线程发出的全局内存加载和存储行为合并为尽可能少的事务（transaction），以最大限度地减少DRAM带宽（在计算能力低于2.0的旧硬件上，事务合并发生在16个线程的半warp内，而不是整个warp内）。为了明确不同CUDA设备架构的发生合并的条件，我们在三个Tesla卡上进行了一些简单的实验：Tesla C870（计算能力1.0）、Tesla C1060（计算能力1.3）和Tesla C2050（计算能力2.0）。

我们进行了两个实验，使用了以下代码中显示的做自增运算的内核的两个变形（也可在GitHub上获得），一个是数组偏移，可能会导致对输入数组的不对齐访问，另一个是对输入数组进行跨步访问。

#include <stdio.h>
#include <assert.h>

// Convenience function for checking CUDA runtime API results
// can be wrapped around any runtime API call. No-op in release builds.
inline
cudaError_t checkCuda(cudaError_t result)
{
#if defined(DEBUG) || defined(_DEBUG)
  if (result != cudaSuccess) {
    fprintf(stderr, "CUDA Runtime Error: %s\n", cudaGetErrorString(result));
    assert(result == cudaSuccess);
  }
#endif
  return result;
}

template <typename T>
__global__ void offset(T* a, int s)
{
  int i = blockDim.x * blockIdx.x + threadIdx.x + s;
  a[i] = a[i] + 1;
}

template <typename T>
__global__ void stride(T* a, int s)
{
  int i = (blockDim.x * blockIdx.x + threadIdx.x) * s;
  a[i] = a[i] + 1;
}

template <typename T>
void runTest(int deviceId, int nMB)
{
  int blockSize = 256;
  float ms;

  T *d_a;
  cudaEvent_t startEvent, stopEvent;
    
  int n = nMB*1024*1024/sizeof(T);

  // NB:  d_a(33*nMB) for stride case
  checkCuda( cudaMalloc(&d_a, n * 33 * sizeof(T)) );

  checkCuda( cudaEventCreate(&startEvent) );
  checkCuda( cudaEventCreate(&stopEvent) );

  printf("Offset, Bandwidth (GB/s):\n");
  
  offset<<<n/blockSize, blockSize>>>(d_a, 0); // warm up

  for (int i = 0; i <= 32; i++) {
    checkCuda( cudaMemset(d_a, 0, n * sizeof(T)) );

    checkCuda( cudaEventRecord(startEvent,0) );
    offset<<<n/blockSize, blockSize>>>(d_a, i);
    checkCuda( cudaEventRecord(stopEvent,0) );
    checkCuda( cudaEventSynchronize(stopEvent) );

    checkCuda( cudaEventElapsedTime(&ms, startEvent, stopEvent) );
    printf("%d, %f\n", i, 2*nMB/ms);
  }

  printf("\n");
  printf("Stride, Bandwidth (GB/s):\n");

  stride<<<n/blockSize, blockSize>>>(d_a, 1); // warm up
  for (int i = 1; i <= 32; i++) {
    checkCuda( cudaMemset(d_a, 0, n * sizeof(T)) );

    checkCuda( cudaEventRecord(startEvent,0) );
    stride<<<n/blockSize, blockSize>>>(d_a, i);
    checkCuda( cudaEventRecord(stopEvent,0) );
    checkCuda( cudaEventSynchronize(stopEvent) );

    checkCuda( cudaEventElapsedTime(&ms, startEvent, stopEvent) );
    printf("%d, %f\n", i, 2*nMB/ms);
  }

  checkCuda( cudaEventDestroy(startEvent) );
  checkCuda( cudaEventDestroy(stopEvent) );
  cudaFree(d_a);
}

int main(int argc, char **argv)
{
  int nMB = 4;
  int deviceId = 0;
  bool bFp64 = false;

  for (int i = 1; i < argc; i++) {    
    if (!strncmp(argv[i], "dev=", 4))
      deviceId = atoi((char*)(&argv[i][4]));
    else if (!strcmp(argv[i], "fp64"))
      bFp64 = true;
  }
  
  cudaDeviceProp prop;
  
  checkCuda( cudaSetDevice(deviceId) )
  ;
  checkCuda( cudaGetDeviceProperties(&prop, deviceId) );
  printf("Device: %s\n", prop.name);
  printf("Transfer size (MB): %d\n", nMB);
  
  printf("%s Precision\n", bFp64 ? "Double" : "Single");
  
  if (bFp64) runTest<double>(deviceId, nMB);
  else       runTest<float>(deviceId, nMB);
}

这段代码可以通过传递fp64命令行选项以单精度（默认值）或双精度运行offset和stride内核。每个内核接受两个参数，一个输入数组和一个整数，整数表示用于访问数组元素的偏移量或步长。内核在一系列偏移和步长上循环调用。

不对齐数据访问

offset内核在Tesla C870、C1060和C2050上的结果如下图所示。

CUDA驱动程序将设备内存中分配的数组与256字节内存段（memory segments）对齐。设备可以通过与其大小对齐的32、64或128字节事务访问全局内存。对于C870或计算能力为1.0的任何其他设备，半个warp的线程的任何不对齐访问（或半个warp的线程不按顺序访问内存的对齐访问）都会导致16个独立的32字节事务。由于每个32字节事务只请求4个字节，因此当offset不是16个元素的倍数（对应于半个warp的线程数）时，预计有效带宽将减少为原来的1/8（半个warp一共请求4*16字节的数据，不对齐访问由于使用了16个独立的32字节事务，需要占用32*16字节的带宽，因此是1/8，译者注），如上图（棕色线）所示。

对于Tesla C1060或其他计算能力为1.2或1.3的设备，未对齐访问的问题较小。基本上，半个warp的线程对连续数据的未对齐访问是由“覆盖（cover）”所请求数据的少量事务来提供服务的。由于传输了未请求的数据以及不同的half-warp所请求的数据的一些重叠，因此相对于对齐的情况，仍然存在性能损失，但损失远小于C870。

计算能力为2.0的设备，如Tesla C2050，在每个line size为128字节的multiprocessor中都有一个L1缓存。这类设备将一个warp的线程的访问合并为尽可能少的缓存行（cache line），从而使线程间的顺序内存访问的吞吐量受是否对齐的影响可以忽略不计。

跨步（stride）内存访问

stride内核的结果如下图所示。

对于跨步式的全局内存访问，其结果与之前的不同。对于大的stride，无论哪种架构，有效带宽都很低。这并不奇怪：当并发线程同时访问物理内存中相距很远的内存地址时，硬件没有机会将访问进行合并。你可以在上图中看到，在Tesla C870上，任何除了1之外的stride都会导致有效带宽大幅降低。这是因为计算能力1.0和1.1的硬件需要线程间进行线性、对齐的访问从而对访问进行合并，因此我们看到了在offset内核结果中熟悉的1/8带宽。计算能力1.2及更高版本的硬件可以将对对齐的段（segment）的访问合并（计算能力1.2/1.3的设备上为32、64或128字节segment，计算能力2.0及更高的设备上为128字节cache line），因此这类硬件的带宽曲线是平滑的。

访问多维数组时，线程通常需要对数组的较高维度进行索引，因此跨步访问是不可避免的。我们可以通过使用一种称为共享内存的CUDA内存来处理这些情况。共享内存是由线程块（block）中的所有线程共享的片上内存（on-chip memory）。共享内存的一个用途是以合并的方式将多维数组的2D数据片（2D tile）从全局内存中提取到共享内存中，然后使连续的线程跨过共享内存数据片。与全局内存不同，跨步访问共享内存不会带来额外开销。我们将在下一篇文章中详细介绍共享内存。

总结

在这篇文章中，我们讨论了如何从CUDA内核代码中有效地访问全局内存的一些方面。设备上的全局内存访问与主机上的数据访问有相同的性能特性；即数据位置非常重要。在早期的CUDA硬件中，内存访问对齐与线程间的位置一样重要，但在最近的硬件中，对齐并不重要。另一方面，跨步式内存访问会损害性能，使用片上共享内存可以缓解这种情况。在下一篇文章中，我们将详细探讨共享内存，再下一篇，将展示如何使用共享内存来避免矩阵转置时的跨步全局内存访问。