济沧浪
2023/12/17
目录

CUDA 入门

# 简介

  • CPU （Central Processing Unit）被设计为擅长执行通用逻辑 —— 采用数据缓存和流程控制提升效率 —— 负责执行控制
  • GPU （Graphic Processing Unit）被设计为擅长高度并行计算 —— 采用大量的线程核提升并行效率 —— 负责数据处理

早期 GPU 的功能仅仅是处理图像，所以基本是由 GPU 厂商负责实现类如 OpenGL 之类的图形库，以供开发者使用。

—— 核心是图像处理

但随着时代发展，人们发现了 GPU 在处理大数据上的潜力，所以出现了 CUDA、OpenCL、Metal等调用 GPU 进行数据处理加速的 API。

—— 核心是数据处理

CUDA:

NVIDIA 在 2006 年 11 月推出的利用 NVIDIA GPU 中的并行计算引擎的通用并行计算平台和编程模型。

CUDA 附带一个 C++ 软件环境，需要 NVIDIA 显卡硬件、 NVIDIA 显卡驱动 和 CUDA 开发工具包。

CUDA 并行编程主要有三个抽象：线程组层次结构、共享内存、屏障同步。

相关链接：

CUDA-C-开发手册 (opens new window)

CUDA开发者下载 (opens new window)

同时 CUDA 也有相应的官方编程库：

# Hello World

#include <cstdio>

__global__ void cuda_hello()
{
    printf("Hello World!");
}

int main()
{
    cuda_hello<<<1,1>>>();
    cudaDeviceSynchronize();
    return 0;
}

nvcc -O3 hello_world.cu -o hello_world
./hello_world

# 编译

CUDA 编译器 —— nvcc。

CUDA 应用程序的源文件由传统 C++ 主机代码和 GPU 设备函数混合组成。所以 CUDA 编译时会区分这两部分，使用专有的 NVIDIA 编译器和汇编器编译 GPU 设备函数，使用可用的 C++ 主机编译器编译传统主机代码，最后再链接 CUDA 运行时库组成完整的程序。

其中，编译 GPU 设备函数时通常会经历两大步：

  • 编译器将设备代码编译为汇编形式（PTX代码）—— 离线编译
  • 设备驱动程序将 PTX 代码进一步编译为二进制代码（cubin）—— 即时编译

即时编译：

  • 增加了应用程序的初次加载时间（设备驱动程序会自动缓存生成的二进制代码的副本，且当驱动程序升级时缓存会失效。）。

  • 允许应用程序从每个新设备驱动程序附带的任何新编译器优化中受益。（编译时不存在，未来可能新增的优化）

  • 同时还可以使程序运行在非编译设备的环境里。

例：

nvcc x.cu
 -gencode arch=compute_50,code=sm_50
 -gencode arch=compute_60,code=sm_60
 -gencode arch=compute_70,code=\"compute_70,sm_70\"

  • arch ：指定 NVIDIA GPU 架构（虚拟）

  • code ：指定生成的代码

    • compute_ ：表示生成的 ptx 版本
    • sm_ ：表示生成的 cubin 版本

当然，也可以直接编译，不使用此功能优化。正如上一节 hello_world 程序那样。

# 基础概念

  • kernel 函数

    调用时将任务指派给 GPU 完成。（如 Hello World 中的 cuda_hello 函数）。

    kernel_function<<<grid_size, block_size, shared_size, stream>>>();

    grid_size 指线程块网格的大小；

    block_size 指线程块的大小；

    shared_size 指除了静态分配的共享内存之外，该核函数动态分配的共享内存的字节数；这个动态分配的内存由任何声明为外部数组的变量使用

    stream 指同步流；

  • kernel 函数执行说明符

    • _global_：在 GPU 设备上执行且可由主机调用。

    • _device_：在 GPU 设备上执行且只能由 GPU 设备调用。

    • _host_：在主机上执行且只能从主机上调用。

      ……

  • 内存空间说明符

    • _device_：驻留在 GPU 设备上。每个 GPU 设备有一个不同的对象

    • _constant_：驻留在常量内存空间。每个 GPU 设备有一个不同的对象

    • _shared_：驻留线程块的共享内存空间。每个线程块有一个不同的对象且只能从块内线程访问。

      ……

  • 线程层次结构

    线程块 ——> 线程块簇（可选） ——> 线程块网格

