CUDA教程一、认识cuda

什么是cuda

统一计算设备架构（Compute Unified Device Architecture, CUDA），是由NVIDIA推出的通用并行计算架构。解决的是用更加廉价的设备资源，实现更高效的并行计算。

点击下面链接就可以下载cuda。我个人使用的是10.2版，截止到目前官方已经发布了11.0版。

有人就问了，std::thread它不香吗，为什么要用cuda？

别忘了，cuda是英伟达（NVIDIA）推出的——这个英伟达是何方神圣？没听说过英伟达，也应该听说过GPU了吧。没错，提出GPU概念的，正是英伟达。和中央处理器（Central Processing Unit, CPU）相对，图形处理器（Graphics Processing Unit, GPU）是显卡的核心芯片。而cuda正是暴露了英伟达开发的GPU的编程接口。

几乎所有的编程语言，不使用特定框架，都只能实现CPU编程——std::thread也是将线程开在CPU中。而不同于每一位程序员都接触过的CPU编程，GPU编程可以使用更多的流处理器、更多的线程数。

举个例子，打开我们的任务管理器，我们可以查看自己电脑CPU的线程数：

而GPU的核数和线程数更是丰富，只要合理使用，几乎没有用得完的情况：

如何合理使用，才最高效、利用率最高呢？

从矩阵相乘说起

举一个简单的例子，计算两个矩阵AB相乘的结果：

我们知道：

如果使用单线程编程，我们可以这么写：

 1 template <typename T>
 2 void MatrixProduct(matrix<T>& C /* Out */, matrix<T> const& A /* In */, matrix<T> const& B /* In */) {
 3     assert(A.Width() == B.Height());
 4     C.Set_size(A.Height(), B.Width());
 5     for(size_t i = 0; i < A.Height(); ++i) {
 6         for(size_t j = 0; j < B.Width(); ++j) {
 7             T ans = 0;
 8             for(size_t k = 0; k < A.Width(); ++k) {
 9                 ans += A[i][k] * B[k][j];
10             }
11             C[i][j] = ans;
12         }
13     }
14 }

这是毋庸置疑的，相信即使是编程初学者也可以写出。当然它在局部性上有待优化，不过优化空间非常非常小。

如果是多线程编程呢？

由于多个线程同时写入同一个变量是很危险的，但同时读取同一个变量不会出现问题，所以如果不想引入非常耗时的互斥锁，我们就必须根据可写入变量的数量来划分线程数——比如上式两个矩阵的乘积是 $2\times2$ 的方阵，因此我们可以分成四个线程来计算，每个线程只计算结果矩阵的一个元素。代码如下：

 1 template <typename T>
 2 void thrd_MatrixProduct(matrix<T>& C /* Out */, matrix<T> const& A /* In */, matrix<T> const& B /* In */,
 3             size_t row, size_t col) {
 4 //每一个线程因row和col的不同而被区分
 5     T ans = 0;
 6     for(size_t k = 0; k < A.Width(); ++k) {
 7         ans += A[row][k] * B[k][col];
 8     }
 9     C[row][col] = ans;
10 }
11 
12 template <typename T>
13 void MatrixProduct(threadPool& tp, matrix<T>& C /* Out */, matrix<T> const& A /* In */, matrix<T> const& B /* In */) {
14     assert(A.Width() == B.Height());
15     C.Set_size(A.Height(), B.Width());
16 //将每个计算单元的任务放进线程池中
17     for(size_t i = 0; i < A.Height(); ++i) {
18         for(size_t j = 0; j < B.Width(); ++j) {
19             tp.addTask(thrd_MatrixProduct, C, A, B, i, j);
20         }
21     }
22 //线程池启动并阻塞等待
23     tp.start();
24     tp.join();
25 }

要注意的是，threadPool和matrix类都不是STL，需要自行实现。其中threadPool可以用STL封装的线程std::thread或者C语言的pthread库来实现。std::thread类的声明中delete了拷贝构造函数和赋值运算符重载，这一点请在实现阶段注意。

事实上，这个操作在数据量较小的时候看不出优势；而数据量过大时，线程调度又是操作系统的难题。我们需要可以容纳更多线程的廉价计算资源。

