Cython与CUDA之Add

技术背景

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

首先我们需要一个CUDA头文件cuda_add.cuh来定义CUDA函数的接口：

1. #include <stdio.h>
3. extern "C" float Add(float *A, float *B, float *res, int N);

复制代码

其他头文件如异常捕获，可以参考这篇文章，CUDA函数计时可以参考这篇文章。
CUDA文件

CUDA文件cuda_add.cu中包含了核心部分的算法：

1. // nvcc -shared ./cuda_add.cu -Xcompiler -fPIC -o ./libcuadd.so
2. #include <stdio.h>
3. #include "cuda_add.cuh"
4. #include "error.cuh"
5. #include "record.cuh"
7. __global__ void AddKernel(float *A, float *B, float *res, int N) {
8.     int tid = blockIdx.x * blockDim.x + threadIdx.x;
9.     // 每个线程处理多个元素
10.     int stride = blockDim.x * gridDim.x;
11.     for (int i = tid; i < N; i += stride) {
12.         res[i] = A[i] + B[i];
13.     }
14. }
16. extern "C" float Add(float *A, float *B, float *res, int N){
17.     float *A_device, *B_device, *res_device;
18.     CHECK(cudaMalloc((void **)&A_device, N * sizeof(float)));
19.     CHECK(cudaMalloc((void **)&B_device, N * sizeof(float)));
20.     CHECK(cudaMalloc((void **)&res_device, N * sizeof(float)));
21.     CHECK(cudaMemcpy(A_device, A, N * sizeof(float), cudaMemcpyHostToDevice));
22.     CHECK(cudaMemcpy(B_device, B, N * sizeof(float), cudaMemcpyHostToDevice));
24.     int block_size, grid_size;
25.     cudaOccupancyMaxPotentialBlockSize(&grid_size, &block_size, AddKernel, 0, N);
26.     grid_size = (N + block_size - 1) / block_size;
28.     float timeTaken = GET_CUDA_TIME((AddKernel<<<grid_size, block_size>>>(A_device, B_device, res_device, N)));
29.     CHECK(cudaGetLastError());
30.     CHECK(cudaDeviceSynchronize());
31.     CHECK(cudaMemcpy(res, res_device, N * sizeof(float), cudaMemcpyDeviceToHost));
32.     CHECK(cudaFree(A_device));
33.     CHECK(cudaFree(B_device));
34.     CHECK(cudaFree(res_device));
35.     return timeTaken;
36. }

复制代码

此处代码是部分经过DeepSeek优化的，例如在核函数中使用for循环对多个数据进行处理，而不是只处理一个数据。另外block_size由cudaOccupancyMaxPotentialBlockSize自动生成，也避免了手动设定带来的一些麻烦。不过这里我们没有使用Stream来优化，只是简单的演示一个功能算法。
Cython接口文件

由于我们的框架是通过Cython来封装CUDA函数，然后在Python中调用，所以这里需要一个Cython接口文件wrapper.pyx。

1. # cythonize -i -f wrapper.pyx
3. import numpy as np
4. cimport numpy as np
5. cimport cython
7. cdef extern from "<dlfcn.h>" nogil:
8.     void *dlopen(const char *, int)
9.     char *dlerror()
10.     void *dlsym(void *, const char *)
11.     int dlclose(void *)
12.     enum:
13.         RTLD_LAZY
15. ctypedef float (*AddFunc)(float *A, float *B, float *res, int N) noexcept nogil
17. cdef void* handle_add = dlopen('/path/to/cuda/libcuadd.so', RTLD_LAZY)
19. @cython.boundscheck(False)
20. @cython.wraparound(False)
21. cpdef float[:] cuda_add(float[:] x, float[:] y):
22.     cdef:
23.         AddFunc Add
24.         float timeTaken
25.         int N = x.shape[0]
26.         float[:] res = np.zeros((N, ), dtype=np.float32)
27.     Add = dlsym(handle_add, "Add")
28.     timeTaken = Add(&x[0], &y[0], &res[0], N)
29.     print (timeTaken)
30.     return res
32. while not True:
33.     dlclose(handle)