注：目前一个线程块最多 1024 个线程。

在内置变量 threadIdx 中可以得到线程的唯一 ID 。（threadIdx 是一个有 3 个元素的向量，可以用来从一维、二维、三维角度标识线程）

一维线程块：线程 ID = threadIdx.x

二维线程块：线程 ID = threadIdx.x + threadIdx.y * blockDim.x

三维线程块：线程 ID = threadIdx.x + threadIdx.y * blockDim.x + threadIdx.z * blockDim.x * blockDim.y

内置变量 blockDim 存储 线程块 的维度，gridDim 存储 线程块网格 的维度。

同理，在内置变量 blockIdx 中可以得到线程块的唯一 ID 。过程同上

  • 内存层次结构

    内存层次结构 与 线程层次结构 保持一致，每一层次的线程结构都有对应的内存空间。

  • 流同步控制

    流是按顺序执行的命令序列。不同的流可能会相互乱序或同时执行命令，但流内的命令时顺序执行的。

    未显示指定流的命令会加入到 NULL 流，即默认同步流。

# 编程范式

  • 申请 GPU 设备内存

    cudaMalloc();

e.g.
    double* d_elements = nullptr;
    cudaMalloc(&d_elements, 100 * sizeof(double));

  • 拷贝数据

    cudaMemcpy();

e.g.
    double* h_elements = new double[100];
    cudaMemcpy(d_elements, h_elements, 100 * sizeof(double), cudaMemcpyHostToDevice);

  • 调用 GPU kernel核函数

    kernel_function<<<grid_size, block_size, shared_size, stream>>>(param1, param2);

  • 释放内存

    cudaFree();

e.g.
    cudaFree(d_elements)

  • 同步：同步线程块内的线程

    __syncthreads();

  • 流

    cudaStream_t stream = NULL;
cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking);

cudaMemcpyAsync(d_elements, h_elements, 100 * sizeof(double), cudaMemcpyHostToDevice, stream);

// 将流作为参数并等待，直到给定流中的所有命令都完成。
cudaStreamSynchronize(stream);
cudaStreamDestroy(stream)

    未指定流的 kernel 函数或 CUDA 调用会被放置在默认流（NULL）中。

# 入门实例——矩阵乘法

基础函数

gemm_test.h：

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>

typedef struct
{
    int width;
    int height;
    int stride;
    double* elements;
} Matrix;

#define GetElement(M, row, col) (M.elements[row * M.stride + col])

void RandomData(Matrix M)
{
    srand(time(0));
    int size = M.width * M.height;
    for (int i = 0; i < size; ++i)
    {
        M.elements[i] = double(rand()) / RAND_MAX;
    }
}

bool InitTestData(Matrix& A, Matrix& B, Matrix& C)
{
    // 计算复杂度大概为 2^30 约 1000000000
    // A 512 * 2048
    // B 2048 * 1024
    // C 512 * 1024
    A.stride = A.width = 2048;
    A.height = 512;
    B.stride = B.width = 1024;
    B.height = 2048;
    C.stride = C.width = 1024;
    C.height = 512;

    A.elements = new double[A.width * A.height]{};
    B.elements = new double[B.width * B.height]{};
    C.elements = new double[C.width * C.height]{};

    if (nullptr == A.elements || nullptr == B.elements || nullptr == C.elements)
    {
        printf("new elements error!\n");
        return false;
    }

    RandomData(A);
    RandomData(B);

    return true;
}

