38、通用GPU编程入门指南

通用GPU编程入门指南

1. 多核GPU概述

多核GPU在如今几乎所有的硬件平台上都有广泛应用，从标准桌面计算机到计算机集群。最初，这些处理器是为图形应用程序设计的，但现在，它们在科学计算和科学模拟等非图形应用领域的重要性日益增加。对于数据并行程序而言，使用GPU可以显著提高效率，这主要得益于GPU的特定硬件设计，它针对图形应用的大量数据和大量线程执行的浮点运算高吞吐量进行了优化。如今，一个GPU可能包含数百个核心来执行这些线程。

2. GPU编程环境

早期，由于图形编程环境（如DirectX或OpenGL）的限制，使用GPU进行通用非图形应用和模拟的难度极大。近年来，出现了更适合GPU编程的环境，如CUDA和OpenCL。
- CUDA ：由NVIDIA自2007年起为新一代GPU提供支持的通用并行编程环境，也可以在CPU上进行模拟。
- OpenCL ：由包括苹果、英特尔、AMD/ATI和NVIDIA等工业伙伴联合开发的GPU标准化编程模型。

3. GPU架构

3.1 架构特点

GPU从图形加速器发展而来，其架构与通用CPU架构的发展相对独立。由于图形应用中存在高度的数据并行性，GPU架构比CPU架构更早地采用了多个专用处理器核心。如今，单个GPU可以包含数百个计算核心，远多于当前CPU技术中的核心数量。

除了图形处理，GPU还可用于通用非图形应用，前提是这些应用程序具有足够的数据并行性，能够充分利用GPU中的大量计算核心。科学模拟通常具备这种特性，许多实例表明，这类应用程序使用多核GPU架构可以获得比CPU更好的计算性能。这一趋势促使GPU制造商（如NVIDIA）开发了CUDA和OpenCL等并行编程环境，以支持GPU上的非图形应用，并提供适合GPU架构的特定编程模型。

3.2 架构组成

GPU由多个多线程SIMD处理器组成，这些处理器是独立的MIMD处理器核心，各自处理独立的计算指令序列。一个GPU中包含的SIMD处理器数量取决于具体的GPU型号。例如，属于Fermi架构家族的NVIDIA GTX480 GPU最多可拥有15个独立的SIMD处理器。

每个SIMD处理器包含多个SIMD功能单元，这些单元可以对不同的数据执行相同的SIMD指令。每个SIMD功能单元都有独立的寄存器组，SIMD指令所需的数据必须存储在这些本地寄存器中。GPU提供了特定的传输操作，用于将数据从内存传输到寄存器，但由于数据可能位于GPU的全局内存（片外）中，实际传输可能需要多个机器周期。较新的GPU架构包含缓存内存层次结构，因此一些数据传输操作可以使用缓存数据，从而比访问非缓存数据更快。

3.3 线程调度

为了隐藏数据传输操作带来的不确定等待时间，SIMD处理器可以同时执行多个独立的SIMD线程。SIMD线程调度器会选择一个准备好执行的SIMD线程，并启动该线程的下一条SIMD指令的执行。每个SIMD线程使用独立的寄存器组，这是一种特殊的多线程形式。

在每个执行步骤中，由于SIMD线程相互独立，SIMD线程调度器可以选择不同的SIMD线程。为了支持线程选择，调度器使用一个计分板，其中包含每个SIMD线程当前要执行的指令以及该指令的操作数是否位于寄存器中的信息。支持的最大SIMD线程数取决于计分板的大小，NVIDIA Fermi架构的典型计分板大小为32。实际的独立SIMD线程数由应用程序决定。

3.4 功能单元