复制代码

Python调用文件

最后，我们写一个Python的案例test_add.py来调用Cython封装后的CUDA函数：

1. import numpy as np
2. np.random.seed(0)
3. from wrapper import cuda_add
5. N = 1024 * 1024 * 100
7. x = np.random.random((N,)).astype(np.float32)
8. y = np.random.random((N,)).astype(np.float32)
10. np_res = x+y
12. res = np.asarray(cuda_add(x, y))
13. print (res.shape)
14. print ((res==np_res).sum())

复制代码

运行python文件即可获得CUDA核函数的耗时，以及相应的返回结果输出。
查看GPU信息

为了更加深刻的理解一下CUDA计算的性能，我们可以查看GPU的一些关键参数，以此来推理CUDA运算的理论运算极限。在一些版本的CUDA里面会自带一个deviceQuery：

1. $ cd /usr/local/cuda-10.1/samples/1_Utilities/deviceQuery

复制代码

里面包含有一些可以查询获取本地GPU配置参数的文件：

1. $ ll
2. 总用量 44
3. drwxr-xr-x 2 root root  4096 7月  13  2021 ./
4. drwxr-xr-x 8 root root  4096 7月  13  2021 ../
5. -rw-r--r-- 1 root root 12473 7月  13  2021 deviceQuery.cpp
6. -rw-r--r-- 1 root root 10812 7月  13  2021 Makefile
7. -rw-r--r-- 1 root root  1789 7月  13  2021 NsightEclipse.xml
8. -rw-r--r-- 1 root root   168 7月  13  2021 readme.txt

复制代码

可以将这些文件进行编译，但是因为这些代码强行指定了nvcc的地址在/usr/local/cuda下，所以如果本地没有这个路径的，可能需要使用ln -s来创建一个路径软链接：

1. $ sudo ln -s /usr/local/cuda-10.1 /usr/local/cuda

复制代码

然后再执行编译指令：

1. $ sudo make
2. /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
3. /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
4. mkdir -p ../../bin/x86_64/linux/release
5. cp deviceQuery ../../bin/x86_64/linux/release

复制代码

编译完成后直接执行编译好的可执行文件：

