第六章 kenel矩阵计算实战篇(下篇)

前言

学习我的教程专栏，你将绝对能实现CUDA工程化，实现环境安装、index计算、kernel核函数编程、内存优化与steam性能优化、原子操作、nms的cuda算子、yolov5的cuda部署等内容，并开源教程源码。

到此为止，之前章节已让我们大致熟悉cuda相关理论与cuda函数运作原理，特别是kernel函数的线程index计算规则。为此，本节我们将正式使用cuda编写基本矩阵运算，主要涉及矩阵加法与乘法，而矩阵减法与除法与本节介绍的方法相似，将不在作为教程。我们选择四则运算也是因其具有一定实用性和代表性，矩阵加减运算是基于对应像素处理过程，矩阵乘除法运算相对较为复杂，是行与列运算逻辑。更重要的是，我们也通过cuda的基本运算帮助大家熟悉cuda编程的线程index使用的相关逻辑。

专栏概括

1、cuda教程目录

第一章 指针篇–>点击这里
第二章 CUDA原理篇–>点击这里
第三章 CUDA编译器环境配置篇–>点击这里
第四章 kernel函数基础篇–>点击这里
第五章 kernel索引(index)篇–>点击这里
第六章 kenel矩阵计算实战篇–>点击这里-上篇–点击这里-下篇

第七章 kenel实战强化篇–>点击这里
第八章 CUDA内存应用与性能优化篇–>点击这里-上篇–点击这里-中篇–点击这里-下篇
第九章 CUDA原子(atomic)实战篇–>点击这里
第十章 CUDA流(stream)实战篇–>点击这里
第十一章 CUDA的NMS算子实战篇–>点击这里-上篇–点击这里-下篇
第十二章 YOLO的部署实战篇–>点击这里-上篇–点击这里-下篇

第十三章 基于CUDA的YOLO部署实战篇–>点击这里

2、cuda教程背景

随着人工智能的发展与人才的内卷，很多企业已将深度学习算法的C++部署能力作为基本技能之一。面对诸多arm相关且资源有限的设备，往往想更好的提速，满足更高时效性，必将更多类似矩阵相关运算交给CUDA处理。同时，面对市场诸多教程与诸多博客岑子不起的教程或高昂教程费用，使读者(特别是小白)容易迷糊，无法快速入手CUDA编程，实现工程化。
因此，我将结合我的工程实战经验，我将在本专栏实现CUDA系列教程，帮助读者(或小白)实现CUDA工程化，掌握CUDA编程能力。学习我的教程专栏，你将绝对能实现CUDA工程化，完全从环境安装到CUDA核函数编程，从核函数到使用相关内存优化与流stream优化，从内存、stream优化到深度学习算子开发(如：nms)，从算子优化到模型(以yolo系列为基准)部署。最重要的是，我的教程将简单明了直切主题，CUDA理论与实战实例应用，并附相关代码，可直接上手实战。我的想法是掌握必要CUDA相关理论，去除非必须繁杂理论，实现CUDA算法应用开发，待进一步提高，将进一步理解更高深理论。

3、cuda教程内容

第一章到第三章探索指针在cuda函数中的作用与cuda相关原理及环境配置；

第四章初步探索cuda相关函数编写(global、device、__host__等)，实现简单入门；

第五章探索不同grid与block配置，如何计算kernel函数的index，以便后续通过index实现各种运算；

第六、七章由浅入深探索核函数矩阵计算，深入探索grid、block与thread索引对kernel函数编写作用与影响，并实战多个应用列子(如：kernel函数实现图像颜色空间转换)；

第八章探索cuda内存纹理内存、常量内存、全局内存等分配机制与内存实战应用(附代码)，通过不同内存的使用来优化cuda计算性能；

第九章探索cuda原子(atomic)相关操作，并实战应用(如:获得某些自加索引等)；

第十章探索cuda流stream相关应用，并给出相关实战列子(如:多流操作等)；