cuda正是给显卡计算这一廉价而高效的并行计算方式提供了接口，同时也不需要线程池的维护。比如上述问题，用cuda实现，或许过程有点复杂，但核心还是相当容易的：

 1 template <typename T>
 2 __global__ void kern_MatrixProduct(T* C_ptr, const T* A_ptr, const T* B_ptr,
 3                    size_t A_Height, size_t A_Width, size_t B_Width) {
 4 //获取当前线程调用者的Id
 5     size_t Idx = blockIdx.x * blockDim.x + threadIdx.x;
 6     if(Idx >= A_Height * B_Width) return;
 7     
 8 //规定计算第i行j列的线程Id为(i * width + j)
 9     size_t row = Idx / B_Width;
10     size_t col = Idx % B_Width;
11 //核心计算
12     T ans = 0;
13     for(size_t k = 0; k < A_Width; ++k) {
14         ans += A_ptr[row * A_Width + k] * B_ptr[k * B_Width + col];
15     }
16     C_ptr[row * B_Width + col] = ans;
17 }
18 
19 template <typename T>
20 void MatrixProduct(matrix<T>& C /* Out */, matrix<T> const& A /* In */, matrix<T> const& B /* In */) {
21     assert(A.Width() == B.Height());
22     C.Set_size(A.Height(), B.Width());
23     
24     T *A_ptr, *B_ptr, *C_ptr;
25 //在显卡中申请内存（显存）
26     assert(cudaSuccess == cudaMalloc(&A_ptr, sizeof(T) * A.Height() * A.Width()));
27     assert(cudaSuccess == cudaMalloc(&B_ptr, sizeof(T) * B.Height() * B.Width()));
28     assert(cudaSuccess == cudaMalloc(&C_ptr, sizeof(T) * A.Height() * B.Width()));
29 //将数据从内存复制到显存
30     assert(cudaSuccess == cudaMemcpy(A_ptr, A.pos_ptr(), sizeof(T) * A.Height() * A.Width(), cudaMemcpyHostToDevice));
31     assert(cudaSuccess == cudaMemcpy(B_ptr, B.pos_ptr(), sizeof(T) * B.Height() * B.Width(), cudaMemcpyHostToDevice));
32 //调用核函数
33     size_t thread_count = 1024;
34     size_t block_count = (A.Height() * B.Width() - 1) / thread_count + 1;
35     kern_MatrixProduct<<<block_count, thread_count>>> (C_ptr, A_ptr, B_ptr, A.Height(), A.Width(), B.Width());
36     cudaDeviceSynchronize();
37     assert(cudaSuccess == cudaGetLastError());
38 //将数据从显存复制到内存
39     assert(cudaSuccess == cudaMemcpy(C.pos_ptr(), C_ptr, sizeof(T) * C.Height() * C.Width(), cudaMemcpyDeviceToHost));
40 //释放临时内存
41     assert(cudaSuccess == cudaFree(A_ptr));
42     assert(cudaSuccess == cudaFree(B_ptr));
43     assert(cudaSuccess == cudaFree(C_ptr));
44 }

当然，如果你在实现matrix的时候就把其元素数组申请在显存中，这个代码就会变得更加容易：

 1 template <typename T>
 2 __global__ void kern_MatrixProduct(T* C_ptr, const T* A_ptr, const T* B_ptr,
 3                    size_t A_Height, size_t A_Width, size_t B_Width) {
 4 //获取当前线程调用者的Id
 5     size_t Idx = blockIdx.x * blockDim.x + threadIdx.x;
 6     if(Idx >= A_Height * B_Width) return;
 7     
 8 //规定计算第i行j列的线程Id为(i * width + j)
 9     size_t row = Idx / B_Width;
10     size_t col = Idx % B_Width;
11 //核心计算
12     T ans = 0;
13     for(size_t k = 0; k < A_Width; ++k) {
14         ans += A_ptr[row * A_Width + k] * B_ptr[k * B_Width + col];
15     }
16     C_ptr[row * B_Width + col] = ans;
17 }
18 
19 template <typename T>
20 void MatrixProduct(matrix<T>& C /* Out */, matrix<T> const& A /* In */, matrix<T> const& B /* In */) {
21     assert(A.Width() == B.Height());
22     C.Set_size(A.Height(), B.Width());
23 //调用核函数
24     size_t thread_count = 1024;
25     size_t block_count = (A.Height() * B.Width() - 1) / thread_count + 1;
26     kern_MatrixProduct<<<block_count, thread_count>>> (C.pos_ptr(), A.pos_ptr(), B.pos_ptr(), A.Height(), A.Width(), B.Width());
27     cudaDeviceSynchronize();
28     assert(cudaSuccess == cudaGetLastError());
29 }