1. $ ./deviceQuery
2. ./deviceQuery Starting...
4. CUDA Device Query (Runtime API) version (CUDART static linking)
6. Detected 2 CUDA Capable device(s)
8. Device 0: "Quadro RTX 4000"
9.   CUDA Driver Version / Runtime Version          12.2 / 10.1
10.   CUDA Capability Major/Minor version number:    7.5
11.   Total amount of global memory:                 7972 MBytes (8358723584 bytes)
12.   (36) Multiprocessors, ( 64) CUDA Cores/MP:     2304 CUDA Cores
13.   GPU Max Clock rate:                            1545 MHz (1.54 GHz)
14.   Memory Clock rate:                             6501 Mhz
15.   Memory Bus Width:                              256-bit
16.   L2 Cache Size:                                 4194304 bytes
17.   Maximum Texture Dimension Size (x,y,z)         1D=(131072), 2D=(131072, 65536), 3D=(16384, 16384, 16384)
18.   Maximum Layered 1D Texture Size, (num) layers  1D=(32768), 2048 layers
19.   Maximum Layered 2D Texture Size, (num) layers  2D=(32768, 32768), 2048 layers
20.   Total amount of constant memory:               65536 bytes
21.   Total amount of shared memory per block:       49152 bytes
22.   Total number of registers available per block: 65536
23.   Warp size:                                     32
24.   Maximum number of threads per multiprocessor:  1024
25.   Maximum number of threads per block:           1024
26.   Max dimension size of a thread block (x,y,z): (1024, 1024, 64)
27.   Max dimension size of a grid size    (x,y,z): (2147483647, 65535, 65535)
28.   Maximum memory pitch:                          2147483647 bytes
29.   Texture alignment:                             512 bytes
30.   Concurrent copy and kernel execution:          Yes with 3 copy engine(s)
31.   Run time limit on kernels:                     Yes
32.   Integrated GPU sharing Host Memory:            No
33.   Support host page-locked memory mapping:       Yes
34.   Alignment requirement for Surfaces:            Yes
35.   Device has ECC support:                        Disabled
36.   Device supports Unified Addressing (UVA):      Yes
37.   Device supports Compute Preemption:            Yes
38.   Supports Cooperative Kernel Launch:            Yes
39.   Supports MultiDevice Co-op Kernel Launch:      Yes
40.   Device PCI Domain ID / Bus ID / location ID:   0 / 3 / 0
41.   Compute Mode:
42.      < Default (multiple host threads can use ::cudaSetDevice() with device simultaneously) >
44. Device 1: "Quadro RTX 4000"
45.   CUDA Driver Version / Runtime Version          12.2 / 10.1
46.   CUDA Capability Major/Minor version number:    7.5
47.   Total amount of global memory:                 7974 MBytes (8361738240 bytes)
48.   (36) Multiprocessors, ( 64) CUDA Cores/MP:     2304 CUDA Cores
49.   GPU Max Clock rate:                            1545 MHz (1.54 GHz)
50.   Memory Clock rate:                             6501 Mhz
51.   Memory Bus Width:                              256-bit
52.   L2 Cache Size:                                 4194304 bytes
53.   Maximum Texture Dimension Size (x,y,z)         1D=(131072), 2D=(131072, 65536), 3D=(16384, 16384, 16384)
54.   Maximum Layered 1D Texture Size, (num) layers  1D=(32768), 2048 layers
55.   Maximum Layered 2D Texture Size, (num) layers  2D=(32768, 32768), 2048 layers
56.   Total amount of constant memory:               65536 bytes
57.   Total amount of shared memory per block:       49152 bytes
58.   Total number of registers available per block: 65536
59.   Warp size:                                     32
60.   Maximum number of threads per multiprocessor:  1024
61.   Maximum number of threads per block:           1024
62.   Max dimension size of a thread block (x,y,z): (1024, 1024, 64)
63.   Max dimension size of a grid size    (x,y,z): (2147483647, 65535, 65535)
64.   Maximum memory pitch:                          2147483647 bytes
65.   Texture alignment:                             512 bytes
66.   Concurrent copy and kernel execution:          Yes with 3 copy engine(s)
67.   Run time limit on kernels:                     Yes
68.   Integrated GPU sharing Host Memory:            No
69.   Support host page-locked memory mapping:       Yes
70.   Alignment requirement for Surfaces:            Yes
71.   Device has ECC support:                        Disabled
72.   Device supports Unified Addressing (UVA):      Yes
73.   Device supports Compute Preemption:            Yes
74.   Supports Cooperative Kernel Launch:            Yes
75.   Supports MultiDevice Co-op Kernel Launch:      Yes
76.   Device PCI Domain ID / Bus ID / location ID:   0 / 166 / 0
77.   Compute Mode:
78.      < Default (multiple host threads can use ::cudaSetDevice() with device simultaneously) >
79. > Peer access from Quadro RTX 4000 (GPU0) -> Quadro RTX 4000 (GPU1) : Yes
80. > Peer access from Quadro RTX 4000 (GPU1) -> Quadro RTX 4000 (GPU0) : Yes
82. deviceQuery, CUDA Driver = CUDART, CUDA Driver Version = 12.2, CUDA Runtime Version = 10.1, NumDevs = 2
83. Result = 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代码输出的结果为：

1. $ python3 test_add.py
2. 3.3193600177764893
3. (104857600,)
4. 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
  

来源：程序园用户自行投稿发布，如果侵权，请联系站长删除
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！