第十一到十三章探索基于tensorrt部署yolo算法，我们首先将给出通用tensorrt的yolo算法部署，该部署的前后处理基于C++语言的host端实现，然后给出基于cuda的前后处理的算子核函数编写，最后数据无需在gpu与host间复制操作，实现gpu处理，提升算法性能。

目前，以上为我们的cuda教学全部内容，若后续读者有想了解知识，可留言，我们将根据实际情况，更新相关教学内容。

大神忽略

源码链接地址点击这里
yolov5部署代码链接点击这里
yolov5的cuda部署代码链接：文中链接源码

基于我的代码实测，cuda部署yolov5加速10倍，只想说cuda太香了！！！

一、矩阵乘除法

矩阵除法实际就是kernel函数符号“*”变成符号“/”，其运算原理基本与矩阵乘法一致，基于此，我们只介绍矩阵乘法相关规则。

假设：矩阵 a[m,n] 与 b[n,k]
目的：实现矩阵a行与矩阵b列逐像素之积

二、矩阵乘除法编码步骤

1、初始化与内存分配

①、条件设置

假设m与n及block_size

    const int  BLOCK_SIZE = 2;
    int m = 8;    //行
    int n = 4;   //中间变量
    int k = 10;    //列

②、host分配内存与赋值

	int* a, * b;
    // 初始化 a与b
    cudaMallocHost((void**)&a, sizeof(int) * m * n);
    cudaMallocHost((void**)&b, sizeof(int) * n * k);
    std::cout << "value of a:" << endl;
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            a[i * n + j] = rand() % 6;
            std::cout << "\t" << a[i * n + j];
        }
        std::cout << "\n";
    }

    std::cout << "value of b:" << endl;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < k; j++) {
            b[i * k + j] = rand() % 10;
            std::cout << "\t" << b[i * k + j];
        }
        std::cout << "\n";
    }

变量a与变量b为host分配内存，for循环为a与b始化赋值。
变量a与b是指针，也就是我为什么在第一章写了一篇指针原由，后面基本host端数据给到device端均需指针传递。

③、device分配内存与赋值

    cudaMalloc((void**)&g_a, sizeof(int) * m * n);
    cudaMalloc((void**)&g_b, sizeof(int) * n * k);
    cudaMemcpy(g_a, a, sizeof(int) * m * n, cudaMemcpyHostToDevice);
    cudaMemcpy(g_b, b, sizeof(int) * k * n, cudaMemcpyHostToDevice);
    

此时,变量g_a与g_b通过cudaMalloc将其在device端开辟内存空间，在gpu为全局变量，在使用cudaMemcpy方式将host端变量a与变量b的值复制给g_a与g_b变量。至于如何进行，可参考我的专栏第二章cuda原理点击这里。

④、设置kernel相关配置

	unsigned int grid_rows = (m + BLOCK_SIZE - 1) / BLOCK_SIZE;  //行写给y
    unsigned int grid_cols = (k + BLOCK_SIZE - 1) / BLOCK_SIZE;  //列写给x
    dim3 dimGrid(grid_cols, grid_rows);
    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);

设置kernel函数的grid与block配置，这些配置设置影响一维索引index寻找(详情点击这里)。

⑤kernel函数调用

gpu_matrix_mult << <dimGrid, dimBlock >> > (g_a,  g_b, g_c,  m,  n,  k);	方法一 
gpu_matrix_multiply_thread << <dimGrid, dimBlock >> > (g_a, g_b, g_c2, m,n, k);	方法二

以上为kernel函数调用，即可实现矩阵乘除运算。在这里，我使用了二种方式编写kernel函数计算逻辑，与其说二种计算逻辑，更换成二种index方法使用更合理， 因其重要，我将在下面详细介绍。

三、矩阵乘除法多种kernel实现方法

1、kernel使用row与col作为index

