12:57

自编教材分享：第十章—CUDA程序优化（一）

4446

0

CUDA编程简介

CUDA是什么

在一个CPU与GPU组成的异构计算系统中，CPU处理器有四个算术逻辑单元ALU，GPU由多个流式多处理器（Streaming Multiprocessor,SM)组成，一个SM相当于一个完整的多核处理单元，每个SM中包含有公用的控制单元和缓存，以及大量的算术逻辑单元。与CPU相比，在同样的芯片面积上GPU将更多的元器件用于计算，因此GPU对逻辑简单、数据量大的计算密集型任务有着天然的计算优势。

统一计算设备架构CUDA是NVIDIA提出的通用并行计算平台和编程模型，为使用GPU的异构计算开发提供了便捷高效的开发环境。异构计算采用并行或分布式计算方式，通过协调地使用性能、结构各异的计算器件以满足不同的计算需求，由CPU处理器与众核GPU可组成一个典型的异构计算架构。如图所示。

CUDA平台支持开发者使用C/C++、FORTRAN、Python等行业标准程序语言的扩展来构建CUDA程序，同时CUBLAS、Thrust等丰富的CUDA加速库也为CUDA程序的开发提供了便利，如图所示。CUDA C是标准ANSI C语言的一个扩展，被广泛应用于各领域CUDA程序的开发，本章后续范例将统一使用CUDA C进行编写。

CUDA提供了两层应用程序接口API来管理GPU设备，分别是CUDA驱动和CUDA运行时，如图所示。CUDA驱动能够细致全面的控制GPU设备的运行状态，但使用驱动API编程的难度较大。CUDA运行时作为更高级的API实现在CUDA驱动API的上层，使用运行时API能够简化管理GPU设备的操作、降低编程难度。运行时与驱动API在使用时是相互排斥的，无法实现混合调用，且合理利用运行时API或CUDA驱动API 都能构建高效的CUDA程序，本章后续范例都将使用CUDA运行时API来实现对GPU设备的管理。

CUDA程序用NVCC编译器进行编译，其编译流程自上而下如图所示，NVCC编译器在编译过程中会将主机代码与设备代码进行分离，经过代码分离后，使用C语言编写的主机端代码将由本地C语言编译器进行编译，使用CUDA C语言编写的设备端代码会通过NVCC编译器进行编译。

CUDA编程模型

CUDA异构编程模型将异构系统分为主机端(host)与设备端(device)，主机端对应于CPU，设备端对应于GPU。如图所示，CUDA程序主机端代码运行在CPU上，设备端代码运行在GPU上，在设备端执行的函数被称为核函数（Kernel）。

一个典型的CUDA程序实现流程如下：

（1）获取GPU设备

（2）开辟GPU上显存空间

（3）发起主机向设备的数据传输

（4）启动核函数

（5）发起设备向主机的数据传输

（6）释放GPU的显存空间，重置设备

与实现流程对应的CUDA程序主要代码如下。在编写CUDA程序时，通过调用CUDA运行时的cudaMalloc、cudaFree等函数能够显式地控制GPU设备进行内存开辟与内存释放；通过调用cudaMemcpy函数能够控制CUDA程序中主机端与设备端的数据传输；使用语句kernel_name <<<grid,block>>>能够实现对核函数的调用；通过调用cudaDeviceResset函数能够对GPU设备进行重置。

cudaSetDevice(0);
cudaMalloc((void**) &d_a, sizeof(float) * n);
cudaMemcpy(d_a, a, size_t count, cudaMemcpyHostToDevice);
kernel<<<blocks,threads>>>;
cudaMemcpy(a, d_a, size_t count, cudaMemcpyDeviceToHost);
cudaFree(d_a);
cudaDeviceReset();

下面概要地对一些常用CUDA运行时函数进行介绍，常用的设备管理类函数有：

cudaError_t cudaGetDeviceCount( int* count )//用来获取当前系统中可用GPU设备的数量
cudaSetDevice (int *device)//用来在系统中选择希望调用的GPU设备
cudaDeviceReset(void)//用来显式销毁和清理当前GPU设备上的所有资源
cudaDeviceSynchronize(void)//用来显式地阻塞主机端进程直至系统中的GPU设备完成其上的计算任务

常用的内存管理类函数有：

cudaMalloc (void** devptr, size_t size)//用来在GPU设备上分配一定字节的线性内存，并以devPtr的形式返回指向所分配内存的指针
cudaMemcpy (void* dst, const void* src, size_t count, cudaMemcpyKind kind )//用来实现主机端与设备端数据传输
cudaMallocPitch(void** devPtr, size_t* pitch, size_t width, size_t height)//用来在GPU设备上分配线性内存，并以*devPtr的形式返回指向所分配内存的指针
cudaMemcpy2D(void* dst, size_t* pitch, const void* src, size_t width, size_t height, cudaMemcpyKind kind )//针对cudaMallocPitch在GPU设备上分配内存后，进行主机端到设备端的数据传输
cudaFree(void *devPtr)//释放devPtr指向的由cudaMalloc()调用开辟的GPU上内存空间

设备核函数：

__global__ void kernel_name (argument list)//设备端执行的代码称为核函数，在程序中使用__global__声明定义，函数返回类型必须为void类型

错误处理函数：

const char* cudaGetErrorString (cudaError_t error)//将CUDA程序运行时产生的错误信息error进行转化为可读的错误信息

CUDA程序编写

