笔记02：CUDA编程模型

m0_46521579

已于 2023-08-28 22:57:55 修改

阅读量446

点赞数

分类专栏： CUDA 文章标签：笔记

于 2023-08-01 17:39:29 首次发布

本文链接：https://blog.csdn.net/m0_46521579/article/details/132022587

版权

CUDA 专栏收录该内容

7 篇文章 0 订阅

订阅专栏

CUDA是一种通用的并行计算平台和编程模型，是在C语言基础上扩展的。

一、CUDA编程模型概述

1. CUDA编程结构

在一个异构环境中包含多个CPU和GPU，每个GPU和CPU的内存都由一条PCI-e总线分隔开，需要注意区分

（1）主机：CPU及其内存（主机内存）

（2）设备：GPU及其内存（设备内存）

从CUDA6.0开始，NVIDIA提出了名为“统一寻址”的编程模型的改进，连接了主机内存和设备内存空间，可使用单个指针访问CPU和GPU内存，无须彼此之间手动拷贝数据。重要的是如何为主机和设备分配内存空间以及如何在CPU和GPU之间拷贝共享数据。

内核（kernel）是CUDA编程模型的一个重要组成部分，其代码在GPU上运行。

多数情况下，主机可以独立地对设备进行操作。内核一旦被启动，管理权立刻返回给主机，释放CPU来执行由设备上运行的并行代码实现的额外的任务。

CUDA编程模型主要是异步的，因此在GPU上进行的运算可以与主机-设备通信重叠。一个典型的CUDA程序包括由并行代码互补的串行代码。串行代码在主机CPU上执行，而并行代码在GPU上执行。主机代码按照ANSI C标准进行编写，而设备代码使用CUDA C进行编程。可以将所有的代码统一放在一个源文件中，也可以使用多个源文件来构建应用程序和库。NVIDIA的C编译器（nvcc）为主机和设备生成可执行代码。

一个典型的CUDA程序实现流程遵循以下模式：

（1）把数据从CPU内存拷贝到GPU内存

（2）调用核函数对存储在GPU内存中的数据进行操作

（3）将数据从GPU内存传送回到CPU内存

2. 内存管理

CUDA编程模型假设系统是由一个主机和一个设备组成的，并且各自拥有独立的内存。核函数是在设备上运行的。为了拥有充分的控制权并使系统达到最佳性能，CUDA运行时负责分配与释放设备内存，并且在主机内存和设备内存之间传输数据。

cudaError_t cudaMalloc(void** devPtr, size_t size)

该函数负责向设备分配一定字节的线性内存，并以devPtr的形式返回指向所分配内存的指针。cudaMalloc与标准C语言中的malloc函数几乎一样，只是此函数在GPU的内存里分配内存。

cudaError_t cudaMemcpy(void* dst, const void* src, size_t count, cudaMemcpyKind kind)

此函数负责主机和设备之间的数据传输，从src指向的源存储区复制一定数量的字节到dst指向的目标存储区。复制方向由kind指定。

3. 线程管理

当核函数在主机端启动时，它的执行会移动到设备上，此时设备中会产生大量的线程并且每个线程都执行由核函数指定的语句。了解如何组织线程是CUDA编程的一个关键部分。

下图中：dim3 grid(3, 2, 1), block(5, 3, 1);

由一个内核启动所产生的所有线程统称为一个网格。同一个网格中的所有线程共享相同的全局内存空间。一个网格由多个线程块构成，一个线程块包含一组线程，同一线程块内的线程协作可以通过以下方式来实现：（1）同步；（2）共享内存。不同块内的线程不能协作。

线程依靠以下两个坐标变量来区分彼此：

blockIdx（线程块在线程网格内的索引）

threadIdx（线程块内的线程索引）

以上变量是核函数中需要初始化的内置变量。当执行一个核函数时，CUDA运行时为每个线程分配坐标变量blockIdx和threadIdx。基于这些坐标，可以将部分数据分配给不同的线程。

CUDA可以组织三维的网格和块，网格和块的维度由下列两个内置变量指定：

（1）blockDim（线程块的维度，用每个线程块中的线程数来表示）