我们分配grid_rows与grid_cols是按照m行n列方式，且grid与block均为2维度。在每一个线程内部，我们也可通过矩阵加法原理，获得该线程所在的行和列。同时，矩阵乘法是行与列的相乘。然而，a是m行n列，b是n行k列，a*b实际为a的每行n个数与b的每列n个数一 一相乘的和。
以下代码的grid与block设置为m行k列设计的，为此row与col对应为m行与k列。
代码如下：

__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; // 行线程 y
    int col = blockIdx.x * blockDim.x + threadIdx.x; // 列线程 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];  //看出row与col不动的方式计算
        }
        c[row * k + col] = sum;
    }
}

其中if (col < k && row < m)表示行和列不能超过此条件，a[row * n + i]表示行row不变，i从0便利到n-1的方法，为每行遍历；b[i * k + col]表示行col不变，i从0遍历到n-1的方法，为每列遍历；

int sum = 0;
   if (col < k && row < m) {
       for (int i = 0; i < n; i++) {
           sum += a[row * n + i] * b[i * k + col];  //看出row与col不动的方式计算
       }
       c[row * k + col] = sum;
   }

以上代码，可理解为每次都需要执行，或理解成有一个for循环外包即可。

2、kernel核使用线程id

我们分配grid_rows与grid_cols是按照m行n列方式，且grid与block均为2维度。在每一个线程内部，我们可以通过上一节索引方法，使用线程方式，实现2个矩阵的乘法。因其线程id索引上节有介绍，也将不在赘述。以下变量与上设置一样，仅仅将id转为行和列的思想计算。
代码如下：

_global__ void gpu_matrix_multiply_thread(int* a, int* b, int* c,int m, int n, int k)
{
   // [m*n]* [n*k]<---m n
   //方法一：通过id方式计算
   //grid为2维度，block为2维度,使用公式id=blocksize * blockid + threadid
   int blocksize = blockDim.x * blockDim.y;
   int blockid = gridDim.x * blockIdx.y + blockIdx.x;
   int threadid = blockDim.x * threadIdx.y + threadIdx.x;
   int id = blocksize * blockid + threadid;

   int row = id / k;
   int col = id % k;
   int sum = 0;

   for (int i = 0; i < n; i++) {
       sum += a[row*n + i] * b[i * k +col ];
   }
   c[row * k + col] = sum;

}

至于解释，我们将不在过多说明，其理解与上相似，仅将线程id转为row与col的逻辑处理。

四、矩阵乘除法完整代码

#include <iostream>
#include <time.h>
#include "opencv2/highgui.hpp"  //实际上在/usr/include下
#include "opencv2/opencv.hpp"
#include "device_launch_parameters.h"
#include <cuda_runtime_api.h>
using namespace cv;
using namespace std;

void Print_2dim(int* ptr, int m,int n) {
    std::cout << "result:\n";
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            std::cout <<"\t" << ptr[i * n + j];
        
        }
        std::cout << "\n";
    }
    
}


__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; // 行线程 y
    int col = blockIdx.x * blockDim.x + threadIdx.x; // 列线程 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];  //看出row与col不动的方式计算
        }
        c[row * k + col] = sum;
    }
}


__global__ void gpu_matrix_multiply_thread(int* a, int* b, int* c,int m, int n, int k)
{
    // [m*n]* [n*k]<---m n
    //方法一：通过id方式计算
    //grid为2维度，block为2维度,使用公式id=blocksize * blockid + threadid
    int blocksize = blockDim.x * blockDim.y;
    int blockid = gridDim.x * blockIdx.y + blockIdx.x;
    int threadid = blockDim.x * threadIdx.y + threadIdx.x;
    int id = blocksize * blockid + threadid;

    int row = id / k;
    int col = id % k;
    int sum = 0;

    for (int i = 0; i < n; i++) {
        sum += a[row*n + i] * b[i * k +col ];

    }
    c[row * k + col] = sum;


}


