Cython与CUDA之Add

rXFJfZxJUBae

技术背景

在前一篇文章中，我们介绍过使用Cython结合CUDA实现了一个Gather算子以及一个BatchGather算子。这里我们继续使用这一套方案，实现一个简单的求和函数，通过CUDA来计算数组求和。由于数组求和对于不同的维度来说都是元素对元素进行求和，因此高维数组跟低维数组没有差别，这里我们全都当做是一维的数组输入来处理，不做Batch处理。
头文件

首先我们需要一个CUDA头文件cuda_add.cuh来定义CUDA函数的接口：
#include <stdio.h>extern "C" float Add(float *A, float *B, float *res, int N);其他头文件如异常捕获，可以参考这篇文章，CUDA函数计时可以参考这篇文章。
CUDA文件

CUDA文件cuda_add.cu中包含了核心部分的算法：
// nvcc -shared ./cuda_add.cu -Xcompiler -fPIC -o ./libcuadd.so#include <stdio.h>#include "cuda_add.cuh"#include "error.cuh"#include "record.cuh"__global__ void AddKernel(float *A, float *B, float *res, int N) {    int tid = blockIdx.x * blockDim.x + threadIdx.x;    // 每个线程处理多个元素    int stride = blockDim.x * gridDim.x;    for (int i = tid; i < N; i += stride) {        res = A + B;    }}extern "C" float Add(float *A, float *B, float *res, int N){    float *A_device, *B_device, *res_device;    CHECK(cudaMalloc((void **)&A_device, N * sizeof(float)));    CHECK(cudaMalloc((void **)&B_device, N * sizeof(float)));    CHECK(cudaMalloc((void **)&res_device, N * sizeof(float)));    CHECK(cudaMemcpy(A_device, A, N * sizeof(float), cudaMemcpyHostToDevice));    CHECK(cudaMemcpy(B_device, B, N * sizeof(float), cudaMemcpyHostToDevice));    int block_size, grid_size;    cudaOccupancyMaxPotentialBlockSize(&grid_size, &block_size, AddKernel, 0, N);    grid_size = (N + block_size - 1) / block_size;    float timeTaken = GET_CUDA_TIME((AddKernel<<<grid_size, block_size>>>(A_device, B_device, res_device, N)));    CHECK(cudaGetLastError());    CHECK(cudaDeviceSynchronize());    CHECK(cudaMemcpy(res, res_device, N * sizeof(float), cudaMemcpyDeviceToHost));    CHECK(cudaFree(A_device));    CHECK(cudaFree(B_device));    CHECK(cudaFree(res_device));    return timeTaken;}此处代码是部分经过DeepSeek优化的，例如在核函数中使用for循环对多个数据进行处理，而不是只处理一个数据。另外block_size由cudaOccupancyMaxPotentialBlockSize自动生成，也避免了手动设定带来的一些麻烦。不过这里我们没有使用Stream来优化，只是简单的演示一个功能算法。
Cython接口文件

由于我们的框架是通过Cython来封装CUDA函数，然后在Python中调用，所以这里需要一个Cython接口文件wrapper.pyx。
# cythonize -i -f wrapper.pyximport numpy as npcimport numpy as npcimport cythoncdef extern from "<dlfcn.h>" nogil:    void *dlopen(const char *, int)    char *dlerror()    void *dlsym(void *, const char *)    int dlclose(void *)    enum:        RTLD_LAZYctypedef float (*AddFunc)(float *A, float *B, float *res, int N) noexcept nogilcdef void* handle_add = dlopen('/path/to/cuda/libcuadd.so', RTLD_LAZY)@cython.boundscheck(False)@cython.wraparound(False)cpdef float[:] cuda_add(float[:] x, float[:] y):    cdef:        AddFunc Add        float timeTaken        int N = x.shape[0]        float[:] res = np.zeros((N, ), dtype=np.float32)    Add = <AddFunc>dlsym(handle_add, "Add")    timeTaken = Add(&x[0], &y[0], &res[0], N)    print (timeTaken)    return reswhile not True:    dlclose(handle)Python调用文件

最后，我们写一个Python的案例test_add.py来调用Cython封装后的CUDA函数：
import numpy as npnp.random.seed(0)from wrapper import cuda_addN = 1024 * 1024 * 100x = np.random.random((N,)).astype(np.float32)y = np.random.random((N,)).astype(np.float32)np_res = x+yres = np.asarray(cuda_add(x, y))print (res.shape)print ((res==np_res).sum())运行python文件即可获得CUDA核函数的耗时，以及相应的返回结果输出。
查看GPU信息