void DestoryData(Matrix& A, Matrix& B, Matrix& C)
{
    delete[] A.elements;
    delete[] B.elements;
    delete[] C.elements;
    A.elements = nullptr;
    B.elements = nullptr;
    C.elements = nullptr;
}

int64_t GetNowTimeNs()
{
    struct timespec time_spec {};
    clock_gettime(CLOCK_REALTIME, &time_spec);
    return time_spec.tv_sec * 1000000000 + time_spec.tv_nsec;
}

void CheckResult(Matrix& A, Matrix& B, Matrix& C)
{
    for (int i = 0; i < C.height; ++i)
    {
        for (int j = 0; j < C.width; ++j)
        {
            double result = 0.0;
            for (int k = 0; k < A.width; ++k)
            {
                result += (GetElement(A, i, k) * GetElement(B, k, j));
            }
            if (fabs(result - GetElement(C, i, j)) > 1e-5)
            {
                printf("result is false:(%d, %d): (%lf, %lf)\n", i, j, result, GetElement(C, i, j));
            }
        }
    }
}


// CUDA 错误处理
#include <cuda_runtime_api.h>

#define CUDA_CHECK(err)                                                             \
    do                                                                              \
    {                                                                               \
        cudaError_t err_ = (err);                                                   \
        if (err_ != cudaSuccess)                                                    \
        {                                                                           \
            printf("CUDA error %d at %s:%d\n", err_, __FILE__, __LINE__);           \
        }                                                                           \
    } while (0)

#define CUBLAS_CHECK(err)                                                           \
    do                                                                              \
    {                                                                               \
        cublasStatus_t err_ = (err);                                                \
        if (err_ != CUBLAS_STATUS_SUCCESS)                                          \
        {                                                                           \
            printf("cublas error %d at %s:%d\n", err_, __FILE__, __LINE__);         \
        }                                                                           \
    } while (0)

# CPU单线程

#include "gemm_test.h"

void MatMul(const Matrix A, const Matrix B, Matrix C)
{
    for (int i = 0; i < C.height; ++i)
    {
        for (int j = 0; j < C.width; ++j)
        {
            for (int k = 0; k < A.width; ++k)
            {
                GetElement(C, i, j) += (GetElement(A, i, k) * GetElement(B, k, j));
            }
        }
    }
}

int main()
{
    Matrix A, B, C;
    if (!InitTestData(A, B, C))
    {
        return -1;
    }

    // 计时
    int64_t start_time = 0, end_time = 0;
    const int size = 100;
    // 开始时间
    start_time = GetNowTimeNs();

    for (int i = 0; i < size; ++i)
    {
        MatMul(A, B, C);
    }

    // 结束时间
    end_time = GetNowTimeNs();

    printf("complete!%lf s\n", double(end_time - start_time) / size / 1000000000);

    DestoryData(A, B, C);

    return 0;
}

# CUDA 共享内存优化

每个线程计算一个元素：

A 矩阵从全局内存被读取了 B.width 次， B 矩阵从全局内存被读取了 A.height 次。

每个线程计算一个元素，但是按照线程块大小切割原矩阵并读到线程块共享内存中：

A 矩阵从全局内存被读取了 B.width/BLOCK_SIZE 次， B 矩阵从全局内存被读取了 A.height/BLOCK_SIZE 次。

#include "gemm_test.h"

// 线程块行列大小
#define BLOCK_SIZE 16

// 将 A 矩阵划分为 BLOCK_SIZExBLOCK_SIZE 大小分割的子矩阵后，
// 获取位于 (row, col) 的子矩阵 Asub
__device__ Matrix GetSubMatrix(Matrix A, int row, int col)
{
    Matrix Asub;
    Asub.width    = BLOCK_SIZE;
    Asub.height   = BLOCK_SIZE;
    Asub.stride   = A.stride;
    Asub.elements = &A.elements[A.stride * BLOCK_SIZE * row + BLOCK_SIZE * col];
    return Asub;
}