（2）gridDim（线程格的维度，用每个线程格中的线程数来表示）

它们是dim3类型的变量，是基于uint3定义的整数型向量，用来表示维度。

当定义一个dim3类型的变量时，所有未指定的元素都被初始化为1。

dim3类型变量中的每个组件可以通过它的x、y、z字段获得。

#include <stdio.h>


__global__ void checkIndex(void)
{
        printf("threadIdx: (%d, %d, %d); "
                "blockIdx:  (%d, %d, %d); "
                "blockDim:  (%d, %d, %d); "
                "gridDim:   (%d, %d, %d)\n",
                threadIdx.x, threadIdx.y, threadIdx.z,
                blockIdx.x, blockIdx.y, blockIdx.z,
                blockDim.x, blockDim.y, blockDim.z,
                gridDim.x, gridDim.y, gridDim.z);
}


int main(void)
{
        int nElem = 6;
        dim3 block(3);
        dim3 grid((nElem+block.x-1)/block.x);

        printf("block: %d, %d, %d\n", block.x, block.y, block.z);
        printf("grid: %d, %d, %d\n", grid.x, grid.y, grid.z);

        checkIndex<<<grid, block>>>();

        cudaDeviceReset();

        return 0;
}

对于一个给定的数据大小，确定网格和块尺寸的一般步骤为：

（1）确定块的大小

（2）在已知数据大小和块大小的基础上计算网格维度

要确定块尺寸，通常需要考虑：（1）内核性能特性；（2）GPU资源的限制

4. CUDA内核函数

kernel_name <<<grid, block>>> (argument list);

核函数的调用与主机线程是异步的。核函数调用结束后，控制权立刻返回给主机端。可以调用以下函数强制主机端程序等待所有的核函数执行结束。

cudaError_t cudaDeviceSynchronize(void);

__global__函数必须具有void返回类型，并且不能是类的成员。

__device__：在设备上执行，只能从设备调用。__global__和__device__执行空间说明符不能一起使用。

__host__：在主机上执行，只能从主机调用，__global__和__host__执行空间说明符不能一起使用。但是，__device__和__host__执行空间说明符可以一起使用，在这种情况下，该函数是为主机和设备编译的。

核函数的相关限制：

（1）只能访问设备内存

（2）必须具有void返回类型

（3）不支持可变数量的参数

（4）不支持静态变量

（5）显示异步行为

二、给核函数计时

1. 用CPU计时器计时

可以使用gettimeofday系统调用来创建一个CPU计时器，以获取系统的时钟时间，其将返回1970年1月1日零点以来，到现在的秒数。程序中需要添加sys/time.h头文件。

#include <sys/time.h>

double cpuSecond(void)
{
    struct timeval tp;
    gettimeofday(&tp, NULL);
    
    return ((double)tp.tv_sec + (double)tp.tv_usec*1.e-6);
}

用cpuSecond函数来测试核函数

double iStart = cpuSecond();
kernel_name <<< grid, block >>> (argument list);
cudaDeviceSynchronize();
double iElaps = cpuSecond() - iStart;

由于核函数调用与主机端程序是异步的，需要用cudaDeviceSynchronize函数来等待所有的GPU线程运行结束。

#include <stdio.h>
#include <cuda_runtime.h>
#include <sys/time.h>

void initialData(float *ip, int size)
{
        time_t t;
        srand((unsigned) time(&t));

        for(int i=0; i<size; i++)
        {
                ip[i] = (float)(rand()&0xFF)/10.0f;
        }
}

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

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

void checkResult(float *hostRef, float *gpuRef, const int N)
{
        double epsilon = 1.0E-8;
        bool match = 1;
        for(int i=0; i<N; i++)
        {
                if(abs(hostRef[i]-gpuRef[i]) > epsilon)
                {
                        match = 0;
                        printf("Arrays do not match!\n");
                        printf("host %5.2f gpu %5.2f at count %d\n", hostRef[i], gpuRef[i], i);
                        break;
                }
        }
        if(match) printf("Arrays match.\n\n");
}