SIMD处理器的功能单元数量也因GPU型号而异。每个功能单元包含一个整数单元和一个浮点单元。以具有16个物理功能单元（FU）的单个SIMD处理器为例，总共提供32768个32位寄存器，每个功能单元拥有2048个物理寄存器。相邻的两个物理寄存器可以组合成一个64位寄存器。在当前程序执行中，一个功能单元的2048个寄存器会分配给可用的SIMD线程。例如，当有32个SIMD线程时，每个线程可以使用每个功能单元的64个独立寄存器来存储数据。寄存器在SIMD线程创建时动态分配。

每个SIMD处理器都有一个本地内存（片上），为运行在该处理器上的SIMD线程提供快速访问。全局内存（片外）由所有SIMD处理器共享，CPU也可以访问，但访问全局内存的速度比访问SIMD处理器的本地内存慢得多。

3.5 不同架构对比

以下是几种不同架构的NVIDIA GPU的重要特性对比：
| GPU型号 | 架构 | 晶体管数量 | SIMD处理器数量 | 每个SIMD处理器的SIMD核心数 | 总SIMD核心数 | L2缓存 | 性能（GFLOPS） | 内存带宽 | 内存时钟频率 | 功耗 |
| ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- |
| GTX 285 | Tesla | 1.4×10⁹ | 30 | 8 | 240 | / | 1063 GF | 159 GB/sec | 2484 MHz | 204 W |
| GTX 480 | Fermi | 3.2×10⁹ | 15 | 32 | 480 | 768 KB | 1344 GF | 177 GB/sec | 3696 MHz | 250 W |
| GTX 580 | Fermi | 3.0×10⁹ | 16 | 32 | 512 | 768 KB | 1581 GF | 192 GB/sec | 4008 MHz | 244 W |
| GTX 680 | Kepler | 3.54×10⁹ | 8 | 192 | 1536 | 512 KB | 3090 GF | 192 GB/sec | 6008 MHz | 195 W |

从表中可以看出，随着架构的发展，SIMD核心数量和计算能力（以GFLOPS为单位）都有显著增加。

3.6 架构设计

为了实现高效执行，程序的编写需要充分利用GPU架构的并行单元。特别是，应用程序需要提供足够数量的SIMD线程，以便每个SIMD处理器的线程调度器有足够的线程可供切换，从而隐藏内存访问的延迟。此外，合理组织应用程序数据的内存布局也非常重要。这可以通过在CUDA或OpenCL中使用适当的命令来实现，并可以借助特定的并行编程技术来支持。

4. CUDA编程入门

4.1 CUDA程序结构

CUDA程序由在主机（传统中央处理器CPU）上执行的程序和在设备（如GPU等大规模并行处理器）上执行的程序组成。主机程序是可以用标准C编译器编译的C程序，而设备代码则是用C语言和CUDA特定扩展编写的，用于指定数据并行执行。要在设备上执行的数据并行函数称为内核函数。

CUDA程序的并行性体现在内核函数中，这些函数通常会生成大量的CUDA线程。与CPU线程相比，CUDA线程是轻量级线程，生成和调度所需的周期较少。

4.2 编译过程

包含主机程序和内核函数的CUDA程序文件（扩展名为*.cu）由NVIDIA C编译器（nvcc）进行编译。编译器会将程序的两部分分开处理：
- 内核函数被翻译成PTX（并行线程执行）汇编代码，PTX是NVIDIA的汇编语言，类似于X86汇编语言，可确保不同代NVIDIA GPU的兼容性。
- 主机程序中的内核函数调用被翻译成CUDA运行时系统调用，用于在GPU上启动相应的函数，主机程序则由标准C编译器进行翻译。

4.3 程序执行

CUDA程序的执行从主机程序开始，主机程序调用内核函数，在内核函数调用时会启动一组CUDA线程，这些线程组成一个线程网格。当网格中的所有线程完成内核函数的执行后，网格终止。调用内核函数后，CPU会继续处理主机程序，可能会调用其他内核函数。