// 前向声明
__global__ void MatMulKernel(const Matrix, const Matrix, Matrix);

// 矩阵乘法 主机c++接口
void MatMul(const Matrix A, const Matrix B, Matrix C)
{
    // 加载 A, B, C 矩阵到 GPU 设备
    Matrix d_A;
    d_A.width = d_A.stride = A.width; d_A.height = A.height;
    size_t Asize = A.width * A.height * sizeof(double);
    Matrix d_B;
    d_B.width = d_B.stride = B.width; d_B.height = B.height;
    size_t Bsize = B.width * B.height * sizeof(double);
    Matrix d_C;
    d_C.width = d_C.stride = C.width; d_C.height = C.height;
    size_t Csize = C.width * C.height * sizeof(double);

    CUDA_CHECK(cudaMalloc(&d_A.elements, Asize));
    CUDA_CHECK(cudaMalloc(&d_B.elements, Bsize));
    CUDA_CHECK(cudaMalloc(&d_C.elements, Csize));

    CUDA_CHECK(cudaMemcpy(d_A.elements, A.elements, Asize, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_B.elements, B.elements, Bsize, cudaMemcpyHostToDevice));

    // 调用 GPU 矩阵计算函数
    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
    dim3 dimGrid(B.width / dimBlock.x, A.height / dimBlock.y);
    MatMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C);
    
    // 从 GPU 拷贝结果到 主机内存
    CUDA_CHECK(cudaMemcpy(C.elements, d_C.elements, Csize, cudaMemcpyDeviceToHost));

    // 释放 GPU 内存
    CUDA_CHECK(cudaFree(d_A.elements));
    CUDA_CHECK(cudaFree(d_B.elements));
    CUDA_CHECK(cudaFree(d_C.elements));

    CUDA_CHECK(cudaDeviceReset());
}
// 矩阵乘法 GPU 设备 调用函数(kernel)
__global__ void MatMulKernel(Matrix A, Matrix B, Matrix C)
{
    // 获取线程块的坐标
    int blockRow = blockIdx.y;
    int blockCol = blockIdx.x;

    // 每个线程块负责计算一个 C 的子矩阵 Csub
    Matrix Csub = GetSubMatrix(C, blockRow, blockCol);
    
    // 每个线程负责计算 子矩阵 中的一个元素
    double Cvalue = 0;
    // 获取线程在线程块中的坐标
    int row = threadIdx.y;
    int col = threadIdx.x;
    // 循环遍历计算 Csub 所需的 A 和 B 的所有子矩阵,将每对子矩阵相乘并累加
    for (int m = 0; m < (A.width / BLOCK_SIZE); ++m) {
        // 获取 A 的对应子矩阵
        Matrix Asub = GetSubMatrix(A, blockRow, m);
        // 获取 B 的对应子矩阵
        Matrix Bsub = GetSubMatrix(B, m, blockCol);
        // 分别用于存储 Asub 和 Bsub 的共享内存
        __shared__ double As[BLOCK_SIZE][BLOCK_SIZE];
        __shared__ double Bs[BLOCK_SIZE][BLOCK_SIZE];
        // 将 Asub 和 Bsub 从设备内存加载到共享内存
 		// 每个线程加载每个子矩阵的一个元素
        As[row][col] = GetElement(Asub, row, col);
        Bs[row][col] = GetElement(Bsub, row, col);
        // 同步以确保子矩阵加载
        __syncthreads();
        // 计算 Cvalue
        for (int e = 0; e < BLOCK_SIZE; ++e)
            Cvalue += As[row][e] * Bs[e][col];
        // 同步，以确保在下一次迭代中逻辑正确
        __syncthreads();
    }
    // 将结果写入 GPU 设备内存
    GetElement(Csub, row, col) = Cvalue;
}