double cpuSecond(void)
{
        struct timeval tp;
        gettimeofday(&tp, NULL);

        return ((double)tp.tv_sec + (double)tp.tv_usec*1.e-6);
}
int main(int argc, char **argv)
{
        double iStart, iElaps;

        printf("%s Starting ...\n", argv[0]);

        // set up device
        int dev = 0;
        cudaDeviceProp deviceProp;
        cudaGetDeviceProperties(&deviceProp, dev);
        printf("Using Device %d : %s\n", dev, deviceProp.name);
        cudaSetDevice(dev);

        // set up data size of vectors
        int nElem = 1<<24;
        printf("Vector size %d\n", nElem);

        // malloc host memory
        size_t nBytes = nElem * sizeof(float);

        float *h_A, *h_B, *hostRef, *gpuRef;
        h_A = (float *)malloc(nBytes);
        h_B = (float *)malloc(nBytes);
        hostRef = (float *)malloc(nBytes);
        gpuRef = (float *)malloc(nBytes);

        // initialize data at host side
        initialData(h_A, nElem);
        initialData(h_B, nElem);

        memset(hostRef, 0, nBytes);
        memset(gpuRef, 0, nBytes);

        // malloc device global memory
        float *d_A, *d_B, *d_C;
        cudaMalloc((float **)&d_A, nBytes);
        cudaMalloc((float **)&d_B, nBytes);
        cudaMalloc((float **)&d_C, nBytes);

        // transfer data from host to device
        cudaMemcpy(d_A, h_A, nBytes, cudaMemcpyHostToDevice);
        cudaMemcpy(d_B, h_B, nBytes, cudaMemcpyHostToDevice);

        // invoke kernel at host side
        int iLen = 1024;
        dim3 block(iLen);
        dim3 grid((nElem+block.x-1)/block.x);
        printf("<<< grid, block >>> = %d, %d \n", grid.x, block.x);


        iStart = cpuSecond();
        for(int i=0; i<100; i++)
        {
                sumArraysOnGPU<<< grid, block>>>(d_A, d_B, d_C, nElem);
                //printf("<<<grid, block>>> = %d, %d\n", grid.x, block.x);
                cudaDeviceSynchronize();
        }
        iElaps = cpuSecond() - iStart;
        printf("GPU : %f sec \n", iElaps);

        // copy kernel result back to host side
        cudaMemcpy(gpuRef, d_C, nBytes, cudaMemcpyDeviceToHost);

        // add vector at host side for result checks
        iStart = cpuSecond();
        for(int i=0; i<100; i++)
        {
                sumArraysOnHost(h_A, h_B, hostRef, nElem);
        }
        iElaps = cpuSecond() - iStart;
        printf("CPU: %f sec\n", iElaps);

        // check device results
        checkResult(hostRef, gpuRef, nElem);

        // free device global memory
        cudaFree(d_A);
        cudaFree(d_B);
        cudaFree(d_C);

        // free host memeory
        free(h_A);
        free(h_B);
        free(hostRef);
        free(gpuRef);

        return 0;
}

2. 用nvprof工具计时

自CUDA5.0以来，NVIDIA提供了一个名为nvprof的命令行分析工具，可以帮助从应用程序的CPU和GPU活动情况中获取时间线信息，其包括内核执行、内存传输以及CUDA API的调用。

nvprof [nvprof_args] <application> [application_args]

三、组织并行线程

Thread -> Thread Block -> Thread Grid

使用合适的网格和块大小来组织线程，可以对内核性能产生很大的影响。

对于矩阵运算，传统的方法是在内核中使用一个包含二维网格与二维块的布局来组织线程。

在矩阵加法中使用以下布局将有助于了解更多关于网格和块的启发性的用法：

（1）由二维线程块构成的二维网格

（2）由一维线程块构成的一维网格

（3）由一维线程块构成的二维网格

1. 使用块和线程建立矩阵索引

对于一个给定的线程，首先可以通过把线程和块索引映射到矩阵坐标上来获取线程块和线程索引的全局内存偏移量，然后将这些矩阵坐标映射到全局内存的存储单元中。

第一步，把线程和块索引映射到矩阵坐标上：