为了更加深刻的理解一下CUDA计算的性能，我们可以查看GPU的一些关键参数，以此来推理CUDA运算的理论运算极限。在一些版本的CUDA里面会自带一个deviceQuery：
$ cd /usr/local/cuda-10.1/samples/1_Utilities/deviceQuery里面包含有一些可以查询获取本地GPU配置参数的文件：
$ ll总用量 44drwxr-xr-x 2 root root  4096 7月  13  2021 ./drwxr-xr-x 8 root root  4096 7月  13  2021 ../-rw-r--r-- 1 root root 12473 7月  13  2021 deviceQuery.cpp-rw-r--r-- 1 root root 10812 7月  13  2021 Makefile-rw-r--r-- 1 root root  1789 7月  13  2021 NsightEclipse.xml-rw-r--r-- 1 root root   168 7月  13  2021 readme.txt可以将这些文件进行编译，但是因为这些代码强行指定了nvcc的地址在/usr/local/cuda下，所以如果本地没有这个路径的，可能需要使用ln -s来创建一个路径软链接：
$ sudo ln -s /usr/local/cuda-10.1 /usr/local/cuda然后再执行编译指令：
$ sudo make/usr/local/cuda/bin/nvcc -ccbin g++ -I../../common/inc  -m64    -gencode arch=compute_30,code=sm_30 -gencode arch=compute_35,code=sm_35 -gencode arch=compute_37,code=sm_37 -gencode arch=compute_50,code=sm_50 -gencode arch=compute_52,code=sm_52 -gencode arch=compute_60,code=sm_60 -gencode arch=compute_61,code=sm_61 -gencode arch=compute_70,code=sm_70 -gencode arch=compute_75,code=sm_75 -gencode arch=compute_75,code=compute_75 -o deviceQuery.o -c deviceQuery.cpp/usr/local/cuda/bin/nvcc -ccbin g++   -m64      -gencode arch=compute_30,code=sm_30 -gencode arch=compute_35,code=sm_35 -gencode arch=compute_37,code=sm_37 -gencode arch=compute_50,code=sm_50 -gencode arch=compute_52,code=sm_52 -gencode arch=compute_60,code=sm_60 -gencode arch=compute_61,code=sm_61 -gencode arch=compute_70,code=sm_70 -gencode arch=compute_75,code=sm_75 -gencode arch=compute_75,code=compute_75 -o deviceQuery deviceQuery.o mkdir -p ../../bin/x86_64/linux/releasecp deviceQuery ../../bin/x86_64/linux/release编译完成后直接执行编译好的可执行文件：
$ ./deviceQuery ./deviceQuery Starting... CUDA Device Query (Runtime API) version (CUDART static linking)Detected 2 CUDA Capable device(s)Device 0: "Quadro RTX 4000"  CUDA Driver Version / Runtime Version          12.2 / 10.1  CUDA Capability Major/Minor version number:    7.5  Total amount of global memory:                 7972 MBytes (8358723584 bytes)  (36) Multiprocessors, ( 64) CUDA Cores/MP:     2304 CUDA Cores  GPU Max Clock rate:                            1545 MHz (1.54 GHz)  Memory Clock rate:                             6501 Mhz  Memory Bus Width:                              256-bit  L2 Cache Size:                                 4194304 bytes  Maximum Texture Dimension Size (x,y,z)         1D=(131072), 2D=(131072, 65536), 3D=(16384, 16384, 16384)  Maximum Layered 1D Texture Size, (num) layers  1D=(32768), 2048 layers  Maximum Layered 2D Texture Size, (num) layers  2D=(32768, 32768), 2048 layers  Total amount of constant memory:               65536 bytes  Total amount of shared memory per block:       49152 bytes  Total number of registers available per block: 65536  Warp size:                                     32  Maximum number of threads per multiprocessor:  1024  Maximum number of threads per block:           1024  Max dimension size of a thread block (x,y,z): (1024, 1024, 64)  Max dimension size of a grid size    (x,y,z): (2147483647, 65535, 65535)  Maximum memory pitch:                          2147483647 bytes  Texture alignment:                             512 bytes  Concurrent copy and kernel execution:          Yes with 3 copy engine(s)  Run time limit on kernels:                     Yes  Integrated GPU sharing Host Memory:            No  Support host page-locked memory mapping:       Yes  Alignment requirement for Surfaces:            Yes  Device has ECC support:                        Disabled  Device supports Unified Addressing (UVA):      Yes  Device supports Compute Preemption:            Yes  Supports Cooperative Kernel Launch:            Yes  Supports MultiDevice Co-op Kernel Launch:      Yes  Device PCI Domain ID / Bus ID / location ID:   0 / 3 / 0  Compute Mode:     < Default (multiple host threads can use ::cudaSetDevice() with device simultaneously) >Device 1: "Quadro RTX 4000"  CUDA Driver Version / Runtime Version          12.2 / 10.1  CUDA Capability Major/Minor version number:    7.5  Total amount of global memory:                 7974 MBytes (8361738240 bytes)  (36) Multiprocessors, ( 64) CUDA Cores/MP:     2304 CUDA Cores  GPU Max Clock rate:                            1545 MHz (1.54 GHz)  Memory Clock rate:                             6501 Mhz  Memory Bus Width:                              256-bit  L2 Cache Size:                                 4194304 bytes  Maximum Texture Dimension Size (x,y,z)         1D=(131072), 2D=(131072, 65536), 3D=(16384, 16384, 16384)  Maximum Layered 1D Texture Size, (num) layers  1D=(32768), 2048 layers  Maximum Layered 2D Texture Size, (num) layers  2D=(32768, 32768), 2048 layers  Total amount of constant memory:               65536 bytes  Total amount of shared memory per block:       49152 bytes  Total number of registers available per block: 65536  Warp size:                                     32  Maximum number of threads per multiprocessor:  1024  Maximum number of threads per block:           1024  Max dimension size of a thread block (x,y,z): (1024, 1024, 64)  Max dimension size of a grid size    (x,y,z): (2147483647, 65535, 65535)  Maximum memory pitch:                          2147483647 bytes  Texture alignment:                             512 bytes  Concurrent copy and kernel execution:          Yes with 3 copy engine(s)  Run time limit on kernels:                     Yes  Integrated GPU sharing Host Memory:            No  Support host page-locked memory mapping:       Yes  Alignment requirement for Surfaces:            Yes  Device has ECC support:                        Disabled  Device supports Unified Addressing (UVA):      Yes  Device supports Compute Preemption:            Yes  Supports Cooperative Kernel Launch:            Yes  Supports MultiDevice Co-op Kernel Launch:      Yes  Device PCI Domain ID / Bus ID / location ID:   0 / 166 / 0  Compute Mode:     < Default (multiple host threads can use ::cudaSetDevice() with device simultaneously) >> Peer access from Quadro RTX 4000 (GPU0) -> Quadro RTX 4000 (GPU1) : Yes> Peer access from Quadro RTX 4000 (GPU1) -> Quadro RTX 4000 (GPU0) : YesdeviceQuery, CUDA Driver = CUDART, CUDA Driver Version = 12.2, CUDA Runtime Version = 10.1, NumDevs = 2Result = PASS这里就输出了两块GPU的相关参数。其中Memory Bus Width:                              256-bit表示总位宽，数值越高越好。Memory Clock rate:                             6501 Mhz表示显存的访问速率，经常被用于估计GPU的性能，因为很多时候GPU的性能瓶颈可能在内存-显存的传输上。GPU Max Clock rate:                            1545 MHz (1.54 GHz)可以用来估计显存操作速率。
性能估算