int main()
{
    Matrix A, B, C;

    if (!InitTestData(A, B, C))
    {
        return -1;
    }

    // 计时
    int64_t start_time = 0, end_time = 0;
    const int size = 100;
    // 开始时间
    start_time = GetNowTimeNs();

    for (int i = 0; i < size; ++i)
    {
        MatMul(A, B, C);
    }

    // 结束时间
    end_time = GetNowTimeNs();

    printf("complete!%lf s\n", double(end_time - start_time) / size / 1000000000);

    CheckResult(A, B, C);

    DestoryData(A, B, C);

    return 0;
}

# CUDA cublas库

#include <cublas_v2.h>

#include "gemm_test.h"

// cublas默认按照列优先存储,c++代码一般按照行优先存储。
// 所以利用 cublas 库得到的结果是 C 的转置 CT。
// 所以我们可以利用 CT = (A*B)T = BT * AT,来计算 CT , 然后经过 cublas 计算直接得到 (CT)T = C
// 所以下述函数计算 BT * AT
void MatMul(const Matrix A, const Matrix B, Matrix C)
{
    cublasHandle_t cublasH = NULL;
    cudaStream_t stream = NULL;

    const double alpha = 1.0;
    const double beta = 0.0;
    double* d_A = nullptr;
    double* d_B = nullptr;
    double* d_C = nullptr;

    cublasOperation_t transa = CUBLAS_OP_N;
    cublasOperation_t transb = CUBLAS_OP_N;

    size_t Asize = A.width * A.height * sizeof(double);
    size_t Bsize = B.width * B.height * sizeof(double);
    size_t Csize = C.width * C.height * sizeof(double);

    // 创建 cublas 句柄，绑定数据流
    CUBLAS_CHECK(cublasCreate(&cublasH));

    CUDA_CHECK(cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking));
    CUBLAS_CHECK(cublasSetStream(cublasH, stream));

    // 拷贝数据到 GPU 设备内存上
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void**>(&d_A), Asize));
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void**>(&d_B), Bsize));
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void**>(&d_C), Csize));

    CUDA_CHECK(cudaMemcpyAsync(d_A, A.elements, Asize, cudaMemcpyHostToDevice, stream));
    CUDA_CHECK(cudaMemcpyAsync(d_B, B.elements, Bsize, cudaMemcpyHostToDevice, stream));

    // 调用 cublas 库计算
    // cublasDgemm 计算 F = αD*E + βF (此处 D = BT, E = AT, F = CT)
    // 注意 cublas 默认列优先，
    CUBLAS_CHECK(cublasDgemm(cublasH, 	// cuBLAS 库上下文的句柄。
                             transa, 	// D 矩阵是否转置
                             transb, 	// E 矩阵是否转置
                             B.width, 	// D 矩阵的行数(F的行数)
                             A.height, 	// E 矩阵的列数(F的列数)
                             A.width, 	// D 矩阵的列数， E 矩阵的行数
                             &alpha, 	// α 的值
                             d_B, 		// D 矩阵地址
                             B.width, 	// D 矩阵存储时第一维维数
                             d_A, 		// E 矩阵地址
                             A.width, 	// E 矩阵存储时第一维维数
                             &beta, 	// β
                             d_C, 		// F 矩阵地址
                             C.width	// F 矩阵存储时第一维维数
                            ));

    // 拷贝结果到 主机内存上
    CUDA_CHECK(cudaMemcpyAsync(C.elements, d_C, Csize, cudaMemcpyDeviceToHost, stream));

    CUDA_CHECK(cudaStreamSynchronize(stream));

    // 释放内存
    CUDA_CHECK(cudaFree(d_A));
    CUDA_CHECK(cudaFree(d_B));
    CUDA_CHECK(cudaFree(d_C));

    CUBLAS_CHECK(cublasDestroy(cublasH));

    CUDA_CHECK(cudaStreamDestroy(stream));

    CUDA_CHECK(cudaDeviceReset());
}

