CUDA 编程笔记：使用 CUDA 加速矩阵乘法

一、前言

        学习 CUDA 编程的第一个课题，就是使用 GPU 来完成矩阵乘法的计算。大致原理，可以概括为：每个坐标相互独立的乘法和加法，从 CPU 的逐个计算，改为使用 GPU 并行计算。

        硬件平台：NVIDIA Jetson nano；

        编译器及：NVCC V10.0.326；

        Jetpack：Jetpack 4.3 [L4T 32.3.1]；

二、CPU 矩阵计算

        先放源码：

void cpu_matrix_mult(int *h_a, int *h_b, int *h_result, int m, int n, int k) {for (int i = 0; i < m; ++i) // a 矩阵按行遍历{for (int j = 0; j < k; ++j) // b 矩阵按列遍历{int tmp = 0.0;for (int h = 0; h < n; ++h) // 相乘相加{// h_a[i][h] = h_a[i * n + h];// h_b[h][j] = h_b[h * k + j];tmp += h_a[i * n + h] * h_b[h * k + j];}// h_result[i][j] = h_result[i * k + j];h_result[i * k + j] = tmp;}}
}

        输入参数分别是：矩阵a，矩阵b，输出矩阵；m是矩阵a行数， n是矩阵a列数（也是矩阵b的行数）， k是矩阵b列数。

        三个 for 循环不难理解，根据矩阵乘法规则梳理一遍就清楚了。

        在传统的 CPU 矩阵计算中，数据的存储方式均使用二维数组，但是在 CUDA 编程中，数据都是以一维的形式保存和读取，所以为了方便迁移和理解，这里我做了一个二维转一维的处理，可以看注释的代码行进行对应。

        时间复杂度可想而知，要是 Block Size 为 100 以内，计算还算不慢，要是增加到 1000，可就太慢了。

三、CUDA 矩阵计算

        首先将对应的矩阵计算部分做一个简述：

__global__ void gpu_matrix_mult(int *a, int *b, int *c, int m, int n, int k) {int row = blockIdx.y * blockDim.y + threadIdx.y; int col = blockIdx.x * blockDim.x + threadIdx.x;int sum = 0;if(col < k && row < m) {for(int i = 0; i < n; i++) {sum += a[row * n + i] * b[i * k + col];}c[row * k + col] = sum;}
}

        这是 CUDA 编程中的核函数，也叫GPU线程（Thread）可运行函数，有点类似 OpenCL 的编程逻辑。

        输入参数相同，参数 row 和 col 用于计算每一个元素在整个 Grid 中的绝对索引，由于该函数由每个线程独立执行，所以每个线程的 row 和 col 都是不一样的。再往后是判断取值越界，然后是独立的一次行列相乘相加操作，结果保存在输出矩阵中。矩阵的索引需要经过一个二维到一维的转换，读者需要注意。

        这里涉及到 thread,block,grid 的概念，详见：

CUDA软件架构—网格（Grid）、线程块（Block）和线程（Thread）的组织关系以及线程索引的计算公式_grid、block 和 thread 的关系-CSDN博客https://blog.csdn.net/dcrmg/article/details/54867507

        有关前缀 __global__ 的详细解释，可以看这篇文章：

cuda 函数前缀 __host__ __device__ __global__ ____noinline__ 和 __forceinline__ 简介_force inline-CSDN博客https://blog.csdn.net/zdlnlhmj/article/details/104896470

四、完整源码

        最关键的就是多线程编程思维和矩阵索引的转换，剩下就是常规的内存申请，数据拷贝等异构计算和数据交换常用步骤。

        由于存在内存动态申请和异常检测，所以我使用了 goto 语句。另外有个问题需要注意，nvcc 编译器下，goto 语句后面不能有变量申明，否则会报错，具体原因我也不是很清楚。下面是 matrix_mul.cu 源文件，可以根据代码中的注释结合来看。

        需要关注一下核函数的调用：

gpu_matrix_mult<<<dimGrid, dimBlock>>>(d_a, d_b, d_c, m, n, k);

dimGrid, dimBlock 分别是线程的组织结构。以矩阵大小 100 为例，dimGrid计算后为 7，dimBlock 为固定值 16，那么这里会使用 16 * 16 * 7 * 7 = 12544 个线程，但矩阵元素有 100 * 100 = 10000个，所以实际参与计算的线程只有 10000 个，剩余的线程会由于核函数内的边界检测而快速跳过。这主要是由于 GPU 硬件的 Grid 和 Block 线程组织方式，需要覆盖所有矩阵元素，所以会有一定冗余。

matrix_mul.cu