CUDA扩展了C函数声明语法，以区分主机函数和内核函数：
- __global__ ：表示该函数是可以从主机函数调用并在GPU上执行的内核函数。
- __device__ ：表示该函数是可以从另一个内核或设备函数调用的内核函数，但这些设备函数不允许递归调用或通过指针进行间接函数调用。
- __host__ ：表示该函数是在主机上执行的传统C函数，只能由另一个主机函数调用。默认情况下，没有任何关键字的函数是主机函数。

4.4 数据传输

为了在设备上执行内核函数，数据必须存储在设备内存中。因此，CUDA程序通常包含从主机到GPU内存以及从GPU到CPU内存的数据传输操作，这些操作使用CUDA特定的函数进行显式编码。

4.4.1 内存分配

在进行数据传输之前，需要在GPU的全局内存中分配适当的内存。这可以通过调用 cudaMalloc 函数来实现：

cudaMalloc(void **, size_t);

该函数有两个参数：第一个参数是指向要分配内存的指针，第二个参数指定分配内存的字节大小。

4.4.2 内存释放

在计算完成后，可以使用 cudaFree 函数释放指定数据对象的存储空间：

cudaFree(void *);

4.4.3 数据复制

数据从主机到GPU或从GPU到主机的传输可以通过调用 cudaMemcpy 函数来实现：

cudaMemcpy(void *, const void *, size_t, enum cudaMemcpyKind);

该函数有四个参数：
- 指向传输操作目标的指针。
- 指向传输操作源的指针。
- 要复制的字节数。
- 预定义的符号常量，指定传输操作中使用的内存操作类型，例如 cudaMemcpyHostToDevice 表示从主机到设备的传输， cudaMemcpyDeviceToHost 表示从设备到主机的传输。

4.5 线程组织

在调用内核函数时，需要指定执行配置，以确定为执行该内核函数而生成的线程网格的组织方式。线程网格中的线程采用两级层次结构进行组织：
- 第一级 ：每个网格由多个线程块组成，所有线程块中的线程数量相同。
- 第二级 ：每个线程块内的线程也有自己的组织方式，同一网格中的所有线程块的组织方式相同。

线程块具有二维（CUDA版本2及更早版本）或三维（CUDA版本3及更高版本）结构，每个线程块有一个唯一的坐标，通过CUDA特定的关键字 blockIdx.x 和 blockIdx.y （二维）或 blockIdx.x 、 blockIdx.y 和 blockIdx.z （三维）表示。线程块内的线程采用三维结构组织，每个线程有一个唯一的三维坐标，通过 threadIdx.x 、 threadIdx.y 和 threadIdx.z 表示。线程块中的最大线程数在CUDA版本2及更早版本中为512，在版本3及更高版本中为1024。

根据这种两级层次结构，网格中的每个线程都可以通过其所属线程块的坐标和线程块内的坐标唯一标识。这些坐标值可用于内核函数中，以区分并行线程，实现数据并行执行。

4.6 执行配置示例

为了指定执行配置，需要声明两个 dim3 类型的结构体变量， dim3 是一种整数向量类型，默认初始值为1。这些结构体参数描述了线程网格和线程块的二维或三维组织方式，并包含在内核调用语法中，用 <<< 和 >>> 包围。例如：

// 执行配置
dim3 gsize(gx, gy);
dim3 bsize(bx, by, bz);
// 调用内核函数
KernelFct <<< gsize, bsize>>>(...);

其中， gsize 指定了内核函数执行的二维网格结构，大小为 gx × gy ； bsize 指定了线程块的三维结构，大小为 bx × by × bz 。当前网格和其线程块的大小存储在CUDA特定的变量 gridDim 和 blockDim 中。

4.7 示例：向量加法

下面是一个使用CUDA实现两个整数向量加法的示例程序：

// 内核函数：向量加法
__global__ void vecadd(int *a, int *b, int *c, int N) {
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    if (tid < N) {
        c[tid] = a[tid] + b[tid];
    }
}

#include <stdio.h>
#include <cuda_runtime.h>

#define N 1000