int main()
{
    Matrix A, B, C;

    if (!InitTestData(A, B, C))
    {
        return -1;
    }

    // 计时
    int64_t start_time = 0, end_time = 0;
    const int size = 100;
    // 开始时间
    start_time = GetNowTimeNs();

    for (int i = 0; i < size; ++i)
    {
        MatMul(A, B, C);
    }

    // 结束时间
    end_time = GetNowTimeNs();

    printf("complete!%lf s\n", double(end_time - start_time) / size / 1000000000);

    CheckResult(A, B, C);

    DestoryData(A, B, C);

    return 0;
}

# 效率对比

运行平台为 Google Colaboratory，不能说明绝对效率，仅仅考察相对效率。

原始数据

矩阵大小 复杂度 CPU 单线程 CUDA 共享内存优化代码 CUDA cublas 库
（16*16)*(16*16） 2^12 0.000006 0.091786 0.644954
（64*64)*(64*64） 2^18 0.000279 0.090339 0.637052
（16*1024)*(1024*16） 2^18 0.000394 0.089789 0.626902
（256*256)*(256*256） 2^24 0.031874 0.089906 0.625845
（1024*16)*(16*1024） 2^24 0.032913 0.087774 0.643814
（256*512)*(512*256） 2^25 0.075586 0.092250 0.622241
（512*256)*(256*512） 2^26 0.206449 0.088731 0.630597
（512*512)*(512*512） 2^27 0.449558 0.096542 0.630839
（512*1024)*(1024*512） 2^28 0.910838 0.098378 0.632863
（1024*1024)*(1024*1024） 2^30 3.726373 0.121930 0.665017
（2048*2048)*(2048*2048） 2^33 93.738654 0.192806 0.718952
（4096*4096)*(4096*4096） 2^36 0.855580 1.305331
（4096*8192)*(8192*4096） 2^37 1.503124 1.902856
（8192*4096)*(4096*8192） 2^38 2.904895 3.094869
（8192*8192)*(8192*8192） 2^39 5.991753 5.442394
（8192*16384)*(16384*8192） 2^40 11.731695 10.189183
（16384*8192)*(8192*16384） 2^41 22.515630 19.574606
（16384*16384)*(16384*16384） 2^42 43.217113 37.017278

​如上图，cpu 的时间开销随着计算复杂度线性上升，而在其区间内 cuda 的时间开销几乎不变。说明在此区间的计算复杂度下，cuda 的主要开销是数据拷贝。且数据量小时，cuda GPU 并行计算的加速效果抵消不了拷贝的时间开销。

​如上图，将横坐标取对数后可以将更多的数据更好的呈现在统计图中。cpu 的时间开销将迅速增长，而 cuda GPU 并行计算仍然没有全负荷运行。

​如上两图，利用了共享内存优化的 cuda 核函数虽然在较小数据量时有更高的效率（主要是共享内存的读取优化），但其时间增长斜率仍然大于 cuda cublas 库函数的斜率。所以实际上，cuda cublas 库函数不仅更加通用，且依旧高效（说明当库中存在符合需求的函数接口时，尽量使用库接口）。

# 总结

  • CUDA 编程的主要流程比较固定，基本如下图：

  • CUDA 编程解决实际问题的两大核心：

    • 哪些问题可以且值得被GPU并行加速计算。
    • 如何把问题分解为互不相关的高度并行化的子问题。
六大原则
Centos6.7 安装 python3

← Centos6.7 安装 python3→

最近更新
01
深度学习简单介绍
04/18
02
VScode-ssh 不兼容旧版本 Linux
02/08
03
Centos6.7 安装 python3
01/25
更多文章>
Theme by Vdoing | Copyright © 2023-2024 苏ICP备2023014117号-1 苏公网安备 32011502011853号