最后一种写法效率是最高的——第一种没有使用并行，第二种的并行优势不明显（而且申请线程太多，消耗大量时间），第三种I/O的耗时过大。

Windows10系统下的实验表明，当矩阵长宽小于10时，1、2、4三种写法的函数都能在0.01秒内完成，第三种也能在0.02秒内计算完毕。而当矩阵长宽为512时，第一种在1.1s左右计算完毕，第二种10s内计算不出结果，第三种则在1.3s左右结束，而第四种只有0.03s。

如果说即使如此也看不出1、3、4三种方式的孰优孰劣的话，试想一下，如果是若干个矩阵连乘，效果会是如何？结果不言而喻：第三种最慢，因为频繁调用I/O；第一种次之，因为没有并行；第四种的优势将会极其明显，因为I/O少、并行效率高。

很多人好奇，具体是如何并行的呢？上文代码中的三层尖括号（<<< >>>）是什么？莫急，让我一一解释。

首先我们来了解显卡的结构。

显卡的简化结构

显卡内部，有三级结构：网格（grid）、块（block）、线程（thread）。

每个显卡只有很少的网格，一个核函数目前只能运行在一个网格中，而一个网格里有多个块，每个块包含了若干线程。具体包含了多少，可以使用上文图中的deviceQuery来查询，通常是256、1024这样的值。

同样是这个矩阵相乘的例子：

上文的核函数在执行时，选取了一个网格中的一个块进行并行计算：

执行核函数时GPU一个执行块的简化示意图，大多数线程都被直接return掉

关于显卡的介绍到此为止，太过深入的知识我将不作介绍。

如何将计算部署在这些块和线程之上呢？这就是核函数（kernel function）及其运行时参数（runtime parameter）的工作。

核函数及其运行时参数

相信大家也猜到了，核函数就是上文代码里有__global__修饰的函数，其运行时参数就是放在三重尖括号<<< >>>之中的值。

每次调用核函数，都要指定其运行时参数，以便在显卡中弹性地部署计算。而其函数参数和我们学过的C/C++写法是一样的，不过要写在运行时参数的后面。

运行时参数有两种写法：

整数的写法：第一个参数是部署的核数 $b$ ，第二个参数是每个核分配的线程数 $t$ ，即在调用时：

kernelFunc<<<b, t>>> (arg1, arg2);

　　注意不能超过每个核的最大线程数量，否则cudaGetLastError会返回一个

　　cudaErrorInvalidConfiguration错误。运行时通过blockIdx.x * blockDim.x + threadIdx.x来获取线程ID，其区间范围为。如果的值大于你的预期，请在不必要的线程中及时return，以免出现数组越界、访问无效内存等情况。有人可能会问，我挨个指定每个核开辟多少线程不行吗？不行。核函数的线程，往往都是要稍微开多一些，然后return掉不必要的。英伟达就是这么霸道，我不要你觉得，我要我觉得。

　　三维向量的写法：第一个参数指定了网格中的维度，即如何分配块；第二个参数指定了每个块的维度，即如何分配线程；第三个参数可选，指定了核函数运行时每个block除了使用共享显存外，还能再动态申请多少空间，默认是0，单位是字节——这和malloc函数参数的意义一样；第四个参数也是可选的，指定了该核函数运行在哪一个流中，默认是0——这个流方便主显异步并行，但相对麻烦，暂时不作介绍。

　　其中的类型是dim3，即三维无符号整数向量：x坐标表示该网格中沿x坐标的核数，即每行的核数；y坐标表示该网格中沿y坐标的核数，即每列的核数；z坐标只能为1，因为网格中每一个核都是平面摆放的。因此核函数部署的核数就是。

　　而的类型也是dim3：x坐标表示一个核中沿x坐标的线程数，即每行的线程数；y坐标表示一个核中沿y坐标的线程数，即每列的线程数；z坐标表示一个核中沿z坐标的线程数，即每个核的线程有多少个水平面。因此核函数在每个核开辟的线程数就是。线程ID的值的取值范围则在区间中。

　　而这些值都是有最大值的，同样也可以通过deviceQuery查询，我的显卡中五个参数的最大值分别为65536、65536、1024、1024、64。

　　而核函数运行在显卡中，这就限制了它访问内存中的变量——但我们可以将内存中的变量传递到显存中。下一篇文章将会介绍主存与显存的交互。

posted on 2022-02-11 15:17 一杯清酒邀明月阅读(2018) 评论(0) 编辑收藏举报