#include <stdio.h>
#include <math.h>
#include "error.cuh"#define BLOCK_SIZE 16/*** @Description: m: 矩阵A行数。 n: 矩阵A列数（也是矩阵B的行数）。 k: 矩阵B列数* @return {*}*/
void cpu_matrix_mult(int *h_a, int *h_b, int *h_result, int m, int n, int k) {// for (int i = 0; i < A_rows; i++) { // A 矩阵按行遍历//     for (int j = 0; j < B_cols; j++) { // B 矩阵按列遍历//         for (int h = 0; h < A_cols; h++) { // 相乘相加//             C[i][j] += A[i][h] * B[h][j];//         }//     }for (int i = 0; i < m; ++i) {for (int j = 0; j < k; ++j) {int tmp = 0.0;for (int h = 0; h < n; ++h) {// h_a[i][h] = h_a[i * n + h];// h_b[h][j] = h_b[h * k + j];tmp += h_a[i * n + h] * h_b[h * k + j];}// h_result[i][j] = h_result[i * k + j];h_result[i * k + j] = tmp;}}
}__global__ void gpu_matrix_mult(int *a, int *b, int *c, int m, int n, int k) {int row = blockIdx.y * blockDim.y + threadIdx.y; int col = blockIdx.x * blockDim.x + threadIdx.x;int sum = 0;if(col < k && row < m) {for(int i = 0; i < n; i++) {sum += a[row * n + i] * b[i * k + col];}c[row * k + col] = sum;}
}int main() {int m = 100;int n = 100;int k = 100;int ret = 0;unsigned int grid_rows;unsigned int grid_cols;dim3 dimGrid;dim3 dimBlock;int *h_a, *h_b, *h_c, *h_c_check;// 分配 Host 固定内存（强制让系统在物理内存中完成内存申请）// 固定内存的分配有可能会失败，所以需要检查返回值cudaError_t status;status = cudaMallocHost((void **) &h_a, sizeof(int)*m*n); // 矩阵A rows * colsif (status != cudaSuccess) {printf("Error allocating pinned h_a memoryn");return 1;}status = cudaMallocHost((void **) &h_b, sizeof(int)*n*k); // 矩阵B rows * colsif (status != cudaSuccess) {printf("Error allocating pinned h_b memoryn");ret = 1;goto free_h_a;}status = cudaMallocHost((void **) &h_c, sizeof(int)*m*k); // 输出矩阵C rows(A) * cols(B)if (status != cudaSuccess) {printf("Error allocating pinned h_c memoryn");ret = 1;goto free_h_b;}status = cudaMallocHost((void **) &h_c_check, sizeof(int)*m*k);if (status != cudaSuccess) {printf("Error allocating pinned h_c_check memoryn");ret = 1;goto free_h_c;}// 随机生成矩阵 A 和 B 的元素for (int i = 0; i < m; ++i) {for (int j = 0; j < n; ++j) {h_a[i * n + j] = rand() % 1024;}}for (int i = 0; i < n; ++i) {for (int j = 0; j < k; ++j) {h_b[i * k + j] = rand() % 1024;}}int *d_a, *d_b, *d_c;// 分配 Device 固定内存status = cudaMalloc((void **) &d_a, sizeof(int)*m*n);if (status != cudaSuccess) {printf("Error allocating pinned d_a memoryn");ret = 1;goto free_h_c_check;}status = cudaMalloc((void **) &d_b, sizeof(int)*n*k);if (status != cudaSuccess) {printf("Error allocating pinned d_b memoryn");ret = 1;goto free_d_a;}status = cudaMalloc((void **) &d_c, sizeof(int)*m*k);if (status != cudaSuccess) {printf("Error allocating pinned d_c memoryn");ret = 1;goto free_d_b;}// 拷贝至 device memoryCHECK(cudaMemcpy(d_a, h_a, sizeof(int)*m*n, cudaMemcpyHostToDevice));CHECK(cudaMemcpy(d_b, h_b, sizeof(int)*n*k, cudaMemcpyHostToDevice));// 计算网格在行和列方向的 block 数量// (m + BLOCK_SIZE - 1) / BLOCK_SIZE 通过向上取整确保所有数据被覆盖// 补足的范围是 ​​[1, BLOCK_SIZE - 1]​​，也就是确保只增加一个 BLOCK_SIZE 的数据grid_rows = (m + BLOCK_SIZE - 1) / BLOCK_SIZE;grid_cols = (k + BLOCK_SIZE - 1) / BLOCK_SIZE;// dim3 dimGrid(grid_cols, grid_rows); // goto 语句后不能申明变量（nvcc 编译器）// dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);dimGrid = dim3(grid_cols, grid_rows); // 顺序为 x, y, z，不填时默认为 1dimBlock = dim3(BLOCK_SIZE, BLOCK_SIZE);gpu_matrix_mult<<<dimGrid, dimBlock>>>(d_a, d_b, d_c, m, n, k);    CHECK(cudaMemcpy(h_c, d_c, sizeof(int)*m*k, cudaMemcpyDeviceToHost));//cudaThreadSynchronize();cpu_matrix_mult(h_a, h_b, h_c_check, m, n, k);for (int i = 0; i < m; ++i) {for (int j = 0; j < k; ++j) {// 相差过大，则认为错误if (fabs(h_c_check[i * k + j] - h_c[i * k + j]) > (1.0e-10)) {ret = 1;}}}if (!ret)printf("Pass!!!\n");elseprintf("Error!!!\n");cudaFree(d_c);
free_d_b:cudaFree(d_b);
free_d_a:cudaFree(d_a);
free_h_c_check:    cudaFreeHost(h_c_check);
free_h_c:cudaFreeHost(h_c);
free_h_b:cudaFreeHost(h_b);
free_h_a:cudaFreeHost(h_a);return ret;
}

error.cuh

        头文件内定义了一个错误检测函数：

#pragma once
#include <stdio.h>#define CHECK(call)                                   \
do                                                    \
{                                                     \const cudaError_t error_code = call;              \if (error_code != cudaSuccess)                    \{                                                 \printf("CUDA Error:\n");                      \printf("    File:       %s\n", __FILE__);     \printf("    Line:       %d\n", __LINE__);     \printf("    Error code: %d\n", error_code);   \printf("    Error text: %s\n",                \cudaGetErrorString(error_code));          \exit(1);                                      \}                                                 \
} while (0)

最后使用 NVCC 编译运行即可：

nvcc matrix_mul.cu -o matrix_mul
./matrix_mul