当核函数在主机端启动时，GPU设备中会产生大量的线程（thread），一定量的线程组成线程块（block），一个核函数启动产生的所有线程统称为一个网格（grid），它由多个相同的线程块构成。CUDA运行时为网格内的每个线程分配了内置坐标变量threadIdx和blockIdx，threadIdx表示线程在线程块内的索引，blockIdx表示线程块在网格内的索引，在编写核函数时使用这两个坐标变量可以将不同线程区分开来，从而控制不同线程完成指定的操作。

下面将以CUDA向量相加为例完整展示CUDA程序的编写过程，一个在主机端执行的向量相加函数的代码如下所示。

void sumArraysOnHost(float *A, float *B, float *C, const int N){
    for (int idx = 0; idx < N; idx++){
        C[idx] = A[idx] + B[idx];
    }
}

__global__限定符表示该函数是在设备上执行的核函数，核函数sumArraysOnGPU通过开启N个线程实现了向量A、B内N个元素相加的并行实现，消除了函数sumArraysOnHost中的循环体，并使用线程坐标变量threadIdx替换了数组索引。该函数将两个大小为N的向量A和B相加，通过N次循环实现计算操作，该函数对应的CUDA核函数代码如下。

__global__ void sumArraysOnGPU(float *A, float *B, float *C, const int N){
    int tx = threadIdx.x;
    C[tx] = A[tx] + B[tx];
}

CUDA向量相加的实现流程如下：

（1）使用cudaGetDeviceProperties 获取GPU设备

（2）完成数组A和数组B的初始化

（3）使用cudaMalloc函数开辟用于存储数组A、B、C元素的GPU内存空间d_A、d_B和d_C

（4）使用cudaMemcpy控制CPU端向GPU端的传输数组A、B的元素

（5）启动核函数sumArraysOnGPU在GPU上进行数组相加运算

（6）使用cudaMemcpy控制GPU端向CPU端传输结果数组C的元素

（7）验证CUDA数组相加的正确性

（8）释放GPU的显存空间，重置设备

CUDA版矩阵乘

矩阵乘法作为一种基本的数学运算，在科学计算领域有着非常广泛的应用，矩阵乘法的快速算法对科学计算有着极为重要的意义。一般矩阵乘法的实现思想如图所示，阵A中的一行和矩阵B中的一列进行向量内积得到矩阵C中的一个元素。

使用C语言对一般矩阵乘法进行实现，其中A、B矩阵与结果矩阵C均为长宽是width的方阵。实现代码中的最内层循环依次取A矩阵的一行元素B矩阵的一列元素，进行累乘加运算从而得到C矩阵的一个元素，通过外层循环的迭代完成对结果矩阵C中所有元素的计算。部分核心代码如下：

void MatrixMulOnHost(float *A, float *B, float *C, int width){
    for(int i=0; i<width; i++){
        for(int j=0; j<width; j++){
            float sum = 0.0;
            for(int k=0; k<width; k++){
                float a = A[i*width+k];
                float b = B[k*width+j];
                sum+=a * b;
            }
            C[i*width+j]=sum;
        }
    }
}

通过上面对CUDA编程的介绍，对于一个结果矩阵规模为width*width的一般矩阵乘法，可以启用width*width个线程，每个线程负责计算结果矩阵中的一个元素。其中Ad、Bd和Cd是在GPU上显式开辟的内存空间，用于存储矩阵A、矩阵B和结果矩阵C中的元素。在核函数中，使用内置坐标变量threadIdx.x确定各线程负责计算矩阵C上位置坐标为（row，col）的元素；核心计算代码sum += Ad[row*width+i]*Bd[i*width+col]表示每个线程依次取A矩阵的第row行与B矩阵第col列上的元素进行乘加操作，计算出结果矩阵C上的一个元素，相较于一般矩阵乘法C语言代码，CUDA矩阵乘法实现的核函数中没有了一般矩阵乘函数中最外的两层循环体，通过多线程的并行操作，从而提升矩阵乘的运算效率。核函数代码如下：

__global__ void MatrixMulKernel(float* Ad, float* Bd, float* Cd, int width){
	int offset = threadIdx.x;
	int row = offset /width;
	int col = offset %(width-1);
	float sum = 0;
	for(int i=0;i<width;i++){
	    sum += Ad[row*width+i]*Bd[i*width+col];   
	}
	Cd[row*width+col] = sum;
}

使用NVCC对CUDA矩阵乘的代码进行编译并运行测试，编译命令为：nvcc matrixmul.cu -o matrixmul，测试环境中GPU设备为NVIDIA RTX 3090，CUDA版本为11.6，使用性能分析工具Nsight System对核函数进行计时用于性能监测，执行命令为：nsys profile --stats=true ./matrixmul。矩阵A、矩阵B以及结果矩阵C是数据规模为32*32的方阵，矩阵内有1024个元素，CUDA矩阵乘核函数MatrixMulKernel在GPU设备上启动了1个线程块，线程块线程数目为1024个。测试结果表明在该数据规模下，MatrixMulKernel核函数的执行时间仅为5.18us，使用GPU设备的CUDA矩阵乘相较于使用CPU设备的一般矩阵乘性能取得了大幅提升。测试结果如表所示：

通过本节内容，优化人员对CUDA的基本概念和编程方法能够形成初步了解，本章后续内容将对CUDA程序的优化方法展开描述，测试使用的GPU设备均为NVIDIA RTX 3090，并主要通过Nsight System记录核函数的运行时间来评估CUDA程序的性能。