ix = threadIdx.x + blockIdx.x * blockDim.x
iy = threadIdx.y + blockIdx.y * blockDim.y

第二步，把矩阵坐标映射到全局内存中的索引/存储单元上：

idx = iy * nx + ix

#include <stdio.h>

void initialInt(int *ip, int size)
{
        for(int i=0; i<size; i++)
        {
                ip[i] = i;
        }
}

void printMatrix(int *C, const int nx, const int ny)
{
        int *ic = C;

        printf("\n Matrix: (%d, %d) \n", nx, ny);
        for(int iy=0; iy<ny; iy++)
        {
                for(int ix=0; ix<nx; ix++)
                {
                        printf("%3d ", ic[ix]);
                }
                ic += nx;
        }
        printf("\n");
}

__global__ void printThreadIndex(int *A, const int nx, const int ny)
{
        int ix = threadIdx.x + blockIdx.x * blockDim.x;
        int iy = threadIdx.y + blockIdx.y * blockDim.y;
        unsigned int idx = iy*nx + ix;

        printf("thread_id (%d, %d)"
                "block_id (%d, %d)"
                "coordinate (%d, %d)"
                "global index %2d ival %2d \n",
                threadIdx.x, threadIdx.y, blockIdx.x,
                blockIdx.y, ix, iy, idx, A[idx]);
}

int main(int argc, char **argv)
{
        printf("%s Starting ...\n", argv[0]);

        // get device information
        int dev = 0;
        cudaDeviceProp deviceProp;
        cudaGetDeviceProperties(&deviceProp, dev);
        printf("Using Device %d : %s\n", dev, deviceProp.name);
        cudaSetDevice(dev);

        // set matrix dimension
        int nx = 8;
        int ny = 6;
        int nxy = nx*ny;
        int nBytes = nxy * sizeof(int);

        int *h_A;
        h_A = (int *)malloc(nBytes);

        // initialize host matrix with integer
        initialInt(h_A, nxy);
        printMatrix(h_A, nx, ny);

        // malloc device memory
        int *d_MatA;
        cudaMalloc((void **)&d_MatA, nBytes);

        // transfer data from host to device
        cudaMemcpy(d_MatA, h_A, nBytes, cudaMemcpyHostToDevice);

        // set up execution configuration
        dim3 block(4, 2);
        dim3 grid((nx+block.x-1)/block.x, (ny+block.y-1)/block.y);

        // invoke the kernel
        printThreadIndex <<< grid, block >>> (d_MatA, nx, ny);
        cudaDeviceSynchronize();

        //free host and Device memory
        cudaFree(d_MatA);
        free(h_A);

        // reset device
        cudaDeviceReset();

        return 0;
}

2. 使用二维网格和二维块对矩阵求和

使用二维网格和二维块来实现一个矩阵加法核函数。

#include <stdio.h>
#include <sys/time.h>

void initialData(float* h_A, int nx, int ny)
{
        float* pointer = h_A;
        for(int i=0; i<ny; i++)
        {
                for(int j=0; j<nx; j++)
                {
                        *pointer = (float)(rand()&0xFF)/10.0f;
                        pointer ++;
                        //printf("%d, %d = %f \n", i, j, *(h_A + i*nx + j));
                }
        }
}

void sumMatrixOnCPU(float* h_A, float* h_B, float* hostRef, int nx, int ny)
{
        int offset = 0;
        for(int i=0; i<ny; i++)
        {
                for(int j=0; j<nx; j++)
                {
                        hostRef[offset] = h_A[offset] + h_B[offset];
                        offset ++;
                }
        }
}

__global__ void sumMatrixOnGPU(float* d_A, float* d_B, float* d_C, const int nx, const int N)
{
        int ix = threadIdx.x + blockIdx.x * blockDim.x;
        int iy = threadIdx.y + blockIdx.y * blockDim.y;
        int idx = iy * nx + ix;
        if(idx < N) d_C[idx] = d_A[idx] + d_B[idx];
}