以普通的CUDA加法为例，有效速率的大致公式为：

\[有效速率(Gbps)=\frac{物理频率\times 2}{1000}\]

进而可以计算带宽：

\[带宽(GB/s)=\frac{有效速率\times 总线宽度}{8}\]

最后，根据带宽估算计算速率的上限，也就等同于估算一个CUDA加法的计算耗时的下限：

\[计算耗时(s)=\frac{总操作数据量(B)}{带宽(B/s)}\]

实际计算的话，单次的加法操作涉及到四个步骤：读取数组A元素，读取数组B元素，加和，写入C数组。也就是说，涉及到3次内存操作和1次加和操作。关于内存部分的耗时（假定N=1024*1024*100）：

\[T_{mem}=\frac{N*4*3}{\frac{\frac{6501}{1000}*256}{8}*10^9}\approx 0.0030243 s\]

耗时估计在3ms左右（这里的4是单精度浮点数到Byte的换算）。至于单次的加法运算耗时，其实可以忽略不计，因为指令吞吐率大概为：

\[指令吞吐率(TFLOPS)=核心数\times 时钟频率=2304\times 1.54e09\approx 3.55\]

那么理论最小耗时为（假定N=1024*1024*100）：

\[T_{theo}(s)=\frac{总计算量}{指令吞吐率}\approx 2.95e-05\]

指令运算部分耗时大约在0.03 ms，跟显存IO部分的耗时3 ms比起来可以忽略的量级。
真实测试

运行Python代码输出的结果为：
$ python3 test_add.py 3.3193600177764893(104857600,)104857600这个数据3.32 ms已经很接近于极限速率3 ms了，应该说在这样的算法框架下已经很难再往下去优化了，更多时候优化点还是在于CPU到GPU的内存传输效率上。
总结概要

本文介绍了使用CUDA和Cython来实现一个CUDA加法算子的方法，并介绍了使用CUDA参数来估算性能极限的算法。经过实际测试，核函数部分的算法性能优化空间已经不是很大了，更多时候可以考虑使用Stream来优化Host和Device之间的数据传输。
版权声明

本文首发链接为：https://www.cnblogs.com/dechinphy/p/cuda-cython-add.html
作者ID：DechinPhy
更多原著文章：https://www.cnblogs.com/dechinphy/
请博主喝咖啡：https://www.cnblogs.com/dechinphy/gallery/image/379634.html
参考链接

• https://blog.csdn.net/sunyuhua_keyboard/article/details/145633805