CUDA实现任意行列大小常数矩阵相乘

大江东第一深情

已于 2023-05-23 11:13:01 修改

阅读量259

点赞数 1

分类专栏： CUDA 文章标签：矩阵算法 c++ linux

于 2023-05-17 22:43:55 首次发布

本文链接：https://blog.csdn.net/Hughie999/article/details/130736463

版权

CUDA 专栏收录该内容

3 篇文章 0 订阅

订阅专栏

如果对您有帮助请不要吝啬您的点赞哦，还有其他相关问题可以直接在下方留言，笔者看到会及时回复。

1. 前言

楼主因科研需求需要实现高性能计算，但C++及CUDA编程基础几乎为0，所以花了很大的精力来学习，为了让后来学习的朋友能够尽快的扫除简单的障碍，我将自己的CUDA学习过程分享出来，无论是Debug还是自己接触到的新知识都会相应写一点博客，希望对大家的学习有所帮助。

2. 构思

2.1 main函数

首先，我们来看main函数：

1. 规定矩阵行列大小，即：n、theta、len；

2. 创造两个数组，数组大小即为数组本身行列数的相乘，数组内容为每个元素对应的行列索引值的和，例如第一行第一列的元素，索引值为0和0（C++下标索引从0开始），所以其值为0，第一行第二列的元素索引为0和1，其值为1；

3. 建立Device上的指针，包括sig、time_delay、res，并对其分配在Device上的内存，内存大小为：指针指向数组的大小乘上数组所含元素类型的字节数，使用cudaMalloc函数实现内存分配；

4. 将在Host端的数组复制到Device上准备进行计算，同样设置内存大小，标识复制方向，这里是cudaMemcpyHostToDevice；

5. 分配进行计算的网格和线程块大小，dimBlock设置的是线程块的维度，指明包含线程数量；dimGrid设置网格大小，指明网格包含线程块数量；

6. 使用核函数<<< >>>进行计算；

7. 将计算结果复制回Host端，这时候和之前的复制反过来，还是要设置复制的内存大小为多少，方向是cudaMemcpyDevcieToHost；

8. 输出结果。

2.2 核函数

然后我们来看核函数：

1. 首先，核函数应该是无返回类型函数，因为我们是直接将结果通过cudaMemcpy函数复制回Host端；

2. 其输入设置需要包含所有我们在计算中需要的数组和变量，比如这里就是sig、time_delay、res，还有行列数theta、len、n；

3. row和col的索引分别代表我们最终获得数组的行和列数，这个是和核函数的网格、线程块大小设置对应的，简单解释一下：blockDim.y * blockIdx.y + threadIdx.y是索引y方向上的线程，blockDim.x * blockIdx.x + threadIdx.x是索引x方向上的线程，不清楚的话同样照抄就行；

4. 然后进行数组计算，可以看到在核函数内的循环是可以用if代替的，正常来说我们在Host端直接进行行列式计算时需要三个for循环，这里是将外部的两个for循环替换为了线程的索引；

5. 这里我是将二维数组用一维来表示的，由于C++基础薄弱，我对二维数组的理解不是太好，就用了最笨的办法，仅为了实现矩阵乘法目的，大家下来可以试一试res[row][col]的方法。

3. 代码及输出

#include <stdio.h>
#include <cuda_runtime.h>
#include <iostream>

using namespace std;

/*  常数的矩阵乘法  */   

__global__ void matrix_mul(int* sig, int* time_delay, int* res,int theta,int len,int n)
{
	int row = blockDim.y*blockIdx.y + threadIdx.y;
	int col = blockDim.x*blockIdx.x + threadIdx.x;

    if(row<theta && col<len)
    {
        for(int i = 0 ; i < n ;i++)
        {
            //theta x len大小的矩阵res = time_delay矩阵的row行i列 X sig矩阵i行col列
            res[row*len+col] += time_delay[row*n+i] * sig[col+i*len];  
        }   
    }
}
 
int main(int argc, char **argv)
{
    int n = 2;
    int theta = 3;
    int len = 3;

    printf("第一个矩阵是%d行%d列的矩阵A:\n",theta,n);
    int t_delay[theta*n];
    for(int i = 0 ;i<theta;i++)
    {
        for(int j = 0;j<n;j++)
        {
            t_delay[i*n+j] = i*n+j;
            cout << t_delay[i*n+j] << '\t';
        }
        cout << endl;
    }

    printf("第二个矩阵是%d行%d列的矩阵B:\n",n,len);
    int sig_[len*n];
    for(int i = 0 ;i<n;i++)
    {
        for(int j = 0;j<len;j++)
        {
            sig_[i*len+j] = i*len+j;
            cout << sig_[i*len+j] << '\t';
        }
        cout << endl;
    }


    
    int sig_out[theta*len];
    int *time_delay;
    int *sig;
    int *res;

	cudaMalloc((void**)(&time_delay), theta*n*sizeof(int));
    cudaMalloc((void**)(&sig), n*len*sizeof(int));
	cudaMalloc((void**)(&res), theta*len*sizeof(int));  //从dc开始分配ROWS*COLS*sizeof(int)大小的地址

 	cudaMemcpy(sig, sig_, n*len*sizeof(int), cudaMemcpyHostToDevice);
	cudaMemcpy(time_delay, t_delay, theta*n*sizeof(int), cudaMemcpyHostToDevice);
	dim3 dimBlock(32,32);
	dim3 dimGrid((len+dimBlock.x-1)/dimBlock.x,(theta+dimBlock.y-1)/dimBlock.y); //分配线程块大小

	matrix_mul<<<dimGrid, dimBlock>>>(sig,time_delay,res,theta,len,n);

	cudaMemcpy(sig_out,res, theta*len*sizeof(int), cudaMemcpyDeviceToHost);

    printf("得到的矩阵是%d行%d列的矩阵C:\n",theta,len);
    for(int i = 0;i<theta;i++)
    {
        for(int j = 0 ;j<len;j++)
        {
            cout << sig_out[i*len+j]<<'\t';
        }
        cout <<endl;
    }
	cudaFree(sig);
	cudaFree(time_delay);
    cudaFree(res);

	return 0;
}

输出：