void checkResult(float* hostRef, float* gpuRef, int nElem)
{
        int offset = 0;
        for(int i=0; i<nElem; i++)
        {
                if((*(hostRef+offset)) != (*(gpuRef+offset)))
                {
                        printf("Error \n");
                        break;
                }
                offset ++;
        }
}

double cpuSecond()
{
        struct timeval tp;
        gettimeofday(&tp, NULL);

        return ((double)tp.tv_sec + (double)tp.tv_usec*1.e-6);
}

int main()
{
        double iStart, iElaps;
        // set up device
        int dev = 0;
        cudaDeviceProp deviceProp;
        cudaGetDeviceProperties(&deviceProp, dev);
        printf("Using Device %d : %s \n", dev, deviceProp.name);
        cudaSetDevice(dev);

        // set up data size of Matrix
        int nx = 1<<14;
        int ny = 1<<14;
        int nElem = nx * ny;
        printf("Matrix size: %d * %d \n", ny, nx);

        // malloc host memory
        int nBytes = nElem * sizeof(float);

        float *h_A, *h_B;
        h_A = (float *)malloc(nBytes);
        h_B = (float *)malloc(nBytes);

        // initialize data at host side
        initialData(h_A, nx, ny);
        initialData(h_B, nx, ny);

        // sum on CPU
        float *hostRef;
        hostRef = (float *)malloc(nBytes);
        memset(hostRef, 0, nBytes);

        iStart = cpuSecond();
        sumMatrixOnCPU(h_A, h_B, hostRef, nx, ny);
        iElaps = cpuSecond() - iStart;
        printf("CPU time: %f sec \n", iElaps);

        // sum on GPU
        dim3 block(16, 16);
        dim3 grid((nx+block.x-1)/block.x, (ny+block.y-1)/block.y);
        printf("block: %d, %d, %d \n", block.x, block.y, block.z);
        printf("grid : %d, %d, %d \n", grid.x,  grid.y,  grid.z);

        float *d_A, *d_B, *d_C;
        cudaMalloc((float **)&d_A, nBytes);
        cudaMalloc((float **)&d_B, nBytes);
        cudaMalloc((float **)&d_C, nBytes);
        float *gpuRef;
        gpuRef = (float *)malloc(nBytes);

        cudaMemcpy(d_A, h_A, nBytes, cudaMemcpyHostToDevice);
        cudaMemcpy(d_B, h_B, nBytes, cudaMemcpyHostToDevice);

        iStart = cpuSecond();
        sumMatrixOnGPU <<< grid, block >>> (d_A, d_B, d_C, nx, nElem);
        cudaDeviceSynchronize();
        iElaps = cpuSecond() - iStart;
        printf("GPU time: %f sec \n", iElaps);

        cudaMemcpy(gpuRef, d_C, nBytes, cudaMemcpyDeviceToHost);

        // check device results
        checkResult(hostRef, gpuRef, nElem);


        free(h_A);
        free(h_B);
        free(hostRef);
        free(gpuRef);

        cudaFree(d_A);
        cudaFree(d_B);
        cudaFree(d_C);

        return 0;
}

四、设备管理

1. 使用运行时API查询GPU信息

在CUDA运行时API中有很多函数可以帮助管理这些设备，可以使用以下函数查询关于GPU设备的所有信息：

cudaError_t    cudaGetDeviceProperties(cudaDeviceProp* prop, int device);

name	NVIDIA GeForce MX150
totalGlobalMem	2098593792 bytes
clockRate	1531500 Hz
memoryClockRate	3004000 Hz
memoryBusWidth	64-bit
l2CacheSize	524288 bytes
maxTexture1D	131072
maxTexture2D	131072 * 65536
maxTexture3D	16384 * 16384 * 16384
maxTexture1DLayered	(32768) * 2048
maxTexture2DLayered	(32768 * 32768) * 2048
totalConstMem	65536 bytes
sharedMemPerBlock	49152 bytes
regsPerBlock	65536
warpSize	32
maxThreadsPerMultiProcessor	2048
maxThreadsPerBlock	1024
maxThreadsDim	1024 * 1024 * 64
maxGridSize	2147483647 * 65535 * 65535
memPitch	2147483647 bytes
multiProcessorCount	3