int main() {
    int a[N], b[N], c[N];
    int *ad, *bd, *cd;

    // 在GPU上分配内存
    cudaMalloc((void**)&ad, N * sizeof(int));
    cudaMalloc((void**)&bd, N * sizeof(int));
    cudaMalloc((void**)&cd, N * sizeof(int));

    // 初始化向量a和b
    for (int i = 0; i < N; i++) {
        a[i] = i;
        b[i] = i * 2;
    }

    // 将数据从主机复制到GPU
    cudaMemcpy(ad, a, N * sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(bd, b, N * sizeof(int), cudaMemcpyHostToDevice);

    // 调用内核函数
    dim3 grid(10);
    dim3 block(16);
    vecadd<<<grid, block>>>(ad, bd, cd, N);

    // 将结果从GPU复制到主机
    cudaMemcpy(c, cd, N * sizeof(int), cudaMemcpyDeviceToHost);

    // 输出结果
    for (int i = 0; i < N; i++) {
        printf("c[%d] = %d\n", i, c[i]);
    }

    // 释放GPU上的内存
    cudaFree(ad);
    cudaFree(bd);
    cudaFree(cd);

    return 0;
}

在这个示例中，主机程序首先在CPU上声明了三个长度为 N 的数组 a 、 b 和 c ，然后使用 cudaMalloc 函数在GPU的全局内存中为对应的向量 ad 、 bd 和 cd 分配内存。接着，使用 cudaMemcpy 函数将输入向量 a 和 b 复制到GPU内存中。

内核函数 vecadd 被调用，使用执行配置 <<<10, 16>>> 生成一个包含10个线程块，每个线程块包含16个线程的线程网格。每个线程计算对应位置的向量元素之和。计算完成后，使用 cudaMemcpy 函数将结果向量 c 从GPU内存复制回CPU内存，并输出结果。最后，使用 cudaFree 函数释放GPU上的数据结构。

4.8 线程执行

内核函数调用会根据给定的执行配置生成一个线程网格，网格和线程块的大小以及标识符存储在 gridDim 、 blockDim 、 blockIdx 和 threadIdx 变量中。生成的线程可以使用这些变量在SIMD编程模型下执行内核函数。

例如，在向量加法示例中，每个线程通过计算 threadIdx.x + blockIdx.x * blockDim.x 得到一个唯一的一维线程标识符 tid ，用于访问向量中的特定元素。

综上所述，通过合理利用GPU的架构特点和CUDA编程环境，我们可以实现高效的数据并行计算。在实际应用中，开发者需要根据具体的问题和硬件平台，合理组织线程和内存，以充分发挥GPU的性能优势。

5. CUDA线程组织与执行配置深入理解

5.1 线程网格与线程块的灵活性

线程网格和线程块的组织方式非常灵活，可以根据具体的应用需求进行调整。例如，对于一维数据处理，我们可以使用一维的线程网格和线程块；对于二维或三维数据处理，则可以使用相应维度的组织方式。

当只需要一维的线程网格和线程块时，可以直接在执行配置中指定大小。例如：

// 一维线程网格和线程块的调用
KernelFct <<<8, 16>>>(...);

这与以下代码是等价的：

dim3 gsize(8, 1);
dim3 bsize(16, 1);
KernelFct <<<gsize, bsize>>>(...);

5.2 线程标识符的计算

在CUDA中，每个线程都有唯一的标识符，这些标识符对于正确地访问数据和执行计算非常重要。以二维线程网格和线程块为例，线程的全局标识符可以通过以下方式计算：

int tid_x = threadIdx.x + blockIdx.x * blockDim.x;
int tid_y = threadIdx.y + blockIdx.y * blockDim.y;

对于三维的情况，只需再加上 z 维度的计算即可：

int tid_x = threadIdx.x + blockIdx.x * blockDim.x;
int tid_y = threadIdx.y + blockIdx.y * blockDim.y;
int tid_z = threadIdx.z + blockIdx.z * blockDim.z;

这样，每个线程就可以根据自己的全局标识符访问数据数组中的特定元素。

5.3 线程执行流程

可以用以下mermaid流程图来表示CUDA线程的执行流程：

graph TD;
    A[主机程序开始] --> B[分配GPU内存];
    B --> C[复制数据到GPU];
    C --> D[调用内核函数];
    D --> E[生成线程网格];
    E --> F[线程执行内核函数];
    F --> G[线程完成执行];
    G --> H[复制结果到主机];
    H --> I[释放GPU内存];
    I --> J[主机程序结束];

6. CUDA内存管理优化

6.1 内存层次结构利用

CUDA的GPU具有不同层次的内存，包括全局内存、共享内存、寄存器等。合理利用这些内存层次可以显著提高程序的性能。
- 全局内存 ：是最常用的内存，但访问速度相对较慢。在使用全局内存时，应尽量减少不必要的内存访问，例如通过合并内存访问来提高带宽利用率。
- 共享内存 ：位于片上，访问速度比全局内存快得多。可以将频繁访问的数据存储在共享内存中，以减少全局内存的访问次数。以下是一个使用共享内存的示例：

__global__ void sharedMemoryExample(int *input, int *output, int N) {
    __shared__ int sharedData[256];
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    if (tid < N) {
        sharedData[threadIdx.x] = input[tid];
        __syncthreads(); // 同步线程，确保所有线程都已将数据写入共享内存
        output[tid] = sharedData[threadIdx.x] * 2;
    }
}

• 寄存器 ：是最快的内存，但数量有限。应尽量让线程的局部变量存储在寄存器中，避免溢出到全局内存。

6.2 内存分配与释放策略

在进行内存分配和释放时，需要注意以下几点：
- 尽量减少内存分配和释放的次数，因为这些操作会带来一定的开销。可以采用内存池的方式，预先分配一定数量的内存，在需要时直接从内存池中获取。
- 确保在不再使用内存时及时释放，避免内存泄漏。可以使用RAII（资源获取即初始化）的思想，在对象的构造函数中分配内存，在析构函数中释放内存。

6.3 内存访问模式优化

优化内存访问模式可以提高内存带宽的利用率。例如，采用合并访问的方式，让相邻的线程访问相邻的内存地址。以下是一个合并访问的示例：

__global__ void coalescedAccess(int *input, int *output, int N) {
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    if (tid < N) {
        output[tid] = input[tid];
    }
}

7. OpenCL简介

7.1 OpenCL概述

OpenCL是一种由多个工业伙伴联合开发的标准化编程模型，可用于GPU等多种并行计算设备。与CUDA不同，OpenCL具有更好的跨平台兼容性，可以在不同厂商的GPU和CPU上运行。

7.2 OpenCL与CUDA的对比

特性	CUDA	OpenCL
兼容性	主要用于NVIDIA GPU	支持多种厂商的GPU和CPU
编程难度	相对较低，针对NVIDIA GPU优化	相对较高，需要考虑不同设备的差异
生态系统	NVIDIA提供丰富的工具和库	有更广泛的社区支持，但工具和库的集成度可能不如CUDA

7.3 OpenCL编程基础

OpenCL编程的基本步骤如下：
1. 平台和设备选择 ：选择要使用的OpenCL平台和设备。
2. 上下文和命令队列创建 ：创建OpenCL上下文和命令队列，用于管理设备和任务。
3. 程序和内核创建 ：编写OpenCL内核代码，并将其编译成程序对象。
4. 内存对象创建 ：在设备上分配内存对象，用于存储数据。
5. 数据传输 ：将数据从主机传输到设备内存。
6. 内核执行 ：设置内核参数，将内核放入命令队列中执行。
7. 结果获取 ：将计算结果从设备内存传输回主机。
8. 资源释放 ：释放所有分配的资源。

以下是一个简单的OpenCL向量加法示例：

#include <CL/cl.h>
#include <stdio.h>

#define N 1000

const char *kernelSource = "__kernel void vecadd(__global const int *a, __global const int *b, __global int *c) {"
                           "    int i = get_global_id(0);"
                           "    if (i < N) {"
                           "        c[i] = a[i] + b[i];}"
                           "}";

int main() {
    cl_platform_id platform;
    cl_device_id device;
    cl_context context;
    cl_command_queue queue;
    cl_program program;
    cl_kernel kernel;
    cl_mem a_mem, b_mem, c_mem;
    int a[N], b[N], c[N];

    // 选择平台和设备
    clGetPlatformIDs(1, &platform, NULL);
    clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL);

    // 创建上下文和命令队列
    context = clCreateContext(NULL, 1, &device, NULL, NULL, NULL);
    queue = clCreateCommandQueue(context, device, 0, NULL);

    // 创建程序对象
    program = clCreateProgramWithSource(context, 1, &kernelSource, NULL, NULL);
    clBuildProgram(program, 1, &device, NULL, NULL, NULL);

    // 创建内核对象
    kernel = clCreateKernel(program, "vecadd", NULL);

    // 初始化向量a和b
    for (int i = 0; i < N; i++) {
        a[i] = i;
        b[i] = i * 2;
    }

    // 创建内存对象
    a_mem = clCreateBuffer(context, CL_MEM_READ_ONLY, N * sizeof(int), NULL, NULL);
    b_mem = clCreateBuffer(context, CL_MEM_READ_ONLY, N * sizeof(int), NULL, NULL);
    c_mem = clCreateBuffer(context, CL_MEM_WRITE_ONLY, N * sizeof(int), NULL, NULL);

    // 将数据写入内存对象
    clEnqueueWriteBuffer(queue, a_mem, CL_TRUE, 0, N * sizeof(int), a, 0, NULL, NULL);
    clEnqueueWriteBuffer(queue, b_mem, CL_TRUE, 0, N * sizeof(int), b, 0, NULL, NULL);

    // 设置内核参数
    clSetKernelArg(kernel, 0, sizeof(cl_mem), &a_mem);
    clSetKernelArg(kernel, 1, sizeof(cl_mem), &b_mem);
    clSetKernelArg(kernel, 2, sizeof(cl_mem), &c_mem);

    // 执行内核
    size_t globalSize = N;
    clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &globalSize, NULL, 0, NULL, NULL);

    // 读取结果
    clEnqueueReadBuffer(queue, c_mem, CL_TRUE, 0, N * sizeof(int), c, 0, NULL, NULL);

    // 输出结果
    for (int i = 0; i < N; i++) {
        printf("c[%d] = %d\n", i, c[i]);
    }

    // 释放资源
    clReleaseMemObject(a_mem);
    clReleaseMemObject(b_mem);
    clReleaseMemObject(c_mem);
    clReleaseKernel(kernel);
    clReleaseProgram(program);
    clReleaseCommandQueue(queue);
    clReleaseContext(context);

    return 0;
}

8. 总结与展望

8.1 总结

GPU编程为数据并行计算提供了强大的支持，CUDA和OpenCL是两种常用的编程环境。通过合理利用GPU的架构特点，如多核并行、内存层次结构等，以及掌握CUDA和OpenCL的编程技巧，如线程组织、内存管理等，可以显著提高程序的性能。

8.2 展望

随着硬件技术的不断发展，GPU的性能将不断提升，同时编程环境也将不断完善。未来，GPU编程将在更多领域得到应用，如人工智能、深度学习、科学计算等。开发者需要不断学习和掌握新的技术，以充分发挥GPU的潜力。同时，跨平台编程的需求也将越来越高，OpenCL等跨平台编程模型将得到更广泛的应用。

在实际应用中，开发者应根据具体的需求和硬件平台选择合适的编程环境和优化策略，以实现高效、稳定的并行计算。