CUDA编程模型详解
本文以vectorAdd为例，通过描述C在CUDA中的使用（vectorAdd这个例子可以在CUDA sample中找到。）来介绍CUDA编程模型的主要概念。CUDA C的进一步描述可以参考《Programming Interface》。
主要内容包括：
1、Kernels（核函数）
2、Thread Hierarchy（线程结构）
3、Memory Hierarchy（存储结构）
4、Heterogeneous Programming（异构编程）
5、Compute Capability（计算能力）
6、实例分析
下面分别进行详细地讲解。
1、Kernels（核函数）
CUDA C通过允许编程人员自定义C功能函数（叫做kernels）来扩展C。此功能可以理解为：以并行方式由N个不同的CUDA线程编译N次的程序相当于由常规C编译一次一样。
核函数用_global_来声明定义，对于一个给定的内核调用，执行该内核的CUDA线程数量用一个<<<...>>>执行配置语法（详见C Language Extensions）来指定。每一个执行该内核的线程会分配一个唯一的线程ID，通过内置的threadIdx变量在内核中访问。
下面给出一个例子进行说明：
//功能：把两个大小为N的向量A和向量B相加，并把结果存储在向量C中。
// 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);
...
}此处，N个线程中的每一个线程都会执行VecAdd（），是一个双向加法的例子。
2、Thread Hierarchy（线程结构）
为了方便，threadIdx是一个“三组成”向量，因此，线程能够被一维、二维、三维的线程索引识别，组成一维、二维、三维的线程块。这提供了一个很自然的方式，即在一个向量、矩阵或册领域调用整个元素进行计算。
线程的索引及其ID以一种简单的方式互相关：对于一维块，它们等价；对于大小为（Dx,Dy）的二维块，一个线程的索引（x,y）的线程ID是（x+yDx）；对于大小为（Dx,Dy,Dz）的三维线程，一个线程的索引（x,y,z）的线程ID是（x+yDx+zDxDy）。
下面给出一个例子：
//功能：将大小为NxN的两个矩阵A和B相加，把结果存储在矩阵C中。
// Kernel definition
__global__ void MatAdd(float A[N][N], float B[N][N],
float C[N][N])
{
int i = threadIdx.x;
int j = threadIdx.y;
C[i][j] = A[i][j] + B[i][j];
}
int main()
{
...
// Kernel invocation with one block of N * N * 1 threads
int numBlocks = 1;
dim3 threadsPerBlock(N, N);
MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
...
}对于每个块，此处对线程的数量有一个限制，因为一个块的所有线程都依赖于同一个处理器的核，它们必须共享此核有限的内存资源。在目前的GPUs，一个线程块最大可包含1024个线程。
然而，一个核函数可以被多个等型的线程块执行，因此线程的总数等于每一个块线程数乘以块数量。
块可以被组织为一维、二维或者三维的线程块网格，如下图所示，在一个网格内线程块的数目经常由将要处理的数据大小或者系统中的处理器核数决定，当然也可以大大超过此数目。
每一个块的线程数目和每一个网格中的块的数目可以由语法<<<...>>>指定，其类型可以是int或者dim3。二维的块或者网格指定方法与上图例子中一致。
在网格中的每一个块可以由一个一维、二维或三维内核索引通过内置变量blockIdx访问确定。线程块的维数可以通过内置变量blockDim得到。
举个例子说明：
若我们把上一个例子扩展为可以处理多个块的函数，其代码如下：
// Kernel definition
__global__ void MatAdd(float A[N][N], float B[N][N],
float C[N][N])
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
int j = blockIdx.y * blockDim.y + threadIdx.y;
if (i < N && j < N)
C[i][j] = A[i][j] + B[i][j];
}
int main()
{
...
// Kernel invocation
dim3 threadsPerBlock(16, 16);
dim3 numBlocks(N / threadsPerBlock.x, N / threadsPerBlock.y);
MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);
...
}
通常，一个线程块的大小为16x16，即256个线程，尽管块大小可以任意，但此种情况是最常见的选择。网格可以像以前一样由足够的块创建，每一个矩阵元素都有一个线程。为简单起见，这个例子假定了在每一个维中每个网格的线程数目，在那个维中被每一个块的线程数目均匀地区分（尽管不需要如此）。
线程块需要独立执行：必须能够以并联、串联或者任何顺序执行。这种独立性的要求允许线程块可以在任意数目的核内以任意顺序被调用，如下图所示，使得程序员能够编写与核数尺寸一致的代码。
一个块中的线程能够通过一些共享内存共享数据和同步它们的执行协调内存访问来配合完成任务。更确切地说，可以通过调用内部函数_syncthreads()在内核中指定同步点；_syncthreads()像一个屏障，使得块中的所有线程必须在得到允许前等待执行，使用共享内存的例子可以参见《Shared Memory》。
为了有效合作，希望共享内存在相邻合作的处理器有较低的延迟缓存，就像L1 cache一样，同时还希望_syncthreads()是低量级的。
3、Memory Hierarchy（存储结构）
如下图所示，CUDA线程在执行过程中可以访问来自多个内存空间的数据。每一个线程都有私有的本地内存，每一个线程块都共享了内存，使得对块中的所有线程都是可见的，并且与该块有相同的生命周期，所有的线程都可以进入相同的全局内存。
对于所有的线程，也有两个额外的只读存储器访问空间：常量和纹理内存空间。全局、常量和纹理存储空间针对不同的内存用途也有被优化（具体见《Device Memory Accesses》）。对于指定的数据格式，纹理内存与数据滤波器一样也提供了不同的寻址方式。
运行相同的应用程序时，全局、常量和纹理内存空间在内核中是持续的（persistent）。
4、Heterogeneous Programming（异构编程）
如下图所示，CUDA编程模型假定CUDA线程执行在独立的物理设备上，就像主机在协处理器上运行C代码。例如在这种情况下，内核执行在一个GPU上，而其余的C程序执行在一个CPU上。
CUDA编程模型也假定主机和设备在DRAM都拥有独立的内存空间，分别称为主机内存和设备内存。因此，一个程序管理着全局、常量和纹理内存空间，通过CUDA runtime调用，使其对内核可见（详见《Programming Interface》），包括在主机和设备内存中的设备内存分配、回收以及数据传输。
5、Compute Capability（计算能力）
设备的计算能力由a major revision number和a minor revision number定义。
相同修订号的设备具有相同的核结构，主修订号是3的设备是基于Kepler架构的，主修订号为2的设备是基于Fermi架构的，主修订号为1的是基于Tesla架构的。
副修订号对应于核架构的升级改进，通常会新添新的特性。
6、实例分析
/**
* Copyright 1993-2014 NVIDIA Corporation. All rights reserved.
*
* Please refer to the NVIDIA end user license agreement (EULA) associated
* with this source code for terms and conditions that govern your use of
* this software. Any use, reproduction, disclosure, or distribution of
* this software and related documentation outside the terms of the EULA
* is strictly prohibited.
*
*/
/**
* Vector addition: C = A + B.
*
* This sample is a very basic sample that implements element by element
* vector addition. It is the same as the sample illustrating Chapter 2
* of the programming guide with some additions like error checking.
*/
#include <stdio.h>
// For the CUDA runtime routines (prefixed with "cuda_")
#include <cuda_runtime.h>
/**
* CUDA Kernel Device code
*
* Computes the vector addition of A and B into C. The 3 vectors have the same
* number of elements numElements.
*/
__global__ void
vectorAdd(const float *A, const float *B, float *C, int numElements)
{
int i = blockDim.x * blockIdx.x + threadIdx.x;
if (i < numElements)
{
C[i] = A[i] + B[i];
}
}
/**
* Host main routine
*/
int main(void)
{
// Error code to check return values for CUDA calls
cudaError_t err = cudaSuccess;
// Print the vector length to be used, and compute its size
int numElements = 50000;
size_t size = numElements * sizeof(float);
printf("[Vector addition of %d elements]\n", numElements);
// Allocate the host input vector A
float *h_A = (float *)malloc(size);
// Allocate the host input vector B
float *h_B = (float *)malloc(size);
// Allocate the host output vector C
float *h_C = (float *)malloc(size);
// Verify that allocations succeeded
if (h_A == NULL || h_B == NULL || h_C == NULL)
{
fprintf(stderr, "Failed to allocate host vectors!\n");
exit(EXIT_FAILURE);
}
// Initialize the host input vectors
for (int i = 0; i < numElements; ++i)
{
h_A[i] = rand()/(float)RAND_MAX;
h_B[i] = rand()/(float)RAND_MAX;
}
// Allocate the device input vector A
float *d_A = NULL;
err = cudaMalloc((void **)&d_A, size);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to allocate device vector A (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
// Allocate the device input vector B
float *d_B = NULL;
err = cudaMalloc((void **)&d_B, size);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to allocate device vector B (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
// Allocate the device output vector C
float *d_C = NULL;
err = cudaMalloc((void **)&d_C, size);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to allocate device vector C (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
// Copy the host input vectors A and B in host memory to the device input vectors in
// device memory
printf("Copy input data from the host memory to the CUDA device\n");
err = cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to copy vector A from host to device (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
err = cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to copy vector B from host to device (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
// Launch the Vector Add CUDA Kernel
int threadsPerBlock = 256;
int blocksPerGrid =(numElements + threadsPerBlock - 1) / threadsPerBlock;
printf("CUDA kernel launch with %d blocks of %d threads\n", blocksPerGrid, threadsPerBlock);
vectorAdd<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, numElements);
err = cudaGetLastError();
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to launch vectorAdd kernel (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
// Copy the device result vector in device memory to the host result vector
// in host memory.
printf("Copy output data from the CUDA device to the host memory\n");
err = cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to copy vector C from device to host (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
// Verify that the result vector is correct
for (int i = 0; i < numElements; ++i)
{
if (fabs(h_A[i] + h_B[i] - h_C[i]) > 1e-5)
{
fprintf(stderr, "Result verification failed at element %d!\n", i);
exit(EXIT_FAILURE);
}
}
printf("Test PASSED\n");
// Free device global memory
err = cudaFree(d_A);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to free device vector A (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
err = cudaFree(d_B);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to free device vector B (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
err = cudaFree(d_C);
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to free device vector C (error code %s)!\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
// Free host memory
free(h_A);
free(h_B);
free(h_C);
// Reset the device and exit
// cudaDeviceReset causes the driver to clean up all state. While
// not mandatory in normal operation, it is good practice. It is also
// needed to ensure correct operation when the application is being
// profiled. Calling cudaDeviceReset causes all profile data to be
// flushed before the application exits
err = cudaDeviceReset();
if (err != cudaSuccess)
{
fprintf(stderr, "Failed to deinitialize the device! error=%s\n", cudaGetErrorString(err));
exit(EXIT_FAILURE);
}
printf("Done\n");
return 0;
}
到这儿，CUDA编程模就基本结束了。