int kernel_multiply()
{   
    /*
    matrix a[m,n], matrix b[n,k]
    a[m,n]*b[n,k]=[m,k]
    */

    const int  BLOCK_SIZE = 2;
    int m = 8;    //行
    int n = 4;   //中间变量
    int k = 10;    //列
    int* a, * b;
   
    // 初始化 a与b
    cudaMallocHost((void**)&a, sizeof(int) * m * n);
    cudaMallocHost((void**)&b, sizeof(int) * n * k);
    
    std::cout << "value of a:" << endl;
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            a[i * n + j] = rand() % 6;
            std::cout << "\t" << a[i * n + j];
        }
        std::cout << "\n";
    }

    std::cout << "value of b:" << endl;
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < k; j++) {
            b[i * k + j] = rand() % 10;
            std::cout << "\t" << b[i * k + j];
        }
        std::cout << "\n";
    }
    
    //a*b相乘
    std::cout << "value of a*b:" << endl;
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < k; j++) {
            int tmp = 0;
            for (int h = 0; h < n; h++) {
                tmp += a[i * n + h] * b[h * k + j];
            }
            //c[i * k + j] = tmp;
            std::cout << "\t" << tmp;
        }
        std::cout << "\n";
    }




    int* g_a, * g_b;
  
    cudaMalloc((void**)&g_a, sizeof(int) * m * n);
    cudaMalloc((void**)&g_b, sizeof(int) * n * k);
    cudaMemcpy(g_a, a, sizeof(int) * m * n, cudaMemcpyHostToDevice);
    cudaMemcpy(g_b, b, sizeof(int) * k * n, cudaMemcpyHostToDevice);
    unsigned int grid_rows = (m + BLOCK_SIZE - 1) / BLOCK_SIZE;  //行写给y
    unsigned int grid_cols = (k + BLOCK_SIZE - 1) / BLOCK_SIZE;  //列写给x
    dim3 dimGrid(grid_cols, grid_rows);
    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
    std::cout << "gridDIM.x:" << grid_cols << "\tgridDIM.y:" << grid_rows << endl;
    std::cout << "blockDIM.x:" << BLOCK_SIZE << "\tblockDIM.y:" << BLOCK_SIZE << endl;

    // 使用row与col计算
    int* c1, * g_c;
    cudaMalloc((void**)&g_c, sizeof(int) * m * k);
    cudaMallocHost((void**)&c1, sizeof(int) * m * k);
    gpu_matrix_mult << <dimGrid, dimBlock >> > (g_a,  g_b, g_c,  m,  n,  k);
    cudaMemcpy(c1, g_c, sizeof(int) * m * k, cudaMemcpyDeviceToHost);
    Print_2dim(c1, m, k);
    cudaFree(g_c);
    cudaFreeHost(c1);
    
    //使用id计算
    int* c2, * g_c2;
    cudaMalloc((void**)&g_c2, sizeof(int) * m * k);
    cudaMallocHost((void**)&c2, sizeof(int) * m * k);
    gpu_matrix_multiply_thread << <dimGrid, dimBlock >> > (g_a, g_b, g_c2, m,n, k);
    cudaMemcpy(c2, g_c2, sizeof(int) * m * k, cudaMemcpyDeviceToHost);
    Print_2dim(c2, m, k);
    cudaFree(g_c2);
    cudaFreeHost(c2);

    //释放内存
    cudaFree(g_a);
    cudaFree(g_b);
    cudaFreeHost(a);
    cudaFreeHost(b);
    

    return 0;
}

五、总结

以上为矩阵乘除法相关实战应用，我个人认觉得矩阵乘除法比矩阵加减法较为复杂些，主要是掌握如何使用kernel编写矩阵a的行与矩阵b的列逐元素相乘，本节重点依然掌握矩阵kernel函数编码的几种索引方式与矩阵行与列逐元素相乘。